數(shù)字信號(hào)處理 吳鎮(zhèn)揚(yáng)-CHP3快速傅立葉變換

第3章

傅立葉變換及其快速算法

之

快速傅里葉變換概述

快速傅里葉變換(FFT)是求解離散傅里葉變換(DFT)的快速算法。問(wèn)題的提出

直接計(jì)算N點(diǎn)DFT需要的計(jì)算量是多少？

計(jì)算一個(gè)X(k)需要N次復(fù)數(shù)乘法和N一1次復(fù)數(shù)加法。算出全部N點(diǎn)X(k)共需N2次復(fù)數(shù)乘法和N(N一1)次復(fù)數(shù)加法.

總運(yùn)算量近似地正比于N2

。當(dāng)N值很大（如2-D圖像處理），運(yùn)算量將非常龐大；同時(shí)，對(duì)于實(shí)時(shí)性很強(qiáng)的信號(hào)處理來(lái)說(shuō)，必將對(duì)計(jì)算速度有十分苛刻的要求。為此，需要改進(jìn)對(duì)DFT的計(jì)算方法，以減少總的運(yùn)算次數(shù)。概述在正交矩陣中，雖然有N2個(gè)元素，但只有N個(gè)不同的值，且有些取值特別簡(jiǎn)單，主要由于旋轉(zhuǎn)因子具有如下的特點(diǎn)：對(duì)稱性周期性下面以四點(diǎn)DFT為例來(lái)說(shuō)明快速算法的思路。如何充分利用這些關(guān)系？概述概述交換矩陣第二列和第三列得從上面的結(jié)果可以看出,利用對(duì)稱性和周期性，求四點(diǎn)DFT只需要一次復(fù)數(shù)乘法，稱為Coolkey-Tukey算法。概述算法分類：N為2的整次冪：按基數(shù)分為基-2FFT算法、基-4FFT算法、混合基FFT算法、分裂基FFT算法;當(dāng)N不是2的整次冪:典型的有Winograd算法.按抽取方法分：時(shí)間抽取(Decimation－in－Time，簡(jiǎn)稱DIT)；頻率抽取(Decimation－in－Frequency，簡(jiǎn)稱DIF)

概述3.3.1時(shí)間抽?。―IT）基2FFT算法

為了將大點(diǎn)數(shù)的DFT分解為小點(diǎn)數(shù)的DFT運(yùn)算，要求序列的長(zhǎng)度N為N＝2M(M為正整數(shù))。該情況下的變換稱為基2FFT。核心思想是N點(diǎn)DFTN/2點(diǎn)

DFTN/4點(diǎn)

DFT2點(diǎn)

DFT

1個(gè)2個(gè)4個(gè)N/2個(gè)問(wèn)題是如何分最有效？可以對(duì)時(shí)間變量分（DIT)，也可對(duì)頻率變量分(DIF）3.3.1時(shí)間抽?。―IT）基2FFT算法基本思路：從時(shí)域?qū)點(diǎn)序列x(n)按奇偶項(xiàng)分解為兩組，分別計(jì)算兩組N/2點(diǎn)DFT，然后再合成一個(gè)N點(diǎn)DFT，按此方法繼續(xù)下去，直到2點(diǎn)DFT，從而減少運(yùn)算量。算法具體步驟：一、算法的推導(dǎo)1、序列x(n)按奇偶項(xiàng)分解為兩組，將DFT運(yùn)算也相應(yīng)分為兩組則3.3.1時(shí)間抽取（DIT）基2FFT算法2、兩個(gè)N/2點(diǎn)的DFT合成一個(gè)N點(diǎn)DFT問(wèn)題：A(k)，B(k)都只有N/2個(gè)點(diǎn)，怎樣得到X(k)的后N/2點(diǎn)？利用周期性和對(duì)稱性得4.2時(shí)間抽?。―IT）基2FFT算法3.3.13.3.1時(shí)間抽?。―IT）基2FFT算法4.2時(shí)間抽?。―IT）基2FFT算法3、繼續(xù)分解（一直分解到兩點(diǎn)DFT變換）A(K)和B(K)仍是高復(fù)合數(shù)(N／2)的DFT，我們可按上述方法繼續(xù)以分解。令r＝2l，r＝2l十1，l＝0，1，…，N／4-1，則A(K)和B(K)可分別表示為3.3.1時(shí)間抽?。―IT）基2FFT算法4.2時(shí)間抽?。―IT）基2FFT算法3.3.1時(shí)間抽?。―IT）基2FFT算法令則4.2時(shí)間抽取（DIT）基2FFT算法3.3.1時(shí)間抽?。―IT）基2FFT算法同理，令則按此方法一直分解下去直到2點(diǎn)DFT，當(dāng)N=8時(shí)，如下：4.2時(shí)間抽?。―IT）基2FFT算法3.3.1時(shí)間抽取（DIT）基2FFT算法3.3.1時(shí)間抽?。―IT）基2FFT算法下面通過(guò)討論尋找FFT的一般規(guī)律。二、算法的討論1、“級(jí)”的概念在分解過(guò)程中，每分一次，稱為一級(jí)運(yùn)算。因?yàn)镸=log2N，所以N點(diǎn)DFT可以分解為M級(jí)，按抽取算法的信號(hào)流圖中來(lái)定義，從左到右分別稱為0級(jí)、1級(jí)到M-1級(jí)。3.3.1時(shí)間抽?。―IT）基2FFT算法2、蝶形單元在算法的信號(hào)流圖中，第m級(jí)存在這種運(yùn)算，這種結(jié)構(gòu)幾何形狀像蝴蝶，稱為蝶形單元p、q是參于本蝶形單元運(yùn)算的上、下節(jié)點(diǎn)的序號(hào)。由于第m級(jí)序號(hào)的兩點(diǎn)只參于這一個(gè)蝶形單元的運(yùn)算，其輸出在第m十l級(jí)。且這一蝶形單元也不再涉及別的點(diǎn)。由于這一特點(diǎn)，在計(jì)算機(jī)編程時(shí)，我們可將蝶形單元的輸出仍放在輸入數(shù)組中，這一特點(diǎn)稱為“同址運(yùn)算”。4.2時(shí)間抽取（DIT）基2FFT算法3.3.1時(shí)間抽?。―IT）基2FFT算法3.3.1時(shí)間抽?。―IT）基2FFT算法

由于一級(jí)都含有N/2個(gè)蝶形單元，每個(gè)蝶形單元需要1次復(fù)數(shù)乘法和兩次復(fù)數(shù)加法，因此完成log2N級(jí)共需要的復(fù)數(shù)乘法和加法分別為

直接計(jì)算DFT時(shí)所需的復(fù)乘數(shù)與復(fù)加數(shù)都是與N2成正比的。所以采用FFT算法使運(yùn)算量大大減少。顯然，N值愈大，節(jié)省的運(yùn)算量愈多。3.3.1時(shí)間抽取（DIT）基2FFT算法3、“組”的概念在分解過(guò)程中，每一級(jí)的N/2個(gè)蝶形單元可以分成若干組，每一組具有相同的結(jié)構(gòu)和W因子分布。第m級(jí)可分成N/2m+1組。例：N=8=23，分3級(jí)。第一級(jí)的分組及Wr因子如下：m=0級(jí),分成四組：因子為m=1級(jí),分成二組,因子為m=2級(jí),分成一組,因子為3.3.1時(shí)間抽取（DIT）基2FFT算法4、Wr因子的分布由上分析可知結(jié)論：每由后向前（m由M-1-->0級(jí)）推進(jìn)一級(jí)，則此系數(shù)為后級(jí)系數(shù)中偶數(shù)序號(hào)的那一半。3.3.1時(shí)間抽取（DIT）基2FFT算法5、碼位倒置在FFT算法中，輸出的頻譜依照正常次序排列，但輸入的序列x(n)是按奇偶分開的，分開的規(guī)律，以N=8為例，按如下方法進(jìn)行排序（1）、將x(n)的序號(hào)寫成二進(jìn)制

x(000),x(001),…,x(110)

，x(111)。（2）將二進(jìn)制的碼進(jìn)行翻轉(zhuǎn)，得

x(000),x(100),…,x(011)

，

x(111)。（3）將二進(jìn)制的翻轉(zhuǎn)碼轉(zhuǎn)換為對(duì)應(yīng)的十進(jìn)制

x(0),x(4),…,x(3)，x(7)。這就是按奇偶抽取得到的順序。3.3.1時(shí)間抽?。―IT）基2FFT算法說(shuō)明：①在上述的基2FFT算法中，由于每一步分解都是按輸入序列x(n)在時(shí)域上的次序是屬于偶數(shù)還是奇數(shù)來(lái)抽取的，所以稱為“按時(shí)間抽取法”或“時(shí)間抽取”。

②上述的基2FFT算法中，抽取也可在頻域進(jìn)行，引出頻率抽取（DIF）基2FFT算法。3.3.2頻率抽?。―IF）基2FFT算法

設(shè)輸入序列長(zhǎng)度為N=2M(M為正整數(shù))，頻率抽取法將輸入序列不是按奇、偶分組，而是按前后對(duì)半分開，這樣可將N點(diǎn)DFT寫成前后兩部分;將該序列的頻域的輸出序列X(k)(也是N點(diǎn)序列，按其頻域順序的奇偶分解為越來(lái)越短的子序列，稱為基2按頻率抽取的FFT算法。也稱為Sander-Tukey算法。3.3.2頻率抽?。―IF）基2FFT算法算法分析

現(xiàn)將輸入x(n)按n的順序分前后兩部分:前半子序列x(n),0≤n≤N/2-1;后半子序列x(n+N/2),0≤n≤N/2-1;例：N=8時(shí)，前半序列為：x(0),x(1),x(2),x(3);后半序列為：x(4),x(5),x(6),x(7);考慮N點(diǎn)的DFT,由DFT定義得3.3.2頻率抽?。―IF）基2FFT算法3.3.2頻率抽取（DIF）基2FFT算法按k的奇偶將X(k)分成奇偶兩部分,k=2r和k=2r+1,考慮k為偶數(shù)情況令3.3.2頻率抽?。―IF）基2FFT算法考慮k為奇數(shù)情況令3.3.2頻率抽取（DIF）基2FFT算法結(jié)論

一個(gè)N點(diǎn)的DFT被分解為兩個(gè)N/2點(diǎn);與時(shí)間抽取法的推演過(guò)程一樣，由于N=2M,因此,N/2仍為偶數(shù)，所以可以將N/2點(diǎn)DFT的輸出X(k)再分為偶數(shù)組和奇數(shù)組，這樣就將一個(gè)N/2點(diǎn)的DFT分成兩個(gè)N/4點(diǎn)DFT的輸入，也是將N/2點(diǎn)的DFT的輸入上、下對(duì)半分后通過(guò)蝶形運(yùn)算而形成，直至最后為2點(diǎn)DFT。8點(diǎn)DIF基2FFT算法流圖3.3.2頻率抽?。―IF）基2FFT算法3.3.2頻率抽取（DIF）基2FFT算法3.3.2頻率抽?。―IF）基2FFT算法DIT與DIF的相同之處：（1）DIF與DIT兩種算法均為原位運(yùn)算。（2）DIF與DIT運(yùn)算量相同。DIT與DIF的不同之處：（1）DIF與DIT兩種算法結(jié)構(gòu)倒過(guò)來(lái)。DIF為輸入順序，輸出亂序。運(yùn)算完畢再運(yùn)行“二進(jìn)制倒讀”程序。DIT為輸入亂序，輸出順序。先運(yùn)行“二進(jìn)制倒讀”程序，再進(jìn)行求DFT。（2）DIF與DIT根本區(qū)別：在于蝶形結(jié)不同。DIT的復(fù)數(shù)相乘出現(xiàn)在減法之前。DIF的復(fù)數(shù)相乘出現(xiàn)在減法之后。3.3.3N為組合數(shù)的FFT

以上討論的都是以2為基數(shù)的FFT算法，即N=2M，這種情況實(shí)際上使用得最多。

優(yōu)點(diǎn)：程序簡(jiǎn)單，效率高，使用方便。實(shí)際應(yīng)用時(shí)，有限長(zhǎng)序列的長(zhǎng)度N很大程度上由人為因素確定，因此多數(shù)場(chǎng)合可取N=2M，從而直接使用以2為基數(shù)的FFT算法。如N不能人為確定,N的數(shù)值也不是以2為基數(shù)的整數(shù)次方,處理方法有兩種:①補(bǔ)零:將x(n)補(bǔ)零,使N=2M.

例如N=30,補(bǔ)上x(30)=x(31)=0兩點(diǎn),使N=32=25,這樣可直接采用以2為基數(shù)M=5的FFT程序。有限長(zhǎng)度序列補(bǔ)零后并不影響其頻譜X(ej)，只是頻譜的采樣點(diǎn)數(shù)增加了,上例中由30點(diǎn)增加到32點(diǎn)，所以在許多場(chǎng)合這種處理是可接受的。3.3.3N為組合數(shù)的FFT②如要求準(zhǔn)確的N點(diǎn)DFT值,可采用任意數(shù)為基數(shù)的FFT算法,

其計(jì)算效率低于以2為基數(shù)FFT算法。如N為復(fù)合數(shù)，可分解為兩個(gè)整數(shù)p與q的乘積,像前面以2為基數(shù)時(shí)一樣,FFT的基本思想是將DFT的運(yùn)算盡量分小,因此,在N=pq情況下,也希望將N點(diǎn)的DFT分解為p個(gè)q點(diǎn)DFT或q個(gè)p點(diǎn)DFT,以減少計(jì)算量。步驟：分別為0,1,…，Q-1;分別為0,1,…，P-1。3.3.3N為組合數(shù)的FFT

N點(diǎn)DFT可以重新寫成為考慮到

再令令(1)先將x(n)通過(guò)x(n1Q+n0)改寫成x(n1,n0)。因?yàn)?/p>

Q=4,n1=0,1,2,n0=0,1,2,3,故輸入是按自然順序的，即x(0,0)=x(0)x(0,1)=x(1)x(0,2)=x(2)x(0,3)=x(3)x(1,0)=x(4)x(1,1)=x(5)x(1,2)=x(6)x(1,3)=x(7)x(2,0)=x(8)x(2,1)=x(9)x(2,2)=x(10)x(2,3)=x(11)以P=3，Q=4,N=12為例

(2)求Ｑ個(gè)Ｐ點(diǎn)的DFT（３）X1(k0,n0)乘以得到X1′(k0,n0)。（４）求P個(gè)Ｑ點(diǎn)的DFT，參變量是k0（５）將X2(k0,k1)通過(guò)X(k0+k1P)恢復(fù)為X(k)3.3.3N為組合數(shù)的FFTN=12為組合數(shù)時(shí)的FFT3.3.3N為組合數(shù)的FFT（１）求Ｑ個(gè)P點(diǎn)DFT需要ＱP2次復(fù)數(shù)乘法和Ｑ·P·(P-1)次復(fù)數(shù)加法；（２）乘Ｎ個(gè)W因子需要Ｎ次復(fù)數(shù)乘法；（３）求P個(gè)Ｑ點(diǎn)DFT需要PQ2次復(fù)數(shù)乘法和P·Ｑ（Ｑ-1)次復(fù)數(shù)加法。

總的復(fù)數(shù)乘法量:QP2+N+PQ2=N(P+Q+1)；總的復(fù)數(shù)加法量:Q·P(P-1)+P·Q·(Q-1)=N(P+Q-2)例:N=23*29=667,N2=444889,N(P+Q+1)=35351當(dāng)組合數(shù)

N=P1P2P3…Pm中所有的Pi均為4時(shí)，就是基四FFT算法

以上提出FFT算法，可以很快地求出全部DFT值。即求出有限長(zhǎng)序列x(n)的z變換X(z)在單位園上N個(gè)等間隔抽樣點(diǎn)zk處的抽樣值。它要求N為高度復(fù)合數(shù)。即N可以分解成一些因子的乘積。例N=2L

實(shí)際上：（1）也許對(duì)其它圍線上z變換取樣發(fā)生興趣。（2）只需要計(jì)算單位圓上某一段的頻譜,即M不等于N。如窄帶信號(hào)，希望在窄帶頻率內(nèi)頻率抽樣能夠非常密集，提高分辨率，帶外則不考慮。（3）若N是大素?cái)?shù)時(shí)，不能加以分解，又如何有效計(jì)算這種序列DFT。例N=311，若用基2則須補(bǔ)N=28=512點(diǎn)，要補(bǔ)211個(gè)零點(diǎn)。3.3.4線性調(diào)頻Z變換3.3.4線性調(diào)頻Z變換問(wèn)題提出為了提高DFT的靈活性，須用新的方法。線性調(diào)頻z變換(CZT)就是適用這種更為一般情況下，由x(n)求X(zk)的快速變換。CZT

來(lái)自于雷達(dá)專業(yè)的專用詞匯。Z變換采用螺線抽樣,可計(jì)算單位圓上任一段曲線的Z變換，適用于更一般情況下（M不等于N）由x(n)求X(zr)的快速算法,達(dá)到頻域細(xì)化的目的，這種變換稱為線性調(diào)頻Z變換(簡(jiǎn)稱CZT)。

為適應(yīng)z可以沿平面內(nèi)更一般的路徑取值，我們沿z平面上的一段螺線作等分角的抽樣，則z的取樣點(diǎn)Zr可表示為：

已知N點(diǎn)序列x(n)，0≤n≤N-1,其z變換為其中M：表示欲分析的復(fù)頻譜的點(diǎn)數(shù)。M不一定等于N。A，W都為任意復(fù)數(shù),令

3.3.4線性調(diào)頻Z變換一、CZT的定義3.3.4線性調(diào)頻Z變換上式即為CZT的定義.現(xiàn)在討論A0,W0,θ0,φ0的含義:為輸出M點(diǎn)的變換域值.r=0時(shí)的A0ejθ0是CZT的起點(diǎn),隨著r的變化,r0,r1,…,rM-1構(gòu)成CZT的變化路徑，對(duì)于M-1點(diǎn)其極坐標(biāo)為3.3.4線性調(diào)頻Z變換3.3.4線性調(diào)頻Z變換CZT在現(xiàn)Z平面上的變換路徑是一條螺旋線。（1）A為起始樣點(diǎn)位置起點(diǎn)半徑，大于1時(shí)，表示螺旋線在單位圓外，反之，在單位圓內(nèi)。起點(diǎn)半相角。（2）當(dāng)W0>1,螺旋線內(nèi)旋，反之外旋。（3）當(dāng)A0=W0=1時(shí)，CZT的變換路徑為單位圓上的一段弧，起點(diǎn)為P，終點(diǎn)為Q，且M不一定等于N。3.3.4線性調(diào)頻Z變換3.3.4線性調(diào)頻Z變換

我們希望做的是頻譜分析，因此考慮A0=W0=1時(shí)，在單位圓上CZT，且M不一定等于N。3.3.4線性調(diào)頻Z變換3.3.4線性調(diào)頻Z變換CZT的線性濾波計(jì)算步驟3.3.4線性調(diào)頻Z變換二、CZT的計(jì)算方法分析：從上面的推導(dǎo)過(guò)程可以看出，計(jì)算CZT關(guān)鍵是計(jì)算一個(gè)線性卷積其中，g(n)應(yīng)為N點(diǎn)序列，h(n)應(yīng)為偶對(duì)稱的無(wú)限長(zhǎng)序列，考慮到輸出M點(diǎn)序列，h(n)的實(shí)際長(zhǎng)度應(yīng)為M點(diǎn)。因此，可用DFT來(lái)實(shí)現(xiàn)兩者的卷積，具體步驟如下：3.3.4線性調(diào)頻Z變換（1）計(jì)算并設(shè)置g(n)3.3.4線性調(diào)頻Z變換（2）計(jì)算并設(shè)置h(n)將h(n)設(shè)置成長(zhǎng)度為L(zhǎng)的序列，考慮到其為偶對(duì)稱序列，且取M點(diǎn)（輸出M點(diǎn)序列），取如下圖所示3.3.4線性調(diào)頻Z變換3.3.4線性調(diào)頻Z變換（3）計(jì)算h‘(n)和g’(n)的DFT，得到L點(diǎn)序列H‘(k)和G’(k)。（4）令Y‘(k)=H‘(k)G’(k)（乘積）后作Y‘(k)

的IDFT（反變換）得到時(shí)域輸出序列y(r)

。（5）取y(r)的前M點(diǎn)，并乘以W-r2/2,則得最后的輸出X(zr)，即

與標(biāo)準(zhǔn)FFT算法相比，CZT算法有以下特點(diǎn)：（1）輸入序列長(zhǎng)度N及輸出序列長(zhǎng)度M不需要相等，且N及M不必是高度合成數(shù)，二者均可為素?cái)?shù)。（2）Zk的角間隔是任意的，說(shuō)明其頻率分辨率也是任意可控的，角間隔小，分辨率高，反之，分辨率低。（3）周線不必是z平面上的圓，在語(yǔ)音分析中螺旋周線具有某些優(yōu)點(diǎn)。（4）由于起始點(diǎn)z0可任意選定，因此可以從任意頻率上開始對(duì)輸入數(shù)據(jù)進(jìn)行窄帶高分辨率的分析?？傊珻ZT算法具有很大的靈活性。3.3.4線性調(diào)頻Z變換2.4FFT應(yīng)用中的幾個(gè)問(wèn)題1)IDFT的運(yùn)算方法IDFT：DFT：以上所討論的FFT算法可用于IDFT運(yùn)算——簡(jiǎn)稱為IFFT比較IDFT的定義式：IDFT與DFT的差別：

1）把DFT中的每一個(gè)系數(shù)改為，

2）再乘以常數(shù)1/N，

第二種方法，完全不需要改動(dòng)FFT程序，而是直接利用它作IFFT?？紤]到故

IFFT計(jì)算分三步：①將X（k）取共軛（虛部乘以-1）

②對(duì)直接作FFT③對(duì)FFT的結(jié)果取共軛并乘以1/N，得x（n）。2）實(shí)數(shù)序列的FFT

以上討論的FFT算法都是復(fù)數(shù)運(yùn)算,包括序列x(n)也認(rèn)為是復(fù)數(shù),但大多數(shù)場(chǎng)合,信號(hào)是實(shí)數(shù)序列,任何實(shí)數(shù)都可看成虛部為零的復(fù)數(shù),例如,求某實(shí)信號(hào)x(n)的復(fù)譜,可認(rèn)為是將實(shí)信號(hào)加上數(shù)值為零的虛部變成復(fù)信號(hào)(x(n)+j0),再用FFT求其離散傅里葉變換。這種作法很不經(jīng)濟(jì),因?yàn)榘褜?shí)序列變成復(fù)序列,存儲(chǔ)器要增加一倍,且計(jì)算機(jī)運(yùn)行時(shí),即使虛部為零,也要進(jìn)行涉及虛部的運(yùn)算,浪費(fèi)了運(yùn)算量。合理的解決方法是利用復(fù)數(shù)據(jù)FFT對(duì)實(shí)數(shù)據(jù)進(jìn)行有效計(jì)算,下面介紹兩種方法。

（1）用

一個(gè)N點(diǎn)FFT同時(shí)計(jì)算兩個(gè)N點(diǎn)實(shí)序列的DFT

設(shè)x

(n)、y

(n)是彼此獨(dú)立的兩個(gè)N點(diǎn)實(shí)序列,且

X

(k)=DFT[x

(n)]，Y

(k)=DFT[y(n)]

則X

(k)、Y(k)可通過(guò)一次FFT運(yùn)算同時(shí)獲得。首先將x

(n)、y(n)分別當(dāng)作一復(fù)序列的實(shí)部及虛部,令

g(n)=x

(n)+jy(n)通過(guò)FFT運(yùn)算可獲得g(n)的DFT值利用離散傅里葉變換的共軛對(duì)稱性通過(guò)g(n)的FFT運(yùn)算結(jié)果G(k),由上式也可得到Y(jié)(k)的值。

（2）用一個(gè)N點(diǎn)的FFT運(yùn)算獲得一個(gè)2N點(diǎn)實(shí)序列的DFT

設(shè)x(n)是2N點(diǎn)的實(shí)序列,現(xiàn)人為地將x(n)分為偶數(shù)組x1

(n)和奇數(shù)組x2

(n)

x1

(n)=x(2n)n=0,1,…,N-1

x2

(n)=x(2n+1)

n=0,1,…,N-1然后將x1

(n)及x2

(n)組成一個(gè)復(fù)序列：

y(n)=x1

(n)+jx2

(n)

通過(guò)N點(diǎn)FFT運(yùn)算可得到：Y(k)=X1(k)+jX2(k)，N點(diǎn)根據(jù)前面的討論，得到

為求2N點(diǎn)x(n)所對(duì)應(yīng)X(k)，需求出X(k)與X1

(k)、X2

(k)的關(guān)系：

而1）由x1(n)及x2(n)組成復(fù)序列，經(jīng)FFT運(yùn)算求得Y(k)，2）利用共軛對(duì)稱性求出X1(k)、X2(k)，3）最后利用上式求出X（k）,達(dá)到用一個(gè)N點(diǎn)的FFT計(jì)算一個(gè)2N點(diǎn)的實(shí)序列的DFT的目的。

X(k)=X1(k)+W2NkX2(k)所以

3)線性卷積的FFT算法

線性卷積是求離散系統(tǒng)響應(yīng)的主要方法之一,許多重要應(yīng)用都建立在這一理論基礎(chǔ)上,如卷積濾波等。以前曾討論了用圓周卷積計(jì)算線性卷積的方法歸納如下:

將長(zhǎng)為N2的序列x(n)延長(zhǎng)到L,補(bǔ)L-N2個(gè)零，將長(zhǎng)為N1的序列h(n)延長(zhǎng)到L,補(bǔ)L-N1個(gè)零，如果L≥N1+N2-1,則圓周卷積與線性卷積相等,此時(shí),可用FFT計(jì)算線性卷積，方法如下:

a.計(jì)算X(k)=FFT[x(n)]b.求H(k)=FFT[h(n)]c.求Y(k)=H(k)X(k)k=0--L-1d.求y(n)=IFFT[Y(k)]n=0--L-1

可見(jiàn)，只要進(jìn)行二次FFT,一次IFFT就可完成線性卷積計(jì)算。計(jì)算表明,L>32時(shí),上述計(jì)算線性卷積的方法比直接計(jì)算線卷積有明顯的優(yōu)越性,因此,也稱上述循環(huán)卷積方法為快速卷積法。

上述結(jié)論適用于x(n)、h(n)兩序列長(zhǎng)度比較接近或相等的情況,如果x(n)、h(n)長(zhǎng)度相差較多,例如,h(n)為某濾波器的單位脈沖響應(yīng),長(zhǎng)度有限,用來(lái)處理一個(gè)很長(zhǎng)的輸入信號(hào)x(n)，或者處理一個(gè)連續(xù)不斷的信號(hào)，按上述方法,有三個(gè)問(wèn)題：（1）h(n)要補(bǔ)許多零再進(jìn)行計(jì)算,計(jì)算量有很大的浪費(fèi)，或者根本不能實(shí)現(xiàn)。（2）系統(tǒng)的存儲(chǔ)量要求極高。（3）帶來(lái)了很大的系統(tǒng)延遲。為了克服上述三個(gè)問(wèn)題，保持快速卷積法的優(yōu)越性,可將x(n)分為許多段,每段的長(zhǎng)度與h(n)接近,處理方法有兩種：(1)

重疊相加法——由分段卷積的各段相加構(gòu)成總的卷積輸出h(n)x(n)則輸入序列可表為：于是輸出可分解為：

其中假定xi(n)表示x（n）序列的第i段：

1）先對(duì)h（n）及xi（n）補(bǔ)零，補(bǔ)到具有N點(diǎn)長(zhǎng)度，N=N1+N2-1。一般選N=2M。由于yi（n）的長(zhǎng)度為N1，而xi（n）的長(zhǎng)度為N2，因此相鄰兩段yi（n）序列必然有N-N2=N1-1點(diǎn)發(fā)生重疊。2）用基2FFT計(jì)算yi（n）=xi（n）*h（n）。3）重疊部分相加構(gòu)成最后的輸出序列。計(jì)算步驟：a.

事先準(zhǔn)備好濾波器參數(shù)H（k）=DFT[h(n)]，N點(diǎn)b.

用N點(diǎn)FFT計(jì)算Xi（k）=DFT[xi(n)]c.

Yi(k)=Xi(k)H(k)d.

用N點(diǎn)IFFT求yi（n）=IDFT[Yi(k)]e.

將重疊部分相加(2)重疊保留

這種方法和第一種方法稍有不同，即將上面分段序列中補(bǔ)零的部分不是補(bǔ)零，而是保留原來(lái)的輸入序列值，這時(shí)，如利用DFT實(shí)現(xiàn)h(n)和xi(n)的循環(huán)卷積，則每段卷積結(jié)果中有N1-1個(gè)點(diǎn)不等于線性卷積值需舍去。

重疊保留法與重疊相加法的計(jì)算量差不多，但省去了重疊相加法最后的相加運(yùn)算。

y0(n)中的[N1-1,L-1]點(diǎn)對(duì)應(yīng)于線性卷積

x(n)*h(n)中的[0,N2-1]點(diǎn)y1(n)中的[N1-1,L-1]點(diǎn)對(duì)應(yīng)于線性卷積

x(n)*h(n)中的[N2,2N2-1]點(diǎn)y2(n)中的[N1-1,L-1]點(diǎn)對(duì)應(yīng)于線性卷積

x(n)*h(n)中的[2N2,3N2-1]點(diǎn)依此類推

，并將yi(n)拼接起來(lái)構(gòu)成y

(n)

4）用FFT計(jì)算相關(guān)函數(shù)

相關(guān)的概念很重要，互相關(guān)運(yùn)算廣泛應(yīng)用于信號(hào)分析與統(tǒng)計(jì)分析，如通過(guò)相關(guān)函數(shù)峰值的檢測(cè)測(cè)量?jī)蓚€(gè)信號(hào)的時(shí)延差等。兩個(gè)長(zhǎng)為N的實(shí)離散時(shí)間序列x（n）與y（n）的互相關(guān)函數(shù)定義為

則可以證明，rxy(τ)的離散傅里葉變換為

Rxy（k）=X*(k)Y(k)

其中X(k)=DFT[x(n)],Y(k)=DFT[y(n)],Rxy(k)=DFT[rxy(τ)]，0≤k≤N-1互相關(guān)函數(shù)定義為

x(n)及y(n)的卷積公式相比較，我們可以得到相關(guān)和卷積的時(shí)域關(guān)系：

因

得

證畢。

當(dāng)x（n）=y（n）時(shí)，得到x（n）的自相關(guān)函數(shù)為：

維納——辛欽定理：自相關(guān)函數(shù)與信號(hào)功率譜互為傅立葉變換對(duì)。

上面的推導(dǎo)表明，互相關(guān)和自相關(guān)函數(shù)的計(jì)算可利用FFT實(shí)現(xiàn)。由于離散傅里葉變換隱含著周期性，所以用FFT計(jì)算離散相關(guān)函數(shù)也是對(duì)周期序列而言的。直接做N點(diǎn)FFT相當(dāng)于對(duì)兩個(gè)N點(diǎn)序列x(n)、y(n)作周期延拓，作相關(guān)后再取主值（類似圓周卷積）。實(shí)際一般要求的是兩個(gè)有限長(zhǎng)序列的線性相關(guān)，為避免混淆，需采用與循環(huán)卷積求線性卷積相類似的方法，先將序列延長(zhǎng)補(bǔ)0后再用上述方法。利用FFT求兩個(gè)有限長(zhǎng)序列的線性相關(guān)：(5)

對(duì)R（k）作IFFT，取后N-1項(xiàng)，得取前N項(xiàng)，得（1）

設(shè)x(n)和y(n)的長(zhǎng)均為N，求線性相關(guān)；（2）為了使兩個(gè)有限長(zhǎng)序列的線性相關(guān)可用其循環(huán)相關(guān)代替而不產(chǎn)生混淆，選擇周期L≥2N-1，且L=2m，以使用FFT，將x(n)，y(n)補(bǔ)零至長(zhǎng)為L(zhǎng)。（3）

用FFT計(jì)算X(k)，Y(k)（k=0，1…，L-1）（4）

R（k）=X*（k）y（k）

x=[13-112331];

y=[21-1120-13];

k=length(x);

xk=fft(x,2*k);

yk=fft(y,2*k);

rm=real(ifft(conj(xk).*yk));

rm=[rm(k+2:2*k)rm(1:k)];

m=(-k+1):(k-1);

stem(m,rm)

xlabel('m');ylabel('幅度');

-8-6-4-202468-6-4-2024681012m幅度兩個(gè)序列的自相關(guān)函數(shù)例8例9-20-1001020-50050100150200250300m-20-1001020-50050100150200250300m幅度幅度

(a)(b)延遲序列的互相關(guān)函數(shù)(a)和自相關(guān)函數(shù)(b)

5）用FFT計(jì)算二維離散的傅里葉變換

二維信號(hào)有圖象信號(hào)、時(shí)空信號(hào)、時(shí)頻信號(hào)等。二維離散傅里葉變換可用于處理二維離散信號(hào)。二維離散傅里葉變換的定義為：

二維離散傅里葉變換可通過(guò)兩次一維離散傅里葉變換來(lái)實(shí)現(xiàn):1）作一維N點(diǎn)DFT（對(duì)每個(gè)m做一次，共M次）

k=0,1,…,N-1，m=0,1,…,M-12）作M點(diǎn)的DFT（對(duì)每個(gè)k做一次，共N次）

k=0,1,…,N-1，