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

《测试信号分析及处理》课程作业快速傅里叶变换一、程序设计思路快速傅里叶变换的目的是减少运算量，其用到的方法是分级进行运算。全部计算分解为M级，其中M=log2N；在输入序列xG）中是按码位倒序排列的，输出序列XO是按顺序排列；每级包含N个蝶形单元，第i级有N个群，每个22i群有2i-1个蝶形单元；每个蝶形单元都包含乘Wr和-Wr系数的运算，每个蝶形NN单元数据的间隔为2i-1，i为第i级；同一级中各个群的系数W分布规律完全相同。将输入序列xG）按码位倒序排列时，用到的是倒序算法一一雷德算法。自然序排列的二进制数，其下面一个数总比上面的数大1，而倒序二进制数的下面一个数是上面一个数在最高位加1并由高位向低位仅为而得到的。若已知某数的倒序数是J，求下一个倒序数，应先判断J的最高位是否为0，与k=N进行比较即可得到结果。如果k〉J，说明最高位为0，应把其变成1，2即J+N，这样就得到倒序数了。如果k<J，即J的最高位为1，将最高位化为20，即J-史，再判断次高位；与k=史进行比较，若为0，将其变位1，即J+N，244即得到倒序数，如果次高位为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;i<n/2;i++){m=i;j=2*j;//得到流程图的共几级//得到流程图的共几级〃如果i<j,即进行变址}n1=n-1;for(j=0,i=0;i<n1;i++){if(i<j){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;i<n;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;i<n;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;}四、计算实例及运行结果设输入序列x(i)为x(i)=sin20SiAt,(i=0,1,2,...,n-1)其离散傅里叶变换为X(k)=N-x(i)Wik,(k=0,1,2,...,n-1)Ni=0・2K这里W=e-Jn。选n=512,计算离散傅里叶变换X(k)。所用软件为Turboc2.0，操作界面如图1所示图1Turboc2.0操作界面程序运行结束后的界面如图2所示鼬C:\Windows\5ystein32\cmdexe6.569065fc+Jll.76652756.4772524+J8.87838975.1?1163E+J7.B9?07454.32BJ632+J5.8?2Q1363.6933914+J5.01872B73.21970e5+J4.35707352.«41?7513+J3.83U(JMbl2.553W3U8+J3.41V5b«72.3099713+J3.07475092.1073772+J2.7654472l.V3bM?7^+J1.阳"WY+J幻3Zf>3¥A1.6627248+J2.14077721.5521563+11.97725401.4549735+J1.83191021.3689822+J1.78172291.292^30+,T1.F；R4^RBf；1.2239970+J1.47772341.162^914+J1.38843271.1070124+J1.29114971.05G798E+J1.2Q981451.0112154+J1.122541G6.9697290+J1.06158470.931B864+J0.9953102G.O??3099«J0.93317G4B.aGGtGSl■JG.07471?G0.836671E+J0.81952B40.3100851+J0.76725906.7856943+J0.71760070.7633147+JG.67023230.7427856+J0.62506370.7239659+J0.56173126.7067320+J0.54B09370.690?751+J0.49997930.6765996+J0.46123290.6635214+J0.4237136tl.bblbbb^+J虬38"戏8M.b4tJV6UV+J虬3-1852(!0.6313724+J0.3172840。上22B266+J队28343630.6152881+J0.25036520.6087193+J0.2178324PL朋即丽H+.TR.1A5HR44R.S9A^A71+.1fl.1F4231^6.5945342+J0.12294900.5915714+J0.0919726S.5394G4e+J0.0G120170.5882Q4GS.Q2Q5&G9图2程序运行后的界面例子的具体程序如下：#include<math.h>#include<stdio.h>#include<stdlib.h>#definepi3oidfft(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;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=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;i<n;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;i<n;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;i<n;i++){x[i]=sin(200*pi*i*dt);}for(i=0;i<n;i++){fprintf(p,"%10.7f",x[i]);fprintf(p,"\n");printf("%10.7f",x[i]);printf("\n");}fft(x,n);fprintf(p,"\nDISCRETEFOURIERTRANSFORM^");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;i<n/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;i<n/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-