fft_r4_real_256.c

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

/*
 *	DATE:			09-20-2006
 *
 *	AUTHOR:			CHEN PENG
 *
 *	VERSION:
 *		ORIGINAL	09-20-2006
 */
/**************************************************************************************\
*	函数名：fft_r4_real_256
*
*	参数  ：x[]为输入的实数数据，类型为short
*			w[]为旋转因子，类型为short
*			ls[]数组中每一个元素都表示对应一级FFT运算及输出分离时的左移位数
*			rs[]数组中每一个元素都表示对应一级FFT运算及输出分离时的右移位数
*			points为输入实数的长度
*
*	返回值：无
*
*	说明  ：1.该函数使用radix4和radix2混合的实输入数据的DIF FFT算法，主要针对256，1024
*			和64这类数据长度的实数。
*			2.输出数据占据输入数据空间，并且输出的复数点数为输入实数点数的一半，按照
*			{实部0，虚部0，实部1，虚部1，...}排列。实部和虚部各占一个short。
*			3.旋转因子为points点的旋转因子，按照{实部0，虚部0，实部1，虚部1，...}
*			排列,以Q来表示其标度为Q15。
*			4.该函数不检查ls[]和rs[]的长度，所以，调用函数时必须保证ls[]和rs[]的长度
*			不小于FFT的级数+1,最后一位有效元素的值为输出分离时的右移位数，这里ls[0]和
*			rs[0]表示radix2的移位数。
*
\**************************************************************************************/
extern int reverseN[256];
void	fft_r4_real_256( int x[], int xTemp[], int w[], int rs[], int points ){
	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);
	n2 = fft_points;	// n2 is used for butterfly counts
	ie = 1;				// ie is used for the steps for the index of twist factors 
	counter = 0;

	/* Complex FFT calculation */
	/* 1. 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] = (short)(((long)real1 * w[rw1] - (long)imag1 * w[iw1]) >> 15 >> rs[counter]);
				x[ix2_f] = (short)(((long)real1 * w[iw1] + (long)imag1 * w[rw1]) >> 15 >> rs[counter]);		
		}
		rw1 += 4;
		iw1 = rw1 + 1;
	}
	counter++;

	/* 2. 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] = (short)((((long)real1 * w[rw2]) - ((long)imag1 * w[iw2])) >> rs[counter] >> 15);//B(n) real part
					x[ix1_f] = (short)((((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] =(short)( (((long)real1 * w[rw1]) - ((long)imag1 * w[iw1])) >> rs[counter] >> 15);//C(n) real part
					x[ix2_f] = (short)((((long)imag1 * w[rw1]) + ((long)real1 * w[iw1])) >> rs[counter] >> 15);//C(n) image part
					x[rx3_f] = (short)((((long)real2 * w[rw3]) - ((long)imag2 * w[iw3])) >> rs[counter] >> 15);//D(n) real part
					x[ix3_f] = (short)((((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] = (short)((((long)real1 * w[rw2]) - ((long)imag1 * w[iw2])) >> rs[counter] >> 15);//B(n) real part
					x[ix1_l] = (short)((((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] = (short)((((long)real1 * w[rw1]) - ((long)imag1 * w[iw1])) >> rs[counter] >> 15);//C(n) real part
					x[ix2_l] = (short)((((long)imag1 * w[rw1]) + ((long)real1 * w[iw1])) >> rs[counter] >> 15);//C(n) image part
					x[rx3_l] = (short)((((long)real2 * w[rw3]) - ((long)imag2 * w[iw3])) >> rs[counter] >> 15);//D(n) real part
					x[ix3_l] = (short)((((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;
	}

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

	/* 4. Output separation */
	/* 
	 * [X^(k) + X^*(-k)] / 2 = Xe(k) 
	 * [X^(k) - x^*(-k)] / 2j = Xo(k)
	 * X(k) = Xe(k) + Xo(k) * e(-j*2*pi*k/N)
	 * N = points / 2
	 * X^*(-k) = X^*(N - k)
	 */
	
	xe_real_f = ( xTemp[0] + xTemp[0] ) >> 1;
	xe_imag_f = ( xTemp[1] + (-xTemp[1]) ) >> 1;
	xo_real_f = ( xTemp[1] - (-xTemp[1]) ) >> 1;
	xo_imag_f = ( -(xTemp[0] - xTemp[0]) ) >> 1;	
		x[0] = (short)(((long)xe_real_f + (((long)xo_real_f * w[0] - (long)xo_imag_f * w[1]) >> 15)) >> rs[counter]);//X(0) real part
		x[1] = (short)(((long)xe_imag_f + (((long)xo_real_f * w[1] + (long)xo_imag_f * w[0]) >> 15)) >> rs[counter]);//X(0) image part
	rx1_f = 2;									//index for k
	ix1_f = 3;
	rx2_f = points - 2;							//index for N-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 = ( xTemp[rx1_f] + xTemp[rx2_f] ) >> 1;		//Xe(k) real part of former points
		xe_imag_f = ( xTemp[ix1_f] + (-xTemp[ix2_f]) ) >> 1;	//Xe(k) image part of former points
		xo_real_f = ( xTemp[ix1_f] - (-xTemp[ix2_f]) ) >> 1;	//Xo(k) real part of former points
		xo_imag_f = (-( xTemp[rx1_f] - xTemp[rx2_f] )) >> 1;	//Xo(k) image part of former points

		xe_real_l = ( xTemp[rx2_f] + xTemp[rx1_f] ) >> 1;		//Xe(k) real part of latter points
		xe_imag_l = ( xTemp[ix2_f] + (-xTemp[ix1_f]) ) >> 1;	//Xe(k) image part of latter points
		xo_real_l = ( xTemp[ix2_f] - (-xTemp[ix1_f]) ) >> 1;	//Xo(k) real part of latter points	
		xo_imag_l = (-( xTemp[rx2_f] - xTemp[rx1_f] )) >> 1;	//Xo(k) image part of latter points
		
			x[rx1_f] =(short)( ((long)xe_real_f + (((long)xo_real_f * w[rw1] - (long)xo_imag_f * w[iw1]) >> 15)) >> rs[counter]);	//X(k) real part
			x[ix1_f] =(short)( ((long)xe_imag_f + (((long)xo_real_f * w[iw1] + (long)xo_imag_f * w[rw1]) >> 15)) >> rs[counter]);	//X(k) image part
			x[rx2_f] =(short)( ((long)xe_real_l + (((long)xo_real_l * w[rw2] - (long)xo_imag_l * w[iw2]) >> 15)) >> rs[counter]);	//X(N/2-k) real part
			x[ix2_f] =(short)( ((long)xe_imag_l + (((long)xo_real_l * w[iw2] + (long)xo_imag_l * w[rw2]) >> 15)) >> rs[counter]);	//X(N/2-k) image part		
		
		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;
	}
	xe_real_l = ( xTemp[rx1_f] + xTemp[rx2_f] ) >> 1;
	xe_imag_l = ( xTemp[ix1_f] + (-xTemp[ix2_f]) ) >> 1;
	xo_real_l = ( xTemp[ix1_f] - (-xTemp[ix2_f]) ) >> 1;
	xo_imag_l = ( -(xTemp[rx1_f] - xTemp[rx2_f]) ) >> 1;

		x[rx1_f] = (short)(((long)xe_real_l + (((long)xo_real_l * w[rw1] - (long)xo_imag_l * w[iw1]) >> 15)) >> rs[counter]);//X(N/4) real part
		x[ix1_f] = (short)(((long)xe_imag_l + (((long)xo_real_l * w[iw1] + (long)xo_imag_l * w[rw1]) >> 15)) >> rs[counter]);//X(N/4) image part	 
  return;
}