FFT-C快速傅里叶变换超级详细的原代码

快速傅立叶变换（FFT）的C++实现收藏标准的离散傅立叶DFT变换形式如：yk=Σj=0n-1ajωn-kj=A(ωn-k).（ωnk为复数1的第k个n次方根，且定义多项式A(x)=Σj=0n-1ajxj）而离散傅立叶逆变换IDFT(InverseDFT)形式如：aj=(Σk=0n-1ykωnkj)/n.yk=Σj=0n-1ajωn-kj=A(ωn-k).（ωnk为复数1的第k个n次方根，且定义多项式A(x)=Σj=0n-1ajxj）而离散傅立叶逆变换IDFT(InverseDFT)形式如：aj=(Σk=0n-1ykωnkj)/n.以下不同颜色内容为引用并加以修正：快速傅立叶变换（FastFourierTransform，FFT）是离散傅立叶变换（DiscreteFouriertransform，DFT）的快速算法，它是根据离散傅立叶变换的奇、偶、虚、实等特性，对离散傅立叶变换的算法进行改进获得的。它对傅立叶变换的理论并没有新的发现，但是对于在计算机系统或者说数字系统中应用离散傅立叶变换，可以说是进了一大步。设Xn为N项的复数序列，由DFT变换，任一Xi的计算都需要N次复数乘法和N-1次复数加法，而一次复数乘法等于四次实数乘法和两次实数加法，一次复数加法等于两次实数加法，即使把一次复数乘法和一次复数加法定义成一次“运算”（四次实数乘法和四次实数加法），那么求出N项复数序列的Xi，即N点DFT变换大约就需要N2次运算。当N=1024点甚至更多的时候，需要N2=1048576次运算，在FFT中，利用ωn的周期性和对称性，把一个N项序列（设N为偶数），分为两个N/2项的子序列，每个N/2点DFT变换需要(N/2)2次运算，再用N次运算把两个N/2点的DFT变换组合成一个N点的DFT变换。这样变换以后，总的运算次数就变成N+2*(N/2)2=N+N2/2。继续上面的例子，N=1024时，总的运算次数就变成了525312次，节省了大约50%的运算量。而如果我们将这种“一分为二”的思想不断进行下去，直到分成两两一组的DFT运算单元，那么N点的DFT变换就只需要N*log2N次的运算，N=1024点时，运算量仅有10240次，是先前的直接算法的1%，点数越多，运算量的节约就越大，这就是FFT的优越性。FFT的实现可以自顶而下，采用递归，但是对于硬件实现成本高，对于软件实现都不够高效，改用迭代较好，自底而上地解决问题。感觉和归并排序的迭代版很类似，不过先要采用“位反转置换”的方法把Xi放到合适的位置，设i和j互为s=log2N位二进制的回文数，假设s=3,i=(110)2=6,j=(011)2=3,如果i≠j,那么Xi和Xj应该互换位置。（关于这个回文数的生成，是很有趣而且是很基本的操作，想当初偶初学C++的时候就有这样的习题。）当“位反转置换”完成后，先将每一个Xi看作是独立的多项式，然后两个两个地将它们合并成一个多项式（每个多项式有2项），合并实际上是“蝶形运算”（ButterflyOperation,参考《算法导论》吧^_^），继续合并（第二次的每个多项式有4项），直到只剩下一个多项式（有N项），这样，合并的层数就是log2N，每层都有N次操作，所以总共有N*log2N次操作。迭代过程如下图所示，自底而上。C++源程序，如下：////快速傅立叶变换FastFourierTransform//[e-mail]//2007-07-20//版本2.0//改进了《算法导论》的算法，旋转因子取ωn-kj(ωnkj的共轭复数)//且只计算n/2次，而未改进前需要计算(n*lgn)/2次。//#includestdio.h#includestdlib.h#includemath.hconstintN=1024;constfloatPI=3.1416;inlinevoidswap(float&a,float&b){floatt;t=a;a=b;b=t;}voidbitrp(floatxreal[],floatximag[],intn){//位反转置换Bit-reversalPermutationinti,j,a,b,p;for(i=1,p=0;in;i*=2){p++;}for(i=0;in;i++){a=i;b=0;for(j=0;jp;j++){b=(b1)+(a&1);//b=b*2+a%2;a=1;//a=a/2;}if(bi){swap(xreal[i],xreal[b]);swap(ximag[i],ximag[b]);}}}voidFFT(floatxreal[],floatximag[],intn){//快速傅立叶变换，将复数x变换后仍保存在x中，xreal,ximag分别是x的实部和虚部floatwreal[N/2],wimag[N/2],treal,timag,ureal,uimag,arg;intm,k,j,t,index1,index2;bitrp(xreal,ximag,n);//计算1的前n/2个n次方根的共轭复数W'j=wreal[j]+i*wimag[j],j=0,1,,n/2-1arg=-2*PI/n;treal=cos(arg);timag=sin(arg);wreal[0]=1.0;wimag[0]=0.0;for(j=1;jn/2;j++){wreal[j]=wreal[j-1]*treal-wimag[j-1]*timag;wimag[j]=wreal[j-1]*timag+wimag[j-1]*treal;}for(m=2;m=n;m*=2){for(k=0;kn;k+=m){for(j=0;jm/2;j++){index1=k+j;index2=index1+m/2;t=n*j/m;//旋转因子w的实部在wreal[]中的下标为ttreal=wreal[t]*xreal[index2]-wimag[t]*ximag[index2];timag=wreal[t]*ximag[index2]+wimag[t]*xreal[index2];ureal=xreal[index1];uimag=ximag[index1];xreal[index1]=ureal+treal;ximag[index1]=uimag+timag;xreal[index2]=ureal-treal;ximag[index2]=uimag-timag;}}}}voidIFFT(floatxreal[],floatximag[],intn){//快速傅立叶逆变换floatwreal[N/2],wimag[N/2],treal,timag,ureal,uimag,arg;intm,k,j,t,index1,index2;bitrp(xreal,ximag,n);//计算1的前n/2个n次方根Wj=wreal[j]+i*wimag[j],j=0,1,,n/2-1arg=2*PI/n;treal=cos(arg);timag=sin(arg);wreal[0]=1.0;wimag[0]=0.0;for(j=1;jn/2;j++){wreal[j]=wreal[j-1]*treal-wimag[j-1]*timag;wimag[j]=wreal[j-1]*timag+wimag[j-1]*treal;}for(m=2;m=n;m*=2){for(k=0;kn;k+=m){for(j=0;jm/2;j++){index1=k+j;index2=index1+m/2;t=n*j/m;//旋转因子w的实部在wreal[]中的下标为ttreal=wreal[t]*xreal[index2]-wimag[t]*ximag[index2];timag=wreal[t]*ximag[index2]+wimag[t]*xreal[index2];ureal=xreal[index1];uimag=ximag[index1];xreal[index1]=ureal+treal;ximag[index1]=uimag+timag;xreal[index2]=ureal-treal;ximag[index2]=uimag-timag;}}}for(j=0;jn;j++){xreal[j]/=n;ximag[j]/=n;}}voidFFT_test(){charinputfile[]=input.txt;//从文件input.txt中读入原始数据charoutputfile[]=output.txt;//将结果输出到文件output.txt中floatxreal[N]={},ximag[N]={};intn,i;FILE*input,*output;if(!(input=fopen(inputfile,r))){printf(Cannotopenfile.);exit(1);}if(!(output=fopen(outputfile,w))){printf(Cannotopenfile.);exit(1);}i=0;while((fscanf(input,%f%f,xreal+i,ximag+i))!=EOF){i++;}n=i;//要求n为2的整数幂while(i1){if(i%2){fprintf(output,%disnotapowerof2!,n);exit(1);}i/=2;}FFT(xreal,ximag,n);fprintf(output,FFT:irealimag);for(i=0;in;i++){fprintf(output,%4d%8.4f%8.4f,i,xreal[i],ximag[i]);}fprintf(output,=================================);IFFT(xreal,ximag,n);fprintf(output,IFFT:irealimag);for(i=0;in;i++){fprintf(output,%4d%8.4f%8.4f,i,xreal[i],ximag[i]);}if(fclose(input)){printf(Filecloseerror.);exit(1);}if(fclose(output)){printf(Filecloseerror.);exit(1);}}intmain(){FFT_test();return0;}//////////////////////////////////////////////////sample:input.txt////////////////////////////////////////////////0.575100.451400.043900.027200.312700.012900.384000.683100.092800.035300.612400.608500.015800.016400.190100.58690//////////////////////////////////////////////////sample:output.txt////////////////////////////////////////////////FFT:irealimag04.64850.000010.01760.312221.11140.042931.6776-0.13534-0.23401.389750.3652-1.25896-0.43250.20737-0.13120.37638-0.19490.00009-0.1312-0.376310-0.4326-0.2073110.36521.258912-0.2340-1.3897131.67760.1353141.1113-0.0429150.0176-0.3122=================================IFFT:irealimag00.5751-0.000010.45