fft.c

三类FFT变换(时域基2

*Create Date:   24/07/2007
*Create By:     YiMin.Pang & Rong.Cui
*Modify Date:   
********************************************************/
void FFT_Radix_t2(complex *ptTD,int tTN,complex *tFD,int tFN)
{
	int i,j,k,l,power;
	int block,block_size;
	int tab_1,tab_2;
	complex *tTD;
	complex temp;

	block_size = 2;
	block = tFN / block_size;
	power = cal_power(tFN,2);
	//expand sequence
	tTD = expand(ptTD,tTN,tFN);
	//bit-reverse the tTD
	bit_reverse(tTD,tFN,power);

	//calculate tFN points FFT
	for(i=0;i<power;i++)
	{
		for(j=0;j<block;j++)
		{
			for(k=0;k<block_size/2;k++)
			{
				//the current pair of elements to calculate
				tab_1 = j * block_size + k;
				tab_2 = tab_1 + block_size / 2;
				//calculate the basic butterfly
				temp = complex_mul(tTD[tab_2],phase_factor(1,k,block_size));
				tFD[tab_1] = complex_add(tTD[tab_1],temp);
				tFD[tab_2] = complex_sub(tTD[tab_1],temp);
			}
		}
		block = block >> 1;
		block_size = block_size << 1;
		for(l=0;l<tFN;l++)
		{tTD[l] = tFD[l];}
	}
	free(tTD);
}

/********************************************************
*Function Name: FFT_Radix_f2
*Function:      Calculate length-N discrete fourier transform
				based on Radix-2 decimation-in-frequence.
*Iuput:         tTD[](time domain data)
				tFD[](array for frequence domain data)
				tFN(length of the FFT)
*Output:        
*Return:        
*Create Date:   09/08/2007
*Create By:     Rong.Cui
*Modify Date:   
********************************************************/
void FFT_Radix_f2(complex *ptTD, int tTN,complex *tFD,int tFN)
{
	int i,j,k,l,power;
	int block,block_size;
	int tab_1,tab_2;
	complex *tTD;
	complex temp;

	block_size = tFN;
	block = tFN / block_size;
	power = cal_power(tFN,2);

	//expand sequence
	tTD = expand(ptTD,tTN,tFN);

	//calculate tFN points FFT
	for(i=0;i<power;i++)
	{
		for(j=0;j<block;j++)
		{
			for(k=0;k<block_size/2;k++)
			{
				//the current pair of elements to calculate
				tab_1 = j * block_size + k;
				tab_2 = tab_1 + block_size / 2;
				//calculate the basic butterfly
				temp = complex_sub(tTD[tab_1],tTD[tab_2]);
				tFD[tab_1] = complex_add(tTD[tab_1],tTD[tab_2]);
				tFD[tab_2] = complex_mul(temp,phase_factor(1,k,block_size));
			}
		}
		block = block << 1;
		block_size = block_size >> 1;
		for(l=0;l<tFN;l++)
		{tTD[l] = tFD[l];}
	}
	//bit-reverse the tFD
	bit_reverse(tFD,tFN,power);
	free(tTD);
}

/********************************************************
*Function Name: SRFFT
*Function:      Calculate length-N discrete fourier transform
				by Duhamel's split-radix algorithms(1986,IEEE)
*Iuput:         tTD[](time domain data)
				tFD[](array for frequence domain)
				tFN(length of the FFT)
*Output:        
*Return:        
*Create Date:   09/08/2007
*Create By:     YiMin.Pang & Rong.Cui
*Modify Date:   
********************************************************/
void SRFFT(complex *ptTD,int tTN,complex *tFD,int tFN)
{
	int *block_size;
	int block=1;
	int power;
	int i,j,k,l;
	int tab_1,tab_2,tab_3,tab_4;
	int process_counter;

	complex *tTD;
	complex temp_1,temp_2;
	power = cal_power(tFN,2);
	block_size = (int *)malloc(sizeof(int));
	*block_size = tFN;

	//expand sequence
	tTD = expand(ptTD,tTN,tFN);

	//calculate tFN points IFFT
	for(i=0;i<power;i++)
	{
		process_counter = 0;
		for(j=0;j<block;j++)
		{
			if(block_size[j]>2)
			{
				//by SRFFT
				for(k=0;k<block_size[j]/4;k++)
				{
					//the current pair of elements to calculate
					tab_1 = process_counter + k;
				    tab_2 = tab_1 + block_size[j] / 4;
				    tab_3 = tab_2 + block_size[j] / 4;
				    tab_4 = tab_3 + block_size[j] / 4;
					//calculate the basic butterfly ("L" shaped)
				    tFD[tab_1] = complex_add(tTD[tab_1],tTD[tab_3]);
				    tFD[tab_2] = complex_add(tTD[tab_2],tTD[tab_4]);
			
				    temp_1  = complex_sub(tTD[tab_1],tTD[tab_3]);
				    temp_2 = complex_mul(complex_sub(tTD[tab_2],tTD[tab_4]),complex_set_value(0,1));
				    tFD[tab_3] = complex_mul(complex_sub(temp_1,temp_2),phase_factor(1,k,block_size[j]));
				    tFD[tab_4] = complex_mul(complex_add(temp_1,temp_2),phase_factor(1,3*k,block_size[j]));
				}
			}
			//by Radix-2 FFT to end last stage
			else if(block_size[j]==2)
			{
				//the current pair of elements to calculate
				tab_1 = process_counter;
				tab_2 = tab_1 + 1;
				//calculate the basic butterfly ("X" shaped)
				tFD[tab_1] = complex_add(tTD[tab_1],tTD[tab_2]);
				tFD[tab_2] = complex_sub(tTD[tab_1],tTD[tab_2]);
			}
			else if(block_size[j]<2)
			{}  //finished,nothing to do

			process_counter = process_counter + block_size[j];
		}
		//block state update for SRFFT
		updata_block(&block_size,&block);
		for(l=0;l<tFN;l++)
		{tTD[l] = tFD[l];}
	}
	//bit-reverse the tFD
	bit_reverse(tFD,tFN,power);
	free(tTD);
}


/*********************IFFT Functions**********************/

/********************************************************
*Function Name: IFFT_Radix_t2
*Function:      Calculate length-N inverse discrete fourier
                transform by Cooly-Turkey FFT algorithms,
				which based on Radix-2 decimation-in-time.
*Iuput:         tFD[](frequence domain data)
				tTD[](array for time domain data)
				tN(length of the IFFT)
*Output:        
*Return:        
*Create Date:   09/08/2007
*Create By:     YiMin.Pang & Rong.Cui
*Modify Date:   
********************************************************/
void IFFT_Radix_t2(complex *ptFD,complex *tTD,int tN)
{
	int i,j,k,l,power;
	int block,block_size;
	int tab_1,tab_2;
	complex *tFD;
	complex temp;

	//copy input sequence to tFD
	tFD = (complex *)malloc(tN * sizeof(complex));
	for(i=0;i<tN;i++)
	{tFD[i] = ptFD[i];}

	block_size = 2;
	block = tN / block_size;
	power = cal_power(tN,2);

	//bit-reverse the tFD
	bit_reverse(tFD,tN,power);

	//calculate tN points IFFT
	for(i=0;i<power;i++)
	{
		for(j=0;j<block;j++)
		{
			for(k=0;k<block_size/2;k++)
			{
				//the current pair of elements to calculate
				tab_1 = j * block_size + k;
				tab_2 = tab_1 + block_size / 2;
				//calculate the basic butterfly
				temp = complex_mul(tFD[tab_2],phase_factor(-1,k,block_size));
				tTD[tab_1] = complex_add(tFD[tab_1],temp);
				tTD[tab_2] = complex_sub(tFD[tab_1],temp);
			}
		}
		block = block >> 1;
		block_size = block_size << 1;
		for(l=0;l<tN;l++)
		{tFD[l] = tTD[l];}
	}

	//result divide by tN
	temp = complex_set_value( (1.0/tN), 0);
	for(l=0;l<tN;l++)
	{tTD[l] = complex_mul(tTD[l],temp);}

	free(tFD);
}
/********************************************************
*Function Name: IFFT_Radix_f2
*Function:      Calculate length-N inverse discrete fourier 
                transform based on Radix-2 decimation-in-frequence.
*Iuput:         tFD[](frequence domain data)
				tTD[](array for time domain data)
				tN(length of the FFT)
*Output:        
*Return:        
*Create Date:   10/08/2007
*Create By:     Rong.Cui
*Modify Date:   
********************************************************/
void IFFT_Radix_f2(complex *ptFD,complex *tTD,int tN)
{
	int i,j,k,l,power;
	int block,block_size;
	int tab_1,tab_2;
	complex *tFD;
	complex temp;

	//copy input sequence to tFD
	tFD = (complex *)malloc(tN * sizeof(complex));
	for(i=0;i<tN;i++)
	{tFD[i] = ptFD[i];}
	
	block_size = tN;
	block = tN / block_size;
	power = cal_power(tN,2);
	
	//calculate tN points IFFT
	for(i=0;i<power;i++)
	{
		for(j=0;j<block;j++)
		{
			for(k=0;k<block_size/2;k++)
			{
				//the current pair of elements to calculate
				tab_1 = j * block_size + k;
				tab_2 = tab_1 + block_size / 2;
				//calculate the basic butterfly
				temp = complex_sub(tFD[tab_1],tFD[tab_2]);
				tTD[tab_1] = complex_add(tFD[tab_1],tFD[tab_2]);
				tTD[tab_2] = complex_mul(temp,phase_factor(-1,k,block_size));
			}
		}
		block = block << 1;
		block_size = block_size >> 1;
		for(l=0;l<tN;l++)
		{tFD[l] = tTD[l];}
	}
	bit_reverse(tTD,tN,power);
	
	//result divide by tN
	temp = complex_set_value( (1.0/tN), 0);
	for(l=0;l<tN;l++)
	{tTD[l] = complex_mul(tTD[l],temp);}

	free(tFD);
}

/********************************************************
*Function Name: SRIFFT
*Function:      Calculate length-N inverse discrete fourier 
                transform by Duhamel's split-radix algorithms
				(1986,IEEE)
*Iuput:         tFD[](frequence domain data)
				tTD[](array for time domain data)
				tN(length of the FFT)
*Output:        
*Return:        
*Create Date:   10/08/2007
*Create By:     YiMin.Pang & Rong.Cui
*Modify Date:   
********************************************************/
void SRIFFT(complex *ptFD,complex *tTD,int tN)
{
	int *block_size;
	int block=1;
	int power;
	int i,j,k,l;
	int tab_1,tab_2,tab_3,tab_4;
	int process_counter;
	complex temp_1,temp_2,temp;
	complex *tFD;

	//copy input sequence to tFD
	tFD = (complex *)malloc(tN * sizeof(complex));
	for(i=0;i<tN;i++)
	{tFD[i] = ptFD[i];}

	power = cal_power(tN,2);
	block_size = (int *)malloc(sizeof(int));
	*block_size = tN;
	
	//calculate tN points IFFT
	for(i=0;i<power;i++)
	{
		process_counter = 0;
		for(j=0;j<block;j++)
		{
			//by SRIFFT
			if(block_size[j]>2)
			{
				for(k=0;k<block_size[j]/4;k++)
				{
					//the current pair of elements to calculate
					tab_1 = process_counter + k;
				    tab_2 = tab_1 + block_size[j] / 4;
				    tab_3 = tab_2 + block_size[j] / 4;
				    tab_4 = tab_3 + block_size[j] / 4;
					//calculate the basic butterfly ("L" shaped)
				    tTD[tab_1] = complex_add(tFD[tab_1],tFD[tab_3]);
				    tTD[tab_2] = complex_add(tFD[tab_2],tFD[tab_4]);
			
				    temp_1  = complex_sub(tFD[tab_1],tFD[tab_3]);
				    temp_2 = complex_mul(complex_sub(tFD[tab_2],tFD[tab_4]),complex_set_value(0,1));
				    tTD[tab_3] = complex_mul(complex_add(temp_1,temp_2),phase_factor(-1,k,block_size[j]));
				    tTD[tab_4] = complex_mul(complex_sub(temp_1,temp_2),phase_factor(-1,3*k,block_size[j]));
				}
			}
			//by Radix-2 IFFT to end last stage
			else if(block_size[j]==2)
			{
				//the current pair of elements to calculate
				tab_1 = process_counter;
				tab_2 = tab_1 + 1;
				//calculate the basic butterfly ("X" shaped)
				tTD[tab_1] = complex_add(tFD[tab_1],tFD[tab_2]);
				tTD[tab_2] = complex_sub(tFD[tab_1],tFD[tab_2]);
			}
			else if(block_size[j]<2)
			{}	//finished,nothing to do
			process_counter = process_counter + block_size[j];
		}
		//block state update for SRIFFT
		updata_block(&block_size,&block);
		for(l=0;l<tN;l++)
		{tFD[l] = tTD[l];}
	}
	//bit-reverse the tTD
	bit_reverse(tTD,tN,power);

	//result divide by tN
	temp = complex_set_value( (1.0/tN), 0);
	for(l=0;l<tN;l++)
	{tTD[l] = complex_mul(tTD[l],temp);}

	free(tFD);
}