數字信號處理第四章-快速傅里葉變換

1FFT:FastFourierTransform1965年，JamesW.Cooley和JohnW.Tukey在《計算數學》（《MathematicsofComputation》）上發(fā)表了“一種用機器計算復序列傅立葉級數的算法（AnalgorithmforthemachinecalculationofcomplexFourierseries）”論文。自此之后，新的算法不斷涌現。一種是對N等于2的整數次冪的算法，如基2算法，基4算法。另一種是N不等于2的整數次冪的算法，例如Winagrad算法，素因子算法。4.1FFT：引言Dr.JamesW.CooleyWorkedatIBMWatsonResearchCenterinYorktownHeights,N.Y..AfterhisretirementfromIBMin1991,hejoinedtheElectricalEngineeringDepartmentattheUniversityofRhodeIsland.24.2FFT：直接計算DFT的問題及改進N點有限長序列x(n)的DFT變換對的定義為：

其中假設x(n)是復序列，同時X(k)一般也是復數。3復數乘法復數加法每一個X(k)NN–1N個X(k)(N點DFT)N2N(N–1)實數乘法實數加法一次復乘42一次復加2每一個X(k)4N2N+2(N–1)=2(2N–1)N個X(k)(N點DFT)4N22N(2N–1)4.2FFT：直接計算DFT的問題及改進復乘的加法次數復加的加法次數4如N=512、1024和8192時，DFT的乘法運算

N2=5122=218=262144（26萬次）

N2=10242=220=1048576（105萬次）

N2=81922=226=67108864（6千7百萬次）對于大N，在實際中是不能接受的，無法“實時”應用DFT。一般地，在計算機中，一次加法比一次乘法所需時間要短;在DSP中，由于乘法用特殊的硬件電路專門完成，因此乘法和加法所需機器周期相同。Cooley與Turkey提出的FFT算法，大大減少了計算次數。如N=512時，FFT的乘法次數約為2000次，提高了約128倍，而且簡化隨N的增加而巨增，因而，用數值方法計算頻譜得到實際應用。4.2FFT：直接計算DFT的問題及改進5FFT：WNkn

的性質x(0)x(2)x(1)x(3)-1-1-1-1W41X(0)X(1)X(2)X(3)6以4點DFT為例：直接計算需要：42=16

次復數乘。而按周期性及對稱性，可以將DFT表示為：只需要

1

次復數乘FFT：WNkn

的性質7FFT算法分類:時間抽選法DIT:Decimation-In-Time頻率抽選法DIF:Decimation-In-FrequencyFFT：基本思想基本思路：雖然存在不同的FFT方法，但其核心思想大致相同，即通過迭代，反復利用低點數的DFT完成高點數的DFT計算，以此達到降低運算量的目的。迭代：利用WNkn

的周期性、特殊點和對稱性，合并DFT計算中很多重復的計算，達到降低運算量的目的。低點數：將傅氏變換DFT分解成相繼小的DFT計算，即N變小，而計算量與N2

成正比。8算法原理設序列點數N=2M，M為整數。若不滿足，則補零。

將序列x(n)按n的奇偶分成兩組：N為2的整數冪的FFT算法稱基-2FFT算法。即一組由偶數序號組成，另一組由奇數序號組成。4.3FFT：基2時間抽選法注意其長度為N/2910偶數取樣點DFT為：

奇數取樣點DFT為：

①k的整個范圍為0～(N-1)，而X1(k)、X2(k)是由N/2個樣點形成的DFT，x(2r)和x(2r＋1）的長度為N/2；②由這兩個偶數和奇數N/2個時域樣值可以計算出前N/2個DFT系數，也可以計算出后N/2個DFT系數。③問題：這前后N/2個DFT有無關系？k取N/2～(N-1)時，X1(k)、X2(k)、WN

情況如何？11觀察則12因此：整個X(k)的計算，可以分解為前、后半部分的運算。而只要求出前一半，就可以由上式求出整個序列。同理而因此13上式表示為信號流程圖：此信號流程圖也稱為蝶形流程圖偶數點的DFT奇數點的DFT14以N=8為例，其信號流程圖：偶數序列奇數序列15N/2仍為偶數，進一步分解：N/2N/416同理17以N=8為例，其信號流程圖為N/2點DFTN/2點DFT18由于N=2M，這樣逐級分解，直到2點DFT當N=8時，即分解到X3(k)，X4(k)，X5(k)，X6(k)，k=0,119每一次分解都是按時間域輸入序列的奇偶性次序，分解為兩個半長的子序列，所以稱為按時間抽取法。仍以N=8為例：注：復數乘法5次，原來64次；復數加法為24次，原來56次。m=0級m=1級m=2級20

2m

點DFT的基2時間抽選計算過程可見：為什么引入FFT？出發(fā)點：考慮W因子的特點和性質，簡化算法。DIT-FFT算法：時間域輸入信號逐級分解為奇偶序列。上式中位序調整第0級蝶形復合第1級蝶形復合第m-1級蝶形復合x(n)……X(k)21例

使用基2時間抽選法FFT流圖，計算下列數據的8點DFT：x=[4,-3,2,0,-1,-2,3,1]4-320-1-23142-13-30-214-123-3-201355-1-5-11-1355j-5-11j85+j-25-j-4-1+j-6-1-j85+j-25-j-4j√26jj√245+j+j√2-2+6j5-j+j√2125+j-j√2-2-6j5-j-j√222蝶形運算在FFT信號流圖中每一級里節(jié)點都是成對存在的，計算時總是下節(jié)點的值乘以WNr，然后與上節(jié)點的值相加減，形成下一列兩個節(jié)點值。這種計算的基本關系是：式中p,q是上下節(jié)點的序號。每一級中每對節(jié)點稱為對偶節(jié)點，對偶節(jié)點的間距為：基2時間抽選法-蝶形運算23同址計算（同位特性）基2時間抽選法-同址計算1)不同級數的節(jié)點序號不變從蝶形運算可以看出第m列的xm(p),xm(q)經蝶形運算后，在第(m+1)列中的節(jié)點序號不變，即xm+1(p)中的p與xm+1(q)中的q仍分別是xm(p),xm(q)中的p和q值。242)蝶形運算是自成獨立單元蝶形運算自成獨立單元，即xm+1(p)、xm+1(q)只與xm(p)、xm(q)有關，而與其他節(jié)點的值無關，同時xm(p)、xm(q)也不參與另外的蝶形運算，后面的運算也不再用到前面的數據。因此，就可把計算結果xm+1(p),xm+1(q)放入計算前xm(p)、xm(q)的存儲單元中，而不用再建新的存儲單元—同址計算。同址計算所需存儲量僅等于給定數據所需的存儲量，這可大大節(jié)省存儲單元?；?時間抽選法-同址計算25輸入位序重排N點DFT分解為兩個N/2點DFT輸入序列按奇偶分組再分解再奇偶重排直到2點DFT。即FFT輸出自然序列輸入序列x(n)位序重排?；?時間抽選法-位序重排010011n0n1n226說明：首先把x(n)中的n表示為二進制。假定N=8，則 x(n)→x(n2n1n0)ni=0,1FFT的時間抽選法按n的奇偶分離而形成位置重排的原理如上圖所示。由此可以看出，時間抽選法FFT的輸入位序重排，輸出分成上下兩部分，輸出結果為自然順序。實際操作輸入序列重排實際上就是完成二進制前后位序的相互交換：基2時間抽選法-位序重排27當N=2M時，共有M=log2N級蝶形；每級N/2個蝶形；每個蝶形有1次復數乘法，2次復數加法。復數乘法：復數加法：比較DFT/FFT

FFT：基2時間抽選法-運算量28FFT算法上面討論的FFT算法假定N為2的整數次冪，即N=2M，是基2FFT算法。當N=Rv時，這樣的算法叫做基RFFT算法，而當N=R1v1R2v2R3v3時，叫做混合基算法。當N是一個高度合數時，可得到最有效的算法。最受歡迎也最易編程的算法是基2FFT算法。對于不能分解的質數，或者當N不是2的整數次冪時，可以在信號的末尾補0，使其成為高度合數或2的整數次冪。MatlabFFT實現Matlab提供了內建的X=fft(x,N)

函數來計算矢量x的DFT。fft函數是機器碼寫成的，而不是以Matlab指令寫成的，即不存在.m文件。因此，它的執(zhí)行速度很快。采用混合算法。若N為2的冪，則得到高速的基2FFT算法；若N不是2的冪，則將N分解成質數，得到較慢的混合基FFT；若N為質數，則fft函數采用原始DFT算法。FFT：基2時間抽選法-Matlab實現29DFT定義表示為矩陣形式為：FFT：基2時間抽選法-矩陣表示30令：矩陣WN

為：DFT矩陣簡單地表示為：特性FFT：基2時間抽選法-矩陣表示31FFT基2時間抽選法信號流圖（N=8）-1W821×8點2×4點4×2點x(0)x(1)x(2)x(3)x(4)x(5)x(6)x(7)x(0)x(4)x(2)x(6)x(1)x(5)x(3)x(7)X(0)X(1)X(2)X(3)X(4)X(5)X(6)X(7)-1-1-1-1-1-1-1-1-1-1-1W80W80W82W80W81W82W83C(0)C(1)D(0)D(1)E(0)E(1)F(0)F(1)A(0)A(1)A(2)A(3)B(0)B(1)B(2)B(3)m=0級m=1級m=2級位序重排FFT：基2時間抽選法-矩陣表示32FFT矩陣形式表示為：1)輸入矢量x與P8

相乘，則只是輸入的位序重排，沒有乘法操作。FFT：基2時間抽選法-矩陣表示332)第0級運算：2點DFTFFT：基2時間抽選法-矩陣表示343)第1級運算：4點DFTFFT：基2時間抽選法-矩陣表示354)第2級運算：8點DFTFFT：基2時間抽選法-矩陣表示36FFT算法看成是DFT矩陣W8

的分解：這個因式分解對應于快速算法，因為矩陣F8(8)，F8(4)

，F8(2)和P8

的大部分元素為0，是稀疏矩陣。因此上述每個矩陣的乘法運算最多只需要8次復乘運算，而P8只是置換操作，沒有乘法操作。不同的FFT算法對應于將DFT矩陣WN

分解成不同的稀疏矩陣。FFT：基2時間抽選法-矩陣表示37算法原理：根據時間-頻率的對偶性時間抽選法：是把輸入x(n)按奇偶分解成兩個子序列，即N點x(n)序列→N/2點子序列，而輸出X(k)是按自然順序排列的。頻率抽選法：是把輸入x(n)按照前后對半分開，而不是奇偶數分開，而輸出X(k)逐項分解成偶數點子序列和奇數點子序列。DFT變換表達式為：4.4FFT：基2頻率抽選法(DIF-FFT)如果將輸入x(n)按前后等分，即將求和分成兩部分，范圍分別為：38基2頻率抽選法–算法原理39

按k的奇偶將X(k)分成兩部分：40令:則X(2r)和X(2r+1)分別是x1(n)和x2(n)的N/2點DFT，記為X1(k)和X2(k)。41x1(0)x1(1)-1x1(2)x1(3)-1x2(0)x2(1)-1x2(2)x2(3)-1N/2點DFTN/2點DFTx(0)x(7)x(1)x(2)x(3)x(4)x(5)x(6)X1(0)=X(0)X2(0)=X(1)X1(1)=X(2)X1(2)=X(4)X1(3)=X(6)X2(1)=X(3)X2(2)=X(5)X2(3)=X(7)42FFT基2頻率抽選法信號流圖（N=8）逐級分解，直到2點DFT43基本蝶形運算不同DIT:先復乘后加減，W因子在上下節(jié)點都有體現DIF:先減后復乘，W因子僅體現在下節(jié)點FFT：基2時間和頻率抽選法的異同44運算量相同同址計算FFT：基2時間和頻率抽選法的異同DIF輸出位序重排，DIT輸入位序重排45DIT和DIF的基本蝶形互為信號流圖轉置DITDIFFFT：基2時間和頻率抽選法的異同46輸入倒位序，輸出自然序（DIT）輸入自然序，輸出倒位序（DIF）輸入輸出均自然序各級具有相同幾何形狀輸入倒位序，輸出自然序輸入自然序，輸出倒位序FFT：基2抽選法-其它形式的流圖47

FFT：基2抽選法-其它形式的流圖48FFT：基2抽選法-其它形式的流圖49FFT：基2抽選法-其它形式的流圖50FFT：基2抽選法-其它形式的流圖51FFT：基2抽選法-其它形式的流圖52DFT和IDFT的定義：4.5FFT：IDFT快速算法(IFFT)DFT和IDFT的區(qū)別：①因子W的指數相差一個負號；②相差一個因子1/N。FFT算法中的分組方式、排序方式以及蝶形運算結構都可用于IFFT算法的設計，而這就是可依據現有的FFT算法直接得出IFFT算法的原因。53FFT和IFFT基本蝶形運算之間的關系設有序列x(n)，其DFT為X(k)，則IDFT為:

在FFT的時間抽選法中:

4.5FFT：IDFT快速算法(IFFT)對于IFFT算法，輸入是X(k)和X(k+N/2)，輸出是X1(k)和X2(k)。解上式可以得到X1(K)、X2(K)：X(k)X(k+N/2)X1(k)X2(k)-1548點DIT-IFFT算法-1x(0)x(1)x(2)x(3)x(4)x(5)x(6)x(7)x(0)x(4)x(2)x(6)x(1)x(5)x(3)x(7)X(0)X(1)X(2)X(3)X(4)X(5)X(6)X(7)-1-1-1-1-1-1-1-1-10.5W8-00.5W8-20.5W8-00.5W8-10.5W8-20.5W8-3-1-10.5W8-00.5W8-20.50.50.50.50.50.50.50.50.50.50.50.5說明：1）分組過程是按時間序列x(n)的奇偶性在時域上展開的，故稱此法為時間抽選算法DIT-IFFT；2）1/N的分解，N=2m

。4.5FFT：IDFT快速算法(IFFT)558點DIF-IFFT算法4.5FFT：IDFT快速算法(IFFT)56用FFT程序求IFFT的方法直接法：利用DIT、DIF的FFT程序，改變參量①把X(k)作為輸入序列，而輸出序列就是x(n)；②把因子WNkn

改為WN-kn

；③輸入序列的每一個元素除以N。4.5FFT：IDFT快速算法(IFFT)57流程如下：

共軛法共軛FFT共軛乘1/N4.5FFT：IDFT快速算法(IFFT)58DFT不適用于：DFT是均勻分布在Z平面單位圓上N點處的頻譜，如果取樣點不均勻時，則很麻煩。只研究信號的某一頻段，要求對該頻段取樣密集，提高分辨率(如窄帶信號的頻譜分析）；研究非單位圓上的取樣值（如頻譜峰值探測）；需要準確計算N點DFT，且N為大的素數；當x(n)是短時間序列時，則得到的頻率分辨率2pi/N是很低的。提高頻譜密度的辦法：用補零的方法增加點數，但DFT的點數又大大增加，使計算工作量增大。CZT算法：對z變換采用螺旋線取樣，chirp-z變換(線性調頻z變換，Chirpz-transform)4.9線性調頻z變換CZT及其快速算法59CZT的定義設x(n)為已知時間序列，其Z變換的形式為：式中復變量z為：這里s為拉普拉斯變量，是個實數CZT：定義60

按照上式計算Z[x(n)]必然是從z的實軸開始，為得到任意起始點和以螺旋線規(guī)律變化的z值，做如下假設：

其中，A、W為任意復數，θ0為A的起始角（第1個取樣點（k=0）），A0為A起始半徑，φ0

為在Z平面中相鄰的zk

（即zk

和zk+1

）之間的夾角，W0

為任意正數值。M為要分析的點數，不一定等于N。61CZT表達式1）取樣點沿螺旋線按角度間隔φ0分布，φ0>0時，取樣軌跡逆時針旋轉；φ0<0時，取樣軌跡逆時針旋轉。2）zk

周線是一條螺旋線：W0

表示螺旋線的伸展率；W0>1，隨著k的增加，向內盤旋，朝向原點；W0<1，隨著k的增加，向外盤旋；3)若W0=1，一段圓弧，若同時A0=1，則為單位圓一部分，此時CZT[x(n)]=DFT[x(n)]（M=N）。62CZT表示為線性卷積形式：利用布魯斯坦(Bluestein)提出的等式

把在CZT[x(n)]中的Wnk

因子展開，得：CZT：快速算法計算量：直接計算M點的CZT，需要MN次復數乘法，M(N-1)次復數加法，需要快速算法。把上述運算轉換為線性卷積形式，從而可以采用FFT算法，提高運算速度。63CZT：快速算法令則式中64CZT：快速算法g(n)和h(n)的取值范圍：65CZT：快速算法h(n)x(n)g(n)X(zk)序列：可以想像為頻率隨時間（n）線性增長的復指數序列，在雷達中這類信號，稱為線性調頻信號（chirpSignal)，因此，這里z變換稱為線性調頻z變換。該信號的瞬時頻率Ωi=2Ωt正比于時間t，因此，該信號的頻率隨時間線性增長。類比：66取值范圍n：0～N-1，k：0～M-1因果序列g(n)的長度為N，而序列h(n)的長度為N+M-1，并位于區(qū)間內-(N-1)~(M-1)，所以卷積后的序列y(n)的長度為2N+M-2，且位于區(qū)間內–(N-1)～（N＋M－2）。實際上只要得到區(qū)間0～(M-1）內的線性卷積結果就能夠完成CZT的計算。CZT：快速算法67用循環(huán)卷積代替線性卷積且不產生混迭失真的條件是：循環(huán)卷積的點數（周期）應大于或等于2N＋M－2，但CZT只需要前0～(M－1)值，對以后的其它值是否混迭并不關心。因此，可將循環(huán)卷積的點數縮減到最小為N+M-1。為了基2FFT運算，取L>=N+M-1，L=2m

。CZT：快速算法68g(n)、h(n)序列的構造為了進行循環(huán)卷積，把h(n)以周期L進行周期拓展，再取主值序列，從而得到圓周卷積的一個序列。從n=M開始補L－（N＋M－1）個零或任意值，補到n=L-N處。從n=L－N＋1到L－1，取h(n)的周期拓展序列h(L-n)。進行循環(huán)卷積的另一個序列g(n)只需要補零到L即可。69h(n)的構造：70CZT快速算法下面說明上述計算方案的具體實現，既然X(zk)可用線性卷積表示，我們就可以用DFT來計算線性卷積部分，只不過DFT換成FFT，從而得到CZT的快速算法。L點DFTH(k)h(n)構造長度M＋N-1長度Ng(n)補零L點DFTG(k)G(k)H(k)L=N+M-1L點IDFTy(k)7172線性卷積求解小結：時域直接求解

補N-N1個零x(n)N點DFT補N-N2個零h(n)N點DFTN點IDFTy(n)=x(n)*h(n)z變換法DFT法4.10FFT的應用：線性卷積求解73上面介紹的是兩個有限長序列x1(n)、x2(n)的線性卷積。但有時其中某個序列會很長或無限長，若等長序列存儲或輸入完后再做卷積運算，將產生問題：

(1)存儲量過大，運算量也太大；(2)等待輸入的時間很長，引起不能忍受的延遲；(3)采用DFT/FFT快速算法的效率降低具體方法：重疊保留法、重疊相加法

解決此問題的思路：

把長序列分段，每一分段分別與短序列進行卷積——分段卷積。FFT的應用：分段卷積74誤差分析在分析分段卷積之前，我們首先分析兩個有限長度序列的循環(huán)卷積的誤差情況。設x(n)為N1

點序列，h(n)為N2

點序列，由前面分析知道，線性卷積y(n)和循環(huán)卷積yC(n)關系為：DFT的應用：分段卷積一般地講，循環(huán)卷積是線性卷積的一種混疊形式，即是線性卷積以周期N延拓后取主值周期。但當N≥N1+N2-1，沒有混迭產生，此時線性卷積y(n)和循環(huán)卷積yN(n)相等。75為了用DFT計算線性卷積，必須選擇適當的N，但實際中卻可能做不到這一點，尤其是當一個序列的長度很大而存儲空間有限時。當選擇比所需要的小的N值進行循環(huán)卷積時，會引入誤差。下面我們計算此誤差的大小。首先N≥max(N1,N2)，假設此時，循環(huán)卷積yc(n)的取值范圍為而線性卷積y(n)的非零范圍為76有上式可知，誤差為：77

因為max(N1,N2)≤N≤（N1＋N2－1），上面的求和式中只剩下對應于r=±1的項，其它的r值超出了max(N1,N2)≤N≤（N1＋N2－1）。因此，

一般來講，如果x(n)和h(n)為因果序列，則y(n)也為因果序列，即：因此78它說明當

時，樣本n處的誤差與線性卷積中第n+N

個樣本相等，即混疊部分的樣值。（N1+N2–1）個樣本后，線性卷積結果為零。這說明循環(huán)卷積中的頭L個樣本存在誤差（混疊），混疊長度L等于線性卷積長度和循環(huán)卷積長度之差，而剩下的則為正確的線性卷積樣本。

0循環(huán)卷積長度N(N1+N2-1)-Ny(n)線性卷積長度N1+N2-1循環(huán)卷積的錯誤樣本（混疊部分）79此結論在用分段處理法實現長卷積時是非常有用的。錯誤樣本數=線性卷積長度–

循環(huán)卷積長度N1N20假設N=N2(N1+N2-1)-N=N1-1e(n)=0e(n)=y(n+N)y(n)yC(n)h(n)x(n)N1+N2-1N1-1點有誤差與線性卷積相同

假設L=N1-180例設x1(n)和x2(n)是兩個四點序列：

x1(n)={1,2,2,1}，x2(n)={1,-1,-1,1}計算當N=6，5，4時的循環(huán)卷積，并在每種情況下，驗證它們的誤差。解：顯然，線性卷積為7點x3(n)={1,1,-1,-2,-1,1,1}

當N=6，得到6點序列的循環(huán)卷積為：因此81當N=5時，得到5點序列的循環(huán)卷積為當N=4時，得到4點序列的循環(huán)卷積為82重疊保留法——對輸入x(n)預處理基本思路假設把序列x(n)分成多段N點序列，一個系統(tǒng)的脈沖響應為M點序列，M<N。由上面的結論可知，輸入分段序列和脈沖響應之間的N點循環(huán)卷積產生該段的輸出序列，其中前M-1個樣本不是正確的輸出值。若將x(n)簡單地分成互不重疊的各段，則所得的輸出序列會有不正確樣本區(qū)間存在。輸入數據如何分段？結果如何處理？DFT的應用：重疊保留法MNy

(n)h(n)x(n)M-1點誤差正確的點MNM-1點誤差正確的點83為了解決這個問題，將x(n)分為長度N1

的分段，然后每段再往前多?。∕-1）個樣本，即第k段的最后M-1點保留下來作為第(k+1)段的前M-1點，各段長變?yōu)椋?/p>

其中最前面的一段x0(n)，在其前面補上M-1個零，使其長度也為N。其中任一段xk(n)與h(n)進行N點卷積運算。

循環(huán)卷積：DFT的應用：重疊保留法相應的周期卷積 的周期為：N=N1+M-1線性卷積：84yk(n)的長度為：N1N’M-1M-1x0(n)nN重疊區(qū)域線性卷積循環(huán)卷積N1N’NM-1M-1n重疊區(qū)域線性卷積x1(n)錯誤樣本數=線性卷積長度–

循環(huán)卷積長度

=[N1+2(M-1)]–[N1+(M-1)]=M-185將每一段所得到的循環(huán)卷積的前M-1個值去掉，保留后面的N1

個值，再將各段保留的N1個值前后拼接起來，就得到所要求的線性卷積y(n)。因此，循環(huán)卷積

yk’(n)的前M-1個值是線性卷積yk(n)前面的M-1個值與yk-1(n)后面M-1個值的混迭，是錯誤的樣本。

循環(huán)卷積yk’(n)的后N1個值才與yk(n)中間的

N1個值相同。86

這種通過將與相重疊的部分舍去，即舍去循環(huán)卷積yk’(n)中前M-1項，保留后面的N1

個數值，來得到所求結果的方法稱為重疊保留法，有的書中又稱為重疊舍去法。8788例設x(n)=(n+1),0≤n≤9,h(n)={1,0,-1}，按N=6用重疊舍去法計算

y(n)=x(n)*h(n)解：由于M=3，必須使每一段與前一段重疊兩個樣本，x(n)為10點序列，需要在開頭加(M-1)=2個零。因為N=