華東理工大學(xué)MATLAB課件

第四章快速(kuàisù)傅立葉變換

FastFourierTransform精品資料1第一節(jié)直接計(jì)算DFT的問題及改進(jìn)(gǎijìn)途徑1、問題(wèntí)的提出設(shè)有限長序列x(n)，非零值長度為N，若對(duì)x(n)進(jìn)行一次DFT運(yùn)算，共需多大的運(yùn)算工作量？計(jì)算成本?計(jì)算速度?精品資料22.DFT的運(yùn)算量回憶(huíyì)DFT和IDFT的變換式：1）x(n)為復(fù)數(shù)，也為復(fù)數(shù)。2）DFT與IDFT的計(jì)算量相當(dāng)。注意：精品資料3計(jì)算機(jī)運(yùn)算(yùnsuàn)時(shí)（編程實(shí)現(xiàn)）：N次復(fù)乘，N-1次復(fù)加

N個(gè)點(diǎn)以DFT為例：精品資料4復(fù)數(shù)乘法復(fù)數(shù)加法一個(gè)X(k)NN–1N個(gè)X(k)(N點(diǎn)DFT)N2N(N–1)實(shí)數(shù)乘法實(shí)數(shù)加法一次復(fù)乘42一次復(fù)加2一個(gè)X(k)4N2N+2(N–1)=2(2N–1)N個(gè)X(k)(N點(diǎn)DFT)4N22N(2N–1)運(yùn)算量(a+jb)(c+jd)=(ac-bd)+j(bc+ad)精品資料5例：計(jì)算(jìsuàn)一個(gè)N點(diǎn)DFT，共需N2次復(fù)乘。以做一次復(fù)乘1μs計(jì)，若N=4096，所需時(shí)間為例：石油勘探，有24個(gè)通道的記錄，每通道波形記錄長度(chángdù)為5秒，若每秒抽樣500點(diǎn)/秒，1）每道總抽樣點(diǎn)數(shù)：500*5=2500點(diǎn)2）24道總抽樣點(diǎn)數(shù)：24*2500=6萬點(diǎn)3）DFT復(fù)乘運(yùn)算時(shí)間：N2=(60000)2=36*108次精品資料6由于計(jì)算量大，且要求相當(dāng)大的內(nèi)存，難以實(shí)現(xiàn)實(shí)時(shí)處理，限制了DFT的應(yīng)用。長期以來，人們(rénmen)一直在尋求一種能提高DFT運(yùn)算速度的方法。FFT便是(biànshì)Cooley&Tukey在1965年提出的的快速算法，它可以使運(yùn)算速度提高幾百倍，從而使數(shù)字信號(hào)處理學(xué)科成為一個(gè)新興的應(yīng)用學(xué)科。精品資料7第二節(jié)改善DFT運(yùn)算效率的基本(jīběn)途徑1、利用DFT運(yùn)算的系數(shù)的固有對(duì)稱性和周期性，改善DFT的運(yùn)算效率。1）對(duì)稱性2）周期性3）可約性精品資料8精品資料92、將長序列(xùliè)DFT利用對(duì)稱性和周期性分解為短序列(xùliè)DFT的思路因?yàn)镈FT的運(yùn)算量與N2成正比的，如果一個(gè)大點(diǎn)數(shù)N的DFT能分解(fēnjiě)為若干小點(diǎn)數(shù)DFT的組合，則顯然可以達(dá)到減少運(yùn)算工作量的效果。精品資料10N點(diǎn)DFTN/2點(diǎn)DFTN/2點(diǎn)DFTN/4點(diǎn)DFTN/4點(diǎn)DFTN/4點(diǎn)DFTN/4點(diǎn)DFT…….復(fù)乘：精品資料11FFT算法的基本思想：利用DFT系數(shù)(xìshù)的特性，合并DFT運(yùn)算中的某些項(xiàng)把長序列DFT→短序列DFT，從而減少運(yùn)算量。FFT算法(suànfǎ)分類:時(shí)間抽選法DIT:Decimation-In-Time頻率抽選法DIF:Decimation-In-Frequency精品資料12第三節(jié)按時(shí)間抽選(chōuxuǎn)的基2-FFT算法1、算法(suànfǎ)原理

設(shè)輸入序列長度為N=2M(M為正整數(shù)，將該序列按時(shí)間順序的奇偶分解為越來越短的子序列，稱為基2按時(shí)間抽取的FFT算法。也稱為Coolkey-Tukey算法。

其中基2表示：N=2M，M為整數(shù).若不滿足這個(gè)條件，可以人為地加上若干零值（加零補(bǔ)長）使其達(dá)到N=2M。精品資料13先將x(n)按n的奇偶分為兩組，作變量置換:當(dāng)n=偶數(shù)(ǒushù)時(shí)，令n=2r;當(dāng)n=奇數(shù)時(shí)，令n=2r+1;分組，變量(biànliàng)置換2、算法步驟得到：精品資料14帶入DFT中精品資料15所以(suǒyǐ)由于(yóuyú)?精品資料16X1(k)、X2(k)只有N/2個(gè)點(diǎn)，以N/2為周期；而X(k)卻有N個(gè)點(diǎn)，以N為周期。要用X1(k)、X2(k)表達(dá)全部(quánbù)的X(k)值，還必須利用WN系數(shù)的周期特性。精品資料17后半部分前半部分又考慮到的對(duì)稱性：有：精品資料18后半部分前半部分蝶形運(yùn)算(yùnsuàn)流圖符號(hào)說明(shuōmíng)：(1)左邊兩路為輸入(2)右邊兩路為輸出(3)中間以一個(gè)小圓表示加、減運(yùn)算（右上路為相加輸出、右下路為相減輸出)1個(gè)蝶形運(yùn)算需要1次復(fù)乘，2次復(fù)加精品資料19復(fù)數(shù)乘法復(fù)數(shù)加法一個(gè)N點(diǎn)DFTN2N(N–1)一個(gè)N/2點(diǎn)DFT(N/2)2N/2(N/2–1)兩個(gè)N/2點(diǎn)DFTN2/2N(N/2–1)一個(gè)蝶形12N/2個(gè)蝶形N/2N總計(jì)N2/2+N/2≈N2/2N(N/2-1)+N=N2/2運(yùn)算量減少(jiǎnshǎo)了近一半分解(fēnjiě)后的運(yùn)算量：精品資料20先將N=8點(diǎn)的DFT分解成2個(gè)4點(diǎn)DFT：可知(kězhī)：時(shí)域上：x(0),x(2),x(4),x(6)為偶子序列x(1),x(3),x(5),x(7)為奇子序列頻域上：X(0)~X(3)，由X(k)給出X(4)~X(7)，由X(k+N/2)給出例子(lìzi)：求N=23=8點(diǎn)FFT變換

按N=8→N/2=4，做4點(diǎn)的DFT：精品資料21N=8點(diǎn)的直接(zhíjiē)DFT的計(jì)算量為：復(fù)乘：N2次=64次復(fù)加：N(N-1)次=8×7=56次此外(cǐwài)，還有4個(gè)蝶形結(jié)，每個(gè)蝶形結(jié)需要1次復(fù)乘，2次復(fù)加。一共是：復(fù)乘4次，復(fù)加8次。得到X1(k)和X2(k)需要：復(fù)乘：(N/2)2+(N/2)2次=32次復(fù)加：N/2(N/2-1)+N/2(N/2-1)=12+12=24次用分解的方法得到X

(k)需要：復(fù)乘：32+4=36次復(fù)加：24+8=32次精品資料22N點(diǎn)DFT的一次時(shí)域抽取(chōuqǔ)分解圖(N=8)4點(diǎn)DFT4點(diǎn)DFTx(0)x(2)x(4)x(6)x(1)x(3)x(5)x(7)X1(0)X1(1)X1(2)X1(3)X2(0)X2(1)X2(2)X2(3)X(0)X(1)X(2)X(3)X(4)X(5)X(6)X(7)精品資料23因?yàn)?點(diǎn)DFT還是比較麻煩，所以再繼續(xù)(jìxù)分解。若將N/2(4點(diǎn))子序列(xùliè)按奇/偶分解成兩個(gè)N/4點(diǎn)(2點(diǎn))子序列(xùliè)。即對(duì)將x1(r)和x2(r)分解成奇、偶兩個(gè)N/4點(diǎn)(2點(diǎn))點(diǎn)的子序列(xùliè)。精品資料24那么(nàme)，X1(k)又可表示為精品資料25X2(k)也可以進(jìn)行相同(xiānɡtónɡ)的分解：注意：通常我們會(huì)把寫成。精品資料26N點(diǎn)DFT的第二次時(shí)域抽取(chōuqǔ)分解圖(N=8)2點(diǎn)DFT2點(diǎn)DFT2點(diǎn)DFT2點(diǎn)DFTx(0)x(4)x(2)x(6)x(1)x(5)x(3)x(7)X3(0)X3(1)X4(0)X4(1)X5(0)X5(1)X6(0)X6(1)X1(0)X1(1)X1(2)X1(3)X2(0)X2(1)X2(2)X2(3)X(0)X(1)X(2)X(3)X(4)X(5)X(6)X(7)4點(diǎn)DFT4點(diǎn)DFTx(0)x(2)x(4)x(6)x(1)x(3)x(5)x(7)X1(0)X1(1)X1(2)X1(3)X2(0)X2(1)X2(2)X2(3)X(0)X(1)X(2)X(3)X(4)X(5)X(6)X(7)精品資料2788X3(0)X3(1)x(0)=x3(0)x(4)=x3(1)精品資料28N點(diǎn)DIT―FFT運(yùn)算(yùnsuàn)流圖(N=8)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)精品資料293、DIT―FFT算法與直接計(jì)算(jìsuàn)DFT運(yùn)算量的比較1)、N=2M的DFT運(yùn)算(yùnsuàn)可分成M級(jí)，每一級(jí)有N/2個(gè)蝶形，每個(gè)蝶形有一次復(fù)乘兩次復(fù)加。2)、所以M級(jí)共有次復(fù)乘和次復(fù)加。3)、若直接計(jì)算DFT，需N2次復(fù)乘和N(N-1)次復(fù)加。顯然，當(dāng)N較大時(shí)，有：例如，N=210=1024時(shí)精品資料30FFT算法與直接計(jì)算(jìsuàn)DFT所需乘法次數(shù)的比較曲線精品資料314、DIT―FFT的運(yùn)算(yùnsuàn)規(guī)律及編程思想FFT的每級(jí)（列）計(jì)算都是由N個(gè)復(fù)數(shù)數(shù)據(jù)（輸入(shūrù)）兩兩構(gòu)成一個(gè)蝶型（共N/2個(gè)蝶形）運(yùn)算而得到另外N個(gè)復(fù)數(shù)數(shù)據(jù)（輸出）。當(dāng)數(shù)據(jù)輸入到存儲(chǔ)器以后，每一組運(yùn)算的結(jié)果，仍然存放在這同一組存儲(chǔ)器中直到最后輸出。例：將x(0)放在單元A(0)中，將x(4)放在單元A(1)中，W80

放在一個(gè)暫存器中。將x(0)+W80x(4)→送回A(0)單元將x(0)-W80x(4)→送回A(1)單元X3(0)X3(1)x(0)x(4)1）

原位運(yùn)算(亦稱同址計(jì)算)精品資料32x(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)回顧(huígù)：N點(diǎn)DIT―FFT運(yùn)算流圖(N=8)精品資料33如上所述，N點(diǎn)DIT―FFT運(yùn)算流圖中，每級(jí)都有N/2個(gè)蝶形。每個(gè)蝶形都要乘以因子WNP，稱其為旋轉(zhuǎn)(xuánzhuǎn)因子，p稱為旋轉(zhuǎn)(xuánzhuǎn)因子的指數(shù)。2)旋轉(zhuǎn)(xuánzhuǎn)因子的變化規(guī)律觀察FFT運(yùn)算流圖發(fā)現(xiàn)，第L級(jí)共有2L-1個(gè)不同的旋轉(zhuǎn)因子。N=23=8時(shí)的各級(jí)旋轉(zhuǎn)因子表示如下：L=1時(shí)，WNp=WN/4J，N/4=21=2L，J=0L=2時(shí)，WNp=WN/2J，N/2=22=2L，J=0，1L=3時(shí)，WNp=WNJ，N=23=2L，J=0，1，2，3精品資料34對(duì)N=2M的一般(yībān)情況，第L級(jí)的旋轉(zhuǎn)因子為：精品資料35設(shè)序列x(n)經(jīng)時(shí)域抽選(倒序)后，存入數(shù)組X中。如果蝶形運(yùn)算的兩個(gè)輸入數(shù)據(jù)相距B個(gè)點(diǎn)(B=2L-1)，應(yīng)用原位計(jì)算(jìsuàn)，則蝶形運(yùn)算可表示成如下形式：下標(biāo)L表示(biǎoshì)第L級(jí)運(yùn)算，XL(J)則表示(biǎoshì)第L級(jí)運(yùn)算后數(shù)組元素X(J)的值。精品資料363)編程思想(sīxiǎng)及流程圖開始送入x(n)和N=2M調(diào)整輸入x(n)的順序for(L=1;L<=M;L++)B=2L-1for(J=0;J<=B-1;J++)p=J·2M-Lfor(k=J;k<=N-1;k=k+2L)輸出結(jié)果結(jié)束精品資料374）碼位倒序(dǎoxù)由N=8蝶形圖看出：原位計(jì)算時(shí)，F(xiàn)FT輸出的X(k)的次序正好是順序(shùnxù)排列的，即X(0)…X(7),但輸入x(n)都不能按自然順序(shùnxù)存入到存儲(chǔ)單元中，而是按x(0),x(4),x(2),x(6),x(1),x(5),x(3),x(7)的順序(shùnxù)存入存儲(chǔ)單元，即為亂序輸入，順序(shùnxù)輸出。這種順序(shùnxù)看起來相當(dāng)雜亂，然而它是有規(guī)律的。即碼位倒讀規(guī)則。精品資料38自然順序n二進(jìn)制碼表示碼位倒讀碼位倒置順序n’以N=8為例：0123456700000101001110010111011100010001011000110101111104261537看出：碼位倒讀后的順序剛好(gānghǎo)是數(shù)據(jù)送入計(jì)算機(jī)內(nèi)的順序。精品資料39倒序(dǎoxù)規(guī)律精品資料40對(duì)于(duìyú)數(shù)N，在其二進(jìn)制最高位加1，等于加N/2。若已知某個(gè)反序號(hào)為J，為求下一個(gè)反序號(hào)，可先判J的最高位：1)若為0，則把該位變成1（即加N/2）就得到下一個(gè)反序號(hào)，2)若為1，則需判斷次高位：①若次高位為0，則把最高位變0（相當(dāng)(xiāngdāng)減去N/2）后，再把次高位變1（即加N/4）。②若次高位為1，則需判斷次次高位……分析：精品資料41倒序排列算法的流程圖正序序列已在數(shù)組A[]中，輸入NLH=N/2,j=LH,N1=N-2j=j-kk=k/2k=LHj<kj=j+kT=A(j)A(i)=A(j)A(j)=Tfor(i=1;i<=N1;i++)i≥jNYYN精品資料42第四節(jié)按頻率(pínlǜ)抽選的基2-FFT算法在基2快速算法中，頻域抽取法FFT也是一種(yīzhǒnɡ)常用的快速算法，簡稱DIF―FFT。設(shè)序列x(n)長度為N=2M，首先將x(n)前后對(duì)半分開，得到兩個(gè)子序列，其DFT可表示為如下形式精品資料43精品資料44精品資料45精品資料46DIF―FFT一次分解(fēnjiě)運(yùn)算流圖(N=8)4點(diǎn)DFT4點(diǎn)DFTx(0)x(1)x(2)x(3)x(4)x(5)x(6)x(7)X(0)X(2)X(4)X(6)X(1)X(3)X(5)X(7)x1(0)x1(1)x1(2)x1(3)x2(0)x2(1)x2(2)x2(3)精品資料47DIF―FFT二次分解(fēnjiě)運(yùn)算流圖(N=8)精品資料48DIF―FFT運(yùn)算(yùnsuàn)流圖(N=8)精品資料49時(shí)間抽取算法(suànfǎ)與頻率抽取算法(suànfǎ)的比較1)頻率(pínlǜ)抽選法和時(shí)間抽選法總的計(jì)算量是相同的復(fù)乘：復(fù)加：2)頻率抽取法和時(shí)間抽取法一樣，都適用于原位運(yùn)算，即蝶形的輸入和輸出占用同一個(gè)存儲(chǔ)單元。3)均存在碼位倒序問題。4)頻率抽選法和時(shí)間抽選法一樣，基本運(yùn)算也是蝶形運(yùn)算。但兩者的蝶形形式略有不同。精品資料50第五節(jié)IDFT的快速(kuàisù)算法－IFFT上述(shàngshù)FFT算法流圖也可以用于離散傅里葉逆變換(InverseDiscreteFourierTransform，簡稱IDFT)。比較DFT和IDFT的運(yùn)算公式：1)旋轉(zhuǎn)因子：2)系數(shù)：精品資料51DIT―IFFT運(yùn)算(yùnsuàn)流圖精品資料52DIT―IFFT運(yùn)算(yùnsuàn)流圖(防止溢出)精品資料53如果(rúguǒ)希望直接調(diào)用FFT子程序計(jì)算IFFT，則可用下面的方法：對(duì)上式兩邊(liǎngbiān)同時(shí)取共軛，得：精品資料54例1、如果通用計(jì)算機(jī)的速度(sùdù)為平均每次復(fù)乘需要5s，每次復(fù)加需要0.5s，用它來計(jì)算512點(diǎn)的DFT[x(n)]，問：1）直接(zhíjiē)計(jì)算需要多少時(shí)間？2）用FFT需要多少時(shí)間？解：1）用DFT進(jìn)行運(yùn)算：復(fù)乘：T1=N2×5×10-6=1.31072秒復(fù)加：T2=N(N-1)×0.5×10-6=0.130816秒總共：T=T1+T2=1.441536秒精品資料552）用FFT進(jìn)行(jìnxíng)運(yùn)算：復(fù)乘：T1’=(N/2)log2N×5×10-6=0.01152秒復(fù)加：T2’=Nlog2N×0.5×10-6=0.002304秒總共(zǒnggòng)：T’=T1’+T2’=0.013824秒精品資料56例2、長度為240點(diǎn)的序列x(n)與長度為N點(diǎn)的h(n)卷積。當(dāng)N=10和240時(shí)，直接進(jìn)行(jìnxíng)卷積x(n)*h(n)和用IFFT[X(K)·H(K)]的方法相比，那種方法求解y(n)的效率更高？x(n)h(n)y(n)=x(n)*h(n)L≥N1+N2-1X(k)補(bǔ)L-N1個(gè)零x(n)L點(diǎn)DFT補(bǔ)L-N2個(gè)零h(n)L點(diǎn)DFTL點(diǎn)IDFTy(n)=x(n)*h(n)H(k)Y(k)精品資料57直接(zhíjiē)進(jìn)行卷積（N=10）：乘法(chéngfǎ)次數(shù)：240×10=2400次用FFT的方法（N=10）：添零到256點(diǎn)，L=256乘法次數(shù)：3×(L/2)log2L＋L=3×128×8＋256=3328次精品資料58直接(zhíjiē)進(jìn)行卷積（N=240）：乘法(chéngfǎ)次數(shù)：240×240=57600次用FFT的方法（N=240）：添零到512點(diǎn)，L=512乘法次數(shù)：3×(L/2)log2L＋L=3×256×9＋512=7424次精品資料59第六節(jié)進(jìn)一步而減少(jiǎnshǎo)運(yùn)算量的措施1、多類蝶形單元(dānyuán)運(yùn)算DIF―FFT運(yùn)算流圖(N=8)精品資料60一個(gè)旋轉(zhuǎn)因子(yīnzǐ)對(duì)應(yīng)一個(gè)蝶形單元。(1)若程序中包含(bāohán)了所有的旋轉(zhuǎn)因子，則稱該對(duì)應(yīng)一個(gè)蝶形單元。

(2)若去掉有關(guān)于WNr=±1旋轉(zhuǎn)因子的乘法運(yùn)算，則稱含二類蝶形單元。

(3)若再去掉有關(guān)于WNr=±j旋轉(zhuǎn)因子的乘法運(yùn)算，則稱含三類蝶形單元。

(4)若再特殊處理關(guān)于旋轉(zhuǎn)因子的乘法運(yùn)算，則稱含四類蝶形單元。精品資料61表4.1基2-FFT在各種(ɡèzhǒnɡ)蝶形單元下所需的實(shí)數(shù)乘法的次數(shù)蝶形單元種類一類蝶形單元實(shí)數(shù)乘法次數(shù)二類蝶形單元實(shí)數(shù)乘法次數(shù)三類蝶形單元實(shí)數(shù)乘法次數(shù)四類蝶形單元實(shí)數(shù)乘法次數(shù)N=24000N=8482084N=32320196136108N=128179212841032908N=5129216717261525644N=204845056368683277630372精品資料622、旋轉(zhuǎn)(xuánzhuǎn)因子的生成將旋轉(zhuǎn)(xuánzhuǎn)因子中正弦和余弦函數(shù)值余弦存放在數(shù)組中，在程序執(zhí)行時(shí)查表得到，比直接讓程序計(jì)算正弦余弦值的效率更高。3、實(shí)序列的FFT算法

(1)利用第三章的方法，用一次N點(diǎn)FFT計(jì)算兩個(gè)N點(diǎn)的實(shí)序列的FFT，一個(gè)實(shí)輸入作為實(shí)部，另一個(gè)實(shí)輸入作為虛部，計(jì)算完成后將結(jié)果進(jìn)行簡單的分解即可得到它們各自的FFT結(jié)果。精品資料63(2)另一種方法是用N/2點(diǎn)的FFT計(jì)