《數(shù)字信號(hào)處理》-第4章

第四章離散傅里葉變換的計(jì)算與應(yīng)用

本章主要內(nèi)容

●離散傅里葉變換的高效計(jì)算思路

●兩種基-2FFT算法：（DIT、DIF）●離散傅里葉反變換（IDFT）的快速計(jì)算方法●幾種FFT算法●用DFT的快速算法（FFT）對(duì)信號(hào)進(jìn)行線性卷積及線性相關(guān)

●用DFT的快速算法（FFT）對(duì)信號(hào)進(jìn)行頻譜分析4.1離散傅里葉變換的高效計(jì)算思路

設(shè)x(n)為N點(diǎn)有限長(zhǎng)序列，其DFT正變換為（）反變換（IDFT）為（）可以認(rèn)為兩式的運(yùn)算量是相同的，因此，我們只討論正變換（）式的運(yùn)算量。直接進(jìn)行DFT運(yùn)算需要次復(fù)數(shù)乘法和次復(fù)數(shù)加法。也可以統(tǒng)計(jì)出用實(shí)數(shù)運(yùn)算完成DFT的運(yùn)算量。（）式可以寫(xiě)為（）（）式表明，一次復(fù)數(shù)乘法需用4次實(shí)數(shù)乘法和2次實(shí)數(shù)加法；一次復(fù)數(shù)加法則需要2次實(shí)數(shù)加法。因此，整個(gè)DFT運(yùn)算共需要4N2次實(shí)數(shù)乘法和(4N-2)N次實(shí)數(shù)加法。由以上分析可知，直接計(jì)算N點(diǎn)DFT的乘法和加法運(yùn)算次數(shù)均與N2成正比，當(dāng)很大時(shí)，運(yùn)算量是很大的的。減少運(yùn)算量的思路把N點(diǎn)DFT分解為幾個(gè)較短的DFT，可使乘法次數(shù)減少。利用系數(shù)的對(duì)稱(chēng)性、周期性和可約性，就可以減小DFT的運(yùn)算量。周期性:對(duì)稱(chēng)性:可約性：由此可以得出利用這些特性，使DFT運(yùn)算中有些項(xiàng)可以合并，并且可以將長(zhǎng)序列的DFT分解為幾個(gè)短序列的DFT，以減少DFT的運(yùn)算次數(shù)。

離散傅里葉變換的高效算法正是基于這樣的基本思路發(fā)展起來(lái)的。所有的這類(lèi)算法被人們稱(chēng)為快速傅里葉變換（FastFourierTransform，F(xiàn)FT）算法。4.2按時(shí)間抽?。―IT）的基-2FFT算法4.2.1算法原理

“按時(shí)間抽取法”（DIT）:逐次分解時(shí)間序列x(n)基-2FFT算法:序列長(zhǎng)度為N=2M，M為正整數(shù)，即N為2的整數(shù)冪的FFT算法將N=2M的序列x(n)按n的偶奇分為兩組，形成兩個(gè)子序列，得（）則由的可約性，有，所以上式可表示成

（）其中X1(k)和X2(k)分別是x1(r)及x2(r)的N/2點(diǎn)DFT。由（）式可知，若取k=0,1,…,N/2-1，用（）式計(jì)算X(k)時(shí)，則只能計(jì)算出X(k)前一半的值，

利用系數(shù)的性質(zhì)以及X1(k)和X2(k)的隱含周期性，可得

（）于是，X(k)可表示為前半部分：（）

后半部分：

（）

(4.2.6)式和(4.2.7)式的運(yùn)算可用下圖的蝶形運(yùn)算流圖符號(hào)（又可稱(chēng)為蝶形運(yùn)算單元）表示。

輸入→→輸出

輸入→→輸出

每個(gè)蝶形運(yùn)算單元需要一次復(fù)數(shù)乘法及兩次復(fù)數(shù)加（減）法。

分解過(guò)程圖示共需要復(fù)數(shù)乘和復(fù)數(shù)加的次數(shù)為2(N/2)2+N/2和N2/2每個(gè)N/2點(diǎn)子序列再按其排列的偶奇順序可進(jìn)一步分解為兩個(gè)N/4點(diǎn)的子序列：（）于是，可得x1(r)的DFT為式中由于及X3(k)和X4(k)的隱含周期性，可得到

同樣，將x2(r)按r取偶、奇可分解成2個(gè)長(zhǎng)N/4的子序列k=0,1,…,N/4-1；

下圖給出了N=8時(shí)，將一個(gè)N/2點(diǎn)DFT分解成兩個(gè)N/4點(diǎn)DFT，由這兩個(gè)N/4點(diǎn)DFT組合成一個(gè)N/2點(diǎn)DFT的流圖。下圖給出了將一個(gè)N=8點(diǎn)的DFT分解成四個(gè)N/4=2點(diǎn)的DFT流圖。兩級(jí)蝶形組合運(yùn)算，比只用一次分解蝶形方式的計(jì)算量又減少了大約一半N=8點(diǎn)DIT-FFT運(yùn)算流圖4.2.2DIT-FFT算法的運(yùn)算量

當(dāng)N=2M時(shí)，共有M級(jí)蝶形,每一級(jí)有N/2個(gè)碟形運(yùn)算，每個(gè)蝶形需要一次復(fù)數(shù)乘法和二次復(fù)數(shù)加法，M級(jí)蝶形運(yùn)算總共需要:復(fù)數(shù)乘法次數(shù)為：

復(fù)數(shù)加法次數(shù)為：

直接計(jì)算N點(diǎn)DFT需要復(fù)數(shù)乘法N2次復(fù)數(shù)加法N(N-1)次可見(jiàn)FFT大大減少了運(yùn)算次數(shù)，提高了運(yùn)算速度。

4.2.3按時(shí)間抽取的基-2FFT算法的運(yùn)算特點(diǎn)及編程思想

1蝶形運(yùn)算

DIT-FFT的運(yùn)算是由大量蝶形運(yùn)算單元組成的。它每一級(jí)（列）的計(jì)算都是由N/2個(gè)蝶形運(yùn)算單元完成的，N=2M點(diǎn)的FFT共有M級(jí)蝶形運(yùn)算。每個(gè)蝶形運(yùn)算單元完成的迭代運(yùn)算可表示為(4.2.19)式中L表示第L級(jí)迭代，1≤K≤M，k和j為數(shù)據(jù)所在的行數(shù)，系數(shù)又稱(chēng)為旋轉(zhuǎn)因子（TwiddleFactor）。（）式的蝶形運(yùn)算如圖所示，一個(gè)蝶形運(yùn)算包括一次復(fù)數(shù)乘法和兩次復(fù)數(shù)加法運(yùn)算。

由圖的運(yùn)算流圖可以看出對(duì)于N=2M點(diǎn)的DIT-FFT算法，其第L級(jí)每個(gè)蝶形運(yùn)算單元的兩節(jié)點(diǎn)間“距離”則為B=2L-1。于是，(4.2.19)式的蝶形運(yùn)算就可表示為

2原位運(yùn)算

定義：利用同一組存儲(chǔ)單元存儲(chǔ)蝶形運(yùn)算輸入、輸出數(shù)據(jù)的方法稱(chēng)為原位運(yùn)算，也可稱(chēng)為同址運(yùn)算。

碟形運(yùn)算流圖同一級(jí)中，每個(gè)蝶形的兩個(gè)輸入數(shù)據(jù)只對(duì)本蝶形有用，每個(gè)蝶形的輸入、輸出數(shù)據(jù)節(jié)點(diǎn)在用一條水平線上。這樣，當(dāng)計(jì)算完一個(gè)蝶形后，所得的輸出數(shù)據(jù)可立即存入原輸入數(shù)據(jù)所占用的存儲(chǔ)單元。經(jīng)過(guò)M級(jí)運(yùn)算后，原來(lái)存放輸入序列數(shù)據(jù)的N個(gè)存儲(chǔ)單元中可依次存放X(k)的N個(gè)值。特點(diǎn)：節(jié)省存儲(chǔ)單元，降低設(shè)備成本

3旋轉(zhuǎn)因子的變化規(guī)律

旋轉(zhuǎn)因子和運(yùn)算級(jí)數(shù)的關(guān)系

N=23=8時(shí)的各級(jí)旋轉(zhuǎn)因子：一般地，當(dāng)N=2M（M為正整數(shù)）時(shí)，第L級(jí)蝶形運(yùn)算共有2L-1種不同的旋轉(zhuǎn)因子，可表示為

這樣，就可以確定第L級(jí)蝶形運(yùn)算的旋轉(zhuǎn)因子（程序中，L為循環(huán)變量）。在第L級(jí)中，同一旋轉(zhuǎn)因子對(duì)應(yīng)著“間隔”為2L點(diǎn)的2M-L個(gè)蝶形。

4倒位序規(guī)律

倒位序：輸入序列經(jīng)過(guò)M級(jí)分解后，其排列順序變?yōu)閤(0)x(4)x(2)x(6)x(1)x(5)x(8)x(7)，服從所謂的“碼位倒置”的規(guī)律，稱(chēng)之為倒位序。造成倒位序的原因：將其按標(biāo)號(hào)的偶奇的不斷分組，每次分解總是將偶序列放在上面，把奇序列放在下面。

首先最低位按0、1分為偶、奇兩組，接著次低位也按0、1分組，依此類(lèi)推右圖為描述倒位序的樹(shù)狀圖(N=8)

5倒位序的實(shí)現(xiàn)對(duì)照表變址功能產(chǎn)生倒序數(shù)的十進(jìn)制運(yùn)算規(guī)律

N=2M，用M位二進(jìn)制數(shù)表示，則從左至

右的十進(jìn)制權(quán)值為：N/2、N/4、N/8,…、2,1對(duì)倒序數(shù)J，其下一個(gè)序數(shù)是在該序數(shù)J的二進(jìn)制首位碼加1，相當(dāng)于十進(jìn)制運(yùn)算J+N/2。如果最高位是0（J<N/2），則(J+N/2）得下一個(gè)倒序值；如果最高位是1（即J≥N/2），則要將最高位變?yōu)?（），次最高位加（）。次最高位加1時(shí)，同樣要判斷0、1值，依次類(lèi)推，直到完成最高位加1，逢2向右進(jìn)位的運(yùn)算，得出下一個(gè)倒序值。

形成倒位序以后，把按自然順序存放在存儲(chǔ)單元中的數(shù)據(jù)，重新按倒位序排列，就可以實(shí)現(xiàn)圖的變址功能，。圖給出了實(shí)現(xiàn)倒序的程序框圖，圖中虛線框內(nèi)就是計(jì)算倒序值的運(yùn)算流程。可以看出，I=J時(shí)不需要調(diào)換；

I<J時(shí)，調(diào)換A(I)和A(J)的內(nèi)容；I>J時(shí)，已調(diào)換過(guò)，不再調(diào)換。

6編程思想及程序框圖從輸入端（第一級(jí)）開(kāi)始，逐級(jí)進(jìn)行，共進(jìn)行M級(jí)運(yùn)算。在進(jìn)行第L級(jí)運(yùn)算時(shí)，依次求出2L-1個(gè)不同的旋轉(zhuǎn)因子，每求出一個(gè)旋轉(zhuǎn)因子，就計(jì)算完它對(duì)應(yīng)的所有2M-L個(gè)蝶形。這樣，就可以用三重循環(huán)程序來(lái)實(shí)現(xiàn)DIT-FFT運(yùn)算。程序框圖：

4.2.4按時(shí)間抽取的基-2FFT算法的其它形式流圖

圖4.2.12按時(shí)間抽取、輸入自然順序、輸出倒位序的FFT運(yùn)算流圖（N=8）

圖4.2.13按時(shí)間抽取、輸入和輸出都是自然順序的FFT運(yùn)算流圖（N=8）

圖4.2.14按時(shí)間抽取、各級(jí)具有相同幾何形狀的的FFT運(yùn)算流圖（N=8）4.3按頻率抽?。―IF）的基-2FFT算法

4.3.1算法原理

設(shè)序列x(n)長(zhǎng)度為N=2M，M為正整數(shù)。在把輸出序列X(k)按k的偶奇分組之前，先把輸入序列x(n)按n的順序分成前后兩半，這樣，可將x(n)的DFT寫(xiě)為前后兩部分，于是有

由于，故，可得

當(dāng)k為偶數(shù)時(shí)，(-1)k=1，當(dāng)k為奇數(shù)時(shí)，(-1)k=-1。因此，可以按k的偶、奇將X(k)分解為偶、奇序號(hào)組。

令則（）（）設(shè)（）則（）式和（）式可以寫(xiě)成以上運(yùn)算關(guān)系用蝶形運(yùn)算單元表示：N=8時(shí)分解過(guò)程

將N點(diǎn)DFT按頻率k的偶奇分解為兩個(gè)N/2點(diǎn)的DFT一個(gè)N/2點(diǎn)DFT分解成兩個(gè)N/4點(diǎn)DFT的分解過(guò)程一個(gè)N=8的完整的按頻率抽取的基2-FFT運(yùn)算流圖

\

4.3.2DIF-FFT算法特點(diǎn)及與DIT-FFT算法的異同

按頻率抽取法的運(yùn)算特點(diǎn)與按時(shí)間抽取法基本相同，它們是兩種等價(jià)的FFT運(yùn)算。整個(gè)DIF-FFT運(yùn)算流圖也全是由蝶形運(yùn)算構(gòu)成，每一個(gè)蝶形運(yùn)算單元完成的迭代運(yùn)算為

,式中L表示第L級(jí)迭代，1≤L≤M，k和j表示數(shù)據(jù)所在的行數(shù)。所對(duì)應(yīng)的蝶形運(yùn)算單元:

不同點(diǎn)：DIF-FFT算法與DIT-FFT算法的蝶形運(yùn)算單元不同，差別是乘旋轉(zhuǎn)因子的位置不同。DIF-FFT算法的蝶形運(yùn)算是先加（減）后相乘，而DIT-FFT算法的蝶形運(yùn)算是先乘后加（減）。DIF-FFT算法的輸入是自然順序，而輸出是倒位序的,運(yùn)算完畢后，要通過(guò)變址運(yùn)算將倒位序轉(zhuǎn)換為自然順序，然后再輸出。DIT-FFT算法的輸入是倒位序的，而輸出是自然順序,序列輸入后，要先進(jìn)行變址運(yùn)算將倒位序轉(zhuǎn)換為自然順序，然后再做碟形運(yùn)算。相同點(diǎn)：采用原位計(jì)算；N=2M時(shí)，共有M級(jí)運(yùn)算，每級(jí)共有N/2個(gè)蝶形，整個(gè)計(jì)算是通過(guò)MN/2個(gè)蝶形運(yùn)算完成的，DIT與DIF算法的運(yùn)算次數(shù)相同，復(fù)數(shù)乘法次，復(fù)數(shù)加法次。都有N/2個(gè)不同的旋轉(zhuǎn)因子。

4.3.3按頻率抽取法與按時(shí)間抽取法運(yùn)算流圖的關(guān)系

這兩個(gè)流圖的形式如同是顛倒了信號(hào)的傳輸方向，其關(guān)系可謂是“互為轉(zhuǎn)置”。利用信號(hào)流圖的轉(zhuǎn)置定理，可以導(dǎo)出其相應(yīng)的轉(zhuǎn)置流圖的。

信號(hào)流圖的轉(zhuǎn)置定理:如果將原信號(hào)流圖（網(wǎng)絡(luò)）中的所有支路方向倒轉(zhuǎn)或反向，保持各支路傳輸系數(shù)（增益）不變，并將輸入和輸出相互交換，則其相應(yīng)網(wǎng)絡(luò)的系統(tǒng)函數(shù)H(z)不改變。

對(duì)于N=2M，按頻率抽取的FFT算法與按時(shí)間抽取的FFT算法都有M級(jí)蝶形運(yùn)算，每一級(jí)都由N/2個(gè)蝶形運(yùn)算組成，因此，把每一種按時(shí)間抽取的FFT運(yùn)算流圖按轉(zhuǎn)置定理加以轉(zhuǎn)置都可以得到相應(yīng)的按頻率抽取的FFT運(yùn)算流圖，反之亦然。將右圖圖就可得到下圖加以轉(zhuǎn)置按頻率抽取，輸入倒位序、輸出自然順序的FFT運(yùn)算流圖（N=8）4.4離散傅里葉反變換（IDFT）的快速計(jì)算方法

上面討論的FFT算法，同樣可以適用于離散傅里葉反變換（IDFT）運(yùn)算，即快速傅里葉反變換，簡(jiǎn)稱(chēng)為IFFT。比較IDFT與DFT的公式：可以看出，只要把DFT運(yùn)算中的每一個(gè)旋轉(zhuǎn)因子換成，最后再乘以常數(shù)1/N，則以上所有按時(shí)間抽取或按頻率抽取的FFT算法都可以用來(lái)計(jì)算IDFT。

用FFT流圖實(shí)現(xiàn)IDFT快速算法將DIT-FFT或DIF-FFT蝶形運(yùn)算流圖中旋轉(zhuǎn)因子WNp改為WN-p在IDFT快速算法的最后結(jié)果輸出前，乘以1/N常數(shù)；如要防止溢出，可在每一級(jí)運(yùn)算中，輸出支路分別乘以1/2，實(shí)現(xiàn)系數(shù)分級(jí)分擔(dān)。在IDFT快速算法中，輸入序列為X(k)，而輸出序列為x(n)。流圖對(duì)應(yīng)關(guān)系：

DIT-FFT對(duì)應(yīng)DIF-IFFTDIF-FFT對(duì)應(yīng)DIT-IFFTDIT-IFFT運(yùn)算流圖（N=8）

直接利用FFT程序來(lái)計(jì)算IFFT的方法:由于DFT與IDFT公式結(jié)構(gòu)的對(duì)稱(chēng)性，對(duì)IDFT公式取共軛可得因此實(shí)現(xiàn)方法：先將X(k)取共軛，得到X*(k)；再直接調(diào)用FFT程序計(jì)算DFT[X*(k)]的值；最后將計(jì)算結(jié)果取一次共軛，并乘以1/N，即得到x(n)值。這一方法雖然多用了兩次對(duì)N點(diǎn)序列取共軛的操作，程序?qū)崿F(xiàn)容易，也不會(huì)明顯增加運(yùn)算量。但是，F(xiàn)FT與IFFT運(yùn)算可以共用一個(gè)子程序，實(shí)現(xiàn)方便。4.5N為復(fù)合數(shù)的FFT算法若不滿足基-2FFT算法(N=2M)，則可采用下面的幾種辦法：1.將x(n)補(bǔ)一些零值點(diǎn)，使N增長(zhǎng)到最鄰近的一個(gè)2M數(shù)值。由DFT的性質(zhì)知道，有限長(zhǎng)序列補(bǔ)零之后，并不影響其頻譜，只不過(guò)是使其頻譜的采樣點(diǎn)數(shù)增加而已，但是，有時(shí)計(jì)算量增加太多，會(huì)造成很大的浪費(fèi)。因此人們希望找到N≠2M時(shí)的FFT算法。2.如果要求準(zhǔn)確的N點(diǎn)DFT，而N又是素?cái)?shù)，則只能采用直接DFT方法，或者用后面將要介紹的線性調(diào)頻變換（Chirp-變換）算法。3.若N是一個(gè)復(fù)合數(shù)，即它可以分解成一些因子的乘積，則可以用FFT的一般算法，即混合基FFT算法，例如N=r1r2的FFT算法，而基-2算法只是這種一般算法的特例。其余略。4.6分裂基FFT算法

4.6.1按頻率抽取的基-4FFT算法

如果N=4M，并取分解式N=4×4×…×4，利用以下的標(biāo)號(hào)變換，就可以設(shè)計(jì)出基-4FFT算法。在每一步分解時(shí)，取N=N/4*4時(shí)，這樣得到的就是按頻率抽取的算法。這時(shí)n和k可分別表示為

于是序列x(n)的N點(diǎn)DFT可表示為

(4.6.3)

可以看出，把N=4M點(diǎn)DFT分解N/4為個(gè)4點(diǎn)DFT，將所得結(jié)果乘以旋轉(zhuǎn)因子1、、、，然后分成4組分別作N/4點(diǎn)DFT，即得所需的N點(diǎn)DFT，見(jiàn)圖。

將k0=0,1,2,3分別代入（）式，并用k表示k1，n表示n0，可以得出

（）

當(dāng)k從0增加到N/4-1時(shí)，（）式中的任一式均為頻域（對(duì)x(n)的N點(diǎn)DFT值X(k)）隔4點(diǎn)取1點(diǎn)的N/4點(diǎn)抽選。

4.6.2分裂基FFT算法的原理

從基-4FFT頻率抽選的（）式中可以看到，X(4k)和X(4k+2)這兩組是全部偶序號(hào)處的值，因而合在一起應(yīng)是隔2點(diǎn)取1點(diǎn)的點(diǎn)N/2點(diǎn)抽選，其表達(dá)式也可直接由基-2FFT頻率抽選分解得到。于是，(4.6.4)式又可寫(xiě)成如下形式

（）令

（）

則（）式可寫(xiě)成如下更簡(jiǎn)明的形式（）由此可見(jiàn)，一個(gè)N點(diǎn)DFT的計(jì)算可以分解為一個(gè)N/2點(diǎn)DFT和兩個(gè)N/4點(diǎn)DFT的計(jì)算。這種分解既有基-2（按二進(jìn)制抽選）部分，又有基-4（按四進(jìn)制抽選）部分?；?2部分X(2k)的奇數(shù)點(diǎn)部分又進(jìn)一步分解為基-4抽選分解，而基-4部分的偶數(shù)點(diǎn)部分又進(jìn)一步分解為基-2抽選分解，所以稱(chēng)（）式為分裂基頻率抽選分解，對(duì)應(yīng)的N點(diǎn)DFT一次分解的流圖如圖所示。圖4.6.3分裂基第一次分解L形流圖

運(yùn)算流圖（即許多L形的排列以及最后的4點(diǎn)和2點(diǎn)DFT流圖）。其排列示意圖見(jiàn)圖（a）、（b）。

（a）（b）圖4.6.4分裂基FFT算法L形排列示意圖與結(jié)構(gòu)示意圖

圖4.6.816點(diǎn)分裂基FFT運(yùn)算流圖

4.6.3分裂基FFT算法的運(yùn)算量

將圖的L形排列圖和圖的流圖一起觀察，可以看出N點(diǎn)分裂基FFT算法的全部復(fù)數(shù)乘法次數(shù)就是L形個(gè)數(shù)的2倍，而全部復(fù)數(shù)加法次數(shù)則與基-2FFT算法一樣為Nlog2N次。因此，復(fù)數(shù)乘法的總次數(shù):

（）與基-2FFT算法的復(fù)數(shù)乘法次數(shù)相比，Nlog2N前面的系數(shù)由1/2變?yōu)?/3，僅這一項(xiàng)就使復(fù)數(shù)乘法運(yùn)算次數(shù)下降33%，與基-4FFT算法的復(fù)數(shù)乘法次數(shù)3/8(Nlog2N)相比，可節(jié)省運(yùn)算量11%。4.7線性調(diào)頻變換（Chirp-變換）算法

4.7.1算法原理

已知序列x(n),(0≤n≤N-1)是有限長(zhǎng)序列，其z變換為,為適應(yīng)z可沿z平面更一般的路徑取值，就沿z平面上的一段螺線作等分角的采樣，z的這些采樣點(diǎn)zk為因此有其中A決定起始采樣點(diǎn)z0的位置A0表示z0的矢量半徑長(zhǎng)度,通常取A0≤1

θ0表示z0的相角

φ0表示兩相鄰采樣點(diǎn)之間的角度差W0（一般為正值）表示螺線的伸展率

圖4.7.1線性調(diào)頻變換在平面的螺線采樣

當(dāng)M=N，，（即，）時(shí)，各采樣點(diǎn)zk就均勻等間隔地分布在單位圓上，這就是求序列的DFT。將（）式的zk代入z變換表達(dá)式（）式中，可得(4.7.6)直接計(jì)算這一公式，與直接計(jì)算DFT相似。當(dāng)N和M數(shù)值很大時(shí)，運(yùn)算量會(huì)很大，因而限制了運(yùn)算速度。若采用布魯斯坦所提出的等式得到(4.7.7)

（）

即這一運(yùn)算過(guò)程如圖所示。圖4.7.2Chirp-變換運(yùn)算流圖

當(dāng)w0=1時(shí)，序列可以想象為頻率隨時(shí)間n成線性增長(zhǎng)的復(fù)指數(shù)序列，在雷達(dá)系統(tǒng)中這種信號(hào)稱(chēng)為線性調(diào)頻信號(hào)（ChirpSignal），因此這種變換又稱(chēng)為線性調(diào)頻z變換，即Chirp-z變換.

4.7.2線性調(diào)頻變換的實(shí)現(xiàn)

CZT運(yùn)算的實(shí)現(xiàn)步驟可由圖加以說(shuō)明(略)圖4.7.3Chirp-z變換的循環(huán)卷積計(jì)算

4.8實(shí)序列的FFT算法

4.8.1利用一次N點(diǎn)復(fù)序列的FFT計(jì)算兩個(gè)N點(diǎn)實(shí)序列的FFT

設(shè)x1(n)和x2(n)是兩個(gè)長(zhǎng)度均為N的有限長(zhǎng)實(shí)序列，將此二序列構(gòu)成一個(gè)復(fù)序列，其長(zhǎng)度仍然為N，得到,則利用有限長(zhǎng)序列的的對(duì)稱(chēng)性得:(4.8.1)

(4.8.2)

這樣，做一次N點(diǎn)復(fù)序列的FFT運(yùn)算，再按（）式和（）式將所得結(jié)果進(jìn)行組合，就同時(shí)得到了X1(k)和X2(k)，從而大大減少了總的計(jì)算量，使計(jì)算速度提高近一倍。

4.8.2利用一次N點(diǎn)復(fù)序列的FFT計(jì)算2N點(diǎn)實(shí)序列的FFT

設(shè)x(n)為2N點(diǎn)的實(shí)序列，將其按偶數(shù)點(diǎn)和奇數(shù)點(diǎn)分組形成兩個(gè)N點(diǎn)實(shí)序列，即且，則而所以利用x(n)實(shí)序列的DFT的對(duì)稱(chēng)性，得X(k)的另外N點(diǎn)的值為4.9用DFT的快速算法（FFT）實(shí)現(xiàn)線性卷積及線性相關(guān)

4.9.1用DFT（FFT）實(shí)現(xiàn)線性卷積

考慮一個(gè)線性時(shí)不變系統(tǒng)，輸入序列x(n)的長(zhǎng)度為M，系統(tǒng)的單位脈沖響應(yīng)h(n)的長(zhǎng)度為N，輸出響應(yīng)為y(n)。顯然有，且其長(zhǎng)度為(M+N-1)。為了用循環(huán)卷積來(lái)代替線性卷積的計(jì)算，循環(huán)卷積的點(diǎn)數(shù)L必須滿足L≥M+N-1,此時(shí)有L兩序列的L點(diǎn)循環(huán)卷積就代表它們的線性卷積。由DFT的時(shí)域循環(huán)卷積定理知道，循環(huán)卷積不僅可以在時(shí)域直接計(jì)算，而且也可以在頻域計(jì)算。因此，將序列x(n)和h(n)補(bǔ)上一定的零值點(diǎn)后，用DFT的快速算法（FFT）來(lái)計(jì)算兩序列的L點(diǎn)循環(huán)卷積，也就計(jì)算了它們的線性卷積，從而得到線性時(shí)不變系統(tǒng)的輸出響應(yīng)。

用DFT（FFT）計(jì)算兩序列x(n)和h(n)的線性卷積的步驟如下：（1）取L≥M+N-1，且使L=2J（J為正整數(shù)）。將輸入序列x(n)補(bǔ)上L-M個(gè)零值點(diǎn)，將線性時(shí)不變系統(tǒng)的單位沖激響應(yīng)h(n)補(bǔ)上L-N個(gè)零值點(diǎn)。（2）計(jì)算L點(diǎn)FFT，求（3）計(jì)算L點(diǎn)FFT，求（4）計(jì)算Y(k)=X(k)（5）計(jì)算L點(diǎn)IFFT，求這里計(jì)算DFT和IDFT都采用FFT和IFFT進(jìn)行計(jì)算，當(dāng)M、N足夠長(zhǎng)時(shí)，與直接計(jì)算線性卷積相比，計(jì)算速度可以大大提高。因而用FFT計(jì)算循環(huán)卷積以代替線性卷積的計(jì)算，常稱(chēng)為快速卷積。

用DFT的快速算法（FFT）計(jì)算線性卷積的運(yùn)算流程圖如圖所示圖4.9.1用DFT的快速算法（FFT）計(jì)算線性卷積的流程圖

4.9.2分段卷積

在實(shí)際應(yīng)用中，經(jīng)常遇到輸入信號(hào)x(n)的長(zhǎng)度與系統(tǒng)的單位脈沖響應(yīng)h(n)的長(zhǎng)度N相差很大的情況，用以上所述的快速卷積算法計(jì)算線性卷積，則要求對(duì)短序列補(bǔ)很多個(gè)零值點(diǎn)，且長(zhǎng)序列必須全部輸入后才能進(jìn)行計(jì)算。因此要求存儲(chǔ)容量大，運(yùn)算時(shí)間長(zhǎng)，很不經(jīng)濟(jì)，并會(huì)使處理延時(shí)很大，很難實(shí)時(shí)處理。解決這個(gè)問(wèn)題的方法是將長(zhǎng)序列分段計(jì)算，也就是將長(zhǎng)序列分成長(zhǎng)度和短序列相仿的段，分別求出每一段與短序列的卷積結(jié)果，然后用一定方式把它們合在一起，就得到總的輸出，這就是所謂的分段卷積。有重疊相加法和重疊保留法兩種。

1.重疊相加法設(shè)序列h(n)的長(zhǎng)度為N，x(n)為無(wú)限長(zhǎng)序列，且當(dāng)n<0時(shí)，x(n)=0。我們將x(n)均勻分段，每段為M點(diǎn)，M選擇成和N的數(shù)量級(jí)相同。序列x(n)可以表示成長(zhǎng)度為M的許多個(gè)平移有限長(zhǎng)序列之和，即

，式中

圖給出了有限長(zhǎng)單位脈沖響應(yīng)h(n)和未定義長(zhǎng)度的信號(hào)x(n)，以及對(duì)x(n)的分段方法的示意圖。要注意的是，每段中的第一個(gè)樣本均在n=0處，而第k段的第零個(gè)樣本xk(n)是序列的第kM個(gè)樣本。

從圖中可以看出，當(dāng)?shù)诙€(gè)分段卷積y1(n)計(jì)算完后，迭加重疊點(diǎn)便可得輸出序列y(n)的前2M個(gè)值，這種由分段卷積的序列段的重疊部分相加形成輸出序列的方法，通常稱(chēng)為重疊相加法。這樣進(jìn)行計(jì)算，要求存儲(chǔ)量小，且運(yùn)算量和延時(shí)也大大減少。

2.重疊保留法(略)

4.9.3快速相關(guān)

利用FFT計(jì)算相關(guān)函數(shù)，即是用循環(huán)相關(guān)代替線性相關(guān)，通常稱(chēng)為快速相關(guān)。循環(huán)相關(guān)也是按照循環(huán)移位的規(guī)則來(lái)定義相關(guān)函數(shù)的，與循環(huán)卷積不同之處在于沒(méi)有“翻褶”這一步驟。如果兩個(gè)序列長(zhǎng)度分別為M和N，若需用兩序列的循環(huán)相關(guān)來(lái)代替它們的線性相關(guān)，同樣需將這兩個(gè)序列補(bǔ)零到序列長(zhǎng)度L≥N+M-1才可避免混疊失真。由循環(huán)相關(guān)定理的（）式及（）式，可得出用DFT（FFT算法）計(jì)算計(jì)算線性相關(guān)的步驟見(jiàn)書(shū)p175：

4.10用DFT的快速算法（FFT）對(duì)信號(hào)進(jìn)行頻譜分析

假設(shè)xa(t)是經(jīng)過(guò)預(yù)濾波和截取處理的有限長(zhǎng)帶限信號(hào)。通常，對(duì)連續(xù)時(shí)間信號(hào)的離散傅里葉分析的處理步驟如圖所示。從圖中可以看到，在計(jì)算x(n)的DFT之前，用了一個(gè)窗函數(shù)W(n)加到x(n)上，這即是對(duì)x(n)的截取過(guò)程。

圖4.10.1連續(xù)時(shí)間信號(hào)離散傅里葉分析的處理步驟

假設(shè)連續(xù)時(shí)間非周期信號(hào)xa(t)持續(xù)時(shí)間為T(mén)p，最高頻率為fh。xa(t)的傅里葉變換對(duì)為(4.10.1)

(4.10.2)

xa(t)及如圖（a）所示。

圖4.10.2用DFT計(jì)算連續(xù)時(shí)間信號(hào)頻譜原理

（a）

（b）

（c）

時(shí)域離散化:對(duì)作如下近似，即則（）式可近似為(4.10.3)由于時(shí)域采樣，fs=1/T，所以引起頻譜以為周期的周期延拓，如果xa(t)是帶限信號(hào)，則有可能不產(chǎn)生混疊，如圖（b）所示。這里用表示采樣信號(hào)x(n)=x(nT)的頻譜。

頻域離散化：即對(duì)頻域函數(shù)的一個(gè)周期進(jìn)行N點(diǎn)等間隔采樣，采樣間隔為F，則fs=NF。于是，時(shí)域和頻域都是離散周期序列，Tp=1/F,因此，將其限制在一個(gè)周期之內(nèi)，分別取其主值序列進(jìn)行分析。各參數(shù)間的關(guān)系為

頻域采樣:由（）式，得令可得同理，由傅里葉反變換式（）可得：Xa(k)及x(n)如圖（c）所示。結(jié)論:(1)連續(xù)信號(hào)的頻譜特性可以通過(guò)對(duì)連續(xù)信號(hào)采樣，并進(jìn)行DFT再乘以T的近似方法得到。(2)連續(xù)信號(hào)的時(shí)域采樣信號(hào)可以通過(guò)對(duì)其頻譜函數(shù)進(jìn)行采樣，并進(jìn)行IDFT再乘以1/T的近似方法得到。4.10.2用DFT（FFT）對(duì)連續(xù)時(shí)間信號(hào)進(jìn)行頻譜分析時(shí)的幾個(gè)問(wèn)題

采樣頻率fs設(shè)連續(xù)時(shí)間信號(hào)xa(t)頻譜中最高頻率分量為fh。若采樣頻率為fs，按照奈奎斯特采樣定理，為了避免在DFT運(yùn)算中發(fā)生頻譜（頻率響應(yīng)）混疊的現(xiàn)象，必須使fs滿足

fs>2

fh→fh<fs/2若不能滿足這一要求，就會(huì)在折疊頻率f=fs/2（在數(shù)字頻率ω=）附近產(chǎn)生頻率響應(yīng)的混疊失真。在實(shí)際應(yīng)用中，一般取fs=(2.5~3)

fs。fs取得太大，會(huì)大大增加計(jì)算量，且要增加計(jì)算機(jī)內(nèi)存。

頻率分辨率F

對(duì)頻率函數(shù)采樣，使之變成離散的序列，其采樣間隔F就表示信號(hào)的頻率分辨率。它表示頻譜分析中能夠分辨的兩個(gè)頻譜分量的最小間隔。顯然，F(xiàn)越小，頻譜分析的結(jié)果就越接近相應(yīng)連續(xù)時(shí)間信號(hào)的頻率特性。所以，F(xiàn)較小時(shí)，稱(chēng)頻率分辨率較高。由及兩個(gè)公式來(lái)看，信號(hào)的最高頻率分量（又稱(chēng)高頻容量）與頻率分辨率之間存在著矛盾，要想增加，則時(shí)域采樣間隔就一定減小，而采樣頻率就增加，由于采樣點(diǎn)數(shù)滿足在采樣點(diǎn)數(shù)保持不變的情況下，就必然要增大，即頻率分辨率下降。相反，要提高頻率分辨率（減小），就要增加Tp，當(dāng)N一定時(shí)，必須增加T，即減小采樣頻率。為了不產(chǎn)生混疊失真，則必然會(huì)減小高頻容量fh，也就是減小了頻譜分析的范圍。

信號(hào)的最高頻率分量fk（又稱(chēng)高頻容量）與頻率分辨率F之間存在著矛盾，如果要想兼顧高頻容量與頻率分辨率，即使一個(gè)性能提高，而另一個(gè)性能不變（或也得以提高）的唯一辦法就是增加記錄長(zhǎng)度的點(diǎn)數(shù)N，因此必須滿足這個(gè)公式是在未采用任何特殊數(shù)據(jù)處理（例如加窗處理）的情況下，為實(shí)現(xiàn)基本DFT(FFT)算法必須滿足的最低條件。如果進(jìn)行加窗處理，相當(dāng)于時(shí)域相乘，對(duì)應(yīng)于頻率函數(shù)的卷積，必然會(huì)加寬頻譜分量，就可能使頻率分辨率下降。為了保證頻率分辨率不變，就必須增加記錄長(zhǎng)度Tp，這也就增加了采樣點(diǎn)數(shù)N。對(duì)于Tp，一般應(yīng)滿足

用DFT對(duì)連續(xù)信號(hào)進(jìn)行譜分析的參數(shù)選擇原則：(1)采樣頻率fs：fs>2fh(2)譜分辯率:F=fs/N(3)采樣點(diǎn)數(shù)N的選擇：N>2fh/F(4)信號(hào)觀察時(shí)間Tp的選擇：Tp1/F例4.10.1有一頻譜分析用的FFT處理器，其采樣點(diǎn)數(shù)必須是2的整數(shù)冪。假定沒(méi)有采用任何特殊的數(shù)據(jù)處理措施，若要求頻率分辨率F≥10Hz，信號(hào)最高頻率fh≤2kHz，試確定信號(hào)的最小記錄時(shí)間Tpmin；最大采樣間隔Tmax；在一個(gè)記錄中的最少點(diǎn)數(shù)Nmin。如果fh不變，要求頻率分辨率增加一倍，最少的采樣點(diǎn)數(shù)和最小的記錄時(shí)間是多少？解：（1）因?yàn)門(mén)p必須滿足,所以Tpmin=

（2）由，得（3）取N為2的整數(shù)冪,Nmin=29=512

若使頻率分辨率提高一倍，即F=5Hz，則

取Nmin=210=1024,

可見(jiàn)，記錄時(shí)間及采樣點(diǎn)數(shù)的增加可以提高頻率分辨率，但b必須滿足時(shí)域采樣定理。

2.頻譜泄漏

由于實(shí)際的需要，往往要把觀測(cè)信號(hào)x1(n)限制在一定的時(shí)間之內(nèi)，也就是要取出信號(hào)的某一時(shí)間段x2(n)，對(duì)于無(wú)限長(zhǎng)的信號(hào)序列，往往要將它截短成若干段有限長(zhǎng)序列來(lái)進(jìn)行處理，這種過(guò)程就是截?cái)鄶?shù)據(jù)的過(guò)程。數(shù)據(jù)截?cái)嘞喈?dāng)于是加窗處理。如果用矩形窗函數(shù)RN(n)，則x2(n)=x1(n)RN(n)，窗內(nèi)的數(shù)據(jù)并不改變。時(shí)域的截?cái)?，在頻域中則相當(dāng)于所研究的信號(hào)的頻譜與矩形窗函數(shù)頻譜的卷積，即其中

矩形窗的頻譜與原信號(hào)頻譜卷積的結(jié)果將造成失真的頻譜。對(duì)比下圖中與，可以看出，頻譜產(chǎn)生了“拖尾”（擴(kuò)展）現(xiàn)象，這種頻譜展寬（“擴(kuò)散”）的現(xiàn)象，稱(chēng)為頻譜泄漏。

圖4.10.3圖

例4.10.2設(shè)余弦序列，求其頻譜，并定性畫(huà)出的圖形及y(n)=x

(n)RN(n)的幅度譜示意圖。解：由于,r取整數(shù)所以

圖4.10.5余弦序列加矩形窗前后的頻譜(a)(