FFT算法

/****************************************************************************//**//*快速傅里叶变换/快速傅里叶逆变换测试*//**//*2014年04月20日*//**//****************************************************************************/#includestdio.h//C语言标准输入输出函数库#includemath.h//C数学函数库#includemathlib.h//DSP数学函数库#includedsplib.h//DSP函数库/****************************************************************************//**//*宏定义*//**//****************************************************************************///软件断点#defineSW_BREAKPOINTasm(SWBP0);//快速傅里叶变换//π及浮点数极小值#definePI3.14159#defineF_TOL(1e-06)/****************************************************************************//**//*全局变量*//**//****************************************************************************///快速傅里叶变换测试//测试快速傅里叶变换点数//注意:TIDSP库最大支持一次性计算128K个点的FFT#defineTn1024//采样频率#defineFs1000.0//信号floatInput[2*Tn+4];//FFT输入信号#pragmaDATA_ALIGN(CFFT_In,8);floatCFFT_In[2*Tn+4];//FFT输入信号副本floatCFFT_InOrig[2*Tn+4];//FFT输出#pragmaDATA_ALIGN(CFFT_Out,8);floatCFFT_Out[2*Tn+4];//IFFT输出#pragmaDATA_ALIGN(CFFT_InvOut,8);floatCFFT_InvOut[2*Tn+4];//中间运算临时变量floatCTemp[2*Tn+4];//存储旋转因子floatCw[2*Tn];//模floatCmo[Tn+2];//二进制位翻转#pragmaDATA_ALIGN(brev,8);unsignedcharbrev[64]={0x0,0x20,0x10,0x30,0x8,0x28,0x18,0x38,0x4,0x24,0x14,0x34,0xc,0x2c,0x1c,0x3c,0x2,0x22,0x12,0x32,0xa,0x2a,0x1a,0x3a,0x6,0x26,0x16,0x36,0xe,0x2e,0x1e,0x3e,0x1,0x21,0x11,0x31,0x9,0x29,0x19,0x39,0x5,0x25,0x15,0x35,0xd,0x2d,0x1d,0x3d,0x3,0x23,0x13,0x33,0xb,0x2b,0x1b,0x3b,0x7,0x27,0x17,0x37,0xf,0x2f,0x1f,0x3f};/****************************************************************************//**//*函数声明*//**//****************************************************************************///产生旋转因子voidtw_gen(float*w,intn);//FFT测试voidFFTTest();/****************************************************************************//**//*主函数*//**//****************************************************************************/intmain(void){//FFT测试FFTTest();//断点SW_BREAKPOINT;}/****************************************************************************//**//*快速傅里叶变换测试*//**//****************************************************************************///产生旋转因子voidtw_gen(float*w,intn){inti,j,k;doublex_t,y_t,theta1,theta2,theta3;for(j=1,k=0;j=n2;j=j2){for(i=0;in2;i+=j){theta1=2*PI*i/n;x_t=cos(theta1);y_t=sin(theta1);w[k]=(float)x_t;w[k+1]=(float)y_t;theta2=4*PI*i/n;x_t=cos(theta2);y_t=sin(theta2);w[k+2]=(float)x_t;w[k+3]=(float)y_t;theta3=6*PI*i/n;x_t=cos(theta3);y_t=sin(theta3);w[k+4]=(float)x_t;w[k+5]=(float)y_t;k+=6;}}}//快速傅里叶变换voidFFTTest(void){//产生待测试信号unsignedinti;for(i=0;iTn;i++)Input[i]=5*sin(2*PI*150*(i/Fs))+15*sin(2*PI*350*(i/Fs));//确定快速傅里叶变换基unsignedcharrad;if(Tn==16||Tn==64||Tn==256||Tn==1024||Tn==4096||Tn==16384||Tn==65536)rad=4;elseif(Tn==8||Tn==32||Tn==128||Tn==512||Tn==2048||Tn==8192||Tn==32768)rad=2;else{printf(不支持计算%d点快速傅里叶变换！\n,Tn);return;}//复数FFTfor(i=0;i2*Tn;i++)CFFT_In[i]=0.0;for(i=0;iTn;i++){CFFT_In[2*i]=Input[i];//实部CFFT_In[2*i+1]=0;//虚部为0}//保留一份输入信号副本memcpy(CFFT_InOrig,CFFT_In,2*Tn*sizeof(float));//产生旋转因子tw_gen(Cw,Tn);//FFT计算DSPF_sp_fftSPxSP(Tn,CFFT_In,Cw,CFFT_Out,brev,rad,0,Tn);//计算振幅for(i=0;iTn;i++)Cmo[i]=0.0;for(i=0;iTn+2;i++){Cmo[i]=sqrtsp(CFFT_Out[2*i]*CFFT_Out[2*i]+CFFT_Out[2*i+1]*CFFT_Out[2*i+1]);Cmo[i]=Cmo[i]*2/Tn;}//保留一份FFT结果副本memcpy(CTemp,CFFT_Out,2*Tn*sizeof(float));//IFFT计算DSPF_sp_ifftSPxSP(Tn,CFFT_Out,Cw,CFFT_InvOut,brev,rad,0,Tn);//恢复FFT结果memcpy(CFFT_Out,CTemp,2*Tn*sizeof(float));printf(\n复数FFT测试结果:);unsignedcharFlag;for(i=0;iTn;i++)if(abs(CFFT_InOrig[i]-CFFT_InvOut[i])F_TOL)Flag=1;if(Flag==1)printf(失败！\n);elseprintf(成功！\n);}