(完整word版)快速傅里叶变换(FFT)原理及源程序.pdf

《(完整word版)快速傅里叶变换(FFT)原理及源程序.pdf》由会员分享，可在线阅读，更多相关《(完整word版)快速傅里叶变换(FFT)原理及源程序.pdf（9页珍藏版）》请在得力文库 - 分享文档赚钱的网站上搜索。

1、测试信号分析及处理课程作业快速傅里叶变换一、程序设计思路快速傅里叶变换的目的是减少运算量，其用到的方法是分级进行运算。全部计算分解为 M 级，其中NM2log；在输入序列ix中是按码位倒序排列的，输出序列kX是按顺序排列；每级包含2N个蝶形单元，第 i 级有iN2个群，每个群有12i个蝶形单元；每个蝶形单元都包含乘rNW 和rNW 系数的运算，每个蝶形单元数据的间隔为12i，i 为第 i 级；同一级中各个群的系数W分布规律完全相同。将输入序列ix按码位倒序排列时，用到的是倒序算法雷德算法。自然序排列的二进制数，其下面一个数总比上面的数大1，而倒序二进制数的下面一个数是上面一个数在最高位加1 并

2、由高位向低位仅为而得到的。若已知某数的倒序数是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 开始进行求解

3、。二、程序设计框图（1）倒序算法雷德算法流程图(2)FFT 算法流程三、FFT 源程序void fft(x,n)int n;double x;int i,j,k,l,m,n1,n2;double c,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=xj;xj=xi;xi=tr;k=n/2;/求 j 的下一个倒位序while(k(j+1)/如果 k(j+1),表示 j 的最高位为 1 j=j-k;/把最高位变成 0

4、k=k/2;/k/2，比较次高位，依次类推，逐个比较，直到某个位为0 j=j+k;/把 0 改为 1 for(i=0;in;i+=2)tr=xi;xi=tr+xi+1;xi+1=tr-xi+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=xi;xi=tr+xi+n2;xi+n2=tr-xi+n2;xi+n2+n4=-xi+n2+n4;a=e;for(j=2;j=(n4-1);j+)/控制计算不同种蝶形结，即计算系数不同的

5、蝶形结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*xi3+ss*xi4;t2=ss*xi3-cc*xi4;xi4=xi2-t2;xi3=-xi2-t2;xi2=xi1-t1;xi1=xi1+t1;四、计算实例及运行结果设输入序列)(ix为)1,.,2,1,0(,200sin)(nitiix其离散傅里叶变换为)1,.,2,1,0(,)()(10nkWixkXNiikN这里njeW2。选 n=512,计算离散傅里叶变换)(kX。所用软件为 Turbo c 2.0，操作界面如图1 所示图 1 Turbo c

6、 2.0 操作界面程序运行结束后的界面如图2 所示图 2 程序运行后的界面例子的具体程序如下：#include#include#include#define pi 3.14159265359 void fft(x,n)int n;double x;int i,j,k,l,i1,i2,i3,i4,n4,m,n1,n2;double a,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=xj;xj=xi;xi=tr;k=n/2;while(k(j+1)j=j-

7、k;k=k/2;j=j+k;for(i=0;in;i+=2)tr=xi;xi=tr+xi+1;xi+1=tr-xi+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=xi;xi=tr+xi+n2;xi+n2=tr-xi+n2;xi+n2+n4=-xi+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*xi3+ss*xi4;t2=ss*

8、xi3-cc*xi4;xi4=xi2-t2;xi3=-xi2-t2;xi2=xi1-t1;xi1=xi1+t1;main()FILE*p;int i,j,n;double dt=0.001;double x512;p=fopen(d:123.c,w);n=512;for(i=0;in;i+)xi=sin(200*pi*i*dt);for(i=0;in;i+)fprintf(p,%10.7f,xi);fprintf(p,n);printf(%10.7f,xi);printf(n);fft(x,n);fprintf(p,n DISCRETE FOURIER TRANSFORMn);printf(n

9、 DISCRETE FOURIER TRANSFORMn);fprintf(p,%10.7f,x0);printf(%10.7f,x0);fprintf(p,%10.7f+J%10.7fn,x1,xn-1);printf(%10.7f+J%10.7fn,x1,xn-1);for(i=2;in/2;i+=2)fprintf(p,%10.7f+J%10.7f,xi,xn-i);fprintf(p,%10.7f+J%10.7f,xi+1,xn-i-1);fprintf(p,n);printf(%10.7f+J%10.7f,xi,xn-i);printf(%10.7f+J%10.7f,xi+1,xn-

10、i-1);printf(n);fprintf(p,%10.7f,xn/2);printf(%10.7f,xn/2);fprintf(p,%10.7f+J%10.7fn,xn/2-1,-xn/2+1);for(i=2;in/2;i+=2)fprintf(p,%10.7f+J%10.7f,xn/2-i,-xn/2+i);fprintf(p,%10.7f+J%10.7f,xn/2-i-1,-xn/2+i+1);fprintf(p,n);printf(%10.7f+J%10.7f,xn/2-i,-xn/2+i);printf(%10.7f+J%10.7f,xn/2-i-1,-xn/2+i+1);printf(n);将程序运行后所得数据绘制成曲线图（其中 FFT 变换的数据要先取绝对值后再画图）如下由上图可知，变换后的图开在频率100Hz 处出现一个峰值，这与理论上的结果一致。