快速傅里葉變換

第四章快速傅里葉變換FastFourierTransform福建農(nóng)林大學金山學院信息與機電工程系(309428110@qq)數(shù)字信號處理DigitalSignalProcessing本章內容介紹傅里葉變換的一些快速算法快速算法的思想根據(jù)原始變換定義的運算規(guī)律，及其中某些算子的特殊性，找出減少乘法和加法運算次數(shù)的有效途徑，實現(xiàn)原始變換的各種高效算法。了解直接計算N點DFT的運算量了解減少運算量的根本途徑理解按時間抽取的基-2FFT算法的算法原理、運算流圖、所需計算量和算法特點了解按時間抽取的基-2FFT算法的編程思想理解按頻率抽取的基-2FFT算法的算法原理、運算流圖、所需計算量和算法特點本章學習目標離散傅里葉變換不僅具有明確的物理意義，更便于用計算機處理。但是，直至上個世紀六十年代，由于數(shù)字計算機的處理速度較低以及離散傅里葉變換的計算量較大，離散傅里葉變換長期得不到真正的應用，快速離散傅里葉變換算法的提出，才得以顯現(xiàn)出離散傅里葉變換的強大功能，并被廣泛地應用于各種數(shù)字信號處理系統(tǒng)中。近年來，計算機的處理速率有了驚人的開展，同時在數(shù)字信號處理領域出現(xiàn)了許多新的方法，但在許多應用中始終無法替代離散傅里葉變換及其快速算法。4.1引言FFT產(chǎn)生故事雖然頻譜分析和DFT運算很重要，但在很長一段時間里，由于DFT運算復雜，并沒有得到真正的運用，而頻譜分析仍大多采用模擬信號濾波的方法解決。當時加文(Garwin)在自已的研究中極需要一個計算傅里葉變換的快速方法。他注意到圖基(J.W.Turkey)正在寫有關傅里葉變換的文章，因此詳細詢問了圖基關于計算傅里葉變換的技術知識。圖基概括地對加文介紹了一種方法，它實質上就是后來的著名的庫利(CooleyJ.W)圖基算法。在加文的迫切要求下，庫利很快設計出一個計算機程序。1965年庫利--圖基在<計算數(shù)學>、MathematicofComputation雜志上發(fā)表了著名的“機器計算傅里級數(shù)的一種算法〞文章，提出一種快速計算DFT的方法和計算機程序--揭開了FFT開展史上的第一頁。促使FFT算法產(chǎn)生原因還有1967年至1968年間FFT的數(shù)字硬件制成，電子數(shù)字計算機的條件，使DFT的運算大簡化了。1984年，法國的杜梅爾(P.Dohamel)和霍爾曼(H.Hollmann)將基2分解和基4分解糅合在一起，提出了分裂基FFT算法。其運算量比前幾種算法都有所減少，運算流圖卻與基2FFT很接近，運算程序也很短。它是目前一種實用的高效新算法。4.2基2FFT算法有限長序列通過離散傅里葉變換(DFT)將其頻域離散化成有限長序列。但其計算量太大(與N的平方成正比〕,很難實時地處理問題，因此引出了快速傅里葉變換(FFT)。FFT并不是一種新的變換形式，它只是DFT的一種快速算法，并且根據(jù)對序列分解與選取方法的不同而產(chǎn)生了FFT的多種算法。FFT在離散傅里葉反變換、線性卷積和線性相關等方面也有重要應用。問題提出：

設有限長序列x(n)，非零值長度為N，計算對x(n)進行一次DFT運算，共需多大的運算工作量?4.2.1直接計算DFT的特點及減少

運算量的根本途徑其中x(n)為復數(shù)，也為復數(shù)，所以DFT與IDFT二者計算量相同。1.比較DFT與IDFT的運算量4.2.1直接計算DFT的特點及減少

運算量的根本途徑計算一個X(k)(一個頻率成分)值，運算量為例k=1那么要進行N次復數(shù)乘法，(N-1)次復數(shù)加法，所以，要完成整個DFT運算，其計算量為：N*N次復數(shù)相乘和N*(N-1)次復數(shù)加法4.2.1直接計算DFT的特點及減少

運算量的根本途徑2.以DFT為例計算復數(shù)運算量一個復數(shù)乘法包括4個實數(shù)乘法和2個實數(shù)加法。(a+jb)(c+jd)=(ac-bd)+j(bc+ad)所以所有X(k)就要4N2次實數(shù)乘法運算,2N2＋2N(N-1)＝N(4N-2)次實數(shù)加法運算。當N很大時,運算量將是驚人的,這樣,難以做到實時處理。

4次實數(shù)乘法2次實數(shù)加法4.2.1直接計算DFT的特點及減少

運算量的根本途徑3.一次復數(shù)乘法換算成實數(shù)運算量例1:當N=1024點時，直接計算DFT需要：N2=220=1048576次，即一百多萬次的復乘運算，這對實時性很強的信號處理(如雷達信號處理)來講，它對計算速度有十分苛刻的要求-->迫切需要改進DFT的計算方法，以減少總的運算次數(shù)。例2：石油勘探，24道記錄，每道波形記錄長度5秒，假設每秒抽樣500點/秒，每道總抽樣點數(shù)=500*5=2500點，24道總抽樣點數(shù)=24*2500=6萬點，DFT運算時間=N2=(60000)2=36*108次4.2.1直接計算DFT的特點及減少

運算量的根本途徑3.一次復數(shù)乘法換算成實數(shù)運算量FFT算法對于N點DFT,僅需(N/2)log2N

次復數(shù)乘法運算和Nlog2N次復數(shù)加法。4.2.1直接計算DFT的特點及減少

運算量的根本途徑4、FFT的計算工作量如果計算機的速度為平均每次復數(shù)乘需要5×10-6秒，每次復加需要1×10-6秒，用來計算N=1024點DFT,問1)直接計算需要多少時間？2〕用FFT算法計算需要多少時間？4.2.1直接計算DFT的特點及減少

運算量的根本途徑4、FFT的計算工作量利用的特性4.2.1直接計算DFT的特點及減少

運算量的根本途徑利用以上特性，得到改進DFT直接算法的方法。4.2.1直接計算DFT的特點及減少

運算量的根本途徑

例快速傅里葉變換(FFT)就是在此特性根底上發(fā)展起來的：〔1〕利用DFT系數(shù)的對稱性和周期性,合并DFT運算中的某些項；〔2〕將長序列分解為短序列,從而減少其運算量。

因合并與分解方法的不同產(chǎn)生了多種DFT的快速算法。4.2.1直接計算DFT的特點及減少

運算量的根本途徑5、DFT的根本思想1.按抽取方法分：時間抽取法(DITDecimation-In-Time);頻率抽取法(DIFDecimation-In-Frequency)2.按“基數(shù)〞分：

基-2FFT算法；基-4FFT算法；混合基FFT算法；分裂基FFT算法3.其它方法：

線性調頻Z變換(Chrip-z法)4.2.1直接計算DFT的特點及減少

運算量的根本途徑6、FFT算法分類：1、時域抽取算法原理設輸入序列長度為N=2M(M為正整數(shù)〕，將該序列按時間順序的奇偶分解為越來越短的子序列，稱為基2按時間抽取的FFT算法。其中基數(shù)2----N=2M，M為整數(shù)。假設不滿足這個條件，可以人為地加上假設干零值〔加零補長〕使其到達N=2M。4.2.2時域抽取法基2FFT根本原理Decimation-in-Time(DIT)設一序列x(n)的長度為L=9，應加零補長為N=24=16

應補7個零值。0123456789101112131415nx(n)1、時域抽取算法原理4.2.2時域抽取法基2FFT根本原理Decimation-in-Time(DIT)設序列點數(shù)N=2M，M為整數(shù)。假設不滿足，那么補零2、算法步驟4.2.2時域抽取法基2FFT根本原理Decimation-in-Time(DIT)4.2.2時域抽取法基2FFT根本原理Decimation-in-Time(DIT)2、算法步驟4.2.2時域抽取法基2FFT根本原理Decimation-in-Time(DIT)3、蝶形運算1111-1(1)N/2點的DFT運算量:復乘次數(shù): 復加次數(shù):(2)兩個N/2點的DFT運算量:復乘次數(shù): 復加次數(shù):(3)N/2個蝶形運算的運算量:復乘次數(shù):

復加次數(shù):

復乘:復加:總共運算量:可見，N點DFT的復乘為N2，復加N(N-1)；與分解后相比可知，計算工作點差不多減少

一半。4.2.2時域抽取法基2FFT根本原理Decimation-in-Time(DIT)例8點FFT的算法首先可以分解為兩個N/2=4點的DFT。具體方法如下:X(5)X(4)~X(7)同學們自已寫〔1〕N=8點分解成2個4點的DFT的信號流圖

x1(0)=x(0)

x1(1)=x(2)N/2點

x1(2)=x(4)DFT

x1(3)=x(6)x2(0)=x(1)

x2(1)=x(3)N/2點

x2(2)=x(5)DFT

x2(3)=x(7)X1(0)X1(1)X1(2)X1(3)X2(0)X2(1)X2(2)X2(3)WN2NW1NW0-1-1-1X(0)X(1)X(2)X(3)X(4)X(6)X(7)WN3-1由于N=2L,所以N/2仍為偶數(shù),可以進一步把每個N/2點的序列再按其奇偶局部分解為兩個N/4的子序列。(2)N/2(4點)-->N/4(2點)FFT(2)N/2(4點)-->N/4(2點)FFT02NN02NN-1-1-1-1-1-1WWWWX(0)X(1)X(2)X(3)X(4)X(5)X(6)X(7)N/4DFTN/4DFTN/4DFTN/4DFTx(0)x(2)x(6)x(1)x(5)x(3)x(7)x(4)WN0WN1WN2WN3X1(0)X1(1)X1(2)X1(3)X2(0)X2(1)X2(2)X2(3)x3(0)x3(1)x4(0)x4(1)x5(0)x5(1)x6(0)x6(1)-1-1最后剩下兩點DFT,它可分解成兩個一點DFT，但一點DFT就等于輸入信號本身，所以兩點DFT可用一個蝶形結表示。取x(0)、x(4)為例。(3)N/4〔2點〕-->2個1點DFT1點DFTx(0)1點DFTx(4)X3(0)X3(1)2個1點的DFT蝶形流圖-1(3)N/4〔2點〕-->2個1點DFT一個完整的按時間抽取的8點FFT1〕輸入倒位序，輸出自然序自然順序n二進制n2n1n0

倒位序二進制n0n1n2倒位順序n

00

0

00

0

00

10

0

11

0

04

20

1

00

1

02

30

1

11

1

06

41

0

00

0

11

51

0

11

0

15

61

1

00

1

13

71

1

11

1

174.2.3DIT-FFT的運算規(guī)律及流程圖畫法對N=2L點FFT，共需L級蝶形運算，每級有N/2個蝶形運算組成，蝶形運算兩節(jié)點的距離:2L-1〔L表示級數(shù)〕每個蝶形運算有一次復乘和再次復加。如1111-14.2.3DIT-FFT的運算規(guī)律及流程圖畫法2〕蝶形運算4.2.3DIT-FFT的運算規(guī)律及流程圖畫法3〕WNr確實定(僅給出方法)-1-1-1-1X(0)X(1)X(2)X(3)X(4)X(5)X(6)X(7)x(0)x(2)x(6)x(1)x(5)x(3)x(7)x(4)000，000001，100010，010011，110100，001101，101110，011111，111x3(0)x3(1)x4(0)

x4(1)x5(0)x5(1)x6(0)x6(1)WN0WN2WN0WN2X1(0)X1(1)X1(2)X1(3)X2(0)X2(1)X2(2)X2(3)WN0WN1WN2WN3-1-1-1-1WN0WN0WN0WN0N=84.2.3DIT-FFT的運算規(guī)律及流程圖畫法-1-1-1-1X(0)X(1)X(2)X(3)X(4)X(5)X(6)X(7)X(8)X(9)X(10)X(11)X(12)X(13)X(14)X(15)0000，00000001，10000010，01000011，11000100，00100101，10100110，01100111，11101000，00011001，10011010，01011011，11011100，00111101，10111110，01111111，1111x(0)x(8)x(4)x(12)x(2)x(10)x(6)x(14)x(1)x(9)x(5)x(13)x(3)x(11)x(7)x(15)X1(0)X1(1)X1(2)X1(3)X1(4)X1(5)X1(6)X1(7)X2(0)X2(1)X2(2)X2(3)X2(4)X2(5)X2(6)X2(7)X3(0)X3(1)X3(2)X3(3)X4(0)X4(1)X4(2)X4(3)X5(0)X5(1)X5(2)X5(3)X6(0)X6(1)X6(2)X6(3)X7(0)X7(1)X8(0)X8(1)X9(0)X9(1)X10(0)X10(1)X11(0)X11(1)X12(0)X12(1)X13(0)X13(1)X14(0)X14(1)WN0WN1WN2WN3WN4WN5WN6WN7WN0WN2WN4WN6WN0WN2WN4WN6WN0WN4WN0WN4WN0WN4WN0WN4WN0WN0WN0WN0WN0WN0WN0WN0N=164.2.3DIT-FFT的運算規(guī)律及流程圖畫法由上述分析可知,N=8需三級蝶形運算N=23=8,由此可知,N=2L共需L級蝶形運算,而且每級都由N/2個蝶形運算組成,每個蝶形運算有一次復乘,兩次復加。因此,N點的FFT的運算量為復乘:mF=〔N/2〕L=〔N/2〕log2N復加:aF=NL=Nlog2N

由于計算機的乘法運算比加法運算所需的時間多得多，故以乘法作為比較基準。4.2.4DIT-FFT算法與直接計算DFT

運算量的比較4.2.4DIT-FFT算法與直接計算DFT

運算量的比較1、算法原理設序列點數(shù)N=2L，L為整數(shù)。將X(k)按k的奇偶分組前，先將輸入x(n)按n的順序分成前后兩半：4.2.5按頻率抽取的FFT1、算法原理

4.2.5按頻率抽取的FFT當k為偶數(shù),即k=2r時,(-1)k=1；當k為奇數(shù),即k=2r+1時(-1)k=-1。這時X(k)可分為兩局部：2、N點DFT按k的奇偶分組可分為兩個N/2的DFT

4.2.5按頻率抽取的FFT

可見,上面兩式均為N/2的DFT。3、蝶形運算

-14.2.5按頻率抽取的FFT4.N=8時,按k的奇偶分解過程先蝶形運算，后DFT:再將N/2點DFT按k的奇偶分解為兩個N/4點的DFT,如此進行下去,直至分解為2點DFT。-1-1-1-1WWWWNNNN01234.2.5按頻率抽取的FFT例如

N=8時頻域抽取的FFT流圖如下：4.2.5按頻率抽取的FFT

x(0)X(0)x(1)X(4)x(2)X(2)x(3)X(6)x(4)X(1)x(5)X(5)x(6)X(3)x(7)X(7)-1-1-1WWWWNNNN0123-1-1-1WWWWNNNN0202-1-1-1WWWWNNNN0000-1-1-11〕根本蝶形不同時域：先復乘后加減頻域：先減后復乘2〕運算量相同3〕都可原位運算4〕兩者根本蝶形互為轉置5、時域抽取和頻域抽取FFT的區(qū)別

4.2.5按頻率抽取的FFT比較兩式可知,只要DFT的每個系數(shù)換成,最后再乘以常數(shù)1/N就可以得到IDFT的快速算法--IFFT。4.2.6用FFT計算IDFT4.2.6用FFT計算IDFT這就是說，先將X(k)取共軛，即將X(k)的虛部乘-1，直接利用FFT程序計算DFT；然后再取一次共軛；最后再乘1/N,即得x(n)。所以FFT，IFFT可用一個子程序。共軛FFT共軛乘1/N直接調用FFT子程序計算IFFT的方法：4.2.6用FFT計算IDFT4.3線性卷積的FFT算法4.3線性卷積的FFT算法4.3線性卷積的FFT算法4.3線性卷積的FFT算法利用FFT計算線性卷積

MATLAB程序N=1024x=[23145ones(1,N)];h=[21745723ones(1,N)];Lenx=length(x);%求序列x的長度Lenh=length(h);%求序列h的長度N=Lenx+Lenh-1;t=cputime(%或tic)Xk=fft(x,N);%計算x序列的DFTHk=fft(h,N);%計算h序列的DFTYk=Xk.*Hk;y=ifft(Yk)%求IDFTt1=cputime-t〔或toc〕stem(y);xlabel('n');4.3線性卷積的FFT算法4.3線性卷積的FFT