快速傅里葉變換(FFT)原理及源程序

1、精選文檔測(cè)試信號(hào)分析及處理課程作業(yè)快速傅里葉變換1、 程序設(shè)計(jì)思路快速傅里葉變換的目的是削減運(yùn)算量，其用到的方法是分級(jí)進(jìn)行運(yùn)算。全部計(jì)算分解為級(jí)，其中；在輸入序列中是按碼位倒序排列的，輸出序列是按挨次排列；每級(jí)包含個(gè)蝶形單元，第級(jí)有個(gè)群，每個(gè)群有個(gè)蝶形單元； 每個(gè)蝶形單元都包含乘和系數(shù)的運(yùn)算，每個(gè)蝶形單元數(shù)據(jù)的間隔為，i為第i級(jí)； 同一級(jí)中各個(gè)群的系數(shù)分布規(guī)律完全相同。將輸入序列按碼位倒序排列時(shí)，用到的是倒序算法雷德算法。 自然序排列的二進(jìn)制數(shù)，其下面一個(gè)數(shù)總比上面的數(shù)大1，而倒序二進(jìn)制數(shù)的下面一個(gè)數(shù)是上面一個(gè)數(shù)在最高位加1并由高位向低位僅為而得到的。 若已知某數(shù)的倒序數(shù)是，求下一個(gè)倒序數(shù)，應(yīng)

2、先推斷的最高位是否為0，與進(jìn)行比較即可得到結(jié)果。假如，說明最高位為0，應(yīng)把其變成1，即，這樣就得到倒序數(shù)了。假如，即的最高位為1，將最高位化為0，即，再推斷次高位；與進(jìn)行比較，若為0，將其變位1，即，即得到倒序數(shù)，假如次高位為1，將其化為0，再推斷下一位即從高位到低位依次推斷其是否為1，為1將其變位0，若這一位為0，將其變位1，即可得到倒序數(shù)。若倒序數(shù)小于挨次數(shù)，進(jìn)行換位，否則不變，防治重復(fù)交換，變回原數(shù)。注：由于0的倒序數(shù)為0，所以可從1開頭進(jìn)行求解。2、 程序設(shè)計(jì)框圖（1）倒序算法雷德算法流程圖(2)FFT算法流程3、 FFT源程序void fft(x,n)int n;double x;i

3、nt i,j,k,l,m,n1,n2;double c,c1,e,s,s1,t,tr;for(j=1,i=1;i<n/2;i+) m=i;j=2*j;if(j=n)break; /得到流程圖的共幾級(jí)n1=n-1;for(j=0,i=0;i<n1;i+)if(i<j) /假如i<j,即進(jìn)行變址tr=xj; xj=xi;xi=tr;k=n/2; /求j的下一個(gè)倒位序while(k<(j+1) /假如k<(j+1),表示j的最高位為1j=j-k; /把最高位變成0k=k/2; /k/2，比較次高位，依次類推，逐個(gè)比較，直到某個(gè)位為0j=j+k; /把0改為1for

4、(i=0;i<n;i+=2)tr=xi;xi=tr+xi+1;xi+1=tr-xi+1;n2=1;for(l=1;l<=m;l+) / 把握蝶形結(jié)級(jí)數(shù)n4=n2;n2=2*n4; n1=2*n2;e=6.28318530718/n1;for(i=0;i<n;i+=n1) /把握同一蝶形結(jié)運(yùn)算，即計(jì)算系數(shù)相同蝶形結(jié)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+) /把握計(jì)算不同種蝶形結(jié)，即計(jì)算系數(shù)不同的蝶形結(jié)i1=i+j;i2=i-j+n2;i3=i+j+n2;i4=

5、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;4、 計(jì)算實(shí)例及運(yùn)行結(jié)果設(shè)輸入序列為其離散傅里葉變換為這里。選n=512,計(jì)算離散傅里葉變換。所用軟件為Turbo c 2.0，操作界面如圖1所示圖1 Turbo c 2.0操作界面程序運(yùn)行結(jié)束后的界面如圖2所示圖2 程序運(yùn)行后的界面例子的具體程序如下：#include<math.h>#include<stdio.h>#include<stdl

6、ib.h>#define pi 3oid 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;i<n/2;i+) m=i;j=2*j;if(j=n)break;n1=n-1;for(j=0,i=0;i<n1;i+)if(i<j)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;i<n;i+=2)tr=xi;xi=tr+

7、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;i<n;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*xi3-cc*xi4;xi4=xi2-t2;xi3=-xi

8、2-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;i<n;i+)xi=sin(200*pi*i*dt);for(i=0;i<n;i+) fprintf(p,"%10.7f",xi);fprintf(p,"n");printf("%10.7f",xi);printf("n");fft(x

9、,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;i<n/2;i+=2)fprintf(p,"%

10、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;i<n/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