《快速傅里叶变换原理及源程序》由会员分享，可在线阅读，更多相关《快速傅里叶变换原理及源程序（8页珍藏版）》请在金锄头文库上搜索。
1、测试信号分析及处理课程作业一、程序设计思路快速傅里叶变换的目的是减少运算量，其用到的方法是分级进行运算。全部计算分解为M级，其中M = log2N ；在输入序列xC）中是按码位倒序排列的， 输出序列X（k）是按顺序排列；每级包含N个蝶形单元，第i级有个群，每个22 i群有2 i-1个蝶形单元；每个蝶形单元都包含乘Wr和Wr系数的运算，每个蝶形NN单元数据的间隔为2i-1，i为第i级；同一级中各个群的系数W分布规律完全相 同。将输入序列xC）按码位倒序排列时，用到的是倒序算法一一雷德算法。自然 序排列的二进制数，其下面一个数总比上面的数大1，而倒序二进制数的下面一 个数是上面一个数在最高位加 1
2、 并由高位向低位仅为而得到的。若已知某数的倒序数是J，求下一个倒序数，应先判断J的最高位是否为0, 与k = N进行比较即可得到结果。如果k J,说明最高位为0,应把其变成1,2即J+N，这样就得到倒序数了。如果kJ,即J的最高位为1,将最高位化为20,即J -，再判断次高位；与k = 进行比较，若为0,将其变位1,即J + -,2 4 4即得到倒序数，如果次高位为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;/得到流程图的共几级/如果ivj,即进行变址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的下一个倒位序/如果kv(j+l),表示j的最高位为1j=j-k;/把最高位
4、变成0k=k/2;k/2,比较次高位，依次类推，逐个比较，直到某个位为0j=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；H押1.6627248+J 2.1407
5、772 1.5521563+1 1.9772540 1.4549735+J 1.8319102 1.3689822+J 1.7817229 1 .29230+,T 1 .F；R4RBf； 1 .2239970+J 1 .4777234 1.162914+J 1.3884327 1.1070124+J 1.2911497 1.05G798E+J 1.2Q98145 1.0112154+J 1.122541G 6.9697290+J 1.0615847 0.931B864+J 0.9953102 G.O?3099 J 0.93317G4 B.aGGtGSl J G.07471?G 0.836671
6、E+J 0.81952B4 0.3100851+J 0.7672590 6.7856943+J 0.7176007 0.7633丄叫7+J G.6702323 0.7427856+J 0.6250637 0.7239659+J 0.5617312 6.7067320+J 0.54B0937 0.690575丄十J 0.4999793 0.6765996+J 0.4612329 0.6635214+J 0.4237136 tl.bblbbb+J 虬3出吃丫跖0.6313724+J 0.3172840 0.&22B266+J 6.234363 0.6152881+J 0.2503652 0.608
7、7193+J 0.2178324R.A3finAn+.T R.1 A5HR44 R.S9AA71 +.1 fl.1F4231 6.5945342+J 0.1229490 0.5915714+J 0.0919726S.5394G4e+J 0.0G12017 0.5882Q4GS.Q2Q5&G9图 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
8、,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-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+n
9、2=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*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
10、=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 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;
11、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-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
