func.h

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

void bit_reverse_cplx(long x[], int num){
	int iBefore, iBehind, tempData, temph, templ, temp, i;
	int hBitRef, lBitRef, hBit, lBit;
	char hPos, lPos, bitMov;
	short nBits;
	nBits = 0;
	i = num;
	/* acquire bit numbers occupied by the size of x */
	while(i > 1){
		i >>= 1;
		nBits++;
	}
	hBitRef = 0x00000001 << (--nBits);
	lBitRef = 0x00000001;
	/*
	 *for each index, search for its corresponding bit-reversed index, then
	 *switch the indexed datas.
	 */
	for( iBefore = 1; iBefore < num; iBefore++ ){
		hPos = nBits;
		lPos = 0;
		hBit = hBitRef;
		lBit = lBitRef;
		bitMov = nBits;
		temp = 0;
		/* bit-reverse */
		while( hPos > lPos ){
			temph = iBefore & hBit;
			templ = iBefore & lBit;
			temph >>= bitMov;
			templ <<= bitMov;
			temp += temph + templ;
			hPos--;
			lPos++;
			bitMov -= 2;
			hBit >>= 1;
			lBit <<= 1;
		}
		if( hPos == lPos ){
			temph = iBefore & hBit;
			temp += temph;
		}
		/* data exchange */
		iBehind = temp;  
		if( iBehind < iBefore )
			continue;
		tempData = x[iBefore];
		x[iBefore] = x[iBehind];
		x[iBehind] = tempData;
	}
	return;	
}




void	fft_r4_real_512( short x[], short w[], short ls[], short rs[], short points ){
	short r_or_l;			//right = '1' , left = '0'
	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, fft_points;
	long xe_real_f, xe_imag_f, xo_real_f, xo_imag_f;
	long xe_real_l, xe_imag_l, xo_real_l, xo_imag_l;
	long real0, real1, real2;
	long imag0, imag1, imag2;
	short max_x,min_x,i;


	/* 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;
	}
	
	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;
	r_or_l = 1;			//Default: right shift


    max_x=abs(x[0]);
	min_x=abs(x[0]);
    for(i=1;i<points;i++)
	{
		if(abs(x[i])>max_x)
				max_x=abs(x[i]);
		if(abs(x[i])<min_x)
				min_x=abs(x[i]);
	}

	/* 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 = fft_points; iStage > 1; iStage >>= 2 ){
		n1 = n2;
		n2 >>= 2;
		n3 = n2 << 1;
		rw1 = 0;
		iw1 = 1;
		if(ls[counter] > rs[counter])
			r_or_l = 0;
		for( iButterfly = 0; iButterfly < n2; iButterfly++ ){
			rw2 = rw1 + rw1;
			rw3 = rw2 + rw1;
			iw2 = rw2 + 1;
			iw3 = rw3 + 1;
			for( rx0 = ( iButterfly<<1 ); rx0 < 2*fft_points; rx0 += 2*n1 ){
				rx1 = rx0 + n3;									//real parts index
				rx2 = rx1 + 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];
				if(r_or_l){
					x[rx0] = (short)((real0 + real1) >> rs[counter]);		//A(n) real part
					x[ix0] = (short)((imag0 + imag1) >> rs[counter]);		//A(n) image part
				}
				else {
					x[rx0] = (short)((real0 + real1) << ls[counter]);		//A(n) real part
					x[ix0] = (short)((imag0 + imag1) << ls[counter]);		//A(n) image part
				}
				real1 = real1 - real0;
				imag1 = imag1 - imag0;  
				
				real0 = x[rx1] - x[rx3];
				imag0 = x[ix1] - x[ix3];
				if(r_or_l){
					x[rx1] = ((real1 * w[rw2]) - (imag1 * w[iw2])) >> rs[counter] >> 15;//B(n) real part
					x[ix1] = ((imag1 * w[rw2]) + (real1 * w[iw2])) >> rs[counter] >> 15;//B(n) image part
				}
				else {
					x[rx1] = ((real1 * w[rw2]) - (imag1 * w[iw2])) << ls[counter] >> 15;//B(n) real part
					x[ix1] = ((imag1 * w[rw2]) + (real1 * w[iw2])) << ls[counter] >> 15;//B(n) image part
				}
				/*************************\       
				 *odd points calculations 
				\*************************/
				real1 = real2 + imag0;
				imag1 = imag2 + (-real0);

				real2 = real2 - imag0;
				imag2 = imag2 - (-real0);
				if(r_or_l){
					x[rx2] = ((real1 * w[rw1]) - (imag1 * w[iw1])) >> rs[counter] >> 15;//C(n) real part
					x[ix2] = ((imag1 * w[rw1]) + (real1 * w[iw1])) >> rs[counter] >> 15;//C(n) image part
					x[rx3] = ((real2 * w[rw3]) - (imag2 * w[iw3])) >> rs[counter] >> 15;//D(n) real part
					x[ix3] = ((imag2 * w[rw3]) + (real2 * w[iw3])) >> rs[counter] >> 15;//D(n) image part
				}
				else {
					x[rx2] = ((real1 * w[rw1]) - (imag1 * w[iw1])) << ls[counter] >> 15;//C(n) real part
					x[ix2] = ((imag1 * w[rw1]) + (real1 * w[iw1])) << ls[counter] >> 15;//C(n) image part
					x[rx3] = ((real2 * w[rw3]) - (imag2 * w[iw3])) << ls[counter] >> 15;//D(n) real part
					x[ix3] = ((imag2 * w[rw3]) + (real2 * w[iw3])) << ls[counter] >> 15;//D(n) image part
				}
			}//end of inner for()
			rw1 += 4*ie;
			iw1 = rw1 + 1;

		}//end of outer for()
		counter++;
		ie <<= 2;
		r_or_l = 1;
	    
		max_x=abs(x[0]);
		min_x=abs(x[0]);
        for(i=1;i<points;i++)
		{
		    if(abs(x[i])>max_x)
				max_x=abs(x[i]);
			if(abs(x[i])<min_x)
				min_x=abs(x[i]);
		}

	}//end of outest for()

	/* Bit-reverse */
	bit_reverse_cplx((long*)x, fft_points);


	/* 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 = ( x[0] + x[0] ) >> 1;
	xe_imag_f = ( x[1] + (-x[1]) ) >> 1;
	xo_real_f = ( x[1] - (-x[1]) ) >> 1;
	xo_imag_f = ( -(x[0] - x[0]) ) >> 1;
	if(ls[counter] > rs[counter])
		r_or_l = 0;
	if(r_or_l){
		x[0] = (xe_real_f + ((xo_real_f * w[0] - xo_imag_f * w[1]) >> 15)) >> rs[counter];//X(0) real part
		x[1] = (xe_imag_f + ((xo_real_f * w[1] + xo_imag_f * w[0]) >> 15)) >> rs[counter];//X(0) image part
	}
	else {
		x[0] = (xe_real_f + ((xo_real_f * w[0] - xo_imag_f * w[1]) >> 15)) << ls[counter];//X(0) real part
		x[1] = (xe_imag_f + ((xo_real_f * w[1] + xo_imag_f * w[0]) >> 15)) << ls[counter];//X(0) image part	
	}
	rx1 = 2;									//index for k
	ix1 = 3;
	rx2 = points - 2;							//index for N-k
	ix2 = points - 1;
	rw1 = 2;									//index for twist factor
	iw1 = 3;
	rw2 = points - 2;
	iw2 = points - 1;

	for( ; rx1 < fft_points; ){
		xe_real_f = ( x[rx1] + x[rx2] ) >> 1;		//Xe(k) real part of former points
		xe_imag_f = ( x[ix1] + (-x[ix2]) ) >> 1;	//Xe(k) image part of former points
		xo_real_f = ( x[ix1] - (-x[ix2]) ) >> 1;	//Xo(k) real part of former points
		xo_imag_f = (-( x[rx1] - x[rx2] )) >> 1;	//Xo(k) image part of former points

		xe_real_l = ( x[rx2] + x[rx1] ) >> 1;		//Xe(k) real part of latter points
		xe_imag_l = ( x[ix2] + (-x[ix1]) ) >> 1;	//Xe(k) image part of latter points
		xo_real_l = ( x[ix2] - (-x[ix1]) ) >> 1;	//Xo(k) real part of latter points	
		xo_imag_l = (-( x[rx2] - x[rx1] )) >> 1;	//Xo(k) image part of latter points
		if(r_or_l){
			x[rx1] = (xe_real_f + ((xo_real_f * w[rw1] - xo_imag_f * w[iw1]) >> 15)) >> rs[counter];	//X(k) real part
			x[ix1] = (xe_imag_f + ((xo_real_f * w[iw1] + xo_imag_f * w[rw1]) >> 15)) >> rs[counter];	//X(k) image part
			x[rx2] = (xe_real_l + ((xo_real_l * w[rw2] - xo_imag_l * w[iw2]) >> 15)) >> rs[counter];	//X(N/2-k) real part
			x[ix2] = (xe_imag_l + ((xo_real_l * w[iw2] + xo_imag_l * w[rw2]) >> 15)) >> rs[counter];	//X(N/2-k) image part		
		}
		else {
			x[rx1] = (xe_real_f + ((xo_real_f * w[rw1] - xo_imag_f * w[iw1]) >> 15)) << ls[counter];	//X(k) real part
			x[ix1] = (xe_imag_f + ((xo_real_f * w[iw1] + xo_imag_f * w[rw1]) >> 15)) << ls[counter];	//X(k) image part
			x[rx2] = (xe_real_l + ((xo_real_l * w[rw2] - xo_imag_l * w[iw2]) >> 15)) << ls[counter];	//X(N/2-k) real part
			x[ix2] = (xe_imag_l + ((xo_real_l * w[iw2] + xo_imag_l * w[rw2]) >> 15)) << ls[counter];	//X(N/2-k) image part		
		}
		rx1 += 2;
		rx2 -= 2;
		rw1 += 2;
		rw2 -= 2;
		ix1 = rx1 + 1;
		ix2 = rx2 + 1;
		iw1 = rw1 + 1;
		iw2 = rw2 + 1;
	}//end for
	xe_real_l = ( x[rx1] + x[rx2] ) >> 1;
	xe_imag_l = ( x[ix1] + (-x[ix2]) ) >> 1;
	xo_real_l = ( x[ix1] - (-x[ix2]) ) >> 1;
	xo_imag_l = ( -(x[rx1] - x[rx2]) ) >> 1;
	if(r_or_l){
		x[rx1] = (xe_real_l + ((xo_real_l * w[rw1] - xo_imag_l * w[iw1]) >> 15)) >> rs[counter];//X(N/4) real part
		x[ix1] = (xe_imag_l + ((xo_real_l * w[iw1] + xo_imag_l * w[rw1]) >> 15)) >> rs[counter];//X(N/4) image part
	}
	else {
		x[rx1] = (xe_real_l + ((xo_real_l * w[rw1] - xo_imag_l * w[iw1]) >> 15)) << ls[counter];//X(N/4) real part
		x[ix1] = (xe_imag_l + ((xo_real_l * w[iw1] + xo_imag_l * w[rw1]) >> 15)) << ls[counter];//X(N/4) image part
	}

    max_x=abs(x[0]);
	min_x=abs(x[0]);
    for(i=1;i<points;i++)
	{
		if(abs(x[i])>max_x)
			max_x=abs(x[i]);
		if(abs(x[i])<min_x)
			min_x=abs(x[i]);
	}

	return;
	
}


/*FUNCTION:To get the maximum from the x[];
 *PARAMETER:The returned value is the maximum of x[].
            the  whole parameter 'point' is the point corresponding to the maximum;*/
short int max(short int *x,short int length)
{
	short int max1;
	short int i;
    max1=x[0];
    for (i=1;i<length;i++)
        if (max1<x[i])
           max1=x[i];
    for (i=0;i<length;i++)
        if (max1==x[i])
           {
              point=i;
              break;
           }
    return(max1);
}

unsigned int max_float(unsigned int *x,short int length)
    {
	///int max1;
	unsigned int max1;
	short int i;
    max1=x[0];
    for (i=1;i<length;i++)
        if (max1<x[i])
           max1=x[i];
    for (i=0;i<length;i++)
        if (max1==x[i])
           {
              point=i;
              break;
           }
    return(max1);
    }

/*FUNCTION:To get the M frequency points corresponding to the maximum peaks of A. */
/*PARAMETER:i0[M]:M frequency points corresponding to the maximum peaks.*/

///void maximum(int *A,int *i0)
void maximum(unsigned int *A,short int *i0,int len)
{   short int n,i,i_left,i_right;
    unsigned int mmax;
    for(n=0;n<len;n++)
    {
		mmax=max_float(A,160);
		i=point;//point is a whole parameter which is changed during the call of the function max;
        i0[n]=i;
        while(i>0)
        {
           if (A[i]<A[i-1])
              break;
           i--;
        }
        i_left=i;
        i=i0[n];
        while (i<160)   
        {  
            if (A[i]<A[i+1])
              break;
            i++;
        }
        i_right=i;
        for (i=i_left;i<=i_right;i++)
            A[i]=0;
	}