小波分析快速傅里叶变换

1、快速傅里叶变换 快速傅里叶变换(FFT)并不是一种新的变换，而是离散博里叶变换(DFT)的一种快速算法。因此，为了很好地理解和掌握快速傅里叶变换，必须对离散傅里叶变换有充分的理解与掌握 。 由于DFT的计算量太大即使采用计算机也很难对问题进行实时处理，所以并没有得到真正的运用。直到1965年库利(J. W. Coo1ey)和图基(J. W. Tukey)在计算数学上发表了著名的“机器计算傅里叶级数的一种算法”的文章提出了DFT的一种快速算法后来又有桑德(G. Sande)和图基的快速算法相继出现，情况才发生了根本的改变。经过人们对算法的改进，发展和完善了一套高速有效的运算方法，使DFT的计算大

2、大简化，运算时间般可缩短一、二个数量级，从而使DFT的运算在实际中真正得到了广泛的加用。 1 直接计算DFT的问题及改进途径( )x n对于N点有现长序列10( )( )( )NnkNnX kDFT x nx n W0,1,1kN101( )( )( )NnkNkx nIDFT X kX k WN0,1,1nN二者的差别只在于 的指数符号不同，以及差一个常数乘因子 ,因而我们只讨论DFT正变换的运算量。IDFT的运算量是完全相同的。 NW1NDFT的运算量分析 一次复数乘法需用四次实数乘法和二次实数加法： 一次复数加法则需两次实数加法。()()()()ajb cjdacbdj adbc()()

3、()()ajbcjdacj bd( )X k10( )NnkNnx n WN1N 4N22(1)NNN( )X k2N(1)N N 24N2(21)NN 一个个复数乘法复数加法实数乘法实数加法改进途径改进途径*nknkNNWW()()nkN n kn N kNNNWWWnkmnkNmNWW/nknk mNN mWW01NW/21NNW (/2)kNkNNWW 2jnknkNNWe利用 的特性 (1)对称性(2)周期性 (3) 可约性 (4) 特殊点 FFT算法的基本思想 使DFT运算中有些项可以合并； 可以将长序列的DFT分解为短序列的DFT FFT算法分类: 时间抽选法 DIT: Decim

4、ation-In-Time 频率抽选法 DIF: Decimation-In-Frequency2 按时间抽选(DIT)的基-2 FFT算法(库利图基算法) 算法原理 先没序列点数为 , 为整数。如果不满足这个条件，可以人为地加上若干零值点，使之达到这一要求。这种 为2的整数幂的FFT也称基-2 FFT。 2LN LN2LN ( )x n(0,1,1)nN将 的序列 先按n的奇偶分成以下两组： 12(2 )( ),0,1,1(21)( )2xrx rNrxrx r111000( )( )( )( )( )NNNnknknkNNNnnnX kDFT x nx n Wx n Wx n Wn为偶数

5、n为奇数 11222(21)00(2 )(21)NNrkrkNNrrxr WxrW1122221200( )( )NNrkrkkNNNrrx rWWx rW 式中 和 分别是 和 的N/2点DFT利用系数 的可约性： nkNW22222/2jNjNNNWeeW11221/22/200( )( )( )NNrkkrkNNNrrX kx r WWx r W12( )( )kNX kW Xk1( )X k2( )Xk1( )x r2( )x r112211/2/200( )( )(2 )NNrkrkNNrrX kx r Wxr W112222/2/200( )( )(21)NNrkrkNNrrXkx

6、 r WxrW可看出，一个N点DFT已分解成两个N2点的DFT，它们又组合成一个N点DFT。 然而 和 以及 和 都是N/2点的序列：1( )X k2( )Xk1( )x r2( )x r,0,1,12Nr k ( )X k 却有N点而用上式计算得到的只是 的前一半项数的结果，要用 和 来表达 全部的值，还必须应用系数的周期性，即： ( )X k1( )X k2( )Xk( )X k()2/2/2Nr krkNNWW1122()211/21/2100()( )( )( )2NNNr krkNNrrNX kx r Wx r WX k1122()222/22/2200()( )( )( )2NNN

7、r krkNNrrNXkx r Wx r WXk将将 表达为前后两部分表达为前后两部分 再考虑到 的以下性质： nkNW()/22NkNkkNNNNWWWW ( )X k前半部分 0,1,12Nk 12( )( )( )kNX kX kW Xk后半部分 212()()()222NkNNNNX kX kWXk12( )( )kNX kW Xk0,1,12Nk 可以用蝶形信号流图符号表示，当支路上没有标出系数时，则该支路的传输系数为1。 2LN 3L 计算量分析：复数乘法复数加法一个N2点DFT 两个N2点DFT 22N(1)22NN22N(1)2NN一个蝶形运算 122NN222NN2(1)22

8、NNNNN点DFT 总计 2N(1)N N N2个蝶形运算 计算量减少一半 可以进一步把每一个N2点子序列再按其奇偶部分分解为两个N4点的子序列 1( )x r2LN N2仍是偶数 1314(2 )( ),0,1,1(21)( )4xlx lNlxlx l34( ),( )x lx l分解为奇偶两个序列11442(21)11/21/200( )(2 )(21)NNlklkNNllX kxl WxlW1144223/2/24/200( )( )NNlklkkNNNllx lWWx lW222224/2/4jjNNNNWeeW由于114413/4/24/400( )( )( )NNlkklkNNN

9、llX kx l WWx l W3/24( )( )kNXkWXk0,1,14Nk 其中， 与 分别是 与 的N4点的DFT： 3( )Xk4( )Xk3( )x l4( )x l114433/41/400( )( )(2 )NNlklkNNllXkx l Wxl W114444/41/400( )( )(21)NNlklkNNllXkx l WxlW同样可以计算 在N/4到N/2-1的值3( )Xk0,1,14Nk ()4/4/4Nl klkNNWW()/44/2/2/2/2NkNkkNNNNWWWW 由于1144()433/43/4300()( )( )( )4NNNl klkNNllNX

10、kx l Wx l WXk1144()444/41/4400()( )( )( )4NNNl klkNNllNXkx l Wx l WXk0,1,14Nk ()413/24()()()444NkNNNNX kXkWXk因此：3/24( )( )kNXkWXk2526(2 )( ),0,1,1(21)( )4xlx lNlxlx l25/26( )( )( )kNXkXkWXk25/26()( )( )4kNNXkXkWXk0,1,14Nk 114455/42/400( )( )(2 )NNlklkNNllXkx l Wxl W114466/42/400( )( )(21)NNlklkNNllX

11、kx l WxlW其中， 与 分别是 与 的N4点的DFT： 5( )Xk6( )Xk5( )x l6( )x lN=8的DFT可以分解为4个 点DFT:24N2/2kkNNWWrNW由于： 将系数统一成 形式 计算量又减少一半计算量又减少一半讨论 按偶数与奇数的分解过程中序列标号的变化。对于一个N8点的DFT的例子，输入序列按偶数点与奇数点第一次分解为两个N2点序列： 再经过一起分解，剩下的是2点DFT，对于此例N8，就是4个 点DFT，其输出为 24N3( )Xk4( )Xk5( )Xk6( )Xk0,1k 可以方便地计算出来，例如：0033232(0)(0)(1)(0)(4)(0)(4)

12、XxW xxW xxx110332322(1)(0)(1)(0)(4)(0)(4)(0)(4)XxW xxW xxW xxx0044242(0)(0)(1)(2)(6)(2)(6)XxW xxW xxx110442422(1)(0)(1)(2)(6)(2)(6)(2)(6)XxW xxW xxW xxx0055252(0)(0)(1)(1)(5)(1)(5)XxW xxW xxx110552522(1)(0)(1)(1)(5)(1)(5)(1)(5)XxW xxW xxW xxx0066262(0)(0)(1)(3)(7)(3)(7)XxW xxW xxx110662622(1)(0)(1)(

13、3)(7)(3)(7)(3)(7)XxW xxW xxW xxx偶序列 奇序列1(2 )( )xrx r2(21)( )xrx rr 0，1，2，3 0，1，2，3 n0，2，4，6 1，3，5，7 偶序列奇序列 偶序列奇序列 13(2 )( )xlx l14(21)( )xlx l25(2 )( )xlx l26(21)( )xlx llrn0，10，1 0，1 0，1 0，21，3 0，2 1，3 0，42，6 1，5 3，7 的排列规律，将 写成二进制，将二进制的数码翻转，即得到正确的输入序列 0,1,1nN( )x n这种方法的每一步分解那是按输入序列在时间上的次序是属于偶数还是属于奇

14、数来分解为两个更短的子序列，所以称为“按时间抽选法”。000,001,010,011,100,101,110,111000,100,010,110,001,101,011,111(0), (1), (2), (3), (4), (5), (6), (7)xxxxxxxx(0), (4), (2), (6), (1), (5), (3), (7)xxxxxxxx以N=8为例称为到位序规律运算量分析运算量分析： 由按时间抽选法FFT的流图可见，当 时，共有L级蝶形，每级都由N/2个蝶形运算组成，每个蝶形有一次复乘、二次复加因而每级运算都需N/2次复乘和N次复加，这样L级运算总共需要:2LN 2lo

15、g22FNNmLN2logFaNLNN复乘数 复加数 实际计算量与这个数字稍有不同，因为 01NW/4NNWj 这几个系数都不用乘法运算，当N较大时，这些特例相对而言就很少。 出于计算机上乘法运算所用时间比加法运算所需时间多得多、故以乘法为例，在直接计算DFT与FFT算法的计算量之比为：222222logloglog22NNNNNNNNN2N2log2NN22log2NNN2414.041644.0864125.416256328.032 10248012.864 409619221.4128 1639444836.6256 65536102464.05122621142304113.8102

16、410485765120204.8算法特点算法特点 1 同址运算 2 到位序规律 3蝶形运算两节点的“距离” 4 的确定rNW11( )( )( )( )( )( )rmmNmrmmNmxpxpW xqxqxpW xqp, q是参与本碟形单元运算的上下节点的序号。第m级序号为p, q的两点只参与这一个碟形单元的运算，其输出在第m+1级，并且这一碟形单元不涉及别的点，即某一列的N个数据送到存储器后，经蝶形运算，其结果为另列数据，它们以蝶形为单位仍存储在这同一组存储器中，直到最后输出，中间无需其他存储器。也就是蝶形的两个输出值仍放回蝶形的两个输人所在的存储器中。每列的N2个蝶形运算全部完成后，再开

17、始下一列的蝶形运算。这样存储数据只需N个存储单元。下一级的运算仍采用这种原位方式，只不过进人蝶形结的组合关系有所不同。这种原位运算结构可以节省存储单元，降低设备成本。称同址运算以8点FFT为例，其输入是倒位序的，输出是自然顺序的，其第一级(第一列)每个蝶形的两节点间“距离”为1，第二级每个蝶形的两节点“距离为2，第二级每个蝶形的两节点“距离”为4，由此类推得，对 点FFT，当输人为倒位序，输出为正常顺序时，其第m级运算，每个蝶形的两节点“距离”为 。蝶形运算两节点的“距离”2LN 12m（从左到右，第一级，第二级，第m级）从最右边一级开始，看 分布规律rNW1111111( )( )(2)(2

18、)( )(2)rmmmNmmrmmmNmxpxpW xpxpxpW xp的确定842(2),0,1,123,0,1,2,32,0,11,0LrNrrrNmLN WrmWrmWrmWr12,0,1,21mrmWr任意第m级4按频率抽选(DIF)的基-2 FFT算法(桑德图基算法)这里讨论另一种FFT算法，称为按频率抽选(DIF)的FFT算法，它是把输出序列 (也是N点序列)按其顺序的奇偶分解为越来越短的序列。( )X k算法原理仍设序列点数为 ，L为整数。在把输出 按k的奇偶分组之前，先把输入按n的顺序分成前后两半(注意，这不是频率抽选)；2LN ( )X k1112002( )( )( )(

19、)NNNnknknkNNNNnnnX kx n Wx n Wx n W1122()2001220( )()2( )()2NNNnknkNNnnNNknkNNnNx n Wx nWNx nx nWW/2nkNWnkNW式中用的 是，不是 ，因而这并不是N2点DFT。/21NNW /2( 1)NkkNW 因为所以120( )( )( 1)()2NknkNnNX kx nx nW ( 1)1k( 1)1k k为偶数时k为奇数时221krkr0,1,12Nr 按照k的奇偶将 分为两部分( )X k11222/200(2 )( )()( )()22NNnrnrNNnnNNXrx nx nWx nx nW

20、11222/200(21)( )()( )()22NNnrnnnrNNNNnnNNXrx nx nWWx nx nW W12( )( )()2( )( )()2nNNx nx nx nNx nx nx nW0,1,12Nn 令121/20122/20(2 )( )(21)( )NnrNnNnrNnXrx n WXrx n W0,1,12Nr 与时间抽选法的推导过程一样，由于 ，N/2仍是个偶数，因而可以将每个N/2点DFT的输出再分解为偶数组与奇数组，这就将N/2点DFT进一步分解为两个N/4点DFT。这两个N4点DFT的输入也是先将N/2点DFT的输入上下对半分开后通过蝶形运算而形成。这样的

21、分解可以一直进行到第L次( )，第L次实际L是做两点DFT，它只有加减运算。但是，为了比较并为了统一运算结构，我们仍然采用系数为 的蝶形运算来表示，这N/2个两点DFT的N个输出就是 的N点DFT的结果 。2LN 2LN nNW( )x n( )X k算法特点 1 同址运算 2蝶形运算两节点的“距离” 3 的确定rNW1111( )( )( )( )( )( )mmmrmmmNXkXkXjXjXkXj W以8点FFT为例，其输入是自然顺序的，输出是倒位序的，其第一级(第一列)每个蝶形的两节点间“距离”为4，第二级每个蝶形的两节点“距离为2，第二级每个蝶形的两节点“距离”为1，由此类推得，对 点

22、FFT，当输人为自然顺序，输出为倒位序时，其第m级运算，每个蝶形的两节点“距离”为 。蝶形运算两节点的“距离”2LN 22L mmN（从左到右，第一级，第二级，第m级）从最左边一级开始，看 分布规律rNW11,0,1,1222,0,2,242,2,0,1,122rNrNrlNllNNmWrNNmWrNNmlWrkk的计算1,2,0,1,12rmNmNWrkk任意第m级1111( )( )()2()( )()22mmmmrmmmNmmNXkXkXkNNXkXkXkW按频率抽选法与按时间抽选法的异同DIF的基本蝶形与DIT的基本蝶形运算顺序有所不问，这才是实质的不同，DIF的复数乘法只出现在减法之

23、后，DIT则是先作复乘后作加减法。DIF与DIT就运算量来说则是相同的，即都有L级(列)运算，每级运算需N2个蝶形运算来完成，总共需要 次复乘与 次复加 DIF法与DIT法都可进行同址运算。2log22FNNmLN2logFaNLNN只要把DFT运算中的每一个系数 换成 ，最后再乘以常数1/N，则以上所有按时间抽选或按频率抽选的FFT都可以拿来运算IDFT。上面的FFT算法同样可以适用于离散傅里叶反变换(1DFT)运算，即快速傅里叶反变换(IFFT)。比较IDFT和DFT公式:离散傅里叶反变换(IDFT)的快速计算方法10101( )( )( )( )NnkNkNnkNnx nX k WNX

24、kx n WnkNWnkNW上面这种IFFT算法虽然编程很方便，但是需要稍稍改动FFT的程序和参数才能实现。下面讨论一种完全不用改变FFT的程序就可以计算IFFT的方法。我们对IDFT公式式取共轭：1*01( )( )NnkNkx nXk WN*1*011( )( )( )NnkNkx nXk WDFT XkNN则这说明，只要先将X(k)取共轭，就可以直接利用FF丁子程序，最后再将运算结果取一次共扼，并乘以1/N，即得到x(n)值。因此，FFT运算和IFFT运算就可以共用一个子程序块。*nknkNNWW基-4FFT算法 时就是基-4FFT算法4LN 3111142430424( )( )( )( )( )NNNNnknknknkNNNNNNNnnnnX kx n Wx n Wx n Wx n W111134444()()()42400003( )()()()424NNNNNNNnknknknkNNNNnnnnNNNx n Wx nWx nWx nW /4/23/43( )()()()424NkNkNknkNNNNNNNx nx nWx nWx nWW频率抽取的基-4FFT算法4 ,41,42,43,0,1,14Nkr krkrkrr/4/23/41,1,kNNNNNNNNWWj WWj 144014/403(4 )( )()()