ifft_r4_real_256.c

同样是浮点型的fft算法一样为vc开发的

/*
 *	DATE:			09-21-2006
 *
 *	AUTHOR:			CHEN PENG
 *
 *	VERSION:
 *		ORIGINAL	09-21-2006
 */
/**************************************************************************************\
*	函数名：ifft_r4_real_256
*
*	参数  ：x[]为输入的复数数据，实部和虚部类型为均为short
*			w[]为FFT时的旋转因子，类型为short
*			ls[]数组中每一个元素都表示对应一级FFT运算及输出分离时的左移位数
*			rs[]数组中每一个元素都表示对应一级FFT运算及输出分离时的右移位数
*			points为输出实数的长度
*
*	返回值：无
*
*	说明  ：1.该函数使用radix4和radix2混合的复数据的DIF FFT算法，主要针对256，1024
*			和64这类数据长度的实数输出数据。
*			2.输出数据占据输入数据空间，并且输出的实数点数为输入复数点数的二倍，输入
*			按照{实部0，虚部0，实部1，虚部1，...}排列，实部和虚部各占一个short。
*			3.旋转因子为points点的FFT旋转因子，按照{实部0，虚部0，实部1，虚部1，...}
*			排列,以Q来表示其标度为Q15。
*			4.该函数不检查ls[]和rs[]的长度，所以，调用函数时必须保证ls[]和rs[]的长度
*			不小于FFT的级数+1,ls[0]和rs[0]有效元素的值为输入分离时的移位数，这里ls[1]
*			和rs[1]表示radix2的移位数。
*
\**************************************************************************************/
//#include "stdio.h";
extern int reverseN[256];
void	ifft_r4_real_256( int x[], int xTemp[],
						  int w[],
						  int rs[],
						  int points,
						  int xHalfN
					    ){
	int n1, n2, n3, ie;
	int iw1, iw2, iw3;
	int rw1, rw2, rw3;
	int rx0_f, rx1_f, rx2_f, rx3_f;
	int ix0_f, ix1_f, ix2_f, ix3_f;
	int rx0_l, rx1_l, rx2_l, rx3_l;
	int ix0_l, ix1_l, ix2_l, ix3_l;
	int iStage, iButterfly, counter, fft_points;
	int xe_real_f, xe_imag_f, xo_real_f, xo_imag_f;
	int xe_real_l, xe_imag_l, xo_real_l, xo_imag_l;
	int real0, real1, real2;
	int imag0, imag1, imag2;
	int i;

	/* Exception handling */
	/* points must be 64, 256 or 1024 */
	if( points != 256 && points != 1024 && points != 64 ){
//		printf("#### Error: Points of input data must be 64, 256 or 1024\n");
		return;
	}

	fft_points = points >> 1;
	counter = 0;
	
	/* 1. Input separation */
	/* 
	 * [X(k) + X*(N/2-k)] / 2 = Xe(k) 
	 * j*[X(k) - X*(N/2-k)] / 2/exp(-j*2*pi*k/N) = Xo(k)
	 * X^(k) = Xe(k) + Xo(k)
	 * N = points
	 * X*(N/2-k) = X(k+N/2)
	 */

	/* X^(0)*/
	xe_real_f = ( x[0] + xHalfN ) >> 1;
	xo_imag_f = ( x[0] - xHalfN ) >> 1;
	
		x[0] =( xe_real_f >> rs[counter]);
		x[1] = (xo_imag_f >> rs[counter]);

	rx1_f = 2;									//index for k
	ix1_f = 3;
	rx2_f = points - 2;							//index for N/2-k
	ix2_f = points - 1;
	rw1 = 2;									//index for twist factor
	iw1 = 3;
	rw2 = points - 2;
	iw2 = points - 1;

	for( ; rx1_f < fft_points; ){
		xe_real_f = (x[rx1_f] + x[rx2_f]) >> 1;
		xe_imag_f = (x[ix1_f] - x[ix2_f]) >> 1;
		real0 = (x[rx1_f] - x[rx2_f]) >> 1;
		imag0 = (x[ix1_f] + x[ix2_f]) >> 1;
		xo_real_f = (int)(((long)real0 * w[iw1] - (long)imag0 * w[rw1]) >> 15);
		xo_imag_f = (int)(((long)real0 * w[rw1] + (long)imag0 * w[iw1]) >> 15);

		xe_real_l = (x[rx2_f] + x[rx1_f]) >> 1;
		xe_imag_l = (x[ix2_f] - x[ix1_f]) >> 1;
		real0 = (x[rx2_f] - x[rx1_f]) >> 1;
		imag0 = (x[ix2_f] + x[ix1_f]) >> 1;
		xo_real_l = (int)(((long)real0 * w[iw2] - (long)imag0 * w[rw2]) >> 15);
		xo_imag_l = (int)(((long)real0 * w[rw2] + (long)imag0 * w[iw2]) >> 15);
				
			x[rx1_f] = ((xe_real_f + xo_real_f) >> rs[counter]);
			x[ix1_f] = ((xe_imag_f + xo_imag_f) >> rs[counter]);
			x[rx2_f] = ((xe_real_l + xo_real_l) >> rs[counter]);
			x[ix2_f] = ((xe_imag_l + xo_imag_l) >> rs[counter]);

		rx1_f += 2;
		rx2_f -= 2;
		rw1 += 2;
		rw2 -= 2;
		ix1_f = rx1_f + 1;
		ix2_f = rx2_f + 1;
		iw1 = rw1 + 1;
		iw2 = rw2 + 1;
	}//end of for(; rx1_f < fft_points;)

	/* X^(N/4)*/
	xe_real_f = (x[rx1_f] + x[rx2_f]) >> 1;
	xe_imag_f = (x[ix1_f] - x[ix2_f]) >> 1;
	real0 = (x[rx1_f] - x[rx2_f]) >> 1;
	imag0 = (x[ix1_f] + x[ix2_f]) >> 1;
	xo_real_f = (int)(((long)real0 * w[iw1] - (long)imag0 * w[rw1]) >> 15);
	xo_imag_f = (int)(((long)real0 * w[rw1] + (long)imag0 * w[iw1]) >> 15);	
		x[rx1_f] = ((xe_real_f + xo_real_f) >> rs[counter]);
		x[ix1_f] = ((xe_imag_f + xo_imag_f) >> rs[counter]);	
	counter++;
	n2 = fft_points;	// n2 is used for butterfly counts
	ie = 1;				// ie is used for the steps for the index of twist factors 

	/* Complex FFT calculation */
	/* 2. Radix2 */
	n1 = n2;
	n2 >>= 1;
	n3 = n2 << 1;
	rw1 = 0;
	iw1 = 1;
	for(iButterfly = 0; iButterfly < n2; iButterfly++){
		for(rx0_f = (iButterfly<<1); rx0_f < 2*fft_points; rx0_f += 2*n1){
			rx2_f = rx0_f + n3;
			ix0_f = rx0_f + 1;
			ix2_f = rx2_f + 1;
			real0 = x[rx0_f] + x[rx2_f];
			imag0 = x[ix0_f] + x[ix2_f];
			real1 = x[rx0_f] - x[rx2_f];
			imag1 = x[ix0_f] - x[ix2_f];
	
				x[rx0_f] =(real0 >> rs[counter]);
				x[ix0_f] =(imag0 >> rs[counter]);
				x[rx2_f] = (int)(((long)imag1 * w[iw1] + (long)real1 * w[rw1]) >> 15 >> rs[counter]);
				x[ix2_f] = (int)(((long)imag1 * w[rw1] - (long)real1 * w[iw1]) >> 15 >> rs[counter]);
		}
		rw1 += 4;
		iw1 = rw1 + 1;
	}
	counter++;

	/* 3. Radix4 */
	/*
	 * X(4r) 	 = FFT[A(n)]	A(n) = x(n) + x(n + N/2) + x(n + N/4) + x(n + 3N/4)
	 * X(4r + 2) = FFT[B(n)]	B(n) = [x(n) + x(n + N/2) - x(n + N/4) - x(n + 3N/4)]*Wn(-2n)
	 * X(4r + 1) = FFT[C(n)]    C(n) = [x(n) - x(n + N/2) + jx(n + N/4) - jx(n + 3N/4)]*Wn(-n)
	 * X(4r + 3) = FFT[D(n)]    D(n) = [x(n) - x(n + N/2) - jx(n + N/4) + jx(n + 3N/4)]*Wn(-3n)
	 */
	fft_points >>= 1;
	for( iStage = fft_points; iStage > 1; iStage >>= 2){
		n1 = n2;
		n2 >>= 2;
		n3 = n2 << 1;
		rw1 = 0;
		iw1 = 1;
		for(iButterfly = 0; iButterfly < n2; iButterfly++){
			rw2 = rw1 + rw1;
			rw3 = rw2 + rw1; 
			iw2 = rw2 + 1;
			iw3 = rw3 + 1;
			for(rx0_f = (iButterfly<<1); rx0_f < 2*fft_points; rx0_f += 2*n1){
				rx1_f = rx0_f + n3;	
				rx2_f = rx1_f + n3;
				rx3_f = rx2_f + n3;
				ix0_f = rx0_f + 1;
				ix1_f = rx1_f + 1;
				ix2_f = rx2_f + 1;
				ix3_f = rx3_f + 1;

				/******************************************************\
				 *even points calculations for former fft_points points
				\******************************************************/
				real0 = x[rx1_f] + x[rx3_f];
				imag0 = x[ix1_f] + x[ix3_f];
				real1 = x[rx0_f] + x[rx2_f];
				imag1 = x[ix0_f] + x[ix2_f];
				real2 = x[rx0_f] - x[rx2_f];
				imag2 = x[ix0_f] - x[ix2_f];
					x[rx0_f] = ((real0 + real1) >> rs[counter]);		//A(n) real part
					x[ix0_f] = ((imag0 + imag1) >> rs[counter]);		//A(n) image part
				real1 = real1 - real0;
				imag1 = imag1 - imag0;  
				
				real0 = x[rx1_f] - x[rx3_f];
				imag0 = x[ix1_f] - x[ix3_f];
					x[rx1_f] = (int)(((long)imag1 * w[iw2] + (long)real1 * w[rw2]) >> rs[counter] >> 15);//B(n) real part
					x[ix1_f] = (int)(((long)imag1 * w[rw2] - (long)real1 * w[iw2]) >> rs[counter] >> 15);//B(n) image part
				/********************************************************\       
				 *odd points calculations for former fft_points points
				\********************************************************/
				real1 = real2 - imag0;
				imag1 = imag2 + real0;

				real2 = real2 + imag0;
				imag2 = imag2 - real0;
					x[rx2_f] = (int)(((long)real1 * w[rw1] + (long)imag1 * w[iw1]) >> rs[counter] >> 15);//C(n) real part
					x[ix2_f] = (int)(((long)imag1 * w[rw1] - (long)real1 * w[iw1]) >> rs[counter] >> 15);//C(n) image part
					x[rx3_f] = (int)(((long)real2 * w[rw3] + (long)imag2 * w[iw3]) >> rs[counter] >> 15);//D(n) real part
					x[ix3_f] = (int)(((long)imag2 * w[rw3] - (long)real2 * w[iw3]) >> rs[counter] >> 15);//D(n) image part
				/*******************************************************\
				 *even points calculations for latter fft_points points
				\*******************************************************/
				rx0_l = rx0_f + (fft_points << 1);
				rx1_l = rx1_f + (fft_points << 1);
				rx2_l = rx2_f + (fft_points << 1);
				rx3_l = rx3_f + (fft_points << 1);
				ix0_l = rx0_l + 1;
				ix1_l = rx1_l + 1;
				ix2_l = rx2_l + 1;
				ix3_l = rx3_l + 1;

				real0 = x[rx1_l] + x[rx3_l];
				imag0 = x[ix1_l] + x[ix3_l];
				real1 = x[rx0_l] + x[rx2_l];
				imag1 = x[ix0_l] + x[ix2_l];
				real2 = x[rx0_l] - x[rx2_l];
				imag2 = x[ix0_l] - x[ix2_l];
					x[rx0_l] = ((real0 + real1) >> rs[counter]);		//A(n) real part
					x[ix0_l] = ((imag0 + imag1) >> rs[counter]);		//A(n) image part
				real1 = real1 - real0;
				imag1 = imag1 - imag0;  
				
				real0 = x[rx1_l] - x[rx3_l];
				imag0 = x[ix1_l] - x[ix3_l];
					x[rx1_l] = (int)(((long)real1 * w[rw2] + (long)imag1 * w[iw2]) >> rs[counter] >> 15);//B(n) real part
					x[ix1_l] = (int)(((long)imag1 * w[rw2] - (long)real1 * w[iw2]) >> rs[counter] >> 15);//B(n) image part
				/********************************************************\       
				 *odd points calculations for latter fft_points points
				\********************************************************/
				real1 = real2 - imag0;
				imag1 = imag2 + real0;

				real2 = real2 + imag0;
				imag2 = imag2 - real0;
					x[rx2_l] = (int)(((long)real1 * w[rw1] + (long)imag1 * w[iw1]) >> rs[counter] >> 15);//C(n) real part
					x[ix2_l] = (int)(((long)imag1 * w[rw1] - (long)real1 * w[iw1]) >> rs[counter] >> 15);//C(n) image part
					x[rx3_l] = (int)(((long)real2 * w[rw3] + (long)imag2 * w[iw3]) >> rs[counter] >> 15);//D(n) real part
					x[ix3_l] = (int)(((long)imag2 * w[rw3] - (long)real2 * w[iw3]) >> rs[counter] >> 15);//D(n) image part
			}//end of inner loop
			rw1 += (ie<<3);
			iw1 = rw1 + 1;
		}//end of outer loop
		counter++;
		ie <<= 2;
	}

	/* 4. Bit-reverse */
//	fft_points <<= 1;
	//bit_reverse_cplx(x, fft_points);
	for(i=0;i<256;i++)
	    xTemp[i]=x[ reverseN[i] ];
	return;
}