第5章 离散傅立叶变换与快速算法2(FFT)

1、15.2DFT的快速计算的快速计算(FFT)n一、一、DFT存在的问题及改进的途径存在的问题及改进的途径n二、戈泽尔算法二、戈泽尔算法n三、时域抽取基三、时域抽取基-2算法（库利算法（库利-图基算法）图基算法）n四、频域抽取基四、频域抽取基-2算法（桑德算法（桑德-图基算法）图基算法）n五、线性调频五、线性调频Z变换变换2nN点有限长序列，其点有限长序列，其DFT变换对为：变换对为：nn 可以看出，可以看出，变换与反变换变换与反变换的差别仅在于的差别仅在于 的指数符号和常数乘因子的指数符号和常数乘因子1/N。实际上：实际上：n 因而我们因而我们只讨论只讨论DFT正变换的运算量正变换的运算量，反

2、，反变换的运算量是完全相同的。变换的运算量是完全相同的。 )()(1)()()()(1010nRWkXNnxIDFTkRWnxkXDFTNNknkNNNnnkN：)(1)(*kXDFTNkXIDFT5.2.1DFT存在的问题及改进途径存在的问题及改进途径3n一般来说一般来说 和和 都是复数，都是复数，n每一个每一个 值的计算，需要值的计算，需要N次复数乘次复数乘法法以及（以及（N-1）次复数加法，次复数加法，n完成整个完成整个DFT运算总需要运算总需要 次复数乘次复数乘法法和和N（N-1）次复数加法次复数加法。)(nxnkNW)(kX2N10( )( )( )NnkNNnDFTX kx n W

3、Rk：DFT直接计算的运算量4n复数运算实际上是由实数运算来完成的，如：复数运算实际上是由实数运算来完成的，如：n一次复数乘法需用：一次复数乘法需用：n四次实数乘法四次实数乘法n二次实数加法；二次实数加法；n一次复数加法则需：一次复数加法则需：n二次实数加法。二次实数加法。n整个整个DFT运算总共需要：运算总共需要：n4N2次实数乘法次实数乘法n 次次实数加法实数加法。)()()()()()()(dcjbajdbjcaBAbcadjcdabjdbjcaAB) 12(2) 1(222NNNNN5n上述统计与实际的运算次数有上述统计与实际的运算次数有少许出入少许出入，某些，某些因子不能按照一般的复

4、数来计算运算量因子不能按照一般的复数来计算运算量,如如:n当当N很大时，这种特例很小。很大时，这种特例很小。 n直接计算直接计算DFT运算量是很可观的运算量是很可观的:nN=8时，时，DFT需需64次复乘，次复乘，nN=1024时，时，DFT所需复乘为所需复乘为1 048 576次，次，n对对实时性很强的信号处理实时性很强的信号处理来说，对来说，对硬件的计算硬件的计算速度要求是太高速度要求是太高了。为了实用的需要，了。为了实用的需要，需要改需要改进进DFT的计算方法，减少运算次数。的计算方法，减少运算次数。1j10NW12/NNWjWNN4/6如何才能减少运算量？仔细观察如何才能减少运算量？仔

5、细观察DFTDFT的运算就可看的运算就可看出，其出，其变换系数变换系数 的具有如下的具有如下特性特性：（1 1） 的的对称对称性性（2 2） 的的周期周期性性 = = =（3 3） 的的可约可约性性 由此可得由此可得: :knNWknNWnkNknNWW*knNWknNWkNnNW)( )(NknNWknNWkmnmNknNWWmnkmNknNWW/nkNknNNkNnNWWW)()(/21NNW kNNkNWW)2/(变换基底的特性7n（1）利用）利用 的的对称性对称性，合并合并DFT运算中某些项，并且可以计算反变换；运算中某些项，并且可以计算反变换；n（2）由于）由于DFT的运算量是与的运

6、算量是与 成正成正比，利用比，利用 周期性周期性和和可约性可约性将长序将长序的的DFT分解为短序列分解为短序列的的DFT。nFFT基本可以分成基本可以分成两大类两大类n按时间抽取按时间抽取(Decimation-in-time,DIT）n按频率抽取按频率抽取(Decimation-in-frequency,DIF）。）。knNW2NknNWFFT的基本的基本出发点出发点85.2.1 戈泽尔算法戈泽尔算法戈泽尔算法是一个利用复变换基W的周期性减少运算量的典型例子可以用于单离散频点或少数任意离散频点场合（如DTMF辨识）的有效计算910)(1010)()()()(NrrNkNNrkrNkNNNrk

7、rNWrxWrxWWrxkX10)()()()()(*)()(NrrnkNnxknNkrnuWrxnuWnxny为有限长10)(10)()()()()(| )(NrrNkNNrrNkNNnkkXWrxrNuWrxny变形10推论：n如果构造冲击响应h(n)为 的系统，则该系统对有限长输入在 时刻的响应即为 。n 下图给出了该系统的计算流图：)()(nuWnhknNNn )(kX11n由上图可以看出计算一个X(k)需要N次复乘和N次复加，即4N次实乘和4N次实加。n比直接计算DFT的运算量稍大一些，但是却避免了计算或存储复系数 。n另外，立足于上面用卷积实现卷积实现DFT的方法，我们是有可能通过

8、某种变形将上面的运算量减小一半。 knNW12)1 (*)/2cos(211)1 (*)1)(1 (111)(1211111zWzzNkzWzWzWzWzHkNkNkNkNkNkn注意到系统函数：13利用迭代实现了复因子的计算如果输入序列是复数，由于系数是实数并且-1不必作乘法运算，所以计算一个X(K)只需要2N次实数乘法和4N次实数加法。戈泽尔算法还有另外一个优点是只需要将前馈环节中的复因子取共轭就可以计算X(N-K)，可以将运算量再次减半。戈泽尔算法虽然比直接运算有效的多，但是无论如何变形，其运算量将仍然正比于N2。戈泽尔算法特点14n一、算法原理一、算法原理n 假设序列点数为假设序列点数

9、为 ,L为整数。如果为整数。如果不满足，可以人为加上若干零值点，使不满足，可以人为加上若干零值点，使之达到这一要求。这种之达到这一要求。这种N为为2的整数幂的的整数幂的FFT也称基也称基-2 FFT。n 将将 的序列按的序列按n的奇偶的奇偶分成两组：分成两组： n n r=0,1,N/2-1LN2LN2)() 12()()2(21rxrxrxrx5.2.2 DIT的的2MFFT(库利库利-图基算法图基算法)151202212021120)12(1202101010)()() 12()2()()()()(NrrkNkNNrrkNNrkrNNrrkNnNnnkNnNnnkNNnnkNWrxWWrx

10、WrxWrxWnxWnxWnxkX为奇数为偶数16n利用利用 的的可约性可约性：n式中式中 与与 均为均为N/2点点DFT：knNW222222NNjNjNWeeW)()()()()(212011202/22/1kXWkXWrxWWrxkXkNNrNrrkNkNrkN kX1 kX21202/1202/221202/1202/11) 12()()()2()()(NrrkNNrrkNNrrkNNrrkNWrxWrxkXWrxWrxkX一分为二一分为二17n但但 以及以及 ， 都是都是N/2点的点的序列，即序列，即n而而 却有却有N点，因而点，因而计算全部计算全部的的 ，必须应用必须应用DFT隐含

11、的周期性扩展：隐含的周期性扩展： rxrx21, kX1 kX212, 1 , 0,Nkr kX kX)(2, )(22211kXkNXkXkNX定义的扩展18n同时考虑到同时考虑到 的以下的以下性质性质n n这样就可将这样就可将 表达为前后两部分：表达为前后两部分：n kNWkNkNNNkNNWWWW2/2)(kX1221212( )( )( )0,1,12( )()( )( )0,1,1222kNNrrNNNX kX kW XkkNNNX kX rWXrX rW Xrr1920 上面的运算可以抽象为下面的蝶形信号流图单元：上面的运算可以抽象为下面的蝶形信号流图单元： 可以看出，每个蝶形运算

12、需要可以看出，每个蝶形运算需要一次复数乘法一次复数乘法 及及两次复数加两次复数加（减）法。（减）法。)(1kX)(2kXkNW1)()()(21kXWkXkXkN)()()2(21kXWkXNkXkNkNWkX)(221nN/2点点DFT复乘复乘:n 复加复加:n蝶形运算蝶形运算 复乘复乘:n 复加复加:n总运算量总运算量 复乘复乘:n 复加复加:n直接直接N点点DFT复乘复乘: N2 复加复加:N(N-1)22222NN12122*2NNNN2/NNN2*22212222NNNNN2122NNNN运算量-减少一半22LN2n既然如此，由于既然如此，由于 ，因此，因此N/2仍是仍是偶数，可以进

13、一步偶数，可以进一步按时间按时间的的奇偶性分解奇偶性分解，并且并且一直进行下去一直进行下去。n由于这种方法每一步都是按输入序列由于这种方法每一步都是按输入序列时时间的奇、偶性间的奇、偶性分解为两个更短的子序列，分解为两个更短的子序列，所以称为所以称为“按时间抽选法按时间抽选法”（DIT）。）。n分解过程如下分解过程如下：231、原位运算，2、倒位序，3、叠型间距24n由按时间抽选法由按时间抽选法FFT的流图可见，当的流图可见，当 时，共有：时，共有：nL级蝶形，每级都由级蝶形，每级都由N/2个蝶形个蝶形运算组成运算组成n每个蝶形有一次复乘、二次复加，因而每级每个蝶形有一次复乘、二次复加，因而每

14、级运算都需运算都需N/2次复乘和次复乘和N次复加次复加nL级级FFT运算总共需要：运算总共需要：n 复乘数：复乘数： (DFT: )n 复加数：复加数： (DFT: )LN22log22FNNmLN2NNNNLaF2log) 1(NN二、运算量二、运算量25n由于计算机上乘法运算所需时间比加法由于计算机上乘法运算所需时间比加法运算所需时间多得多，所以乘法是主要运算所需时间多得多，所以乘法是主要的运算量。的运算量。n直接计算直接计算DFT与与FFT算法计算量之比为算法计算量之比为:NNNNNLNN2222log2log222627)(1kXm)(1jXmrNW1)()()(11jXWkXkXmr

15、Nmm)()()(11jXWkXjXmrNmm三、三、DIT算法的特点算法的特点n1、原位运算（同址运算）原位运算（同址运算）n每级（每列）都由每级（每列）都由N/2个蝶形运算构成个蝶形运算构成n每一个蝶形结构完成下述基本迭代运算：每一个蝶形结构完成下述基本迭代运算： nm表示第表示第m级迭代，级迭代， k、j为数据所在行为数据所在行28n按原位计算时，按原位计算时，FFT的的输出输出 是按是按正常顺正常顺序排列序排列在存储单元中在存储单元中:n X（0），），X（1），），X（7），n但是但是输入输入 却不是按自然顺序存储的却不是按自然顺序存储的:n看起来好像是看起来好像是“混乱无序混乱无序

16、”，实际是有规律，实际是有规律的，称之为的，称之为倒位序倒位序。n造成倒位序的原因是输入造成倒位序的原因是输入 按标号按标号n的的奇、奇、偶性的不断分组偶性的不断分组造成。造成。n倒位序的实现倒位序的实现)(kX)(nx)7(,),4(),0(xxx)(nx012nnn210nnn)3(x)011(x)110(x2、倒位序规律、倒位序规律2912m3、蝶形运算两结点的、蝶形运算两结点的“距离距离”n由由8点点FFT的信号流图可以看出，其输入的信号流图可以看出，其输入是倒位序的，输出是自然顺序。是倒位序的，输出是自然顺序。n第一级每个蝶形的两节点间第一级每个蝶形的两节点间“距离距离”为为1，n第

17、二级每个蝶形的两节点第二级每个蝶形的两节点“距离距离”为为2，n第三级每个蝶形的两节点第三级每个蝶形的两节点“距离距离”为为4，n由此类推得，对由此类推得，对N2L ，当输入为倒位当输入为倒位序，输出为正常顺序时，其序，输出为正常顺序时，其第第m级运算，级运算，每个蝶形的两节点每个蝶形的两节点“距离距离”为为 。30四、四、DIT算法的其他形式流图算法的其他形式流图n对于对于任何流图，只要保持各节点所连任何流图，只要保持各节点所连的各支路及其传输系数不变的各支路及其传输系数不变，则不论，则不论节点位置怎么排列所得流图总是等效节点位置怎么排列所得流图总是等效n所得最后结果都是所得最后结果都是DF

18、T的正确结果，的正确结果，只是只是数据的提取和存放的次序不同而数据的提取和存放的次序不同而已已。3132335.2.3 DIF的的2M FFT(桑德桑德-图基算法图基算法)n频率抽选（频率抽选（DIF）的的FFT算法，是把频域算法，是把频域序列（也是序列（也是N点）按点）按K的奇偶性分解为越的奇偶性分解为越来越短的序列。来越短的序列。 n一、算法原理一、算法原理n假设序列点数为假设序列点数为 ，L为整数。为整数。n在将在将 按按k的奇偶分组之前，先把输入的奇偶分组之前，先把输入按按n的顺序分成前后两半的顺序分成前后两半LN2 kX34nkNNkNNnNnknNNNnnkNNNnnkNNnnkN

19、NnnkNWWNnxnxWNnxWnxWnxWnxWnxkX)2()()2()()()()()(2/120120)2(120121201012/NNW35n由于由于 ，故，故 ，可得：，可得：n n k=0,1,N-1n当当k为偶数时，为偶数时， ；n k为奇数时，为奇数时， 。n因此，按因此，按 k的奇偶可将的奇偶可将 分为两部分。分为两部分。 12/NNWkNkNW) 1(2/nkNkNnWNnxnxkX).2() 1()()(1201) 1(k1) 1(k)(kXnkNNkNNnWWNnxnxkX)2()()(2/12036n令令 n n则：则：120,1, ,122NrrkrknrNn

20、NNnrnNNnnrNNnrnNNnWWNnxnxWNnxnxrXWNnxnxWNnxnxrX2/120)12(1202/1202120).2()().2()()12().2()().2()()2(nkNkNnWNnxnxkX).2() 1()()(12037n令令 n则：则：12,.,1 , 0)2()()()2()()(21NnWNnxnxnxNnxnxnxnN12,.,1 , 0)() 12()()2(1202/21202/1NrWnxrXWnxrXNnnrNNnnrN38)(nx)2(NnxnNW1)2()()(1NnxnxnxnNWNnxnxnx)2()()(2DIF运算关系的基本蝶

21、形运算关系的基本蝶形39N=8点点DIF FFT结构结构40二、二、DIT与与DIF的异同的异同n形式上的区别是：形式上的区别是：DIF输入是自然顺序，输入是自然顺序，输出是倒位序的，与输出是倒位序的，与DIT正好相反。正好相反。n但这不是实质性的区别，因为但这不是实质性的区别，因为DIF与与DIT都可将输入或输出按照要求进行重排。都可将输入或输出按照要求进行重排。n实质性的区别是：实质性的区别是：nDIF的基本蝶形与的基本蝶形与DIT的基本蝶形不同。的基本蝶形不同。414243作业：n9.6n9.7 n9.14n9.17n9.26n9.27五、chirp Z变换n对对非单位圆上的抽样感非单位

22、圆上的抽样感兴趣兴趣n语音信号处理中往往需要知道极语音信号处理中往往需要知道极点所在处的复频率，如果极点位点所在处的复频率，如果极点位置离单位圆较远，只利用单位圆置离单位圆较远，只利用单位圆上的频谱，就很难知道极点所在上的频谱，就很难知道极点所在处的复频率。处的复频率。nz变换采用变换采用螺线抽样螺线抽样就适就适应于这些需要，称为线应于这些需要，称为线性调频性调频z变换变换(CZT)4445n已知已知 是有限长序列，其是有限长序列，其z变换为：变换为：n为适应为适应z可以沿可以沿Z平面更一般的路径取值，沿平面更一般的路径取值，沿z平面上的一段平面上的一段螺线螺线作作等分角等分角的的抽样抽样，z

23、的这的这些抽样点些抽样点zk为为: 10Nnnx nNnznxzX101, 1 , 0,MkAWzkk一、算法的基本原理一、算法的基本原理46n M为所要分析的复频谱的点数，不一为所要分析的复频谱的点数，不一定等于定等于N，A和和W都是任意复数都是任意复数:00jeAA 00jeWW00000000kjkjkkjkeWAeWeAz1, 1 , 0,MkAWzkk47n当当 各各zk就均匀等间隔地就均匀等间隔地分布在单位圆上，这时就是求序列的分布在单位圆上，这时就是求序列的DFT。NjeWANM2, 1,48nCZT的快速算法的快速算法n直接计算这一公式，与直接计算直接计算这一公式，与直接计算D

24、FT相相似，总共算出似，总共算出M个抽样点，需要个抽样点，需要NM次复次复数乘法数乘法与（与（N-1）M次复数加法，当次复数加法，当N，M较大时，运算量将很大。较大时，运算量将很大。 10,1010MkWAnxznxzXnknNnnkNnk49n采用采用布鲁斯坦布鲁斯坦（Bluestein）提出的等式，）提出的等式，可以将以上运算可以将以上运算转换为卷积和转换为卷积和形式，进而形式，进而采用采用FFT算法，提高运算速度。算法，提高运算速度。n布鲁斯坦所提出的等式为：布鲁斯坦所提出的等式为：n 由此可得：由此可得：22221nkknnk 2)(102222210222222nkNnnnkknknnNnkWWAnxWWWWAnxzX50 1, 1 , 0,22NnWAnxngnn 1, 1 , 0,1022MknkhngWzXNnkk 22nWnh 2)(102222210222222nkNnnnkknknnNnkWWAnxWWWWAnxzX51nCZT变换可以通过求变换可以通过求g(n)与与h(n)的线性的线性卷积，然后乘上卷积，然后乘上1/h(n)而得到，即：而得到，即：nh(n)为角为角频率频率随时间随时间线性增长线性增长的复指数序的复指数序列，称为线性调频信号（列，称为线性调频信号（Chirp signal），），因此，称为线调频因此，称为线调频z变换。变