/**********************fft programe*********************/
/*FFT-基2DIT-FFT算法*/
/*xr：=信号序列实部，xi：=信号序列虚部，N=2^M*/
/*计算结果X(k)的实部和虚部分别存储在数组xr和xi中*/
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
const float pp=3.14159;

void FFT(double xr[],double xi[],int N,int M)
{   
	int L,B,J,P,k,i;               /*L级数，B=2^L-1为第L级B个不同旋转因子,P=j*2^M-L*/
	float pi;
	double ps;
	double rx_kb,ix_kb;       /*分别代表X(K+B)的实部和虚部*/
	double rcf[64],icf[64];   /*rcf存储旋转因子实部，icf存储旋转因子虚部*/
	                                  /*旋转因子数组长根据需要调整*/
	reverse_order(xr,xi,N)	;   /*调用倒序子程序*/
		
     /*计算各级蝶形*/
	for(L=1;L<=M;L++)
	{   
		pi=3.14159;
		B=(int)(pow(2,(L-1))+0.5);
		for(J=0;J<=B-1;J++)     /*第L级B个不同的旋转因子*/
		{
			P=J*((int)(pow(2,(M-L))+0.5));   /*2^M-L个群体乘以Ｊ，共Ｎ/2个蝶形*/
			ps=(2*pi/N)*P;        /*计算旋转因子*/
               rcf[P]=cos(ps);
               icf[P]=-sin(ps);
			
			for(k=J;k<=N-1;k+=(int)(pow(2,L)+0.5))  /*2^L为蝶距*/
			{
				rx_kb=xr[k+B]*rcf[P]-xi[k+B]*icf[P];
				ix_kb=xi[k+B]*rcf[P]+xr[k+B]*icf[P];
				xr[k+B]=xr[k]-rx_kb;
				xi[k+B]=xi[k]-ix_kb;
				xr[k]=xr[k]+rx_kb;
				xi[k]=xi[k]+ix_kb;

				}
			}
		}
	} 
	/*倒序子程序*/
	 reverse_order(double xr[],double xi[],int N)
	{
		int LH,N1,I,J,K;
		double T;
		LH=N/2;
		J=LH;
		N1=N-2;
		for(I=1;I<=N1;I++)
		{
			if(I<J)
			{
				T=xr[I]; xr[I]=xr[J]; xr[J]=T;
				T=xi[I]; xi[I]=xi[J]; xi[J]=T;
				}
			K=LH;
               while(J>=K)
               {
				J=J-K;
				K=(int)(K/2+0.5);
				}
                  J=J+K;				
			}
		
		}
		
	/**********************fft programe*********************/
	/********************main programe**********************/
    void main()    /*以3cos(2pi*50t+pi/2)为例*/
		{   
			int m,n;
			double real[32],imag[32];
			double result[32];
			for(m=0;m<16;m++)
			{
				real[m]=3*cos(100*pp+pp/2+m*pp/8);
				imag[m]=0;
				}
				
				FFT(real,imag,32,5);
				clrscr();
			for(m=0;m<16;m++)
              {
				/*   printf("%8f",real[m]);
                         printf("     %f\n",imag[m]);*/
                  result[m]=sqrt(pow(real[m],2)+pow(imag[m],2));
                  printf("%8f\n",result[m]); 
                 }
			
			}
	/********************main programe**********************/