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

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

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

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

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

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

6、進(jìn)一 步把每個(gè)N/2點(diǎn)DFT的計(jì)算再各分解成兩個(gè)N/4點(diǎn)DFT的計(jì)算。具體說(shuō)來(lái)，是把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)。這樣，原信號(hào)序列被分成x(0)，x(4) | x(2)，x(6) I x(1)，x(5) I x(3)，x(7)4個(gè)2項(xiàng)信號(hào)。G(k)和H(k)分別計(jì)算如下:(3.67)(3.68)(3.69)(3.70) 這樣，用式(3.67)(3.70)4個(gè)公式就可計(jì)算圖3.15中的兩組N/2點(diǎn)DFT。圖3.16所示的是其中一組G(k)的計(jì)算。 將圖3.1

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

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

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

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

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

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

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

14、倒置 排列的，所以輸出要進(jìn)行變址計(jì)算。圖3. 22所示的流程圖相當(dāng)于最初由庫(kù)利和圖基給出的時(shí) 間抽選算法。 另一種形式的流程圖是將節(jié)點(diǎn)排列成輸入 和輸出兩者都是正序排列，但這類(lèi)流程圖不能進(jìn)行同址計(jì)算，因而需要兩列 長(zhǎng)度為N的復(fù)數(shù)存儲(chǔ)器。 頻率抽選基2FFT算法簡(jiǎn)稱為頻率抽選，它的推導(dǎo)過(guò)程遵循兩個(gè)規(guī)則：對(duì)時(shí)間前后分；對(duì)頻率偶奇分。 設(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)雖然包含兩個(gè)N/2點(diǎn)

15、求和公式，但是在每個(gè)求和公式中出現(xiàn)的是WNkn，而不是WN/2kn，因此這兩個(gè)求和公式都不是N/2點(diǎn)的DFT。如果合并這兩個(gè)求和公式，并利用則得：其次，按規(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ù)項(xiàng)和奇數(shù)項(xiàng)，則有令得到上面兩式所表示的是N/2的DFT。在實(shí)際計(jì)算中，首先形成序列g(shù)(n)和h(n)，然后計(jì)算h(n)WNn，最后分別計(jì)算g(n) 和h(n)WNn的N/2點(diǎn)DFT，便得到偶數(shù)輸出點(diǎn)和奇數(shù)輸出點(diǎn)的DFT。計(jì)算流程圖如圖3. 24所示。 由于N是2

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

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

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

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

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

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

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

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

24、k)，得到X(k)的實(shí)部、虛部、振幅譜和相位譜。頻譜=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信號(hào)x(n)通過(guò)FIR數(shù)字濾波器得到的輸出等于x(n)與h(n)的線性卷積，即設(shè)信號(hào)x(n)的長(zhǎng)度為N1，F(xiàn)IR數(shù)字濾波器的單位取樣響應(yīng)h(n)的長(zhǎng)度為N2，即它們的線性卷積其線性卷積的結(jié)果y(n)是一個(gè)有限長(zhǎng)序列，其非零值長(zhǎng)度為

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

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

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

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

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

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

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

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