05-第五章快速傅里葉變換(蝶形運(yùn)算)

1、第五章快速傅里葉變換本章目錄直接計(jì)算DFT的問題及改進(jìn)的途徑 按時(shí)間抽取的基2-FFT算法 按頻率挾取的基2-FFT算法 快速傅里葉逆變換(IFFT)算法 Matlab 實(shí)現(xiàn)5J引言 DFT在實(shí)際應(yīng)用中很重要:可以計(jì)算信號的頻譜、功率譜和線性卷積等。直接按DFT變換進(jìn)行計(jì)算，當(dāng)序列長度汕艮大時(shí)，計(jì)算量非常大，所需時(shí)間會(huì)很長。 FFT并不是一種與DFT不同的變換，而是DFT的一種快速計(jì)算的算法。竺亙斟算DFT的問題及改進(jìn)的途徑DFT的運(yùn)算量5#設(shè)復(fù)序列MG)長度為N點(diǎn)，其DFT為N1X(Q = 2>(曲n=0k=0,，AM(1)計(jì)算一個(gè)X(k)值的運(yùn)算量復(fù)數(shù)乘法次數(shù)：N復(fù)數(shù)加法次數(shù)：NT5

2、.2J DFT的運(yùn)算量(2)計(jì)算全部N個(gè)X(k)值的運(yùn)算量復(fù)數(shù)乘法次數(shù)：/V2復(fù)數(shù)加法次數(shù)：N(N1)(3)對應(yīng)的實(shí)數(shù)運(yùn)算量NlN7X 伙)=工兀()比篇=Re x(n) + jlmx(n)RcWk + jImWkn=0n=0N-l= Rex(n)ReWk Imx(n)ImWkw=0+yRe x(n) Im Wk + Im x(n) Re7一次復(fù)數(shù)乘法：4次實(shí)數(shù)乘法+ 2次實(shí)數(shù)加法一個(gè)X(k):4N次實(shí)數(shù)乘法+2A/+2(AM)= 2(2 AM)次實(shí)數(shù)加法所以整個(gè)N點(diǎn)DFT運(yùn)算共需要：實(shí)數(shù)乘法次數(shù)：4/V2實(shí)數(shù)加法次數(shù)：NX2(2AM)=2N(2AM)9*2!算量的結(jié)論N點(diǎn)DFT的復(fù)數(shù)乘法次數(shù)

3、舉例NN22441686416256321028NN26440491281638425665 536512262 14410241 048 576結(jié)論：當(dāng)M艮大時(shí)，其運(yùn)算量很大，對實(shí)時(shí)性很強(qiáng)的信號 處理來說，要求計(jì)算速度快，因此需要改進(jìn)DFT的計(jì)算 方法，以大大減少運(yùn)算次數(shù)。115.2.2減少運(yùn)算工作量的途徑主要原理是利用系數(shù)W護(hù)的以下特性對DFT進(jìn)行分解:(1)對稱性(wky = wnk =(2)周期性(3)可約性VVNwW+N)一 VVN另外,w;k 呼=w-nktmN / tn</2=-l 呼+ N/2)=肌5.3按時(shí)間抽取的基ZFFT算法算法原理 按時(shí)間抽取基-2FFT算法與直接

4、計(jì)算D FT運(yùn)算量的比較按時(shí)間抽取的FFT算法的特點(diǎn)按時(shí)間抽取FFT算法的其它形式流程圖5.3J算法原理設(shè)N=2J將x（門）按n的奇偶分為兩組:x(2r)=召(廠) x(2r + 1) = x2(r)N-iX伙)=DFTxM=工兀昭*7?=0NN-=工兀（）*+工兀）”解/?=0n=0幾為偶數(shù)料為奇數(shù)15n=0n為奇數(shù)2=2何呼+ 2（）呼A?=0n為偶數(shù)2=S 雙2門叭"+ 工 X2r + l）W；2r+1U r=0r=07_1=2心）畔+觀工勺（廠）",=X 伙）+呎X2伙）4佇1r=02r=022式中，E（Q和蜀仇）分別是re）和工2何的N/2的DFT。另外，式中氐的

5、取值范圍是：0, 1,，N/2-1 o此，X(Q = X") +說/伙)只能計(jì)算出x(k)的前一半值。后一半X(k)1¥N/2 , N/2 +1, N ?利用 可得到J 2亠JyV/2-1N/21X(- + k)= £ 兀1(門叭；尸)=£ x“)W篇2 = X”)2r=0r=0同理可得X2(y+ £) = X2(k)考慮到W嚴(yán)k) = W閉2 ,Wk =_Wk及前半部分X X(k) = X(k) + WX2(k)k=0, 1,N/2-1因此可得后半部分X(QX(k + -) = Xl(k+) + W+N/2X2 (k + y)= Xl(k)-

6、WX2(k)k=Qf 1,N/2-119以8點(diǎn)為例第_次按奇偶分解蝶形運(yùn)算X伙)=/伙)+叱X?伙)蝶形運(yùn)算式X(k) = x(切-吩2 伙)蝶形運(yùn)算信 號流圖符號#以8點(diǎn)為例第_次按奇偶分解#以8點(diǎn)為例第_次按奇偶分解因此，只要求出2個(gè)N/2點(diǎn)的DFT,即E(k)和X2(k),再 經(jīng)過蝶形運(yùn)算就可求出全部X(k)的值，運(yùn)算量大大減少。#以8點(diǎn)為例第_次按奇偶分解21以8點(diǎn)為例第_次按奇偶分解兀(0) =x(0) >xi(l)=x(2)>Xj(2) =x(4)>X|(3) =x(6) <x2(0) =x(1) x2(1)=x(3)<X2(2) =x(5).以3)

7、=x(7) 葺點(diǎn)罟點(diǎn)b(0)M)川2)DFT加3)總3)2(0)乂X(5)DFT06)X(7)X|(l)基X(0)i/2)以N=8為例,分解為2個(gè)4點(diǎn) 的DFT,然后 做8/2=4次蝶形運(yùn)算即可求出 所有8點(diǎn)X(k)的 值。#以8點(diǎn)為例第_次按奇偶分解蝶形運(yùn)算量比較觸DFT的運(yùn)算量復(fù)數(shù)乘法次數(shù)：N2復(fù)數(shù)加法次數(shù)：N(N1)分解一次后所需的運(yùn)算量=2個(gè)N/2的DFT+23以8點(diǎn)為例第_次按奇偶分解N/2蝶形:復(fù)數(shù)乘法次數(shù)：2*(M2尸+N/2二N2/2+M2復(fù)數(shù)加法次數(shù)：1)+2*"/2=2/2因此通過一次分解后，運(yùn)算工作量減少了差不多一半。#進(jìn)_步按奇偶分解由于N=2S因而N/2仍是

8、偶數(shù)，可以進(jìn)一步把每個(gè)N/2點(diǎn) 子序列再按其奇偶部分分解為兩個(gè)N/4點(diǎn)的子序列。2 = 0,1,，理14以N/2點(diǎn)序列可(廠)為例再)=召x1(2Z + 1) = x4(Z)則有N/2-N/4-1N/4-1X(k)=血=此扇 + + 1)肌$譏r=0/=0/=0N/4-1N/41=陀4+W怎陀41=01=0= X3(k) + W/2X4(k) B0丄,J-125以8點(diǎn)為例第二次按奇偶分解(N A,N 5X - + k =X3(k)-W爲(wèi)X4伙)f ,-1I 4 丿由此可見，一個(gè)N/2點(diǎn)DFT可分解成兩個(gè)N/4點(diǎn)DFT。同理，也可對兀2(n)進(jìn)行同樣的分解，求出X2(k)027以8點(diǎn)為例第二次按

9、奇偶分解#以8點(diǎn)為例第二次按奇偶分解應(yīng)(0)=兀1(1) =x(2) x4(l) =X|(3) =x(6) %3(0)=%|(0)=兀(0) 如1)x3(l)=xi(2)=x(4)e無3（°）益(0)W)孚點(diǎn)DFTML2WnX1Xi-1冷3)掄(0)E占4小、DFT疋X(0)x5(0) =x2(0) =x(l)x5(1)=x2(2) =x(5)/占4小、兀6(0) =X2(1)=X(3) 兀6(1)=七(3) =x(7) DFT屁花(0)總1)無(2)X(3)X(4)X(5)X(6)-1WN 赦1)WN XA2(3)x（0）X29以8點(diǎn)為例第二次按奇偶分解算法原理#以8點(diǎn)為例第二次按

10、奇偶分解對此例N=8,最后剩下的是4個(gè)N/4=2點(diǎn)的DFT, 2點(diǎn) DFT也可以由蝶形運(yùn)算來完成。以禺仇)為例。N/41N/411心伙)二工兀3W/4二工律/4上=°，1/=01=0即X3 (0) = x3 (0) + 必® (!)=班°)+ 必雙4) = x(0) + W：x(4)X3=%3(°)+昭兀3=x(0) +岡兀(4) =%(0)-叱汝這說明，N=2M的DFT可全部由蝶形運(yùn)算來完成。31N=8按時(shí)間抽取法FFT信號流圖5.3.2按時(shí)間抽販基ZFFT算法與H接計(jì)算DFT運(yùn)算豪的比較2次復(fù)加，因此每級運(yùn)算都需M2次復(fù)乘和由按時(shí)間抽取法FFT的信號

11、流圖可知，當(dāng)N=2厶時(shí)，共有L級 蝶形運(yùn)算；每級都由 M2個(gè)蝶形運(yùn)算組成,而每個(gè)蝶形有 1次復(fù)乘、N次復(fù)加。33這樣丄級運(yùn)算總共需要：復(fù)數(shù)乘法：-L = -log22 2復(fù)數(shù)加法：N直接DFT算法運(yùn)算量復(fù)數(shù)乘法：N1復(fù)數(shù)加法：N(N1)直接計(jì)算DFT與FFT算法的計(jì)算量之比為M“ N?2N10g2”=g2N+ FFT算法與直接DFTX法運(yùn)算量的比較-NN2.N、,計(jì)算量 之比MNN2N、. 三呃N計(jì)算量 之比M2414.012816 38444836.641644.025665 5361 02464.0864125.4512262 1442 304113.816256328.010241 04

12、8 5765 120204.83210288012.820484 194 30411 264372.464404919221.4355.3.3按時(shí)間抽取的FFT算法的特點(diǎn)序列的逆序排列同址運(yùn)算（原位運(yùn)算）蝶形運(yùn)算兩節(jié)點(diǎn)間的距離 C的確定#I序列的逆序排列序列的逆序排列由于歡）被反復(fù)地按奇、偶分組，所以流圖輸入端的排列不再是順序的，但仍有規(guī)律可循:37#因?yàn)镹=2M ,對于任意n(0<n<A-l),可以用M個(gè)#二進(jìn)制碼表示為:n(DEC) = (“M-inM-2 A 2“l(fā)“0)(B/N) fo1反復(fù)按奇、偶分解時(shí)，即按二進(jìn)制碼的分解。#倒位序的樹狀圖(N-8)#01x(切屮

13、6;)1"2-一x(OOO)-一 x(100)-一一 x(010)-x(110) 0-一一 x(001)班101)-一兀(011) x(lll)#碼位的倒位序(Nm8)自然順序n二鯉數(shù)倒位序二進(jìn)制數(shù) 倒位序順序數(shù)分000000001100010010011110100001101101110Oil11111139倒位序的變址處理(N8)41#存儲(chǔ)單元自然順序單元力兀(0)變址Tx(0)倒位序/(8)x(7)Tx(7)A(2) A(3) A(4) A(5) A(6) A(7) x(l) x(2) x(3) x(4) x(5) x(6)x(4)兀(2) x(6) x(l) x(5) x(

14、3)同址運(yùn)算（原位運(yùn)算）同址運(yùn)算（原位運(yùn)算）某一列任何兩個(gè)節(jié)點(diǎn)氐和/的節(jié)點(diǎn)變量進(jìn)行蝶形運(yùn)算 后，得到結(jié)果為下一列肛/兩節(jié)點(diǎn)的節(jié)點(diǎn)變量，而和其他 節(jié)點(diǎn)變量無關(guān)。這種原位運(yùn)算結(jié)構(gòu)可以節(jié)省存儲(chǔ)單元， 降低設(shè)備成本。43#運(yùn)算后X伙）一4伙）#觀察原位運(yùn)算規(guī)律雙0)X4)-1雙2)兀XI)雙5)x(3)兀wl1-1總0)M)X(3)X(4)十 X(6)X(7)45#蝶形運(yùn)算兩節(jié)點(diǎn)間的距離蝶形運(yùn)算兩節(jié)點(diǎn)間的距離以N=8為例:第一級蝶形,距離為:第二級蝶形,距離為:第三級蝶形,距離為:規(guī)律：對于共L級的蝶形而言，其m級蝶形運(yùn)算的節(jié)點(diǎn)間的距離為2心忙的確定 C的確定以N=8為例：加=10寸，Wr=Wj/4=

15、WJm =<,j = 0加=20寸，W=W/2=Wm =Wj = O,l加=3時(shí)，W=W =Wm =Wj = O丄2,3N = 2mL 級：肌二昭，丿=0丄2,A ,2i-l&2L =2M x 2lm =Nx 2lm5.4按頻率抽取的基2-FFT算法算法原理先把輸入按門的順序分成前后兩半再把輸出X(k)按幺的奇偶分組設(shè)序列長度為N=2J L為整數(shù)47#f前半子序列只門)NI后半子序列血+亍N 1ox -10<n< 鄉(xiāng) T2495.4J算法原理由DFT定義得7V-1X伙） = 0（曲n=0N/2-1N-=£ x（）W/ + 工 x（）Wf/?=0n=N 12皆

16、 1, Ng N 5+3=£班砒叭+ 2兀（+虧）必2 n=0mn厶N(yùn)/21/?=0n=0- N"TN kE心）+兀（十3）%2 叱yk=0 9 1,，N由于N/2-x(k)= £/?=0x(n) + x(n + )WjkN 2 =嚴(yán)=_51#所以N 12-X閃=Sn=0x(n) + (-1)A x(n + )Kk#k=0，1,，N#然后按氐的奇偶可將X(Q分為兩部分則式k = 2rN .r=0, 19，一 1k = 2廠+ 12nkN/2-1Nx(k)= 2 x(n) + (-l)kx(n + -) W；53/?=0可轉(zhuǎn)化為#N oX(2"二工 x(n

17、)x(n + )吒N/2-171=0nrN/21X"=0Nx(n) + xn + ) W.-nrN/2#n(2r+)NN/2-1NX(2r + 1)=工x(n)-x(n + )-h=o2N/2-1NX 兀-兀+ )%' 昭;2/?=o2#N(n) = x(n) + x(n + )N x2(n) = x(n) -xn + )代入N/2-1X(2r)=工打=0NgNX(2r + 1)=工Wn)-x(H + -W；.C；2x(h) + x(n + )必;2n=0可得心X(2r) = 2>(叫277=0 N/2-1X(2r + 1)=工也睥77=0r=0, 1為2個(gè)N/2點(diǎn)的D

18、FT,合起來正好是N點(diǎn)X(Q的值。55蝶形運(yùn)算N£ (/7)= x(n) + x(n + )稱為蝶形運(yùn)算57#與時(shí)間抽選基2FFT算法中的蝶形運(yùn)算符號略有不同。右例按頻率抽取(N8)例 按頻率抽取，將N點(diǎn)DFT分解為兩個(gè)N/2點(diǎn)DFT的組合(N=8)59#N占-2小、DFTI/O)i/2)i/4)i/6)多點(diǎn)DFTi(l)f(3)i(5)i/7)#與時(shí)間抽取法的推導(dǎo)過程一樣，由于N=2S N/2仍然是一個(gè)偶數(shù), 與奇數(shù)組,因而可以將每個(gè)N/2點(diǎn)DFT的輸出再分解為偶數(shù)組這就將N/2點(diǎn)DFT進(jìn)一步分解為兩個(gè)N/4點(diǎn)DFT。抽取法與時(shí)間抽取法的異同頻率抽取法輸入是自然順序，輸出是倒位序

19、的；時(shí)間抽取法正好相反。頻率抽取法的基本蝶形與時(shí)間抽取法的基本蝶形有所不同。頻率抽取法運(yùn)算量與時(shí)間抽取法相同。頻率抽取法與時(shí)間抽取法的基本蝶形是互為轉(zhuǎn)置的。右55快速傅里葉逆變換（IFFT）算法IDFT公式1 N1兀（町=IDFT X（£）=帚工 X （k）WtjkN k=oDFT公式NlX伙）=DFT 兀（訓(xùn)二工心）呼71=0比較可以看出,IDFT多出 N 二F =>-=-M個(gè)1/2可分解到M級蝶形運(yùn)算中。例頻率抽取IFFT流圖(28)快速傅里葉逆變換另一種算法NIDFTX（k） = XW1 "T/=-工對伙）（"廣）N L-Ck=0廠NlN1nkNk=0

20、k=0IFFTX 伙）=X伙）_杏丿割一x *伙）FFTX*伙）求弊 FFFX*伙）除以何今 x（n）5.8 Matlab實(shí)現(xiàn)用FFT進(jìn)行譜分析的Matlab實(shí)現(xiàn)用CZT進(jìn)行譜分析的Matlab實(shí)現(xiàn)在Matlab中使用的線性調(diào)頻z變換函數(shù)為czt, 其調(diào)用格式為 »X= czt(x3 M, W, A)其中，x是待變換的時(shí)域信號x(n),其長度為N, M是變換的長度，W確定變換的步長，A確定變換的起點(diǎn)。M=N, A=1,貝!|CZT變成DFT。655.8J用FFT進(jìn)行譜分析的Matlab實(shí)現(xiàn)例51設(shè)模擬信號X0 = 2sin(4奴)+ 5cos(8奴)，以t= O.Oln («

21、;=0: -1)進(jìn)行取樣，試用fft函數(shù)對其做頻譜分析。N分別 為:(1)N=45； (2)N=50； (3)N=55； (2)N=60。程序清單如下1=#%計(jì)算N=45的F FT并繪出其幅頻曲線N=45;n=0:N-1 ;t=O.O1*n;q=n*2*pi/N; x=2*sin(4*pi*t)+5*cos(8*pi*t);y=fft(x,N); figure(1)subplot(2,2,1) plot(q5abs(y) title(TFT N=45J例51程序清單%計(jì)算N=50的FFT并繪出其幅頻曲線N=50;n=0:N-1 ;t=0.01*n; q=n*2*pi/N;x=2*sin(4*pi*t)+5*cos(8*pi*t);y=fft(x5N); figure subplot(2,2