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變換。變