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

1、分治法：快速傅里叶变换算法学院：网研院姓名：xxx学号：xxx分治法原理分治法是把一个复杂的问题分成两个或更多的相同或相似的子问题， 再把子 问题分成更小的子问题,直到最后子问题可以简单的直接求解， 原问题的解即 子问题的解的合并。分治法的精髓：分-将问题分解为规模更小的子问题；治-将这些规模更小的子问题逐个击破；合-将已解决的子问题合并，最终得出“母”问题的解；二、快速傅里叶变换(FFTT简介快速傅里叶变换(Fast Fourier Transform, FFT ),是离散傅里叶变换的 快速算法，也可用于计算离散傅里叶变换的逆变换。它是根据离散傅里叶变换的奇、偶、虚、实等特性，对离散傅立叶变

2、换的算法进行改进获得的。快速傅里叶 变换有广泛的应用，如数字信号处理、计算大整数乘法、求解偏微分方程等等。序列xi, i =0,，N -1的离散傅里叶变换公式为：Nxk：=xne-jnk2 二 n三n =0令WN =e寸则上式可写为：N 1Xk八 xnW；nn 0从算法分析角度：NN，一12N 122Xk- xnWN；n =、x2n WN2nk 八 x2n 1WN2n1)kn =0n =0n =0于是：分别考虑对其奇数项和偶数项作傅氏变换：NN1 22XEk- x2n WN21k 八 x2nWNnk2 n =0n=0N dN d1 12 2XOk =、x2n 1 wN2nk =、 x2n 1W

3、Nk；2n=Sn=SN 4则 Xk=" xnwNkn =XEk W： XOk n z0从而，可以将对N个量的傅氏变换变成为对两个规模更小的序列的变换。这 样，将变换的量大大减小。快速傅里叶变换是分治法的一种具体实现。快速傅里叶变换算法FFT1 .算法输入采样值；对采用值进行补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重排后的

4、序列用三层嵌套的循环来完成 M=log2N (阶数)次迭 代，三层循环的功能是：最里的一层循环完成相同W的蝶形运算，中间的一层循环完成因子 WP的变化，而最外的一层循环则是完成M次迭代过程。2 .算法复杂度T(N)=2T(?)£T(2) =0令 N=2m ,则得 T(2m) = (m-1)2m，从而快速傅氏变换的复杂性度为O( N log N)3 .可能的改进本次实验使用的是最基本的基 2FFT算法，根据维基百科，有如下比较 出名的FFT算法，在性能上应该都优于本实现。Cooley-Tukey算法是最常见的FFT算法。这一方法以分治法为策略递归 地将长度为N =的DFT分解为长度为

5、M的A为个较短序列的DFT,以及与0(八)个旋转因子的复数乘法。在Cooley-Tukey算法之外还有其他DFT的快速算法。对于长度N = MM且M与冲互质的序列,可以采用基于中国剩余定理的互质因子算法将N长序列的DFT分解为两个子序列的DFT1与Cooley-Tukey 算法 不同的是，互质因子算法不需要旋转因子。Rader-Brenner算法是类似于Cooley-Tukey 算法，但是采用的旋转因 子都是纯虚数，以增加加法运算和降低了数值稳定性为代价减少了乘法运算。 这方法之后被 split-radix variant of Cooley-Tukey所取代，与Rader-Brenner算法

6、相比，有一样多的乘法量，却有较少的加法量，且不牺 牲数值的准确性。Bruun以及QFT算法是不断的把DFT分成许多较小的DFT运算。(Rader-Brenner以及QFT算法是为了 2的指数所设计的算法，但依然可以适 用在可分解的整数上。Bruun算法则可以运用在可被分成偶数个运算的数字 )尤其是Bruun算法，把FFT看成是七网一 1,并把它分解成zMl与/打+ azM +1的形式。另一个从多项式观点的快速傅立叶变换法是Winograd算法。此算法把了" - 1分 解成cyclotomic多项式，而这些多项式的系数通常为1,0, -1。 这样只需要很少的乘法量(如果有需要的话)，所

7、以winograd是可以得到最 少乘法量的快速傅立叶算法，对于较小的数字，可以找出有效率的算方式。 更精确地说，winograd算法让DFT可以用2人点的DFT来简化，但减少乘法 量的同时，也增加了非常多的加法量。Winograd也可以利用剩余值定理来简 化DFTRader算法提出了利用点数为N(N为质数)的DFT进行长度为N-1的回旋 折积来表示原本的DFT如此就可利用折积用一对基本的 FFT来计算DFT。 另一个prime-size的FFT算法为chirp-Z算法。此法也是将 DFT用折积来 表示，此法与Rader算法相比，能运用在更一般的转换上，其转换的基础为Z转换。四、算法实现框架本次

8、实验使用的语言是java ,只有一个类FFT.java ,类有一个属性xConv 用来存储补0和重排后的序列；FFT()构造方法对序列进行补0操作，然后调用 i2Sort()方法对补0后序列进行重排，最后调用myFFT()方法进行快速傅里叶变 换；i2Sort()法对补0后序列进行重排；myFFT()方法对补0重排后序列进行快 速傅里叶变换操作；main()方法是程序的入口，用于输入数据并调用构造方法 FFT()。如下图所示。1 import java.util-Scanner; 23 public class Pft private doublef xtonv；/' f ?打：父=二

9、爰二- =56+ public FFT(dotjbl& x) thr<>w5 Exemption |.private vaid 12Sort(double-f int m) |&4&5+ private void myFFT(double xtonv2 int m) 113lL23+ public static void main(trirlg args) throws Exception 1| 13G i2Sort()方法对补 0 后序列进行重排,如 1,2,3,4,5,6,7->1,5,3,7,2,6,4,835 3637 383g 4 41 42

10、 4344 4546 47 阴49 5051 52S3 5455 5657S8 59 60 61 6263private void i2Sort(double xSmv?int m) int index 工 new intxConv2. length; '''/ index不上里子.1. 生品之int bits = newdouble temp = new doubl&xConv2.length;for (int i =1 < xConvJ. length; i-i-K)/ xCorwN非定子班言temp(i' = xConv2i-for (int

11、 i = 0; i < index,length; i+) indexi = i; 7 常i十二 .-7 1-ifor (int j = 0; j < m; j) bits j = ±ndexi - indexi / 2 * 2; / 费.in<k)cL的箫-z indexi /= 2;)indfixi - Bj / 遭针i 也一for (int j - (n, power lj j > 0; j -) indexi +- bitsj - 1 * power; 'r二.三：二宣 power *= 2;)1System.out. printlnf, )s

12、for (int i = 0; 1 < xCofiv2.length; f+) /傻字士翟xCanv2i = temp-indeKiSystem . out. print(Stringyntjrm£rt( "6. 2f ", xConv2i ) + " " ) System.Qwt, printlnt1'r');)myFFT()方法又t补0重排后序列进行快速傅里叶变换操作，方法的原理是用 三层嵌套的循环来完成M=lo%N (阶数)次迭代，三层循环的功能是：最里的一层循环完成相同 W的蝶形运算，中间的一层循环完成因子WP的变

13、化，而最外的一层循环则是完成M次迭代过程。35- privite wid nyFFTfdoublel xconv即 int n) 66int div3 $f diwBydouble j匕,Xi J, hr毋；八 士三大于FFT_：3;二，叱二r父音至宁笠二/七工彳.ydwbltij 妙呻即.tempxij/耳E9int n - WCcmvN Jmcth;7fidfMjble pi - Math 4 Pjj71divBy - 1-?2Xr m wwdouble|科；73旭.newdoublcn|;741tMpJGr -ne* doublefn75tefapMi -ne*i dcMjblE ni

14、jM ww ck>ublefi / 2J;?7hl new doublcfn $ 2 E;78*OT Jint i =1 « n; in") / 扪一比也，妁-之之冶丁±SSFI-ff黑二二盘Mri »C0flv2lJ；SIXii 7iB21S3Fm (Int i - ej i < i”J / 冷若至防d九到“ZJ的f<»r flux k t t; k < dlv&> f 2i k+> / 云物引二Wr k | = Math .£-i?3(L x 2 pi / d五vB*'jB&am

15、p;Wllk - -hath P sin t»c 2 pidivfiy)酊)箝for (iiiit j = % j ¥ n： j+> / 贰，.要士力 f K七- ?Jrempxij xij：,a95-Gfw Lilt lit k < n / dlvfiy- k+) ( / 电，之荐.初一耗侦1tnf)二界=?电学= 口之士若im，dLwS»至.墀4divBy/史中二iM wlrwiejf -电；”梵受渡for (int j * k divey; j < k * divBy + diBy / zj jw) 孑9double XX = CenpXr

16、fj + dlvBy / 2j * Ur1&0- tempiXif j # divly / 1 * Wi wIndl?Midouble H2 - tpn iLj I divBy，21aL Mr hindex102t terpxrj + dlvey / 2 * Hiwindextlesxrj = t««xrj- xi;164Xxj teiiMi|jJ4- W2;1B5Mrj * dlvBy / 2- te«pMrj-XI; /a.j¥-t= E»-isdivBy/2Ibxij * divB-'y / 2* tsnnpxlj-JUt

17、107ulFidex+u16S)1电。110)111 1L25#5 ten., out. print jn( " = FTq. H );113For (znt i = i < n; i+4-JIM/ FFT此戛E罩115Sysfei*4Kgrintlnfstrinc +fpmwt(Xrfi)1U (Xifl，B - " ： F +117+ St ring. fojf (16 IS - 2 f * j Math .a6s(Xi i )> + * * j*J;1LB FFT()构造方法对序列进行补0操作，然后调用i2Sort()方法对补0后序列 进行重排，最后调用m

18、yFFT()方法进行快速傅里叶变换；E public FFT(double x) tlirows Exception int ft = m. length;Sif (n < 1)9throw new Ejcception("55要存到的个教不能少于工"j10int(n=(ini) Math*ceit(Math* lQ?(n) / H司thJoglZ)：11intj-(int) Math.pow(2T m)：1215 )tSrw=new doublejj /jf工皂比xConv16 4fOr(dnt i. = G; 1 <1+4-J17 xtonvji = x i

19、；18 Tor (dnt 1 = n; i < j; i-H-19 xCa(ivi = 0;L8 士ystum.口匕已print:】n(”n莫河"+ length + 1 ( ' + m + *：ij"20 + ,粗蜃+麦先” + Cj - n);21 System.(Jut.priritlnC* P.xn£t22 for (int 1 = 0; i / w.lengthj i-H-J 23System. cut. print St rin g. frmat 1,3nS. 2-F" j m 'i J + "*)25Syst

20、em,out.println("")j26System.out.println("nfiii： ,r);27For (int 1. * ©j, i < xConv*length; i+) 28Sj/stem.out.print(String./omwt(''6.2f", xConvi) + "2g39System, out. prirtJn("1h) j31i2&art(xtonvjr >>:篦myFFT(xConVj m);331main()方法是程序的入口 ,用于输入数据并调用构

21、造方法FFT()。1201121122123124125 12fi 127 12S 129 工即131 1及1331M 135public static void fnain(StringJ args) throws Exception Scanner in new Scanner(System*in);hhilt (trut) System, out. println(""*"*'*52!7：=；1!?FFU.i=si"*'*r*)；Syst em.out.println(0司联上累=也个二(二安泰)，亘毡注关 0)；int n = i

22、n.rextlnt()Syst em. out. printing 0声较入“+ n +'亲"三(克羚浮与今.主格君子)/叵车=成double x = new doublefn;for (int i = 0; i < nj i+j k l = in . nexrtDouble ();)Syst em L ou t, printl n (" 亘2士求三裁FFT斑清之近诅亲 M) jnew FFT(x);Syst em ,out ,println("11);本科的时候曾经系统学习过傅里叶变换的知识,考研的时候也复习了该方面知识，在理解题目上没有遇到太大的

23、困难。但由于从未尝试编程实现FFT,因止匕在本次实验的过程中还是遇到了许多问题。 经过查看课件及网上查找资料，我对 如何编码实现FFT有了一个较好的理解，知道了算法主要应该包括补 0、排序和 蝶形迭代3步操作。通过查找资料参考别人的代码，并结合自己的理解，我使用 了 java语言完成了该实验。由于时间较为紧张，加上期末考临近等原因，该实 验完成的比较仓促，仍有许多不足的地方，但总的来说，通过本次试验，我再次 复习了 FFT的原理及了解了分治法的原理，同时锻炼了自己动手编写代码的能力， 收获多多。附录程序运行示例:行 Problems a Javadoc Declaration 臣 Outlin

24、e |f Task Lit Q ConsoleFFT Java Application D:Program FilesXJavadkl.7.0 45binjavavr.ete 2014-* 昌？不融 FFT */。走*平-圣轮工自占本:（工卖之），岂生是东 r、一-三盘入0个柒EH（三手浮式已.空第三二）,三主沼家备2甘£疝贰-FT*港工豆开果=5Kn fi.it：1.102.203.3中:一率之苣前xn系或, Ip 10乙阴3*狗二君反坤平i_F咕班口后这5.505,50 号，翎目.触1.10 FFT工果1.59S%3.3G2s43.3Q-2.S43.30-5.96屯."

25、j7.97 = 2.2。章：1.J7 * e.QQ * 1.S7 * ' 2.20 * ： 7.97 " j2.2Q 3第 4.4的总结源程序：import java.util.Scanner;public class FFT private double xConv;/ 对 xn 进行二进制倒序排列的结果public FFT(double x) throws Exception int n = x.length;if (n < 1)throw new Exception(" 采样序列的个数不能小于1！ ");int m = (int) Math.c

26、eil(Math.log(n) / Math.log(2);int j = (int) Math.pow(2, m);xConv = new doublej; / 初始化 xConvfor (int i = 0; i < n; i+)xConvi = xi;for (int i = n; i < j; i+)xConvi = 0;System.out.println("xn 共有 " + x.length + " 个采样值" + '(' + m + " 阶 )"+ ", 补零个数为" +

27、 (j - n);System.out.println(" 原 xn 数组： ");for (int i = 0; i < x.length; i+) System.out.print(String.format("%6.2f", xi) + " ");System.out.println("");System.out.println(" 补零之后的xn 数组： ");for (int i = 0; i < xConv.length; i+) System.out.print(Str

28、ing.format("%6.2f", xConvi) + " ");System.out.println("");i2Sort(xConv, m);myFFT(xConv, m);private void i2Sort(double xConv2, int m) int index = new intxConv2.length; / index数组用于，倒序索引int bits = new intm;double temp = new doublexConv2.length;for (int i = 0; i < xConv2

29、.length; i+)/ xConv2 的原序映像tempi = xConv2i;for (int i = 0; i < index.length; i+) indexi = i; / 第 i 个位置，倒序前的值为ifor (int j = 0; j < m; j+) bitsj = indexi - indexi / 2 * 2; /提取 indexi 的第 j 位二进制的值indexi /= 2;indexi = 0; / 清零第 i 个位置的值for (int j = m, power = 1; j > 0; j-) indexi += bitsj - 1 * pow

30、er; /第 i 个位置，倒序后的位置power *= 2;System.out.println(" 二进制排序之后的xn 数组： ");for (int i = 0; i < xConv2.length; i+) / 倒序实现xConv2i = tempindexi;System.out.print(String.format("%6.2f", xConv2i) + " ");System.out.println("");private void myFFT(double xConv2, int m) 分别

31、表示：FFT结果的实部和虚部、旋转因子的实部和虚部蝶形结果暂存器初始化Xr、 Xi ，之所以这样初始化，是为了方便下面的蝶int divBy; / divBy 等分 double Xr, Xi, Wr, Wi; / double tempXr, tempXi; / int n = xConv2.length;double pi = Math.PI;divBy = 1;Xr = new doublen;Xi = new doublen; tempXr = new doublen; tempXi = new doublen;Wr = new doublen / 2;Wi = new doublen

32、 / 2;for (int i = 0; i < n; i+) / 形结果暂存Xri = xConv2i;Xii = 0;for (int i = 0; i < m; i+) /共需要进行m次蝶形计算divBy *= 2;旋转因子赋值for (int k = 0; k < divBy / 2; k+) /Wrk = Math.cos(k * 2 * pi / divBy);Wik = -Math.sin(k * 2 * pi / divBy);for (int j = 0; j < n; j+) /蝶形结果暂存tempXrj = Xrj;tempXij = Xij;for (int k = 0; k < n / divBy; k+) /蝶形运算：每一轮蝶形运算，都有n/2 对的蝴蝶参与；n/2 分为 n/divBy 组，每组divBy/2 个。int wIndex = 0;