快速傅里叶变换(FFT)原理及源程序

《测试信号分析及处理》课程作业快速傅里叶变换一、程序设计思路快速傅里叶变换的目的是减少运算量，其用到的方法是分级进行运算。全部计算分解为M级，其中NM2log；在输入序列ix中是按码位倒序排列的，输出序列kX是按顺序排列；每级包含2N个蝶形单元，第i级有iN2个群，每个群有12i个蝶形单元；每个蝶形单元都包含乘rNW和rNW系数的运算，每个蝶形单元数据的间隔为12i，i为第i级；同一级中各个群的系数W分布规律完全相同。将输入序列ix按码位倒序排列时，用到的是倒序算法——雷德算法。自然序排列的二进制数，其下面一个数总比上面的数大1，而倒序二进制数的下面一个数是上面一个数在最高位加1并由高位向低位仅为而得到的。若已知某数的倒序数是J，求下一个倒序数，应先判断J的最高位是否为0，与2Nk进行比较即可得到结果。如果Jk，说明最高位为0，应把其变成1，即2NJ，这样就得到倒序数了。如果Jk，即J的最高位为1，将最高位化为0，即2NJ，再判断次高位；与4Nk进行比较，若为0，将其变位1，即4NJ，即得到倒序数，如果次高位为1，将其化为0，再判断下一位……即从高位到低位依次判断其是否为1，为1将其变位0，若这一位为0，将其变位1，即可得到倒序数。若倒序数小于顺序数，进行换位，否则不变，防治重复交换，变回原数。注：因为0的倒序数为0，所以可从1开始进行求解。二、程序设计框图（1）倒序算法——雷德算法流程图(2)FFT算法流程三、FFT源程序voidfft(x,n)intn;doublex[];{inti,j,k,l,m,n1,n2;doublec,c1,e,s,s1,t,tr;for(j=1,i=1;in/2;i++){m=i;j=2*j;if(j==n)break;}//得到流程图的共几级n1=n-1;for(j=0,i=0;in1;i++){if(ij)//如果ij,即进行变址{tr=x[j];x[j]=x[i];x[i]=tr;}k=n/2;//求j的下一个倒位序while(k(j+1))//如果k(j+1),表示j的最高位为1{j=j-k;//把最高位变成0k=k/2;//k/2，比较次高位，依次类推，逐个比较，直到某个位为0}j=j+k;//把0改为1}for(i=0;in;i+=2){tr=x[i];x[i]=tr+x[i+1];x[i+1]=tr-x[i+1];}n2=1;for(l=1;l=m;l++)//控制蝶形结级数{n4=n2;n2=2*n4;n1=2*n2;e=6.28318530718/n1;for(i=0;in;i+=n1)//控制同一蝶形结运算，即计算系数相同蝶形结{tr=x[i];x[i]=tr+x[i+n2];x[i+n2]=tr-x[i+n2];x[i+n2+n4]=-x[i+n2+n4];a=e;for(j=2;j=(n4-1);j++)//控制计算不同种蝶形结，即计算系数不同的蝶形结{i1=i+j;i2=i-j+n2;i3=i+j+n2;i4=i-j+n1;cc=cos(a);ss=sin(a);a=a+e;t1=cc*x[i3]+ss*x[i4];t2=ss*x[i3]-cc*x[i4];x[i4]=x[i2]-t2;x[i3]=-x[i2]-t2;x[i2]=x[i1]-t1;x[i1]=x[i1]+t1;}}}}四、计算实例及运行结果设输入序列)(ix为)1,...,2,1,0(,200sin)(nitiix其离散傅里叶变换为)1,...,2,1,0(,)()(10nkWixkXNiikN这里njeW2。选n=512,计算离散傅里叶变换)(kX。所用软件为Turboc2.0，操作界面如图1所示图1Turboc2.0操作界面程序运行结束后的界面如图2所示图2程序运行后的界面例子的具体程序如下：#includemath.h#includestdio.h#includestdlib.h#definepi3.14159265359voidfft(x,n)intn;doublex[];{inti,j,k,l,i1,i2,i3,i4,n4,m,n1,n2;doublea,e,cc,ss,tr,t1,t2;for(j=1,i=1;in/2;i++){m=i;j=2*j;if(j==n)break;}n1=n-1;for(j=0,i=0;in1;i++){if(ij){tr=x[j];x[j]=x[i];x[i]=tr;}k=n/2;while(k(j+1)){j=j-k;k=k/2;}j=j+k;}for(i=0;in;i+=2){tr=x[i];x[i]=tr+x[i+1];x[i+1]=tr-x[i+1];}n2=1;for(l=1;l=m;l++){n4=n2;n2=2*n4;n1=2*n2;e=6.28318530718/n1;for(i=0;in;i+=n1){tr=x[i];x[i]=tr+x[i+n2];x[i+n2]=tr-x[i+n2];x[i+n2+n4]=-x[i+n2+n4];a=e;for(j=2;j=(n4-1);j++){i1=i+j;i2=i-j+n2;i3=i+j+n2;i4=i-j+n1;cc=cos(a);ss=sin(a);a=a+e;t1=cc*x[i3]+ss*x[i4];t2=ss*x[i3]-cc*x[i4];x[i4]=x[i2]-t2;x[i3]=-x[i2]-t2;x[i2]=x[i1]-t1;x[i1]=x[i1]+t1;}}}}main(){FILE*p;inti,j,n;doubledt=0.001;doublex[512];p=fopen(d:\123.c,w);n=512;for(i=0;in;i++){x[i]=sin(200*pi*i*dt);}for(i=0;in;i++){fprintf(p,%10.7f,x[i]);fprintf(p,\n);printf(%10.7f,x[i]);printf(\n);}fft(x,n);fprintf(p,\nDISCRETEFOURIERTRANSFORM\n);printf(\nDISCRETEFOURIERTRANSFORM\n);fprintf(p,%10.7f,x[0]);printf(%10.7f,x[0]);fprintf(p,%10.7f+J%10.7f\n,x[1],x[n-1]);printf(%10.7f+J%10.7f\n,x[1],x[n-1]);for(i=2;in/2;i+=2){fprintf(p,%10.7f+J%10.7f,x[i],x[n-i]);fprintf(p,%10.7f+J%10.7f,x[i+1],x[n-i-1]);fprintf(p,\n);printf(%10.7f+J%10.7f,x[i],x[n-i]);printf(%10.7f+J%10.7f,x[i+1],x[n-i-1]);printf(\n);}fprintf(p,%10.7f,x[n/2]);printf(%10.7f,x[n/2]);fprintf(p,%10.7f+J%10.7f\n,x[n/2-1],-x[n/2+1]);for(i=2;in/2;i+=2){fprintf(p,%10.7f+J%10.7f,x[n/2-i],-x[n/2+i]);fprintf(p,%10.7f+J%10.7f,x[n/2-i-1],-x[n/2+i+1]);fprintf(p,\n);printf(%10.7f+J%10.7f,x[n/2-i],-x[n/2+i]);printf(%10.7f+J%10.7f,x[n/2-i-1],-x[n/2+i+1]);printf(\n);}}将程序运行后所得数据绘制成曲线图（其中FFT变换的数据要先取绝对值后再画图）如下由上图可知，变换后的图开在频率100Hz处出现一个峰值，这与理论上的结果一致。