fft_r4_real_512.c

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

/*
 *	DATE:09/15/2006
 *
 *	AUTHOR:CHEN PENG
 *
 *	VERSION:ORIGINAL	09/15/2006
 */
/*********************************************************************************\
*	函数名：fft_r4_real_512( short x[], short w[], short rs[], short points )
*	参数  ：x[]为输入的实数数据，类型为short
*			w[]为旋转因子，类型为short
*			rs[]数组中每一个元素都表示对应一级FFT运算的标度，即右移位数
*			points为输入实数的长度
*	说明  ：1.该函数使用radix4的实输入数据的FFT算法，主要针对128，512和1024这类数据
*			长度的实数，进行FFT运算。
*			2.输出数据占据输入数据空间，并且只输出输入数据一半点数的复数值，按照
*			{实部0，虚部0，实部1，虚部1，...}排列。实部和虚部各占一个short。
*			3.旋转因子为points点的旋转因子，按照{实部0，虚部0，实部1，虚部1，...}
*			排列。
*			4.该函数不检查rs[]的长度，所以，调用函数时必须保证rs[]的长度不小于FFT
*			的级数。
*
\*********************************************************************************/


void	fft_r4_real_512( short x[], short w[], short rs[], short points ){
	int n1, n2, n3, ie;
	int iw1, iw2, iw3;
	int rw1, rw2, rw3;
	int rx0, rx1, rx2, rx3;
	int ix0, ix1, ix2, ix3;
	int iStage, iButterfly, counter;
	int xe_real, xe_imag, xo_real, xo_imag;
	long real0, real1, real2;
	long imag0, imag1, imag2;
	/* Exception handling */
	/* points must be 128, 512 or 2048 */
	if( points != 512 || points != 2048 || points != 128 ){
		printf("#### Error: Points of input data must be 128, 512 or 2048\n");
		return;
	}
	
	n2 = (points >> 1);	// 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 */
	/*
	 * 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)
	 */
	for( iStage = (points >> 1); iStage > 1; iStage >>= 2 ){
		n1 = n2;
		n2 >>= 2;
		n3 = n2 << 1;
		rw1 = 0;
		iw1 = 1;
		for( iButterfly = 0; iButterfly < n2; iButterfly++ ){
			rw2 = (rw1 + rw1) << 1;
			rw3 = (rw2 + rw1) << 1;
			iw2 = rw2 + 1;
			iw3 = rw3 + 1;
			for( rx0 = ( iButterfly<<1 ); ix0 < points; ix0 += n1 ){
				rx1 = rx0 + n3;									//real parts index
				rx2 = rx2 + n3;
				rx3 = rx2 + n3;
				ix0 = rx0 + 1;									//image parts index
				ix1 = rx1 + 1;
				ix2 = rx2 + 1;
				ix3 = rx3 + 1;
				/*************************\
				 *even points calculations 
				\*************************/
				real0 = x[rx1] + x[rx3];
				imag0 = x[ix1] + x[ix3];
				real1 = x[rx0] + x[rx2];
				imag1 = x[ix0] + x[ix2];
				real2 = x[rx0] - x[rx2];
				imag2 = x[ix0] - x[ix2];

				x[rx0] = (real0 + real1) >> rs[counter];		//A(n) real part
				x[ix0] = (imag0 + imag1) >> rs[counter];		//A(n) image part

				real1 = real1 - real0;
				imag1 = imag1 - imag0;  
				
				real0 = x[rx1] - x[rx3];
				imag0 = x[ix1] - x[ix3];

				x[rx1] = ((real1 * w[rw2]) - (imag1 * w[iw2])) >> rs[counter];//B(n) real part
				x[ix1] = ((imag1 * w[rw2]) + (real1 * w[iw2])) >> rs[counter];//B(n) image part
				/*************************\       
				 *odd points calculations 
				\*************************/
				real1 = real2 + imag0;
				imag1 = imag2 + (-real0);

				real2 = real2 - imag0;
				imag2 = imag2 - (-real0);

				x[rx2] = ((real1 * w[rw1]) - (imag1 * w[iw1])) >> rs[counter];//C(n) real part
				x[ix2] = ((imag1 * w[rw1]) + (real1 * w[iw1])) >> rs[counter];//C(n) image part

				x[rx3] = ((real2 * w[rw3]) - (imag2 * w[iw3])) >> rs[counter];//D(n) real part
				x[ix3] = ((imag2 * w[rw3]) + (real2 * w[iw3])) >> rs[counter];//D(n) image part

			}//end of inner for()
			rw1 += ie;
			rw1 <<= 1;
			iw1 = rw1 + 1;

		}//end of outer for()
		counter++;
		ie <<= 2;
	
	}//end of outest for()

	/*
	 * Bit-reverse	
	 */

	//Add the code here

	/* 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)
	 */
	rx1 = 0;									//index for k
	ix1 = 1;
	rx2 = points - 2;							//index for N-k
	ix2 = points - 1;
	rw1 = 0;									//index for twist factor
	iw1 = 1;
	for( ; rx1 < points; ){
		xe_real = ( x[rx1] + x[rx2] ) >> 1;		//Xe(k) real part
		xe_imag = ( x[ix1] + (-x[ix2]) ) >> 1;	//Xe(k) image part
		xo_real = ( x[ix1] - (-x[ix2]) ) >> 1;	//Xo(k) real part	
		xo_imag = (-( x[rx1] - x[rx2] )) >> 1;	//Xo(k) image part
		x[rx1] = xe_real + (xo_real * w[rw1] - xo_imag * w[iw1]);	//X(k) real part
		x[ix1] = xe_imag + (xo_real * w[iw1] + xo_imag * w[rw1]);	//X(k) image part
		rx1 += 2;
		rx2 -= 2;
		ix1 = rx1 + 1;
		ix2 = rx2 + 1;
		rw1++;
		iw1 = rw1 + 1;
	}
	return;
	
}