合工大 數(shù)字信號處理 第四章

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

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

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

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

其中基2表示：N=2M，M為整數(shù).若不滿足這個條件，可以人為地加上若干零值（加零補長）使其達到N=2M。精品資料先將x(n)按n的奇偶分為兩組，作變量(biànliàng)置換:當n=偶數(shù)時，令n=2r;當n=奇數(shù)時，令n=2r+1;分組，變量(biànliàng)置換2、算法步驟得到：精品資料帶入DFT中精品資料所以(suǒyǐ)由于(yóuyú)?精品資料X1(k)、X2(k)只有(zhǐyǒu)N/2個點，以N/2為周期；而X(k)卻有N個點，以N為周期。要用X1(k)、X2(k)表達全部的X(k)值，還必須利用WN系數(shù)的周期特性。精品資料后半部分前半部分又考慮到的對稱性：有：精品資料后半部分前半部分蝶形運算(yùnsuàn)流圖符號說明：(1)左邊兩路為輸入(2)右邊兩路為輸出(3)中間以一個小圓表示加、減運算(yùnsuàn)（右上路為相加輸出、右下路為相減輸出)1個蝶形運算需要1次復乘，2次復加精品資料復數(shù)乘法復數(shù)加法一個N點DFTN2N(N–1)一個N/2點DFT(N/2)2N/2(N/2–1)兩個N/2點DFTN2/2N(N/2–1)一個蝶形12N/2個蝶形N/2N總計N2/2+N/2≈N2/2N(N/2-1)+N≈N2/2運算量減少(jiǎnshǎo)了近一半分解(fēnjiě)后的運算量：精品資料先將N=8點的DFT分解成2個4點DFT：可知(kězhī)：時域上：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點FFT變換

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

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

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

原位運算(亦稱同址計算)精品資料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)回顧(huígù)：N點DIT―FFT運算流圖(N=8)精品資料如上所述，N點DIT―FFT運算(yùnsuàn)流圖中，每級都有N/2個蝶形。每個蝶形都要乘以因子WNP，稱其為旋轉因子，p稱為旋轉因子的指數(shù)。2)旋轉(xuánzhuǎn)因子的變化規(guī)律觀察FFT運算流圖發(fā)現(xiàn)，第L級共有2L-1個不同的旋轉因子。N=23=8時的各級旋轉因子表示如下：L=1時，WNp=WN/4J，N/4=21=2L，J=0L=2時，WNp=WN/2J，N/2=22=2L，J=0，1L=3時，WNp=WNJ，N=23=2L，J=0，1，2，3精品資料對N=2M的一般情況(qíngkuàng)，第L級的旋轉因子為：精品資料設序列x(n)經(jīng)時域抽選(倒序)后，存入(cúnrù)數(shù)組X中。如果蝶形運算的兩個輸入數(shù)據(jù)相距B個點，應用原位計算，則蝶形運算可表示成如下形式：下標(xiàbiāo)L表示第L級運算，XL(J)則表示第L級運算后數(shù)組元素X(J)的值。精品資料3)編程思想(sīxiǎng)及流程圖開始送入x(n)和N=2M調整輸入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)輸出結果結束精品資料4）碼位倒序(dǎoxù)由N=8蝶形圖看出：原位計算時，F(xiàn)FT輸出的X(k)的次序正好是順序排列的，即X(0)…X(7),但輸入x(n)都不能按自然順序存入到存儲單元(cúnchǔdānyuán)中，而是按x(0),x(4),x(2),x(6),x(1),x(5),x(3),x(7)的順序存入存儲單元(cúnchǔdānyuán)，即為亂序輸入，順序輸出。這種順序看起來相當雜亂，然而它是有規(guī)律的。即碼位倒讀規(guī)則。精品資料自然順序n二進制碼表示碼位倒讀碼位倒置順序n’以N=8為例：0123456700000101001110010111011100010001011000110101111104261537看出：碼位倒讀后的順序剛好(gānghǎo)是數(shù)據(jù)送入計算機內的順序。精品資料倒序(dǎoxù)規(guī)律精品資料對于(duìyú)數(shù)N，在其二進制最高位加1，等于加N/2。若已知某個反序號為J，為求下一個反序號，可先判J的最高位：1)若為0，則把該位變成1（即加N/2）就得到下一個反序號，2)若為1，則需判斷次高位：①若次高位為0，則把最高位變0（相當(xiāngdāng)減去N/2）后，再把次高位變1（即加N/4）。②若次高位為1，則需判斷次次高位……分析：精品資料倒序排列算法的流程圖正序序列已在數(shù)組A[]中，輸入NLH=N/2,J=LH,N1=N-2J=J-kk=k/2k=LHJ<kJ=J+kT=A(I)A(I)=A(J)A(J)=Tfor(i=1;i<=N1;i++)i≥JNYYN精品資料第四節(jié)按頻率(pínlǜ)抽選的基2-FFT算法在基2快速算法中，頻域抽取法FFT也是一種(yīzhǒnɡ)常用的快速算法，簡稱DIF―FFT。設序列x(n)長度為N=2M，首先將x(n)前后對半分開，得到兩個子序列，其DFT可表示為如下形式精品資料精品資料精品資料精品資料DIF―FFT一次分解(fēnjiě)運算流圖(N=8)4點DFT4點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)精品資料DIF―FFT二次分解(fēnjiě)運算流圖(N=8)精品資料DIF―FFT運算(yùnsuàn)流圖(N=8)精品資料時間抽取(chōuqǔ)算法與頻率抽取(chōuqǔ)算法的比較1)頻率抽選法和時間抽選法總的計算(jìsuàn)量是相同的復乘：復加：2)頻率抽取法和時間抽取法一樣，都適用于原位運算，即蝶形的輸入和輸出占用同一個存儲單元。3)均存在碼位倒序問題。4)頻率抽選法和時間抽選法一樣，基本運算也是蝶形運算。但兩者的蝶形形式略有不同。精品資料第五節(jié)IDFT的快速(kuàisù)算法－IFFT上述FFT算法流圖也可以用于離散(lísàn)傅里葉逆變換(InverseDiscreteFourierTransform，簡稱IDFT)。比較DFT和IDFT的運算公式：1)旋轉因子：2)系數(shù)：精品資料DIT―IFFT運算(yùnsuàn)流圖精品資料DIT―IFFT運算(yùnsuàn)流圖(防止溢出)精品資料如果希望直接調用FFT子程序計算IFFT，則可用下面(xiàmian)的方法：對上式兩邊(liǎngbiān)同時取共軛，得：精品資料例1、如果通用計算機的速度(sùdù)為平均每次復乘需要5s，每次復加需要0.5s，用它來計算512點的DFT[x(n)]，問：1）直接計算(jìsuàn)需要多少時間？2）用FFT需要多少時間？解：1）用DFT進行運算：復乘：T1=N2×5×10-6=1.31072秒復加：T2=N(N-1)×0.5×10-6=0.130816秒總共：T=T1+T2=1.441536秒精品資料2）用FFT進行(jìnxíng)運算：復乘：T1’=(N/2)log2N×5×10-6=0.01152秒復加：T2’=Nlog2N×0.5×10-6=0.002304秒總共(zǒnggòng)：T’=T1’+T2’=0.013824秒精品資料例2、對一個連續(xù)(liánxù)時間信號xa(t)采樣1秒得到4096個采樣點的序列，求：1）若采樣(cǎiyànɡ)后沒有發(fā)生混疊現(xiàn)象，xa(t)的最高頻率是多少？解：1秒內采樣4096個點，說明采樣頻率是4096Hz。精品資料2）若計算(jìsuàn)采樣信號的4096點DFT，DFT系數(shù)之間的頻率間隔是多少？解：(要求解的是頻譜分辨(fēnbiàn)的間隔F)精品資料例3、長度為240點的序列x(n)與長度為N點的h(n)卷積。當N=10和240時，直接進行(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)補L-N1個零x(n)L點DFT補L-N2個零h(n)L點DFTL點IDFTy(n)=x(n)*h(n)H(k)Y(k)精品資料直接(zhíjiē)進行卷積（N=10）：乘法(chéngfǎ)次數(shù)：240×10=2400次用FFT的方法（N=10）：添零到256點，L=256乘法次數(shù)：3×(L/2)log2L＋L=3×128×8＋256=3328次精品資料直接(zhíjiē)進行卷積（N=240）：乘法(chéngfǎ)次數(shù)：240×240=57600次用FFT的方法（N=240）：添零到512點，L=512乘法次數(shù)：3×(L/2)log2L＋L=3×256×9＋512=7424次精品資料直接DFT方法(fāngfǎ)/CZT方法(fāngfǎ)：當要求準確的N點DFT，且N是素數(shù)時第六節(jié)N為復合(fùhé)數(shù)的FFT算法

——混合基算法基-2FFT算法：補零使?jié)M足混合基FFT算法：N是復合數(shù)精品資料1、整數(shù)的多基多(jīduō)進制表示形式（1）二進制：

精品資料（2）r進制：

精品資料（3）多基多(jīduō)進制（混合基）：精品資料例：精品資料精品資料2、的快速算法行列行序號列序號行列行變量列變量精品資料精品資料的DFT算法(1)改寫(gǎixiě)成(2)做個點DFT，得為參量，輸入變量(biànliàng)，輸出變量(biànliàng)的點DFT(3)N個（旋轉因子）(4)做個點DFT，得為參量，輸入變量，輸出變量的點DFT(5)整序精品資料例精品資料精品資料當N為高組合素數(shù)時：個

點DFT，乘以旋轉因子

個

點DFT個

點DFT，乘以旋轉因子個

點DFT，乘以旋轉因子L級r點DFT稱基

算法，

基

算法

混合基算法（基算法）基算法精品資料混合(hùnhé)基算法的運算量不計譯序、整序工作量（2）乘N個旋轉(xuánzhuǎn)因子復乘N總計：（1）個點DFT復乘復加（3）個點DFT復乘復加精品資料混合(hùnhé)基節(jié)省的運算量次乘N個旋轉因子個

點DFT個

點DFT

個

點DFT精品資料第七節(jié)基-4FFT算法(suànfǎ)

當混合基FFT算法中

時，

即為基-4FFT算法，n、k都為4進制數(shù)個

點DFT

乘N個旋轉因子個

點DFT

乘N個旋轉因子個

點DFT精品資料精品資料1)的4點DFT精品資料精品資料的四進制數(shù)按二進制倒位序排列成精品資料3)的4點DFT精品資料精品資料一個4點FFT不需乘法，只需3次乘旋轉因子（除外）而基-2FFT基-4FFT運算量：每級有N/4個4點FFT，共L級（L-1級要乘旋轉(xuánzhuǎn)因子）精品資料第八節(jié)分裂(fēnliè)基FFT算法對偶序號使用(shǐyòng)基-2FFT算法，對奇序號使用(shǐyòng)基-4FFT算法，稱分裂基FFT算法針對

的算法中具有最少乘法次數(shù)，且同址運算。將分成三個序列。精品資料偶序號的點DFT

奇序號的點DFT

精品資料利用周期性

分成四段：精品資料精品資料的第一級分解得4個分裂基同樣(tóngyàng)對、、作進一步分解。和的第二級分解(fēnjiě)分別是基-4的4點DFT。的第二級分解得2個分裂基。一個基-4的4點DFT和2個基-2的4點DFT。精品資料精品資料基