《數(shù)字信號處理》課件第3章

第3章快速傅里葉變換3.1引言3.2直接計算DFT的問題及改進的途徑3.3按時間抽?。―IT）的基2-FFT算法3.4按頻率抽?。―IF）的基2-FFT算法3.5N為復(fù)合數(shù)的FFT算法3.6線性調(diào)頻Z變換（Chirp-Z變換）算法3.7利用FFT分析時域連續(xù)信號頻譜3.8FFT的其他應(yīng)用3.1引言

快速傅里葉變換（FFT）并不是一種新的變換，而是離散傅里葉變換（DFT）的一種快速算法。由于有限長序列在其頻域也可離散化為有限長序列（DFT），因此離散傅里葉變換（DFT）在數(shù)字信號處理中是非常有用的。例如，在信號的頻譜分析、系統(tǒng)的分析、設(shè)計和實現(xiàn)中都會用到DFT的計算。但是，在相當(dāng)長的時間里，由于DFT的計算量太大，即使采用計算機也很難對問題進行實時處理，所以并沒有得到真正的運用。直到1965年首次發(fā)現(xiàn)了DFT運算的一種快速算法以后，情況才發(fā)生了根本的變化。人們開始認(rèn)識到DFT運算的一些內(nèi)在規(guī)律，從而很快地發(fā)展和完善了一套高速有效的運算方法，這就是現(xiàn)在人們普遍稱之為快速傅里葉變換（FFT）的算法。FFT出現(xiàn)后使DFT的運算大大簡化，運算時間一般可縮短一二個數(shù)量級之多，從而使DFT的運算在實際中真正得到了廣泛的應(yīng)用。3.2直接計算DFT的問題及改進的途徑3.2.1直接計算DFT的運算量問題

設(shè)x(n)為N點有限長序列，其DFT為k=0,1,…,N-1（3-1）反變換（IDFT）為n=0,1,…,N-1（3-2）

二者的差別只在于WN的指數(shù)符號不同，以及差一個常數(shù)乘因子1/N，所以IDFT與DFT具有相同的運算工作量。下面我們只討論DFT的運算量。一般來說，x(n)和WNnk都是復(fù)數(shù)，X(k)也是復(fù)數(shù)，因此每計算一個X(k)值，需要N次復(fù)數(shù)乘法和N-1次復(fù)數(shù)加法。而X(k)一共有N個點（k從0取到N-1），所以完成整個DFT運算總共需要N2次復(fù)數(shù)乘法及N(N-1)次復(fù)數(shù)加法。在這些運算中乘法運算要比加法運算復(fù)雜，需要的運算時間也多一些。因為復(fù)數(shù)運算實際上是由實數(shù)運算來完成的，這時DFT運算式可寫成（3-3）由此可見，一次復(fù)數(shù)乘法需用四次實數(shù)乘法和二次實數(shù)加法；一次復(fù)數(shù)加法需二次實數(shù)加法。因而每運算一個X(k)需4N次實數(shù)乘法和2N+2(N-1)=2(2N-1)次實數(shù)加法。所以，整個DFT運算總共需要4N2次實數(shù)乘法和2N(2N-1)次實數(shù)加法。

當(dāng)然，上述統(tǒng)計與實際需要的運算次數(shù)稍有出入，因為某些WNnk可能是1或j，就不必相乘了，例如W0N=1，WNＮ／２=-1，WNN/4=-j等就不需乘法。但是為了便于和其他運算方法作比較，一般都不考慮這些特殊情況，而是把WNnk都看成復(fù)數(shù)，當(dāng)N很大時，這種特例的影響很小。從上面的統(tǒng)計可以看到，直接計算DFT，乘法次數(shù)和加法次數(shù)都是和N2成正比的，當(dāng)N很大時，運算量是很可觀的，有時是無法忍受的。

例3-1

根據(jù)式（3-1），對一幅N×N點的二維圖像進行DFT變換，如用每秒可做10萬次復(fù)數(shù)乘法的計算機，當(dāng)N=1024時，問需要多少時間（不考慮加法運算時間）？

解直接計算DFT所需復(fù)乘次數(shù)為（N2)2≈1012次，因此用每秒可做10萬次復(fù)數(shù)乘法的計算機，則需要近3000小時。這對實時性很強的信號處理來說，要么提高計算速度，而這樣，對計算速度的要求太高了。另外，只能通過改進對DFT的計算方法，以大大減少運算次數(shù)。3.2.2改善途徑能否減少運算量，從而縮短計算時間呢？仔細(xì)觀察DFT的運算就可看出，利用系數(shù)WＮnk的以下固有特性，就可減少運算量：（1）WNnk的對稱性（2）WNnk的周期性（3）WNnk的可約性另外

這樣，利用這些特性，使DFT運算中有些項可以合并，并能使DFT分解為更少點數(shù)的DFT運算。而前面已經(jīng)說到，DFT的運算量是與N2成正比的，所以N越小越有利，因而小點數(shù)序列的DFT比大點數(shù)序列的DFT的運算量要小?？焖俑道锶~變換算法正是基于這樣的基本思想而發(fā)展起來的。它的算法形式有很多種，但基本上可以分成兩大類，即按時間抽取(ＤecimationinTime，縮寫為DIT）法和按頻率抽取（DecimationinＦrequency，縮寫為DIF）法。3.3按時間抽?。―IT）的基-2FFT算法3.3.1算法原理設(shè)序列x(n)長度為N，且滿足N=2M，M為正整數(shù)。按n的奇偶把x(n)分解為兩個N/2點的子序列：（3-4）則可將DFT化為由于 ，故上式可表示成（3-5）式中，X1(k)與X2(k)分別是x1(r)及x2(r)的N/2點DFT：（3-6）（3-7）由此，我們可以看到，一個N點DFT已分解成兩個N/2點的DFT。這兩個N/2點的DFT再按照式（3-5）組合成一個N點DFT。這里應(yīng)該看到X1(k)，X2(k)只有N/2個點，即k=0,1,…,N/2-1。而X(k)卻有N個點，即k=0,1,…,N-1，故用式（3-5）計算得到的只是X(k)的前一半的結(jié)果，要用X1(k)，X2(k)來表達全部的X(k)值，還必須應(yīng)用系數(shù)的周期性，即這樣可得到（3-8）同理可得（3-9）式（3-8）、式（3-9）說明了后半部分k值（N/2≤k≤N-1）所對應(yīng)的X1(k)，Ｘ2(k)分別等于前半部分k值（0≤k≤N/2-1）所對應(yīng)的X1(k)，X2(k)。再考慮到WkN的以下性質(zhì):這樣，把式（3-8）、式（3-9）、式（3-10）代入式（3-5），就可將X(k)表達為前后兩部分：（3-10）（3-11）（3-12）圖3-1時間抽取法蝶形運算流圖符號圖3-2按時間抽取將一個N點DFT分解為兩個N/2點DFT（N=8）

一個N點DFT分解為兩個N/2點DFT，每一個N/2點DFT只需(N/2)2=N2/4次復(fù)數(shù)乘法，N/2(N/2-1)次復(fù)數(shù)加法。兩個N/2點DFT共需2×(N/2)2=N2/2次復(fù)數(shù)乘法和N(N/2-1)次復(fù)數(shù)加法。此外，把兩個N/2點DFT合成為N點DFT時，有N/2個蝶形運算，還需要N/2次復(fù)數(shù)乘法及2×N/2=N次復(fù)數(shù)加法。因而通過第一步分解后，總共需要(N2/2)+(N/2)=N(N+1)/2≈N2/2次復(fù)數(shù)乘法和N(N/2-1)+N=N2/2次復(fù)數(shù)加法。由此可見，通過這樣分解后運算工作量差不多節(jié)省了一半。既然這樣分解是有效的，由于N=2M，因而N/2仍是偶數(shù)，可以進一步把每個N/2點子序列再按其奇偶部分分解為兩個N/4點的子序列。（3-13）且式中：（3-14）（3-15）

圖3-3給出N=8時，將一個N/2點DFT分解成兩個N/4點DFT，由這兩個N/4點DFT組合成一個N/2點DFT的流圖。圖3-3N/2點DFT分解為兩個N/4點DFTX2(k)也可進行同樣的分解：式中：(3-16)(3-17)圖3-4一個N=8點DFT分解為四個N/4點DFT

根據(jù)上面同樣的分析知道，利用四個N/4點的DFT及兩級蝶形組合運算來計算N點DFT，比只用一次分解蝶形組合方式的計算量又減少了大約一半。最后，剩下的是2點DFT，對于此例N=8，就是四個N/4=2點DFT，其輸出為X3(k)，X4(k)，X5(k)，X6(k)，k=0,1，這由式（3-14）～式（3-17）四個式子可以計算出來。例如，由式（3-14）可得：式中，,故上式不需要乘法。類似地,可求出X4(k)，X5(k)，X6(k)。

這種方法的每一步分解，都是按輸入序列在時間上的次序是屬于偶數(shù)還是屬于奇數(shù)來分解為兩個更短的子序列，所以稱為“按時間抽取法”。圖3-5N=8按時間抽取的FFT運算流圖3.3.2按時間抽取的FFT算法與直接計算DFT運算量的比較由按時間抽取法FFT的流圖可見，當(dāng)N=2M時，共有M級蝶形，每級都由N/2個蝶形運算組成，每個蝶形需要一次復(fù)乘、二次復(fù)加，因而每級運算都需N/2次復(fù)乘和N次復(fù)加，這樣M級運算總共需要復(fù)乘數(shù)復(fù)加數(shù)式中，數(shù)學(xué)符號lb=log2。(3-18)(3-19)

實際計算量與這個數(shù)字稍有不同，因為W0N=1，ＷNN/2=-1，ＷNN/2=-j，這幾個系數(shù)都不用乘法運算，但是這些情況在直接計算DFT中也是存在的。此外，當(dāng)N較大時，這些特例相對而言就很少。所以，為了統(tǒng)一作比較起見，下面都不考慮這些特例。由于計算機上乘法運算所需的時間比加法運算所需的時間多得多，故以乘法為例，直接DFT復(fù)數(shù)乘法次數(shù)是N2，F(xiàn)FT復(fù)數(shù)乘法次數(shù)是(N/2)lbN。直接計算DFT與FFT算法的計算量之比為當(dāng)N=2048時，這一比值為372.4，即直接計算DFT的運算量是FFT運算量的372.4倍。當(dāng)點數(shù)N越大時，F(xiàn)FT的優(yōu)點更為明顯。（3-20）

例3-2

用FFT算法處理一幅N×N點的二維圖像，如用每秒可做10萬次復(fù)數(shù)乘法的計算機，當(dāng)N=1024時，問需要多少時間（不考慮加法運算時間）？

解當(dāng)N=1024點時，F(xiàn)FT算法處理一幅二維圖像所需復(fù)數(shù)乘法約為 次，僅為直接計算DFT所需時間的10萬分之一。即原需要3000小時，現(xiàn)在只需要2分鐘。3.3.3按時間抽取的FFT算法的特點及DIT-FFT程序框圖

1.原位運算（同址運算）從圖3-5可以看出這種運算是很有規(guī)律的，其每級（每列）計算都是由N/2個蝶形運算構(gòu)成的，每一個蝶形結(jié)構(gòu)完成下述基本迭代運算：（3-21a）（3-21b）式中，m表示第m列迭代，k,j為數(shù)據(jù)所在行數(shù)。式（3-21）的蝶形運算如圖3-6所示，由一次復(fù)乘和兩次復(fù)加（減）組成。圖3-6蝶形運算單元

由圖3-5的流圖看出，某一列的任何兩個節(jié)點k和j的節(jié)點變量進行蝶形運算后，得到結(jié)果為下一列k,j兩節(jié)點的節(jié)點變量，而和其他節(jié)點變量無關(guān)，因而可以采用原位運算，即某一列的N個數(shù)據(jù)送到存儲器后，經(jīng)蝶形運算，其結(jié)果為下一列數(shù)據(jù)，它們以蝶形為單位仍存儲在這同一組存儲器中，直到最后輸出，中間無需其他存儲器。也就是蝶形的兩個輸出值仍放回蝶形的兩個輸入所在的存儲器中。每列的N/2個蝶形運算全部完成后，再開始下一列的蝶形運算。這樣存儲器數(shù)據(jù)只需N個存儲單元。下一級的運算仍采用這種原位方式，只不過進入蝶形結(jié)的組合關(guān)系有所不同。這種原位運算結(jié)構(gòu)可以節(jié)省存儲單元，降低設(shè)備成本。

2.倒位序規(guī)律觀察圖3-5的同址計算結(jié)構(gòu)，發(fā)現(xiàn)當(dāng)運算完成后，F(xiàn)FT的輸出X(k)按正常順序排列在存儲單元中，即按X(0)，X(1)，…，X(7)的順序排列，但是這時輸入x(n)卻不是按自然順序存儲的，而是按x(0)，x(4)，…，x(7)的順序存入存儲單元，看起來好像是“混亂無序”的，實際上是有規(guī)律的，我們稱之為倒位序。

造成倒位序的原因是輸入x(n)按標(biāo)號n的偶奇的不斷分組。如果n用二進制數(shù)表示為(n2n1n0)2（當(dāng)N=8=23時，二進制為三位），第一次分組，由圖3-2看出，n為偶數(shù)（相當(dāng)于n的二進制數(shù)的最低位n0=0）在上半部分，n為奇數(shù)（相當(dāng)于n的二進制數(shù)的最低位n0=1）在下半部分。下一次則根據(jù)次最低位n1為“0”或是“1”來分偶奇（而不管原來的子序列是偶序列還是奇序列），如此繼續(xù)分下去，直到最后N個長度為1的子序列。圖3-7的樹狀圖描述了這種分成偶數(shù)子序列和奇數(shù)子序列的過程。圖3-7倒位序的形成

一般實際運算中，總是先按自然順序?qū)⑤斎胄蛄写嫒氪鎯卧瑸榱说玫降刮恍虻呐帕校覀兺ㄟ^變址運算來完成。如果輸入序列的自然順序號I用二進制數(shù)（例如n2n1n0）表示，則其倒位序J對應(yīng)的二進制數(shù)就是（n0n1n2），這樣，在原來自然順序時應(yīng)該放x(I)的單元，現(xiàn)在倒位序后應(yīng)放x(J)。例如,N=8時，x(3)的標(biāo)號是I=3，它的二進制數(shù)是011，倒位序的二進制數(shù)是110，即J=6，所以原來存放在x(011)單元的數(shù)據(jù)現(xiàn)在應(yīng)該存放在x(110)內(nèi)。表3-1列出了N=8時的自然順序二進制數(shù)以及相應(yīng)的倒位序二進制數(shù)。表3-1N=8時的自然順序二進制數(shù)和相應(yīng)的倒位序二進制數(shù)自然順序（I）二進制數(shù)倒位序二進制數(shù)倒位序（J）0123456700000101001110010111011100010001011000110101111104261537

由表3-1可見，自然順序數(shù)I增加1，是在順序數(shù)的二進制數(shù)最低位加1，向左進位。而倒序數(shù)J則是在二進制數(shù)最高位加1，逢2向右進位。例如，在（000）最高位加1，則得（100），再在（100）最高位加1，向右進位，則得（010）。因（100）最高位為1，所以最高位加1要向次高位進位，其實質(zhì)是將最高位變?yōu)?，再在次高位加1。用這種算法，可以從當(dāng)前任一倒序值求得下一個倒序值。

對于N=2M，M為二進制數(shù)最高位，其權(quán)值為N/2，且從左向右二進制位的權(quán)值依次為N/4，N/8，…，2，1。因此，最高位加1相當(dāng)于十進制運算J+N/2。如果最高位是0（J<N/2），則直接由J+N/2得下一個倒序值；如果最高位是1（J≥N/2），則要將最高位變?yōu)?（J

J-N/2），次高位加1（J+N/4）。但次高位加1時，同樣要判斷次高位0、1值，如果為0（J<N/4），則直接加1（J

J+N/4）;否則，將次高位變?yōu)?（JJ-N/4），再判斷下一位；以此類推，直到完成最高位加1，逢2向右進位的運算。

把按自然順序存放在存儲單元中的數(shù)據(jù)，換成FFT原位運算流圖所要求的倒位序的變址功能如圖3-8所示，當(dāng)I=J時，不必調(diào)換，當(dāng)I≠J時，必須將原來存放數(shù)據(jù)x(I)的存儲單元內(nèi)調(diào)入數(shù)據(jù)x(J)，而將存放x(J)的存儲單元內(nèi)調(diào)入x(I)。為了避免把已調(diào)換過的數(shù)據(jù)再次調(diào)換，保證只調(diào)換一次（否則又回到原狀），我們只需看J是否比I小。若J比I小，則意味著此x(I)在前邊已和x(J)互相調(diào)換過，不必再調(diào)換了;只有當(dāng)J>I時，才將原存放x(I)及存放x(J)的存儲單元內(nèi)的內(nèi)容互換。這樣就得到輸入所需的倒位序列的順序?？梢钥闯?，其結(jié)果與圖3-5的要求是一致的。圖3-8N=8倒位序的變址處理

3.蝶形運算兩節(jié)點的“距離”以圖3-5的8點FFT為例，其輸入是倒位序的，輸出是自然順序的。其第一級（第一列）每個蝶形的兩節(jié)點間“距離”為1，第二級每個蝶形的兩節(jié)點“距離”為2，第三級每個蝶形的兩節(jié)點“距離”為4。由此類推得，對N=2M點FFT，當(dāng)輸入為倒位序，輸出為正常順序時，其第m級運算，每個蝶形的兩節(jié)點“距離”為2m-1。

4.WNr的確定由于對第m級運算，一個DFT蝶形運算的兩節(jié)點“距離”為2m-1，因而式（3-21）可寫成:（3-22a）（3-22b）

為了完成上兩式運算，還必須知道系數(shù)ＷNr的變換規(guī)律。仔細(xì)觀察圖3-5的流圖可以發(fā)現(xiàn)r的變換規(guī)律是：①把式（3-22）中蝶形運算兩節(jié)點中的第一個節(jié)點標(biāo)號值，即k值，表示成M位（注意N=2M）二進制數(shù)；②把此二進制數(shù)乘上2M-m，即將此M位二進制數(shù)左移M-m位（注意m是第m級運算），把右邊空出的位置補零，此數(shù)即為所求r的二進制數(shù)。從圖3-5看出，WNr因子最后一列有N/2種，順序為WN0，WN1，…，，其余可類推。

5.DIT-FFT程序框圖

總結(jié)上述運算規(guī)律，便可采用下述運算方法。先從輸入端（第一級）開始，逐級進行，共進行M級運算。在進行m級運算時，依次求出2m-1個不同的系數(shù)WrN，每求出一個系數(shù)，就計算完它對應(yīng)的所有2M-m個蝶形。這樣，我們可用三重循環(huán)程序?qū)崿F(xiàn)DIT－FFT運算，程序框圖如圖3－9所示。圖3－9中倒序運算程序框圖見圖3－10。圖3-9DIT-FFT運算程序框圖圖3-10倒序程序框圖3.3.4按時間抽取的FFT算法的其他形式流圖顯然，對于任何流圖，只要保持各節(jié)點所連的支路及傳輸系數(shù)不變，則不論節(jié)點位置怎么排列所得流圖總是等效的，所得最后結(jié)果都是x(n)的DFT的正確結(jié)果，只是數(shù)據(jù)的提取和存放的次序不同而已。這樣就可得到按時間抽取的FFT算法的若干其他形式流圖。將圖3-5中和x(4)水平相連的所有節(jié)點和x(1)水平相連的所有節(jié)點位置對調(diào)，再將和x(6)水平相連的所有節(jié)點與和x(3)水平相連的所有節(jié)點對調(diào)，其余諸節(jié)點保持不變，可得圖3-11的流圖。圖3-11與圖3-5的蝶形相同，運算量也一樣，不同點是：

（1）數(shù)據(jù)存放的方式不同，圖3-5是輸入倒位序、輸出自然順序，圖3-11是輸入自然順序、輸出倒位序；（2）取用系數(shù)的順序不同，圖3-5的最后一列是按 的順序取用系數(shù)，且其前一列所用系數(shù)是后一列所用系數(shù)中具有偶數(shù)冪的那些系數(shù)（例如 ）；圖3-11是最后一列是按的順序取用系數(shù)，且其前一列所用系數(shù)是后一列所用系數(shù)的前一半，這種流圖是最初由庫利和圖基給出的時間抽取法。經(jīng)過簡單變換，也可得輸入與輸出都是按自然順序排列的流圖以及其他各種形式的流圖。圖3-11時間抽取、輸入自然順序、輸出倒位序的FFT流圖3.4按頻率抽?。―IF）的基-2FFT算法3.4.1算法原理仍設(shè)序列點數(shù)為N=2M，M為正整數(shù)。在把輸出X(k)按k的奇偶分組之前，先把輸入序列按前一半、后一半分開（不是按偶數(shù)、奇數(shù)分開），把N點DFT寫成兩部分。k=0,1,…,N-1式中，用的是，而不是，因而這并不是N/2點DFT。由于，故 ，可得

k=0,1,…,N-1（3-23）

當(dāng)k為偶數(shù)時，(-1)k=1；k為奇數(shù)時，(-1)k=-1。因此，按k的奇偶可將X(k)分為兩部分:（3-24）（3-25）

式（3-24）為前一半輸入與后一半輸入之和的N/2點DFT，式（3-25）為前一半輸入與后一半輸入之差再與WNn之積的N/2點DFT。令:（3-26）則有（3-27）圖3-12頻率抽取法蝶形運算單元

這樣，我們就把一個N點DFT按k的奇偶分解為兩個N/2點的DFT了。

N=8時，上述分解過程示于圖3-13。與時間抽取法的推導(dǎo)過程一樣，由于N=2M，N/2仍是一個偶數(shù)，因而可以將每個N/2點DFT的輸出再分解為偶數(shù)組與奇數(shù)組，這就將N/2點DFT進一步分解為兩個N/4點DFT。這兩個N/4點DFT的輸入也是先將N/2點DFT的輸入上下對半分開后通過蝶形運算而形成的，圖3-14示出了這一步分解的過程。圖3-13按頻率抽取的第一次分解圖3-14按頻率抽取的第二次分解(N=8)

這樣的分解可以一直進行到第M次（N=2M），第M次實際上是做兩點DFT，它只有加減運算。然而，為了有統(tǒng)一的運算結(jié)構(gòu)，仍然用一個系數(shù)為WN0的蝶形運算來表示，這N/2個兩點DFT的N個輸出就是x(n)的N點DFT的結(jié)果X(k)。圖3-15表示一個N=8完整的按頻率抽取的基-2FFT運算結(jié)構(gòu)。圖3-15按頻率抽取的FFT（N=8）信號流圖3.4.2按頻率抽取法的運算特點按頻率抽取法的運算特點與時間抽取法基本相同。從圖3-15可以看出，它也是通過(N/2)M個蝶形運算完成的。每一個蝶形結(jié)構(gòu)完成下述基本迭代運算：式中，m表示m列迭代，k,j為整數(shù)所在行數(shù)，此式的蝶形運算如圖3-16所示，也需要一次復(fù)乘和兩次復(fù)加。圖3-16頻率抽取法蝶形運算單元

由圖3-15與圖3-5相比較，初看起來，DIF法與DIT法的區(qū)別是：圖3-15的DIF輸入是自然順序，輸出是倒位序的，這與圖3-5的DIT法正好相反。但這不是實質(zhì)性的區(qū)別，因為DIF法與DIT法一樣，都可將輸入或輸出進行重排，使二者的輸入或輸出順序變成自然順序或倒位序順序。DIF的基本蝶形（圖3-16）與DIT的基本蝶形（圖3-6）則有所不同，這才是實質(zhì)的不同，DIF的復(fù)數(shù)乘法只出現(xiàn)在減法之后，DIT則是先作復(fù)乘后再作加減法。

但是，DIF與DIT就運算量來說則是相同的，即都有M級（列）運算，每級運算需N/2個蝶形運算來完成，總共需要mF=(N/2)lbN次復(fù)乘與aF=NlbN次復(fù)加，DIF法與DIT法都可進行原位運算。頻率抽取FFT算法的輸入是自然順序，輸出是倒位序的。因此運算完畢后，要通過變址計算將倒位序轉(zhuǎn)換成自然位序，然后再輸出。轉(zhuǎn)換方法與時間抽取法相同。

由按時間抽取法與按頻率抽取法的基本蝶形（圖3-6與圖3-16）運算看出，如果將DIT的基本蝶形加以轉(zhuǎn)置，就得到DIF的基本蝶形;反過來,將DIF的基本蝶形加以轉(zhuǎn)置，就得到DIT的基本蝶形，因而DIT法與DIF法的基本蝶形是互為轉(zhuǎn)置的。按照轉(zhuǎn)置定理，兩個流圖的輸入-輸出特性必然相同。轉(zhuǎn)置就是將流圖的所有支路方向都反向，并且交換輸入與輸出，但節(jié)點變量值不交換，這樣即可從圖3-6得到圖3-16或者從圖3-16得到圖3-6，因而對每一種按時間抽取的FFT流圖都存在一個按頻率抽取的FFT流圖。這樣把圖3-5，圖3-11的流圖分別加以轉(zhuǎn)置，就可得到不同DIF的FFT流圖。因此可以說，有多少種按時間抽取的FFT流圖就存在多少種按頻率抽取的FFT流圖。頻率抽取法與時間抽取法是兩種等價的FFT運算。3.5N為復(fù)合數(shù)的FFT算法

上面討論的是序列的點數(shù)N為2的冪次(即N=2M)情況下，按時間抽取和按頻率抽取的基-2FFT算法的基本原理。這種基-2FFT算法在實際中使用得最多，因為它的程序簡單，效率高，使用方便。但實際上無法保證總是處理長度為2的整數(shù)冪次的序列。若不滿足N=2M，可將x(n)增補一些零值點，以使N增長到最鄰近的一個2M數(shù)值。有限長序列補零之后，并不影響其頻譜X(ejω)，只不過其頻譜的抽樣點數(shù)增加了，所造成的結(jié)果是增加了計算量而已。但是，有時計算量增加太多，浪費較大。

例如,x(n)的點數(shù)N=300，則須補到N=29=512，要補212個零值點，因而人們才研究N≠2M時的FFT算法。若N是一個復(fù)合數(shù)，即它可以分解成一些因子的乘積，則可以用FFT的一般算法，即混合基FFT算法，如庫利-圖基（CooleyTukey）算法，而基-2算法只是這種一般算法的特例。這里，我們不作詳細(xì)介紹，感興趣的話可參考文獻［1］。總之，不管采用什么方法，計算DFT的高效算法是把計算長度為N的序列的DFT逐次分解成計算長度較短序列的DFT。這是很多高效算法的標(biāo)準(zhǔn)方法和基本原理。3.6線性調(diào)頻Z變換（Chirp-Z變換）算法前面已講過，采用FFT算法可以很快算出全部DFT值，也就是算出有限長序列x(n)的Z變換X(z)在Z平面單位圓上N個等間隔采樣點zk處的采樣值。實際上，常常只對信號的某一頻段感興趣，也就是只需要計算單位圓上某一段的頻譜值。例如，對窄帶信號就是這樣，希望在窄帶頻帶內(nèi)頻率的采樣能夠非常密集，分辨率要高，帶外則不予考慮；如果用DFT方法，則需增加頻域采樣點數(shù)，增加了窄帶之外不需要的計算量。此外，有時也對非單位圓上的采樣感興趣。例如，在語音信號處理中，常常需要知道其Z變換的極點所在頻率，如果極點位置離單位圓較遠，如圖3－17(a)上圖所示，則其單位圓上的頻譜就很平滑，如圖3－17(a)下圖，這時很難從中識別出極點所在的頻率;要是采樣不是沿單位圓而是沿一條接近這些極點的弧線進行，如圖3－17(b)上圖所示，則所得的結(jié)果將會在極點所在頻率上出現(xiàn)明顯的尖峰，如圖3－17(b)下圖所示，這樣，極點頻率的測定就要準(zhǔn)確得多。再有，如果N是大素數(shù)時，不能加以分解，如何有效計算這種序列的DFT。從以上三方面看，螺線采樣就是一種適應(yīng)于這些需要的變換，并且可以采用FFT來快速計算。這種變換稱為線性調(diào)頻Z變換（ChirpZ變換，CZT），它是適用于這種更為一般情況下由x(n)求X(zk)的快速變換算法。圖3-17單位圓與非單位圓采樣（a）沿單位圓采樣;(b)沿AB弧采樣3.6.1算法基本原理已知x(n)（0≤n≤N-1）是有限長序列，其Z變換為（3-28）為適應(yīng)z可以沿Z平面更一般的路徑取值，故沿Z平面上的一段螺線作等分角的采樣，z的這些采樣點zk為zk=AW-k

k=0,1,…,M-1（3-29）M為所要分析的復(fù)頻率的點數(shù)，即采樣點的總數(shù)，不一定等于N；A和W都是任意復(fù)數(shù)，可表示為:將式（3-30）與式（3-31）代入式（3-29），可得（3-30）（3-31）（3-32）因此有：圖3-18螺線采樣

采樣點在Z平面上所沿的周線如圖3-18所示。由以上討論和圖3-18可以看出：（1）A0表示起始采樣點z0的矢量半徑長度，通常A0≤1;否則z0將處于單位圓|z|=1的外部。（2）θ0表示起始采樣點z0的相角，它可以是正值或負(fù)值。（3）φ0表示兩相鄰采樣點之間的角度差。φ0為正時，表示zk的路徑是逆時針旋轉(zhuǎn)的；φ0為負(fù)時，表示zk的路徑是順時針旋轉(zhuǎn)的。

（4）W0的大小表示螺線的伸展率。W0>1時，隨著k的增加螺線內(nèi)縮;W0<1時，則隨k的增加螺線外伸;W0=1時，表示是半徑為A0的一段圓弧。若又有A0=1,則這段圓弧是單位圓的一部分。當(dāng)M=N，A=A0ejθ0=1，W=W0·e-jφ0=（W0=1，φ0=2π/N）這一特殊情況時，各zk就均勻等間隔地分布在單位圓上，這就是求序列的DFT。將式（3-29）的zk代入變換表達式（3-28），可得0≤k≤M-1(3-33)

直接計算這一公式，與直接計算DFT相似，總共算出M個采樣點，需要NM次復(fù)數(shù)乘法與(N-1)M次復(fù)數(shù)加法。當(dāng)N，M很大時，這個量很大，這就限制了運算速度。但是，下面我們將看到，通過一定的變換，以上運算可以轉(zhuǎn)換為卷積形式，從而可以采用FFT算法，這樣就可以大大提高運算速度。nk可以用以下表達式來替換將式（3-34）代入式（3-33），可得如果定義:n=0,1,…,N-1(3-34)(3-35)(3-36)(3-37)則它們的卷積為式中，k=0,1,…,M-1。式（3-38）正好是式（3-35）的一部分，因此式（3-35）又可以用卷積的形式表示為k=0,1,…,M-1(3-38)(3-39)由式（3-39）可看出，如果我們對信號按式（3-36）先進行一次加權(quán)處理，加權(quán)系數(shù)為;然后，通過一個單位脈沖響應(yīng)為h(n)的線性系統(tǒng)即求g(n)與h(n)的線性卷積;最后，對該系統(tǒng)的前M點輸出再做一次加權(quán)，這樣就得到了全部M點螺線采樣值X(zn)(n=0,1,…,M-1)。這個過程可以用圖3-19表示。從圖中我們看到，運算的主要部分是由線性系統(tǒng)來完成的。由于系統(tǒng)的單位脈沖響應(yīng) 可以想象為頻率隨時間(n)呈線性增長的復(fù)指數(shù)序列。在雷達系統(tǒng)中,這種信號稱為線性調(diào)頻信號（ChirpSignal），因此，這里的變換稱為線性調(diào)頻Z變換。圖3-1９Chirp-Z變換的線性系統(tǒng)表示3.6.2Chirp-Z變換（CZT）的實現(xiàn)步驟由式（3-37）可看出，線性系統(tǒng)h(n)是非因果的，當(dāng)n的取值為0到N-1，k的取值為0,1,…,M-1時，h(n)是在n=-(N-1)到n=M-1取值。也就是說，h(n)是一個有限長序列，點數(shù)為N+M-1，見圖3-20（a）。輸入信號g(n)也是有限長序列，點數(shù)為N。g(n)*h(n)的點數(shù)為2N+M-2，因而用圓周卷積代替線性卷積且不產(chǎn)生混疊失真條件是圓周卷積的點數(shù)應(yīng)大于或等于2N+M-2。但是，由于我們只需要前M個值X（zk）(k=0,1,…,M-1)，對以后的其他值是否有混疊失真并不感興趣，這樣可將圓周卷積的點數(shù)縮減到最小為N+M-1。當(dāng)然，為了進行基-2FFT運算，圓周卷積的點數(shù)應(yīng)取為L≥N+M-1，同時又滿足L=2m的最小L。這樣可將h(n)先補零值點，補到點數(shù)等于L，也就是從n=M開始補L-(N+M-1)個零值點，補到n=L-N處，或補L-(N+M-1)個任意序列值，然后將此序列以L為周期而進行周期延拓，再取主值序列，從而得到進行圓周卷積的一個序列，如圖3-20（b）所示。進行圓周卷積的另一個序列只需要將g(n)補上零值點，使之成為L點序列即可，如圖3-20（c）所示。圖3-20Chirp-Z變換的圓周卷積圖（M≤n≤L-1時，h(n)和g(n)的圓周卷積不代表線性卷積）

這樣，我們可以列出CZT運算的實現(xiàn)步驟：（1）選擇一個最小的整數(shù)L，使其滿足L≥N+M-1，同時滿足L=2m，以便采用基-2FFT算法。（2）將 （見圖3-20（c））補上零值點，變?yōu)長點序列，因而有：0≤n≤N-10≤n≤L-1

（3-40）并利用FFT法求此序列的L點DFT0≤r≤L-1

（3-41）0≤n≤M-1

L-N+1≤n≤L-10≤n≤M-1（3-42）

此h(n)可見圖3-20（b）。實際上它就是圖3-20（a）的序列 以L為周期的周期延拓序列的主值序列。

（3）形成L點序列h(n)，如上所述，在n=0到M-1一段 ，n=M~L-N段取h(n)為任意值（一般為零），在L=N+M-1到L-1段取h(n)為 的周期延拓序列 ，即有:對式（3-42）定義的h(n)序列，用FFT法求其L點DFT0≤r≤L-1（3-43）

（4）將H(r)和G(r)相乘，得Q(r)=H(r)G(r)，Q(r)為L點頻域離散序列。（5）用FFT法求Q(r)的L點IDFT，得h(n)和g(n)的圓周卷積L（3-44）式中，前M個值等于h(n)和g(n)的線性卷積結(jié)果［h(n)*g(n)］;n≥M的值沒有意義，不必去求。g(n)*h(n)即g(n)與h(n)圓周卷積的前M個值可見圖3-20（d）。（6）最后求X(zk)：（3-45）3.6.3運算量的估計

CZT的算法求X(zk)比直接求X(zk)的算法有效得多。CZT所需的乘法如下：（1）形成L點序列 ，但只有其中N點序列值，需要N次復(fù)乘，而系數(shù) 可事先準(zhǔn)備好，不必在實時分析時計算。

（2）形成L點序列h(n)。由于它是由 在-(N-1)≤n≤M-1的序列值構(gòu)成的，而 是偶對稱序列。如果設(shè)N>M，則只需要求得0≤n≤M-1一段N點序列值即可;h(n)也可事先準(zhǔn)備好，不必實時分析時計算，因此，可不用考慮其計算量。同時，h(n)的L點FFT即H(r)也可預(yù)先計算好。

（3）計算G(r)，q(n),需要二次L點FFT（或IFFT），共需要LlbL次復(fù)乘。（4）計算Q(r)=G(r)H(r),需要L次復(fù)乘。（5）計算X(zk)=（0≤k≤M-1），需要M次復(fù)乘。綜上所述，CZT總的復(fù)數(shù)乘法次數(shù)為

mc=LlbL+N+M+L

（3-46）前面說過，直接計算式（3-33）的X(zk)需要NM次復(fù)數(shù)乘法?？梢钥闯?，當(dāng)N,M都較大時（例如N,M都大于50時），CZT的FFT算法比直接算法的運算量要小得多。3.7利用FFT分析時域連續(xù)信號頻譜3.7.1基本步驟圖3-21時域連續(xù)信號離散傅里葉分析的處理步驟

在圖3-21中，前置低通濾波器LPF（預(yù)濾波器）的引入，是為了消除或減少時域連續(xù)信號轉(zhuǎn)換成序列時可能出現(xiàn)的頻譜混疊的影響。在實際工作中，時域離散信號x(n)的時寬是很長的，甚至是無限長的（例如語音或音樂信號）。由于DFT的需要(實際應(yīng)用FFT計算)，必須把x(n)限制在一定的時間區(qū)間之內(nèi)，即進行數(shù)據(jù)截斷。數(shù)據(jù)的截斷相當(dāng)于加窗處理（窗函數(shù)見6.2節(jié)）。因此，在計算FFT之前，用一個時域有限的窗函數(shù)w(n)加到x(n)上是非常必要的。

xc(t)通過A/D變換器轉(zhuǎn)換(忽略其幅度量化誤差)成采樣序列x(n)，其頻譜用X(ejω)表示，它是頻率ω的周期函數(shù)，即(3-47)式中，Xc(jΩ)或為xc(t)的頻譜。

在實際應(yīng)用中，前置低通濾波器的阻帶不可能是無限衰減的，故由Xc(jΩ)周期延拓得到的X(ejω)有非零重疊，即出現(xiàn)頻譜混疊現(xiàn)象。由于進行FFT的需要，必須對序列x(n)進行加窗處理，即v(n)=x(n)w(n)，加窗對頻域的影響，用卷積表示如下:(3-48)最后是進行FFT運算。加窗后的DFT是0≤k≤N-1(3-49)式中，假設(shè)窗函數(shù)長度L小于或等于DFT長度N，為進行FFT運算，這里選擇N為2的整數(shù)冪次即N=2m。

有限長序列v(n)=x(n)w(n)的DFT相當(dāng)于v(n)傅里葉變換的等間隔采樣。(3-50)V(k)便是sc(t)的離散頻率函數(shù)。因為DFT對應(yīng)的數(shù)字域頻率間隔為Δω=2π/N，且模擬頻率Ω和數(shù)字頻率ω間的關(guān)系為ω=ΩΤ，其中Ω=2πf。所以，離散的頻率函數(shù)第k點對應(yīng)的模擬頻率為:(3-51)(3-52)

由式（3-52）很明顯可看出，數(shù)字域頻率間隔Δω=2π/N對應(yīng)的模擬域譜線間距為(3-53)

譜線間距，又稱頻譜分辨率（單位：Hz）。所謂頻譜分辨率是指可分辨兩頻率的最小間距。它的意思是，如設(shè)某頻譜分析的F=5Hz，那么信號中頻率相差小于5Hz的兩個頻率分量在此頻譜圖中就分辨不出來。

長度N=16的時間信號v(n)=(1.1)nR16(n)的圖形如圖3-22(a)所示,其16點的DFTV(k)的示例如圖3-22(b)所示。其中T為采樣時間間隔（單位：s）；fs為采樣頻率（單位：Hz）；tp為截取連續(xù)時間信號的樣本長度（又稱記錄長度，單位：s）；F為譜線間距，又稱頻譜分辨率（單位：Hz）。注意：V(k)示例圖給出的頻率間距F及N個頻率點之間的頻率fs為對應(yīng)的模擬域頻率(單位：Hz)。圖3-22v(n)=(1.1)nR16(n)及DFTV（k)示例圖由圖可知：(3-54)(3-55)

在實際應(yīng)用中,要根據(jù)信號最高頻率fh和頻譜分辨率F的要求，來確定T、tp和N的大小。（1）首先，由采樣定理，為保證采樣信號不失真，fs≥2fh(fh為信號頻率的最高頻率分量，也就是前置低通濾波器阻帶的截止頻率)，即應(yīng)使采樣周期T滿足(3-56)（2）由頻譜分辨率F和T確定N,(3-57)

為了使用FFT運算，這里選擇N為2的冪次即N=2m，由式(3-55)可知，N大，分辨率好，但會增加樣本記錄時間tp。（3）最后由N,T確定最小記錄長度，tp=NT。[例3-3]

有一頻譜分析用的FFT處理器，其采樣點數(shù)必須是2的整數(shù)冪，假定沒有采用任何特殊的數(shù)據(jù)處理措施，已給條件為:①頻率分辨率≤10Hz;②信號最高頻率≤4kHz。試確定以下參量：

(1)最小記錄長度tp；

(2)最大采樣間隔T（即最小采樣頻率）；

(3)在一個記錄中的最少點數(shù)N。解(1)由分辨率的要求確定最小長度tp

所以記錄長度為(2)從信號的最高頻率確定最大可能的采樣間隔T（即最小采樣頻率fs=1/T）。按采樣定理即(3)最小記錄點數(shù)N應(yīng)滿足取

如果我們事先不知道信號的最高頻率，可以根據(jù)信號的時域波形圖來估計它。例如，某信號的波形如圖3-23所示。先找出相鄰的波峰與波谷之間的距離，如圖中t1，t2，t3，t4。然后，選出其中最小的一個如t4。這里,t4可能就是由信號的最高頻率分量形成的。峰與谷之間的距離就是周期的一半。因此，最高頻率為知道fh后就能確定采樣頻率圖3-23估算信號最高頻率fh

3.7.2可能出現(xiàn)的誤差

1.頻譜混疊失真在圖3-21畫出的基本步驟中，A/D變換前利用前置低通濾波器進行預(yù)濾波，使xc(t)頻譜中最高頻率分量不超過fh。假設(shè)A/D變換器的采樣頻率為fs，按照奈奎斯特采樣定理，為了不產(chǎn)生混疊，必須滿足也就是采樣間隔T滿足一般應(yīng)取fs=(2.5～3.0)fh

如果不滿足fs≥2fh，就會產(chǎn)生頻譜混疊失真。對于FFT來說，頻率函數(shù)也要采樣，變成離散的序列，其采樣間隔為F(即頻率分辨率)。由式（3-55）可得（3-58）（3-59）

從以上Ｔ和tp兩個公式來看，信號的最高頻率分量fh與頻率分辨率F存在矛盾關(guān)系，要想fh增加，則時域采樣間隔T就一定減小，而fs就增加，由式（3-57）可知，此時若是固定N，必然要增加F，即分辨率下降。

反之，要提高分辨率（減小F），就要增加tp，當(dāng)N給定時，必然導(dǎo)致T的增加（fs減?。Ｒ划a(chǎn)生混疊失真，則必然會減小高頻容量（信號的最高頻率分量）fh。要想兼顧高頻容量fh與頻率分辨率F，即一個性能提高而另一個性能不變（或也得以提高）的惟一辦法就是增加記錄長度的點數(shù)N，即要滿足（3-60）

這個公式是未采用任何特殊數(shù)據(jù)處理（例如加窗處理）的情況下，為實現(xiàn)基本FFT算法所必須滿足的最低條件。如果加窗處理，相當(dāng)于時域相乘，則頻域周期卷積，必然加寬頻譜分量，頻率分辨率就可能變差，為了保證頻率分辨率不變，則須增加數(shù)據(jù)長度tp。

2.柵欄效應(yīng)

利用FFT計算頻譜，只給出離散點ωk=2πk/N或Ωk=2πk/(NT)上的頻譜采樣值，而不可能得到連續(xù)頻譜函數(shù)，這就像通過一個“柵欄”觀看信號頻譜，只能在離散點上看到信號頻譜，稱之為“柵欄效應(yīng)”,如圖3-22所示。這時，如果在兩個離散的譜線之間有一個特別大的頻譜分量，就無法檢測出來了。

減小柵欄效應(yīng)的一個方法就是要使頻域采樣更密，即增加頻域采樣點數(shù)N，在不改變時域數(shù)據(jù)的情況下，必然是在數(shù)據(jù)末端添加一些零值點，使一個周期內(nèi)的點數(shù)增加，但并不改變原有的記錄數(shù)據(jù)。頻譜采樣為2πk/N，N增加，必然使樣點間距更近（單位圓上樣點更多），譜線更密，譜線變密后原來看不到的譜分量就有可能看到了。必須指出，補零以改變計算FFT的周期時，所用窗函數(shù)的寬度不能改變。換句話說，必須按照數(shù)據(jù)記錄的原來的實際長度選擇窗函數(shù)，而不能按照補了零值點后的長度來選擇窗函數(shù)。補零不能提高頻率分辨率，這是因為數(shù)據(jù)的實際長度仍為補零前的數(shù)據(jù)長度。

3.頻譜泄漏與譜間干擾對信號進行FFT計算，首先必須使其變成有限時寬的信號，這就相當(dāng)于信號在時域乘一個窗函數(shù)如矩形窗，窗內(nèi)數(shù)據(jù)并不改變。時域相乘即v(n)=x(n)·w(n)，加窗對頻域的影響，可用式（3-48）卷積公式表示卷積的結(jié)果，造成所得到的頻譜V(ejω)與原來的頻譜X(ejω)不相同，有失真。這種失真最主要的是造成頻譜的“擴散”（拖尾、變寬），這就是所謂的“頻譜泄漏”。

由上可知，泄漏是由于我們截取有限長信號所造成的。對具有單一譜線的正弦波來說，它必須是無限長的。也就是說，如果我們輸入信號是無限長的，那么FFT就能計算出完全正確的單一線頻譜?？墒俏覀儾豢赡苓@么做，而只能取有限長記錄樣本。如果在該有限長記錄樣本中，正弦信號又不是整數(shù)個周期時，就會產(chǎn)生泄漏。例如，一個周期為N=16的余弦信號x(n)=cos(6πn／16）截取一個周期長度的信號即x1(n)=cos(6πn/16)R16(n),其16點FFT的譜圖見3-24上圖，若截取的長度為13，則其16點FFT的譜圖見3-24下圖。由此可見，頻譜不再是單一的譜線，它的能量散布到整個頻譜的各處。這種能量散布到其他譜線位置的現(xiàn)象即為“頻譜泄漏”。

應(yīng)該說明，泄漏也會造成混疊，因為泄漏將會導(dǎo)致頻譜的擴展，從而使最高頻率有可能超過折疊頻率（fs/2），造成頻率響應(yīng)的混疊失真。泄漏造成的后果是降低頻譜的分辨率。此外，由于在主譜線兩邊形成很多旁瓣，引起不同頻率分量間的干擾(簡稱譜間干擾)，特別是強信號譜的旁瓣可能湮沒弱信號的主譜線，或者把強信號譜的旁瓣誤認(rèn)為是另一信號的譜線，從而造成假信號，這樣就會使譜分析產(chǎn)生較大偏差。

在進行FFT運算時，時域截斷是必然的，從而頻譜泄漏和譜間干擾也是不可避免的。為盡量減小泄漏和譜間干擾的影響，需增加窗的時域?qū)挾?頻域主瓣變窄)，但這又導(dǎo)致運算量及存儲量的增加；其次，數(shù)據(jù)不要突然截斷，也就是不要加矩形窗，而是加各種緩變的窗（例如：三角形窗、升余弦窗、改進的升余弦窗等），使得窗譜的旁瓣能量更小，卷積后造成的泄漏減小，這個問題在FIR濾波器設(shè)計一章（第6章）中將會討論到。圖3-24余弦信號頻譜泄漏示例圖3.7.3應(yīng)用實例

一個長度為N的時域離散序列x(n)，其離散傅里葉變換X(k)（離散頻譜）是由實部和虛部組成的復(fù)數(shù)，即（3-61）對實信號x(n)，其頻譜是共軛偶對稱的，故只要求出k在0，1，2，…,N/2上的X(k)即可。將X(k)寫成極坐標(biāo)形式（3-62）式中，|X(k)|稱為幅頻譜，arg［X(k)］稱為相頻譜。

實際中也常常用信號的功率譜表示，功率譜是幅頻譜的平方，功率譜PSD定義為（3-63）功率譜具有突出主頻率的特性，在分析帶有噪聲干擾的信號時特別有用。將式（3-62）繪成的圖形稱為頻譜圖。由頻譜圖可以知道信號存在哪些頻率分量，它們就是譜圖中峰值對應(yīng)的點。譜圖中最低頻率為k=0，對應(yīng)實際頻率為0（即直流）；最高頻率為k=N/2，對應(yīng)實際頻率為f=fs/2；

對處于0，1，2，…

N/2上的任意點k，對應(yīng)的實際頻率為f=kF=kfs/N。由于所取單位不同，頻率軸有幾種定標(biāo)方式。圖3-25列出頻率軸幾種定標(biāo)方式的對應(yīng)關(guān)系。上圖中f′為歸一化頻率，定義為f′無量綱，在歸一化頻率譜圖中，最高頻率為0.5。專用頻譜分析儀器常用歸一化頻率表示。(3-64)圖3-25模擬頻率與數(shù)字頻率之間的定標(biāo)關(guān)系

［例3-4］已知信號x(t)=0.15sin(2πf1t)+sin(2πf2t)-0.1sin(2πf3t)，其中,f1=1Hz，f2=2Hz，f3=3Hz。x(t)包含三個正弦波，但從時域波形圖3-26（a）來看，似乎是一個正弦信號，很難看到小信號的存在，因為它被大信號所掩蓋。取fs=32Hz作頻譜分析。

解因fs=32Hz，故

該信號為周期信號，其周期為N=32?，F(xiàn)對x(n)作32點的離散傅里葉變換(DFT)，其幅度特性|X(k)|如圖3-26(b)所示。圖中僅給出了k=0,1,…,15的結(jié)果。k=16,17,…,32的結(jié)果可由|X(N-k)|=|X(k)|得出。因N=32,故頻率分辨率F=fs/N=1Hz。由圖3-26(b)可知,k=1,2,3所對應(yīng)的頻譜即為頻率f1=1Hz，f2=2Hz,f3=3Hz的正弦波所對應(yīng)的頻譜。幅頻圖如圖3-26（b）,該圖中小信號成分清楚顯示出來?？梢娦⌒盘柍煞衷跁r域中很難辨識而在頻域中容易識別。圖3-26已知信號的時域圖和幅頻圖圖3-27已知信號的時域圖和幅頻圖

MATLAB中計算序列的離散傅里葉變換和逆變換是采用快速算法，利用fft和ifft函數(shù)實現(xiàn)。函數(shù)fft用來求序列的DFT,調(diào)用格式為：

[XK]=fft(xn,N)其中，xn為有限長序列，N為序列xn的長度，XK為序列xn的DFT.函數(shù)ifft用來求IDFT,調(diào)用格式為：

[xn]=ifft(XK,N)其中，XK為有限長序列，N為序列XK的長度，xn為序列XK的IDFT.3.8FFT的其他應(yīng)用3.8.1線性卷積的FFT算法——快速卷積

1.FFT的快速卷積法以FIR濾波器為例，因為它的輸出等于有限長單位脈沖響應(yīng)h(n)與有限長輸入信號x(n)的離散線性卷積。設(shè)x(n)為L點，h(n)為M點，輸出y(n)為

y(n)也是有限長序列，其點數(shù)為L+M-1點。下面首先討論直接計算線性卷積的運算量。由于每一個x(n)的輸入值都必須和全部的h(n)值相乘一次，因而總共需要LM次乘法，這就是直接計算的乘法次數(shù)，以md表示為md=LM

用FFT算法也就是用圓周卷積來代替這一線性卷積時，為了不產(chǎn)生混疊，其必要條件是使x(n)，h(n)都補零值點，補到至少N=M+L-1，即:（3-65）0≤n≤L-1L≤n≤N-10≤n≤M-1M≤n≤N-1然后計算圓周卷積N這時，y(n)就能代表線性卷積的結(jié)果。用FFT計算y(n)的步驟如下：(1)求H(k)=DFT［h(n)］，N點；(2)求X(k)=DFT［x(n)］,N點；(3)計算Y(k)=X(k)H(k);(4)求y(n)=IDFT［Y(k)］，N點。

步驟(1)、(2)、(4)都可用FFT來完成。此時的工作量如下：三次FFT運算共需要 次相乘，還有步驟③的N次相乘，因此共需要相乘次數(shù)為（3-66）

比較直接計算線性卷積(簡稱直接法)和FFT法計算線性卷積(簡稱FFT法)這兩種方法的乘法次數(shù)。設(shè)式（3-65）與式（3-66）的比值為Km，則（3-67）分兩種情況討論如下：（1）x(n)與h(n)點數(shù)差不多。例如，M=L，則N=2M-1≈2M，則這樣可得下表：M=L8163264128256512102420484096Km0.5720.9411.62.785.928.821629.2453.999.9

當(dāng)M=8時，F(xiàn)FT法的運算量大于直接法；當(dāng)M=32時，二者相當(dāng)；當(dāng)M=512時，F(xiàn)FT法運算速度可快16倍；當(dāng)M=4096時，F(xiàn)FT法約快100倍?？梢钥闯?，當(dāng)M=L且M超過32以后，M越長，F(xiàn)FT法的好處越明顯。因而將圓周卷積稱為快速卷積。

（2）當(dāng)x(n)的點數(shù)很多時，即當(dāng)L>>M。通常不允許等x(n)全部采集齊后再進行卷積;否則，使輸出相對于輸入有較長的延時。此外，若N=L+M-1太大，h(n)必須補很多個零值點，很不經(jīng)濟，且FFT的計算時間也要很長。這時FFT法的優(yōu)點就表現(xiàn)不出來了，因此需要采用分段卷積或稱分段過濾的辦法。即將x(n)分成點數(shù)和h(n)相仿的段，分別求出每段的卷積結(jié)果，然后用一定方式把它們合在一起，便得到總的輸出，其中每一段的卷積均采用FFT方法處理。有兩種分段卷積的辦法:重疊相加法和重疊保留法。2．重疊相加法設(shè)h(n)的點數(shù)為M，信號x(n)為很長的序列。我們將x(n)分解為很多段，每段為L點，L選擇成和M的數(shù)量級相同，用xi(n)表示x(n)的第i段：iL≤n≤(i+1)L-1其他n

i=0,1,…（3-68）則輸入序列可表示成（3-69）這樣，x(n)和h(n)的線性卷積等于各xi(n)與h(n)的線性卷積之和，即（3-70）

每一個xi(n)*h(n)都可用上面討論的快速卷積辦法來運算。由于xi(n)*h(n)為L+M-1點，故先對xi(n)及h(n)補零值點，補到N點。為便于利用基-2FFT算法，一般取N=2m≥L+M-1，然后作N點的圓周卷積：N

由于xi(n)為L點，而yi(n)為（L+M-1）點（設(shè)N=L+M-1），故相鄰兩段輸出序列必然有（M-1）個點發(fā)生重疊，即前一段的后（M-1）個點和后一段的前（M-1）個點相重疊，如圖3-27所示。按照式（3-70），應(yīng)該將重疊部分相加再和不重疊的部分共同組成輸出y(n)。圖3-28重疊相加法圖形和上面的討論一樣，用FFT法實現(xiàn)重疊相加法的步驟如下:

(1)計算N點FFT，H(k)=DFT［h(n)］；(2)計算N點FFT，Xi(k)=DFT［xi(n)］；(3)相乘，Yi(k)=Xi(k)H(k)；(4)計算N點IFFT，yi(n)=IDFT［Yi(k)］；(5)將各段yi(n)（包括重疊部分）相加， 。重疊相加的名稱是由于各輸出段的重疊部分相加而得名的。3．重疊保留法此方法與上述方法稍有不同。先將x(n)分段，每段L=N-M+1個點，這是相同的。不同之處是，序列中補零處不補零，而在每一段的前邊補上前一段保留下來的（M-1）個輸入序列值，組成L+M-1點序列xi(n)，如圖3-29(a）所示。如果L+M-1<2m，則可在每段序列末端補零值點，補到長度為2m，這時如果用DFT實現(xiàn)h(n)和xi(n)圓周卷積，則其每段圓周卷積結(jié)果的前（M-1）個點的值不等于線性卷積值，必須舍去。圖3-29重疊保留法示意圖

為了說明以上說法的