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

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

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

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

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

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

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

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

次復(fù)數(shù)乘法運(yùn)算和Nlog2N次復(fù)數(shù)加法。4.2.1直接計(jì)算DFT的特點(diǎn)及減少

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

運(yùn)算量的根本途徑4、FFT的計(jì)算工作量利用的特性4.2.1直接計(jì)算DFT的特點(diǎn)及減少

運(yùn)算量的根本途徑利用以上特性，得到改進(jìn)DFT直接算法的方法。4.2.1直接計(jì)算DFT的特點(diǎn)及減少

運(yùn)算量的根本途徑

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

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

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

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

線性調(diào)頻Z變換(Chrip-z法)4.2.1直接計(jì)算DFT的特點(diǎn)及減少

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

應(yīng)補(bǔ)7個零值。0123456789101112131415nx(n)1、時域抽取算法原理4.2.2時域抽取法基2FFT根本原理Decimation-in-Time(DIT)設(shè)序列點(diǎn)數(shù)N=2M，M為整數(shù)。假設(shè)不滿足，那么補(bǔ)零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、蝶形運(yùn)算1111-1(1)N/2點(diǎn)的DFT運(yùn)算量:復(fù)乘次數(shù): 復(fù)加次數(shù):(2)兩個N/2點(diǎn)的DFT運(yùn)算量:復(fù)乘次數(shù): 復(fù)加次數(shù):(3)N/2個蝶形運(yùn)算的運(yùn)算量:復(fù)乘次數(shù):

復(fù)加次數(shù):

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

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

x1(0)=x(0)

x1(1)=x(2)N/2點(diǎn)

x1(2)=x(4)DFT

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

x2(1)=x(3)N/2點(diǎn)

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ù),可以進(jìn)一步把每個N/2點(diǎn)的序列再按其奇偶局部分解為兩個N/4的子序列。(2)N/2(4點(diǎn))-->N/4(2點(diǎn))FFT(2)N/2(4點(diǎn))-->N/4(2點(diǎn))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最后剩下兩點(diǎn)DFT,它可分解成兩個一點(diǎn)DFT，但一點(diǎn)DFT就等于輸入信號本身，所以兩點(diǎn)DFT可用一個蝶形結(jié)表示。取x(0)、x(4)為例。(3)N/4〔2點(diǎn)〕-->2個1點(diǎn)DFT1點(diǎn)DFTx(0)1點(diǎn)DFTx(4)X3(0)X3(1)2個1點(diǎn)的DFT蝶形流圖-1(3)N/4〔2點(diǎn)〕-->2個1點(diǎn)DFT一個完整的按時間抽取的8點(diǎn)FFT1〕輸入倒位序，輸出自然序自然順序n二進(jìn)制n2n1n0

倒位序二進(jìn)制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的運(yùn)算規(guī)律及流程圖畫法對N=2L點(diǎn)FFT，共需L級蝶形運(yùn)算，每級有N/2個蝶形運(yùn)算組成，蝶形運(yùn)算兩節(jié)點(diǎn)的距離:2L-1〔L表示級數(shù)〕每個蝶形運(yùn)算有一次復(fù)乘和再次復(fù)加。如1111-14.2.3DIT-FFT的運(yùn)算規(guī)律及流程圖畫法2〕蝶形運(yùn)算4.2.3DIT-FFT的運(yùn)算規(guī)律及流程圖畫法3〕WNr確實(shí)定(僅給出方法)-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的運(yùn)算規(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的運(yùn)算規(guī)律及流程圖畫法由上述分析可知,N=8需三級蝶形運(yùn)算N=23=8,由此可知,N=2L共需L級蝶形運(yùn)算,而且每級都由N/2個蝶形運(yùn)算組成,每個蝶形運(yùn)算有一次復(fù)乘,兩次復(fù)加。因此,N點(diǎn)的FFT的運(yùn)算量為復(fù)乘:mF=〔N/2〕L=〔N/2〕log2N復(fù)加:aF=NL=Nlog2N

由于計(jì)算機(jī)的乘法運(yùn)算比加法運(yùn)算所需的時間多得多，故以乘法作為比較基準(zhǔn)。4.2.4DIT-FFT算法與直接計(jì)算DFT

運(yùn)算量的比較4.2.4DIT-FFT算法與直接計(jì)算DFT

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

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

4.2.5按頻率抽取的FFT

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

-14.2.5按頻率抽取的FFT4.N=8時,按k的奇偶分解過程先蝶形運(yùn)算，后DFT:再將N/2點(diǎn)DFT按k的奇偶分解為兩個N/4點(diǎn)的DFT,如此進(jìn)行下去,直至分解為2點(diǎn)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〕根本蝶形不同時域：先復(fù)乘后加減頻域：先減后復(fù)乘2〕運(yùn)算量相同3〕都可原位運(yùn)算4〕兩者根本蝶形互為轉(zhuǎn)置5、時域抽取和頻域抽取FFT的區(qū)別

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

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);%計(jì)算x序列的DFTHk=fft(h,N);%計(jì)算h序列的DFTYk=Xk.*Hk;y=ifft(Yk)%求IDFTt1=cputime-t〔或toc〕stem(y);xlabel('n');4.3線性卷積的FFT算法4.3線性卷積的FFT