分治法：快速傅里叶变换算法

分治法：快速傅里叶变换算法学院：网研院 姓名：xxx 学号：xxx分治法原理分治法是把一个复杂的问题分成两个或更多的相同或相似的子问题，再把子问题分成更小的子问题……直到最后子问题可以简单的直接求解，原问题的解即子问题的解的合并。分治法的精髓：分--将问题分解为规模更小的子问题；治--将这些规模更小的子问题逐个击破；合--将已解决的子问题合并，最终得出“母”问题的解；快速傅里叶变换（FFT）简介快速傅里叶变换（FastFourierTransform,FFT），是\o"离散傅里叶变换"离散傅里叶变换的快速\o"算法"算法，也可用于计算离散傅里叶变换的逆变换。它是根据离散傅里叶变换的奇、偶、虚、实等特性，对离散傅立叶变换的算法进行改进获得的。快速傅里叶变换有广泛的应用，如\o"数字信号处理"数字信号处理、计算\o"大整数乘法（页面不存在）"大整数乘法、求解\o"偏微分方程"偏微分方程等等。序列的离散傅里叶变换公式为：令则上式可写为：从算法分析角度：于是：分别考虑对其奇数项和偶数项作傅氏变换：则从而，可以将对N个量的傅氏变换变成为对两个规模更小的序列的变换。这样，将变换的量大大减小。快速傅里叶变换是分治法的一种具体实现。快速傅里叶变换算法FFT算法输入采样值；对采用值进行补0操作，使采样值的个数是2的幂次；对补0后的序列进行重排，重排的原则是按照序列的下标奇偶性排序，先偶后奇，分成两个子序列，然后对子序列继续重排。如1,2,3,4,5,6,7->1,3,5,7;2,4,6,8->1,5;3,7;2,6;4,8->1,5,3,7,2,6,4,8；对补0重排后的序列用三层嵌套的循环来完成M=log2N（阶数）次迭代，三层循环的功能是：最里的一层循环完成相同WNP的蝶形运算，中间的一层循环完成因子WNP的变化，而最外的一层循环则是完成M次迭代过程。算法复杂度令 ，则得从而快速傅氏变换的复杂性度为可能的改进本次实验使用的是最基本的基2FFT算法，根据维基百科，有如下比较出名的FFT算法，在性能上应该都优于本实现。Cooley-Tukey算法是最常见的FFT算法。这一方法以\o"分治法"分治法为策略\o"递归"递归地将长度为的\o"DFT"DFT分解为长度为的个较短序列的DFT，以及与个旋转因子的复数乘法。在\o"Cooley-Tukey快速傅里叶变换算法"Cooley-Tukey算法之外还有其他DFT的快速算法。对于长度且与互质的序列，可以采用基于\o"中国剩余定理"中国剩余定理的\o"互质因子算法"互质因子算法将N长序列的DFT分解为两个子序列的DFT。与Cooley-Tukey算法不同的是，\o"互质因子算法"互质因子算法不需要旋转因子。Rader-Brenner算法是类似于Cooley-Tukey算法，但是采用的旋转因子都是纯虚数，以增加加法运算和降低了数值稳定性为代价减少了乘法运算。这方法之后被split-radixvariantofCooley-Tukey所取代，与Rader-Brenner算法相比，有一样多的乘法量，却有较少的加法量，且不牺牲数值的准确性。\o"Bruun（页面不存在）"Bruun以及\o"QFT（页面不存在）"QFT算法是不断的把DFT分成许多较小的DFT运算。(Rader-Brenner以及QFT算法是为了2的指数所设计的算法，但依然可以适用在可分解的整数上。Bruun算法则可以运用在可被分成偶数个运算的数字)。尤其是Bruun算法，把FFT看成是，并把它分解成与的形式。另一个从多项式观点的快速傅立叶变换法是\o"威诺格拉德快速傅里叶变换算法"Winograd算法。此算法把分解成cyclotomic多项式，而这些多项式的系数通常为1，0，-1。这样只需要很少的乘法量(如果有需要的话)，所以winograd是可以得到最少乘法量的快速傅立叶算法，对于较小的数字，可以找出有效率的算方式。更精确地说，winograd算法让DFT可以用点的DFT来简化，但减少乘法量的同时，也增加了非常多的加法量。Winograd也可以利用剩余值定理来简化DFT。Rader算法提出了利用点数为N(N为质数)的DFT进行长度为N-1的回旋折积来表示原本的DFT，如此就可利用折积用一对基本的FFT来计算DFT。另一个prime-size的FFT算法为chirp-Z算法。此法也是将DFT用折积来表示，此法与Rader算法相比，能运用在更一般的转换上，其转换的基础为Z转换。算法实现框架本次实验使用的语言是java，只有一个类FFT.java，类有一个属性xConv用来存储补0和重排后的序列；FFT()构造方法对序列进行补0操作，然后调用i2Sort()方法对补0后序列进行重排，最后调用myFFT()方法进行快速傅里叶变换；i2Sort()法对补0后序列进行重排；myFFT()方法对补0重排后序列进行快速傅里叶变换操作；main()方法是程序的入口，用于输入数据并调用构造方法FFT()。如下图所示。i2Sort()方法对补0后序列进行重排，如1,2,3,4,5,6,7->1,5,3,7,2,6,4,8。myFFT()方法对补0重排后序列进行快速傅里叶变换操作，方法的原理是用三层嵌套的循环来完成M=log2N（阶数）次迭代，三层循环的功能是：最里的一层循环完成相同WNP的蝶形运算，中间的一层循环完成因子WNP的变化，而最外的一层循环则是完成M次迭代过程。FFT()构造方法对序列进行补0操作，然后调用i2Sort()方法对补0后序列进行重排，最后调用myFFT()方法进行快速傅里叶变换；main()方法是程序的入口，用于输入数据并调用构造方法FFT()。总结本科的时候曾经系统学习过傅里叶变换的知识，考研的时候也复习了该方面知识，在理解题目上没有遇到太大的困难。但由于从未尝试编程实现FFT，因此在本次实验的过程中还是遇到了许多问题。经过查看课件及网上查找资料，我对如何编码实现FFT有了一个较好的理解，知道了算法主要应该包括补0、排序和蝶形迭代3步操作。通过查找资料参考别人的代码，并结合自己的理解，我使用了java语言完成了该实验。由于时间较为紧张，加上期末考临近等原因，该实验完成的比较仓促，仍有许多不足的地方，但总的来说，通过本次试验，我再次复习了FFT的原理及了解了分治法的原理，同时锻炼了自己动手编写代码的能力，收获多多。附录程序运行示例：源程序：importjava.util.Scanner;publicclassFFT{ privatedouble[]xConv;//对x[n]进行二进制倒序排列的结果 publicFFT(double[]x)throwsException{ intn=x.length; if(n<1) thrownewException("采样序列的个数不能小于1！"); intm=(int)Math.ceil(Math.log(n)/Math.log(2)); intj=(int)Math.pow(2,m); xConv=newdouble[j];//初始化xConv for(inti=0;i<n;i++) xConv[i]=x[i]; for(inti=n;i<j;i++) xConv[i]=0; System.out.println("x[n]共有"+x.length+"个采样值"+'('+m+"阶)" +",补零个数为"+(j-n)); System.out.println("原x[n]数组："); for(inti=0;i<x.length;i++){ System.out.print(String.format("%6.2f",x[i])+""); } System.out.println(""); System.out.println("补零之后的x[n]数组："); for(inti=0;i<xConv.length;i++){ System.out.print(String.format("%6.2f",xConv[i])+""); } System.out.println(""); i2Sort(xConv,m); myFFT(xConv,m); } privatevoidi2Sort(double[]xConv2,intm){ int[]index=newint[xConv2.length];//index数组用于，倒序索引 int[]bits=newint[m]; double[]temp=newdouble[xConv2.length]; for(inti=0;i<xConv2.length;i++) //xConv2的原序映像 temp[i]=xConv2[i]; for(inti=0;i<index.length;i++){ index[i]=i;//第i个位置，倒序前的值为i for(intj=0;j<m;j++){ bits[j]=index[i]-index[i]/2*2;//提取index[i]的第j位二进制的值 index[i]/=2; } index[i]=0;//清零第i个位置的值 for(intj=m,power=1;j>0;j--){ index[i]+=bits[j-1]*power;//第i个位置，倒序后的位置 power*=2; } } System.out.println("二进制排序之后的x[n]数组："); for(inti=0;i<xConv2.length;i++){ //倒序实现 xConv2[i]=temp[index[i]]; System.out.print(String.format("%6.2f",xConv2[i])+""); } System.out.println(""); } privatevoidmyFFT(double[]xConv2,intm){ intdivBy;//divBy等分 double[]Xr,Xi,Wr,Wi;//分别表示：FFT结果的实部和虚部、旋转因子的实部和虚部 double[]tempXr,tempXi;//蝶形结果暂存器 intn=xConv2.length; doublepi=Math.PI; divBy=1; Xr=newdouble[n]; Xi=newdouble[n]; tempXr=newdouble[n]; tempXi=newdouble[n]; Wr=newdouble[n/2]; Wi=newdouble[n/2]; for(inti=0;i<n;i++){//初始化Xr、Xi，之所以这样初始化，是为了方便下面的蝶形结果暂存 Xr[i]=xConv2[i]; Xi[i]=0; } for(inti=0;i<m;i++){//共需要进行m次蝶形计算 divBy*=2; for(intk=0;k<divBy/2;k++){//旋转因子赋值 Wr[k]=Math.cos(k*2*pi/divBy); Wi[k]=-Math.sin(k*2*pi/divBy); } for(intj=0;j<n;j++){//蝶形结果暂存 tempXr[j]=Xr[j]; tempXi[j]=Xi[j]; } for(intk=0;k<n/divBy;k++){//蝶形运算：每一轮蝶形运算，都有n/2对的蝴蝶参与；n/2分为n/divBy组，每组divBy/2个。 intwIndex=0;//旋转因子下标索引 for(intj=k*divBy;j<k*divBy+divBy/2;j++){ doubleX1=tempXr[j+divBy/2]*Wr[wIndex] -tempXi[j+divBy/2]*Wi[wIndex]; doubleX2=tempXi[j+divBy/2]*Wr[wIndex] +tempXr[j+divBy/2]*Wi[wIndex]; Xr[j]=tempXr[j]+X1; Xi[j]=tempXi[j]+X2; Xr[j+divBy/2]=tempXr[j]-X1;//蝶形对两成员距离相差divBy/2 Xi[j+divBy/2]=