第3章5-8 離散傅里葉變換及其快速算法

1、Discrete Fourier Transform and Fast Algorithm 離散傅里葉變換在實際應(yīng)用中是非常重要的，利用它可以計算信號的頻譜、功率譜和線性卷積等。但是，如果使用定義式(3.20)來直接計算DFT，當(dāng)N很大時，即使使用高速計算機(jī)，所花的時間也太多。因此，如何提高計算DFT的速度，便成了重要的研究課題。1965年庫利 (Cooley)和圖基(Tukey)在前人的研究成果的基礎(chǔ)上提出了快速計算DFT的算法，之后，又出現(xiàn)了各種各樣快速計算DFT的方法，這些方法統(tǒng)稱為快速傅里葉變換(Fast Fourier Transform)，簡稱為FFT。FFT的出現(xiàn)，使計算DFT的

2、計算量減少了兩個數(shù)量級，從而成為數(shù)字信號處理強(qiáng)有力的工具。 快速傅里葉變換(FFT)是離散傅里葉變換(DFT)的快速算法。它是DSP領(lǐng)域中的一項重大突破，它考慮了計算機(jī)和數(shù)字硬件實現(xiàn)的約束條件、研究了有利于機(jī)器操作的運(yùn)算結(jié)構(gòu)，使DFT的計算時間縮短了12個數(shù)量級，還有效地減少了計算所需的存儲容量，F(xiàn)FT技術(shù)的應(yīng)用極大地推動了DSP的理論和技術(shù)的發(fā)展。DFT的定義在導(dǎo)出FFT算法之前，首先來估計一下直接計算DFT所需的計算量。其中將DFT定義式展開成方程組將方程組寫成矩陣形式用向量表示 從矩陣形式表示可以看出，由于計算一個X(k)值需要N次復(fù)乘法和(N-1)次復(fù)數(shù)加法，因而計算N個X(k)值，共

3、需N2次復(fù)乘法和N(N-1)次復(fù)加法。每次復(fù)乘法包括4次實數(shù)乘法和2次實數(shù)加法，每次復(fù)加法包括2次實數(shù)加法，因此計算N點的DFT共需要4N2次實數(shù)乘法和(2N2+2N(N-1)次實數(shù)加法。當(dāng)N很大時，這是一個非常大的計算量。 FFT算法主要利用了WNk的兩個性質(zhì)：(1)對稱性，即(2)周期性，即用復(fù)數(shù)表示：r為任意整數(shù)。 FFT算法是基于可以將一個長度為N的序列的離散傅里葉變換逐次分解為較短的離散傅里葉變換來計算這一基本原理的。這一原理產(chǎn)生了許多不同的算法，但它們在計算速度上均取得了大致相當(dāng)?shù)母纳啤?在本章中我們集中討論兩類基本的FFT算法。 第一類 稱為按時間抽取(Decimation-in

4、-Time)的基2FFT算法，它的命名來自如下事實：在把原計算安排成較短變換的過程中，序列x(n)(通?？醋魇且粋€時間序列)可逐次分解為較短的子序列。 第二類稱為按頻率抽取(Decimation-in-Frequency)的基2FFT算法，在這類算法中是將離散傅里葉變換系數(shù)序列X(k)分解為較短的子序列。 前面兩種算法特別適用于N等于2的冪的情況。 對于N為合數(shù)的情況，本章也將介紹兩種處理方法。 這種算法簡稱為時間抽選FFT算法，其基本出發(fā)點是，利用旋轉(zhuǎn)因子WNk的對稱性和周期性，將一個大的DFT分解成一些逐次變小的DFT來計算。 分解過程遵循兩條規(guī)則：對時間進(jìn)行偶奇分解；對頻率進(jìn)行前后分解。

5、 設(shè)N2M，M為正整數(shù)。為了推導(dǎo)方便，取N238，即離散時間信號為 按照規(guī)則(1)，將序列x(n)分為奇偶兩組，一組序號為偶數(shù)，另一組序號為奇數(shù)，即分別表示為：根據(jù)DFT的定義 因為 WN2=WN/21，所以上式變?yōu)樯鲜街械腉(k)和H(k)都是N/2點的DFT。(3.64)按照規(guī)則(2)，將X(k)分成前后兩組，即由(3.64)表示的是N/2點DFT，前4個k值的DFT可表示為后4個k值的X(k)表示為：因為所以(3.66)(3.65)按照式(3.65)和式(3.66)可畫出圖315所示的信號流程圖。 式(3.65)和式(3.66)把原來N點DFT的計算分解成兩個N/2點DFT的計算。照此可

6、進(jìn)一 步把每個N/2點DFT的計算再各分解成兩個N/4點DFT的計算。具體說來，是把x(0)，x(2)，x(4)，x(6)和x(1)，x(3)，x(5)，x(7)分為x(0)，x(4) | x(2)，x(6)和x(1)，x(5) | x(3)，x(7)。這樣，原信號序列被分成x(0)，x(4) | x(2)，x(6) I x(1)，x(5) I x(3)，x(7)4個2項信號。G(k)和H(k)分別計算如下:(3.67)(3.68)(3.69)(3.70) 這樣，用式(3.67)(3.70)4個公式就可計算圖3.15中的兩組N/2點DFT。圖3.16所示的是其中一組G(k)的計算。 將圖3.1

7、6與圖3.15所示的信號流程圖合并，便得到圖3.17所示的信號流程圖。因為N=8，所以上圖中N/4點的DFT就是2點的DFT，不能再分解了。 將2點DFT的信號流程圖移入圖3.17，得到圖3.19所示的8點時間抽選的完整的FFT流程圖。從圖3.19中可看出幾個特點：(1)流程圖的每一級的基本計算單元都是一個蝶形； (2)輸入x(n)不按自然順序排列，稱為“混序”排列，而輸出 X(k)按自然順序排列，稱為“正序”排列，因而要對輸入進(jìn)行“變址”；(3)由于流程圖的基本運(yùn)算單元為蝶形，所以可以進(jìn)行“同址”或“原位”計算。 任何一個N為2的整數(shù)冪(即N=2M)的DFT，都可以通過M次分解，最后成為2點

8、的 DFT來計算。M次分解構(gòu)成了從x(n)到X(k)的M級迭代計算，每級由N/2個蝶形組成。圖3.20表示了蝶形的一般形式表示。其輸入和輸出之間滿足下列關(guān)系： 從上式可以看出完成一個蝶形計算需一次復(fù)數(shù)乘法和兩次復(fù)數(shù)加法。因此，完成N點的時間抽選FFT計算的總運(yùn)算量為 大多數(shù)情況下復(fù)數(shù)乘法所花的時間最多，因此下面僅以復(fù)數(shù)乘法的計算次數(shù)為例來與直接計算進(jìn)行比較。直接計算DFT需要的乘法次數(shù)為D=N2，于是有例如，當(dāng)N=1024時，則： 205，即直接計算DFT所需復(fù)數(shù)乘法次數(shù)約為FFT的205倍。顯然，N越大，F(xiàn)FT的速度優(yōu)勢越大。 表3. 2列出了不同N值所對應(yīng)的兩種計算方法的復(fù)數(shù)乘法次數(shù)和它們

9、的比值。 圖3. 19包含log2N級迭代運(yùn)算，每級由N/2個蝶形計算構(gòu)成。蝶形計算的優(yōu)點是可以進(jìn)行所謂同址或原位計算。 現(xiàn)在來考察第一級的計算規(guī)律。設(shè)將輸入x(0)，x(4)，x(2)，x(6)， x(1)，x(5)，x(3)，x(7)分別存入計算機(jī)的存儲單元M(1), M(2), M(3),，M(7)和M(8)中。首先，存儲單元M(1)和M(2)中的數(shù)據(jù)x(0)和x(4)進(jìn)入運(yùn)算器并進(jìn)行蝶形運(yùn)算，流圖中各蝶形的輸入量或輸出量是互不相重的，任何一個蝶形的二個輸入量經(jīng)蝶形運(yùn)算后，便失去了利用價值，不再需要保存。這樣，蝶形運(yùn)算后的結(jié)果便可以送到M(1)和M(2)存儲起來。類似地，M(3)和M(4

10、)中的x(2)和x(6)進(jìn)入運(yùn)算器進(jìn)行蝶形運(yùn)算后的結(jié)果也被送回 到M(3)和M(4)保存，等等。第二級運(yùn)算與第一級類似，不過，M(1)和M(3)存儲單元中的數(shù) 據(jù)進(jìn)行蝶形運(yùn)算后的結(jié)果被送回M(1)和M(3)保存，M(2)和M(4)中的數(shù)據(jù)進(jìn)行蝶形運(yùn)算 后送回M(2)和M(4)保存，等等。這樣一直到最后一級的最后一個蝶形運(yùn)算完成。 可以看出， 每一級的蝶形的輸入與輸出在運(yùn)算前后可以存儲在同一地址(原來位置上)的存儲單元中，這種同址運(yùn)算的優(yōu)點是可以節(jié)省存儲單元，從而降低對計算機(jī)存儲量的要求或降低硬件實現(xiàn)的成本。 蝶形運(yùn)算的特點是，首先每一個蝶形運(yùn)算都需要兩個輸入數(shù)據(jù)，計算結(jié)果也是兩個數(shù)據(jù)，與其它結(jié)

11、點的數(shù)據(jù)無關(guān)，其它蝶形運(yùn)算也與這兩結(jié)點的數(shù)據(jù)無關(guān)、因此，一個蝶形 運(yùn)算一旦計算完畢，原輸入數(shù)據(jù)便失效了。這就意味著輸出數(shù)據(jù)可以立即使用原輸 人數(shù)據(jù)結(jié)點所占用的內(nèi)存。原來的數(shù)據(jù)也就消失了。輸出、輸人數(shù)據(jù)利用同一內(nèi)存單 元的這種蝶形計算稱為同位(址)計算。 從圖3. 19所示的流程圖看出，同址計算要求輸入x(n)是“混序”排列的。所謂輸入為“混 序”，并不是說輸入是雜亂無章的，實際上它是有規(guī)律的。如果輸入x(n)的序號用二進(jìn)制碼來 表示，就可以發(fā)現(xiàn)輸入的順序恰好是正序輸入的“碼位倒置”，表3. 3列出了這種規(guī)律。 在實際運(yùn)算中，按碼位倒置順序輸入數(shù)據(jù)x(n)，特別當(dāng)N較大時，是很不方便的。因此，數(shù)

12、據(jù)總是按自然順序輸入存儲，然后通過“變址”運(yùn)算將自然順序轉(zhuǎn)換成碼位倒置順序存儲。實現(xiàn)這種轉(zhuǎn)換的程序可用圖3. 21來說明。 圖中用n表示自然順序的標(biāo)號，用l表示碼位倒置的標(biāo)號。當(dāng)l=n時，x(n)和x(l)不必互相調(diào)換。當(dāng)ln時， 必須將x(l)和x(n)互相調(diào)換，但只能調(diào)換一次，為此必須規(guī)定每當(dāng)ln時，要將x(l)和x(n)相互調(diào)換，即把原來存放x(n)的存儲單元中的數(shù)據(jù)調(diào)入存儲x(l)的存儲單元中，而把原來存儲x(l)的存儲單元中的數(shù)據(jù)調(diào)入到存儲x(n)的存儲單元中。 這樣，按自然序輸入的數(shù)據(jù)x(n)經(jīng)過變址計算后變成了碼位倒置的排列順序，便可進(jìn)入第一級的蝶形運(yùn)算。 最后介紹一下時間抽選F

13、FT算法的另外一些形式的流程圖。對于任何流程圖，只要保持 各節(jié)點所連支路及其傳輸系數(shù)不變，則不論節(jié)點位置怎樣排列，所得到的流程圖總是等效的，因而都能得到DFT的正確結(jié)果，只是數(shù)據(jù)的提取和存儲次序不同而已。 把圖3. 19中與x(4)水平相鄰的所有節(jié)點和與x(1)水平相鄰的所有節(jié)點交換，把與x(6)水平相鄰的所有節(jié)點和與 x(3)水平相鄰的所有節(jié)點交換，而與x(2)、x(5)和x(7)水平相鄰各節(jié)點位置不變，就可以從圖3. 19得到圖3.22。圖3.22與圖3.19的區(qū)別只是節(jié)點的排列不同，而支路傳輸比，即WN的 各次冪保持不變。顯然圖3.22所示流程圖的輸入是正序(自然順序)排列的，輸出是碼位

14、倒置 排列的，所以輸出要進(jìn)行變址計算。圖3. 22所示的流程圖相當(dāng)于最初由庫利和圖基給出的時 間抽選算法。 另一種形式的流程圖是將節(jié)點排列成輸入 和輸出兩者都是正序排列，但這類流程圖不能進(jìn)行同址計算，因而需要兩列 長度為N的復(fù)數(shù)存儲器。 頻率抽選基2FFT算法簡稱為頻率抽選，它的推導(dǎo)過程遵循兩個規(guī)則：對時間前后分；對頻率偶奇分。 設(shè)N2M，M為正整數(shù)。為推導(dǎo)方便，取N=238。首先，根據(jù)規(guī)則(1)，將x(n)分成前一半和后一半，即 x(n)x(0), x(1), x(2), x(3), I x(4), x(5), x(6), x(7) 這樣 (3. 72) 式(3. 72)雖然包含兩個N/2點

15、求和公式，但是在每個求和公式中出現(xiàn)的是WNkn，而不是WN/2kn，因此這兩個求和公式都不是N/2點的DFT。如果合并這兩個求和公式，并利用則得：其次，按規(guī)則(2)，將頻率偶奇分，即X(k)=X(0), X(2), X(4), X(6), | X(1), X(3), X(5), X(7)如果用X(2r)和X(2r+1)分別表示頻率的偶數(shù)項和奇數(shù)項，則有令得到上面兩式所表示的是N/2的DFT。在實際計算中，首先形成序列g(shù)(n)和h(n)，然后計算h(n)WNn，最后分別計算g(n) 和h(n)WNn的N/2點DFT，便得到偶數(shù)輸出點和奇數(shù)輸出點的DFT。計算流程圖如圖3. 24所示。 由于N是2

16、的整數(shù)冪，所以N/2仍然是偶數(shù)。這樣，可以將N/2點DFT的輸出再分為偶數(shù)組和奇數(shù)組，也就是將N/2點的DFT計算進(jìn)一步分解為兩個N/4點的DFT計算，其推導(dǎo)過程如下。 將g(n)分為前后兩組，得到因為所以對頻率再進(jìn)行偶奇分，則得頻率的偶數(shù)項為頻率的奇數(shù)項為 通過類似的推導(dǎo)可得 上面4式所表示的計算都是N/4點的DFT計算，從而得到圖3.25所示的形式。 與時間抽選的FFT算法一樣，圖3.27所示的流程圖的基本運(yùn)算也是蝶形運(yùn)算，但是與時間抽選中的蝶形單元(圖3.20)有所不同。圖3. 28所示的是頻率抽選FFT算法的蝶形單元的一般化的流程圖 顯然，一個蝶形的計算包括1次復(fù)數(shù)乘法和2次復(fù)數(shù)加法。

17、圖3. 27所示的整個流程圖共有l(wèi)og2N級迭代運(yùn)算，每級有N/2個蝶形。因此，總計算量為頻率抽選的FFT算法的計算量與時間抽選FFT算法的計算量是一樣。 與時間抽選算法一樣，頻率抽選FFT算法也具有同址(原位)計算的優(yōu)點。但是，與時間抽選不同的是，頻率抽選FFT算法的信號輸入為正序排列，輸出為碼位倒置排列，因此輸出要進(jìn)行變址計算。 FFT算法同樣可以應(yīng)用于IDFT的計算，稱為快速傅里葉反變換，簡寫為IFFT。前述DFT和IDFT公式為 比較上面兩式，可以看出，只要把DFT公式中的系數(shù) 改為 ，并乘以系數(shù)1/N，就可用FFT算法來計算IDFT，這就得到了IFFT的算法。 當(dāng)把時間抽選FFT算法

18、用于 IFFT計算時，由于原來輸入的時間序列x(n)現(xiàn)在變?yōu)轭l率序列X(k)，原來是將x(n)偶奇分的，而現(xiàn)在變成對X(k)進(jìn)行偶奇分了，因此這種算法改稱為頻率抽選IFFT算法。類似地，當(dāng)把頻率抽選FFT算法用于計算IFFT時，應(yīng)該稱為時間抽選IFFT算法。 在IFFT計算中經(jīng)常把常量1/N分解成M個1/2連乘，即1/N(1/2)M，并且在M級的迭代運(yùn)算中，每級的運(yùn)算都分別乘 上一個1/2因子。圖3.29表示的是時間抽選IFFT流程圖。 上面討論的以2為基(即N2M)的時間抽選和頻率抽選FFT算法，由于具有程序簡單、 計算效率高、對存儲量要求不很高等優(yōu)點，因而在實際中得到了最廣泛的應(yīng)用。如果N

19、不等于 2的冪2M，通常有兩種處理辦法：(1)用補(bǔ)零的辦法將x(n)延長為2M。例如N=60，可在序列x(n)的末尾填補(bǔ)4個0，即 令x(60)=x(61) =x(62)=x(63)=0，使N達(dá)到2664，這樣就可使用基2FFT算法。有限長序列補(bǔ)零以后，只是頻譜的取樣點有所增加而不會影響它的頻譜X(ej)的形狀。 (2)采用以任意數(shù)為基數(shù)的FFT算法。 設(shè)N等于兩個整數(shù)p和q 的乘積，即N=pq，則可將N點DFT分解成p個q點DFT或q個p點DFT來計算。為此，首先將x(n) 分為p組，每組長為q，即 例如，N=18=36，即p=3，q=6；將x(n)分成3組，每組各有6個序列值，即然后，將N

20、點DFT也分解為p組來計算，即 由于WNprk=WN/prk=Wqrk，因此是一個q點DFT，這樣上式可寫成從而說明：一個N=pq點的DFT可以用p個q點DFT來組成，如下圖所示。 在最一般的情況下，設(shè) N=p1p2pm，其中p1pm是m個素因子。首先把N分解為兩個因子，即N=p1q1，其中q1p2p3pm，并用以上討論的方法將DFT分解為p1個q1點DFT； 然后，將q1分解為q1=p2q2，其中q2=p3p4pm，即將每一個q1點DFT分解為p2個q2 點DFT；這樣，通過m次分解，最后達(dá)到pm點 DFT。這種算法可以使DFT的運(yùn)算獲得最高效率?？焖俑道锶~變換在數(shù)字信號處理中有著廣泛的應(yīng)用

21、。下面簡要地介紹它在譜分析和線性卷積計算中的應(yīng)用。所謂譜分析就是計算信號的頻譜，包括振幅譜、相位譜和功率譜。1. 譜分析中參數(shù)的選擇假設(shè)所處理的離散時間信號x(n)是從連續(xù)時間信號xa(t)中取樣得到的。下面的討論采用如下符號：T(取樣周期)，單位為s；fs(取樣頻率)，單位為Hz，fs=1/T；f0(連續(xù)時間信號xa(t)的最高頻率)，單位為Hz；F(xa(t)的頻率分辨率)，單位為Hz。所謂頻率分辨率是指頻域取樣中兩相鄰點間的頻率間隔。tp(信號xa(t)的最小記錄長度)，tp=1/F，單位為s；N(一個記錄長度中的取樣數(shù))。基帶信號的頻譜主要集中在低頻段。根據(jù)取樣定理，為了避免混疊失真，

22、要求最小記錄長度必須按所需的頻率分辨率來選擇，即(3.91)(3.92)在保持分辨率不變的情況下，若希望增加所分析的信號的最高頻率，或在保持信號最高頻率不變的情況下，提高分辨率，唯一的辦法是增加在記錄長度內(nèi)的取樣取樣點數(shù)N。如果f0和F都給定，那么N必須滿足條件(3.93)(1) 數(shù)據(jù)準(zhǔn)備設(shè)待分析的信號為任意長的連續(xù)時間信號xa(t) (t0)。若已知信號的最高頻率為f0，頻率分辨率為F，那么根據(jù)式(3.91)、式(3.92)和式(3.93)可分別求出取樣周期T、最小記錄長度tp和取樣點數(shù)N。如果由式(3.93)計算得到的N不是2的整數(shù)冪，則應(yīng)補(bǔ)充零取樣值，使N等于2的整數(shù)冪，以便利用基2FF

23、T算法。 在記錄長度中對xa(t)進(jìn)行N點取樣，得到(2)用FFT計算信號的頻譜。用FFT計算x(n)的頻譜，X(k)一般是由實部XR(k)和虛部XI(k)組成的復(fù)數(shù)，即X(k) = XR(k) + jXI(k) (3.95)(3)由X(k)求振幅、相位譜和功率譜。由式(3. 95)可求出振幅譜|X(k)|、相位譜(k)和功率譜S(k)，它們分別為解: 由頻率分辨率確定最小記錄長度為由信號最高頻率計算取樣周期記錄長度內(nèi)的取樣點數(shù)為取N=16，即在使用FFT計算e-nT的頻譜時，T不是一個重要參量，因而可令T=1，于是得運(yùn)行Matlab程序，便可計算出信號x(n)=e-n (0n15)的頻譜X(

24、k)，得到X(k)的實部、虛部、振幅譜和相位譜。頻譜=1.5820 1.4490-0.3090i 1.2029-0.4229i 1.0064-0.3981i 0.8808-0.3240i 0.8051-0.2399i 0.7611-0.1571i 0.7382-0.0776i 0.7311+0.0000i 0.7382+0.0776i 0.7611+0.1571i 0.8051+0.2399i 0.8808+0.3240i 1.0064+0.3981i 1.2029+0.4229i 1.4490+0.3090i振幅譜=1.5820 1.4816 1.2751 1.0823 0.9385 0.8

25、401 0.7772 0.7423 0.7311 0.7423 0.7772 0.8401 0.9385 1.0823 1.2751 1.4816相位譜=0 -0.2101 -0.3381 -0.3767 -0.3525 -0.2896 -0.2036 -0.1047 0 0.1047 0.2036 0.2896 0.3525 0.3767 0.3381 0.2101信號x(n)通過FIR數(shù)字濾波器得到的輸出等于x(n)與h(n)的線性卷積，即設(shè)信號x(n)的長度為N1，F(xiàn)IR數(shù)字濾波器的單位取樣響應(yīng)h(n)的長度為N2，即它們的線性卷積其線性卷積的結(jié)果y(n)是一個有限長序列，其非零值長度為

26、N1+N2-1。如果直接計算線性卷積，根據(jù)線性卷積的計算方法，其計算量為乘法次數(shù)：PDN1N2加法次數(shù)：QD(N1-1)(N2-1)兩個有限長序列的線性卷積可以用循環(huán)卷積來代替，而循環(huán)卷積可使用FFT來計算。為了使循環(huán)卷積與線性卷積計算結(jié)果相等，必須把x(n)和h(n)都延長到N點(NNl+N2-1)，延長的部分均為零取樣值。這樣，y(n)的計算由以下5步來完成。(1)將x(n)和h(n)都延長到N點，NN1+N2-1；(2)計算x(n)的N點FFT，即X(k)FFTx(n)；(3)計算h(n)的N點FFT，即H(k)FFTh(n)；(4)計算Y(k)X(k)H(k)；(5)計算Y(k)的反變

27、換，即y(n)IFFTX(k)H(k) 。實際中常使用基2FFT算法，因此，當(dāng)NN1+N2-1不為2的整數(shù)冪時，應(yīng)該用補(bǔ)零取樣值的方法將序列x(n)和h(n)都延長到最鄰近的2的整數(shù)冪的值。第2、第3和第5步各需要做一次FFT，第4步要做N次乘法，因此總計算量為現(xiàn)以乘法運(yùn)算次數(shù)為例，對線性卷積的直接計算方法和FFT計算方法進(jìn)行比較。令r表示PD與PF的比值，考慮到N=N1+N2-1，有首先討論x(n)和h(n)的長度相等的情況，這時N=2N1-12N1，于是對于不同的N1值，按上式計算得到的r值如下：由此可見，N1越大，用FFT或循環(huán)卷積計算線性卷積的優(yōu)越性就越大。因此，通常把循環(huán)卷積成為快速

28、卷積，而把直接（線性）卷積稱為慢速卷積。其次，討論信號x(n)相對于單位取樣響應(yīng)很長，即N1N2的情況，這時NN1，于是顯然，r值將下降，從而使循環(huán)卷積算法的優(yōu)點不能發(fā)揮。此時可采用分段卷積來實現(xiàn)即將x(n)分成許多小段，每段長度與h(n)的長度相近，然后用FFT算法進(jìn)行分段計算。具體見下一節(jié)。例3.3 設(shè)x(n)和h(n)分別為用循環(huán)卷積計算它們的線性卷積。解：(1)用補(bǔ)充零取樣值的方法將x(n)和h(n)都延長到16點，即(2)分別計算x(n)和h(n)的N=16點的FFT, 得X(k)=FFTx(n)和H(k)=FFTh(n)。(3)計算Y(k)=X(k)H(k)。(4)計算Y(k)的1

29、6點反變換，得y(n)=IFFTX(k)H(k)。從上面的討論知道，當(dāng)信號x(n)的長度和濾波器的單位取樣響應(yīng)h(n)的長度相差不大時，用循環(huán)卷積計算線性卷積比直接計算線性卷積的速度要快。但是，當(dāng)x(n)是一個很長的序列時，由于h(n)必須補(bǔ)很多零，因而循環(huán)卷積計算方法的效率將降低。同時，如果一次輸入點數(shù)太多，則要求計算時占用很大的存儲量。為了克服這些困難，可采用分段卷積的方法。分段卷積方法的計算過程是：首先把x(n)分成長度為L的若干段，這里L(fēng)比h(n)的長度M略長，如左圖所示；然后用循環(huán)卷積方法計算每段與h(n)的卷積；最后把各段的卷積結(jié)果以適當(dāng)?shù)姆绞綌M合在一起。分段卷積方法有重疊相加法和

30、重疊保留法兩種。將x(n)分成長為L的段，每段表示為x(n)可表示為xk(n)之和，即于是其中，yk(n)的長度為L+M-1。對h(n)和xk(n)都增添零取樣值，使它們的長度都為N=L+M-1。計算h(n)與xk(n)的N點循環(huán)卷積，得到由于yk(n)的長度為L+M-1，而xk(n)長度為L，所以相鄰兩段yk(n)序列必然有M-1個點發(fā)生重疊，如圖3.35(b)所示。這些重疊部分應(yīng)該相加起來才能構(gòu)成最后的輸出序列，這就是“重疊相加法”這一名稱的由來。重疊相加法用FFT處理的步驟歸納如下：(1)計算h(n)的N點DFT H(k)，N=L+M-1；(2)計算xk(n)的N點DFT X(k)，N=

31、L+M-1；(3)計算Yk(k)=X(k)H(k)；(4)求Yk(k)的N點反變換，即yk(n)IFFTXk(k)H(k)；(5)將yk(n)的重疊部分相加起來，最后輸出為如果將重疊相加法中分段序列中補(bǔ)零的部分改為保留原來輸入序列的值，如圖3.36(b)中的虛線間的部分所示，且用yk(n)表示每段與h(n)的循環(huán)卷積，那么yk(n)將出現(xiàn)混疊，即每一段開始的M-1個點的序列樣本是前一段最后M-1個點的序列樣本。長度混疊發(fā)生在yk(n)的起始部分，因此，yk(n)的開始部分有M-1個點是不正確的，yk(n)可表示為現(xiàn)在，來分析圖3.36中yk(n)的局部混疊是怎樣產(chǎn)生的。由于h(n)長為M，當(dāng)0

32、nM-2時，h(n-m)NRN(n)在xk(n)的尾部出現(xiàn)非零值，所以yk(n)中混有xk(n)與h(n-m)NRN(n)的卷積值，于是如圖3.36(c)和(d)所示。當(dāng)M-1nN-1時，h(n-m)NRN(n)在xk(n)尾部無非零值，這樣，如圖3.36(e)所示。因此，yk(n)的前M-1個點，即圖3.36(f)打叉部分不是xk(n)與h(n)的真正卷積值，應(yīng)把它去掉，yk(n)中只有后面L個點才是正確的。理解為什么產(chǎn)生混疊：1.由于將補(bǔ)零的部分改為保留原來輸入序列的值，導(dǎo)致此部分值在相鄰兩段中計算了兩次與h(n)的卷積，這兩次計算就會導(dǎo)致混疊。2.圖c為h(n)折疊后移位前的情況。每次移

33、位就循環(huán)右移一位，可以看出，如果在移位次數(shù)n小于M-1時，均會導(dǎo)致xk(m)與h(n-m)NRN(n)的乘積相加在L-1與N-1之間出現(xiàn)兩次，從而出現(xiàn)混疊?，F(xiàn)將x(n)分解為用FFT計算循環(huán)卷積令因此在實際應(yīng)用中，有時只對信號的某一頻段感興趣，或只需計算單位圓上某一段的頻譜值。例如，在對窄帶信號進(jìn)行分析時，常希望在窄頻帶內(nèi)對頻率的取樣很密集，以便提高頻率分辨率，而在窄頻帶外不予以考慮。對于這種情況，如果采用DFT方法，則需要在窄頻帶內(nèi)外都增加頻域取樣點，而窄頻帶外的計算量是浪費的。此外，有時對非單位圓上的取樣感興趣，例如在語音信號處理中，常常需要知道其Z變換的極點所在處的復(fù)頻率，這時就需要在這些極點附近的曲線上進(jìn)行頻域取樣，這樣，就要沿著螺旋線對Z變換取樣。這種沿螺旋線上取樣點計算的Z變換，稱為線性調(diào)頻Z變換(Chirp Z Transform，簡稱CZT)。下面就來討論這種變換的原理及其計算方法。一個長度為N的序列x(n)，其Z變換為(3.110)為了使z可以沿著z平面更一般的路徑(不只是單位圓)取值，可以沿一段螺旋線對z作等分角取樣，這些取樣點上的zk表示為(3.111)其中M為所要分析的復(fù)頻譜的點數(shù)，不一定等于N。W和A為任意復(fù)數(shù)，可表示為(3.112)得到(3.114)取樣點zk所在的路徑如圖3.38所示，圖中：(1)A0表示起始取樣點z0的矢量長度，通常A01，否