1、第四章 快速傅立葉變換 Fast Fourier Transform,第一節(jié) 直接計(jì)算DFT的問題及改進(jìn)途徑,1、問題的提出,設(shè)有限長序列x(n)，非零值長度為N，若對x(n)進(jìn)行一次DFT運(yùn)算，共需多大的運(yùn)算工作量？,計(jì)算成本? 計(jì)算速度?,2. DFT的運(yùn)算量,回憶DFT和IDFT的變換式：,計(jì)算機(jī)運(yùn)算時(shí)（編程實(shí)現(xiàn)）：,以DFT為例：,運(yùn)算量,(a+jb)(c+jd)=(ac-bd)+j(bc+ad),例：計(jì)算一個(gè) N點(diǎn)DFT ，共需N2次復(fù)乘。以做一次 復(fù)乘1s計(jì)，若N =4096，所需時(shí)間為,例：石油勘探，有24個(gè)通道的記錄，每通道波形記 錄長度為5秒，若每秒抽樣500點(diǎn)/秒， 1）每

2、道總抽樣點(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次,由于計(jì)算量大，且要求相當(dāng)大的內(nèi)存，難以實(shí)現(xiàn)實(shí)時(shí)處理，限制了DFT的應(yīng)用。長期以來，人們一直在尋求一種能提高DFT運(yùn)算速度的方法。,FFT便是 Cooley 當(dāng)n=奇數(shù)時(shí)，令n=2r+1;,分組，變量置換,2、算法步驟,得到：,帶入DFT中,所以,由于,?,X1(k)、X2(k)只有N/2個(gè)點(diǎn)，以N/2為周期；而X (k)卻有N個(gè)點(diǎn)，以N為周期。要用X1(k)、X2(k)表達(dá)全部的X (k) 值，還必須利用WN系數(shù)的周期特性。,又考慮到 的對稱性：,

3、有：,蝶形運(yùn)算流圖符號(hào),說明： (1) 左邊兩路為輸入 (2) 右邊兩路為輸出 (3) 中間以一個(gè)小圓表示加、 減運(yùn)算（右上路為相加 輸出、右下路為相減輸 出),1個(gè)蝶形運(yùn)算需要1次復(fù)乘，2次復(fù)加,運(yùn)算量減少了近一半,分解后的運(yùn)算量：,先將N=8點(diǎn)的DFT分解成2個(gè)4點(diǎn)DFT： 可知：時(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)給出,例子：求 N=23=8點(diǎn)FFT變換 按N=8N/2=4，做4點(diǎn)的DFT：,N=8點(diǎn)的直接DFT的計(jì)算量為： 復(fù)乘：N2次

4、= 64次 復(fù)加：N(N-1)次 = 87=56次,此外，還有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次,N點(diǎn)DFT的一次時(shí)域抽取分解圖(N=8),因?yàn)?點(diǎn)DFT還是比較麻煩，所以再繼續(xù)分解。,若將N/2(4點(diǎn))子序列按奇/偶分解成兩個(gè)N/4點(diǎn)(2點(diǎn))子序列。即對將x1(r)和x2(r)分解成奇、偶兩個(gè)N/4點(diǎn)(2點(diǎn)

5、)點(diǎn)的子序列。,那么，X1(k)又可表示為,X2(k)也可以進(jìn)行相同的分解：,注意：通常我們會(huì)把 寫成 。,N點(diǎn)DFT的第二次時(shí)域抽取分解圖(N=8),8,8,N點(diǎn)DITFFT運(yùn)算流圖(N=8),3、DITFFT算法與直接計(jì)算DFT運(yùn)算量的比較,1)、N=2M的DFT運(yù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í)，有：,FFT算法與直接計(jì)算DFT所需乘法次數(shù)的比較曲線,4、DITFFT的運(yùn)算規(guī)律及編程思想,FFT的每級(jí)（列）計(jì)算都是由N個(gè)復(fù)數(shù)數(shù)據(jù)（輸入）兩

6、兩構(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) + W80 x(4) 送回A(0)單元,將x(0) - W80 x(4) 送回A(1)單元,1） 原位運(yùn)算 (亦稱同址計(jì)算),回顧：N點(diǎn)DITFFT運(yùn)算流圖(N=8),如上所述，N點(diǎn)DITFFT運(yùn)算流圖中，每級(jí)都有N/2個(gè)蝶形。每個(gè)蝶形都要乘以因子WNP，稱其為旋轉(zhuǎn)因子，p稱為旋轉(zhuǎn)因子的指數(shù)。,2)旋轉(zhuǎn)因子的變化規(guī)律,觀察FFT運(yùn)算流圖發(fā)現(xiàn)

7、，第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=0 L=2時(shí)， WNp =WN/2J， N/2 =22 =2L， J=0，1 L=3時(shí)， WNp =WNJ， N =23 =2L， J=0，1，2，3,對N=2M的一般情況，第L級(jí)的旋轉(zhuǎn)因子為：,設(shè)序列x(n)經(jīng)時(shí)域抽選(倒序)后，存入數(shù)組X中。如果蝶形運(yùn)算的兩個(gè)輸入數(shù)據(jù)相距B個(gè)點(diǎn)(B=2L-1)，應(yīng)用原位計(jì)算，則蝶形運(yùn)算可表示成如下形式：,下標(biāo)L表示第L級(jí)運(yùn)算，XL(J)則表示第L級(jí)運(yùn)算后數(shù)組元素X(J)的值。,3) 編程思想及流程圖,4）碼位倒序,由N

8、=8蝶形圖看出：原位計(jì)算時(shí)，F(xiàn)FT輸出的X(k)的次序正好是順序排列的，即X(0)X(7),但輸入x(n)都不能按自然順序存入到存儲(chǔ)單元中，而是按x(0),x(4),x(2),x(6) ,x(1),x(5),x(3),x(7)的順序存入存儲(chǔ)單元，即為亂序輸入，順序輸出。 這種順序看起來相當(dāng)雜亂，然而它是有規(guī)律的。即碼位倒讀規(guī)則。,以N=8為例：,看出：碼位倒讀后的順序剛好是數(shù)據(jù)送入計(jì)算機(jī)內(nèi)的順序。,倒序規(guī)律,對于數(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，

9、則需判斷次高位： 若次高位為0，則把最高位變0（相當(dāng)減去 N/2）后，再把次高位變1（即加N/4）。 若次高位為1，則需判斷次次高位,分析：,倒 序 排 列 算 法 的 流 程 圖,第四節(jié) 按頻率抽選的基2-FFT算法,在基2快速算法中，頻域抽取法FFT也是一種常用的快速算法，簡稱DIFFFT。,設(shè)序列x(n)長度為N=2M，首先將x(n)前后對半分開，得到兩個(gè)子序列，其DFT可表示為如下形式,DIFFFT一次分解運(yùn)算流圖(N=8),DIFFFT二次分解運(yùn)算流圖(N=8),DIFFFT運(yùn)算流圖(N=8),時(shí)間抽取算法與頻率抽取算法的比較,1) 頻率抽選法和時(shí)間抽選法總的計(jì)算量是相同的,2) 頻

10、率抽取法和時(shí)間抽取法一樣，都適用于原位運(yùn) 算， 即蝶形的輸入和輸出占用同一個(gè)存儲(chǔ)單元。,3) 均存在碼位倒序問題。,4) 頻率抽選法和時(shí)間抽選法一樣，基本運(yùn)算也是蝶形 運(yùn)算。但兩者的蝶形形式略有不同。,第五節(jié) IDFT的快速算法IFFT,上述FFT算法流圖也可以用于離散傅里葉逆變換(Inverse Discrete Fourier Transform，簡稱IDFT)。,比較DFT和IDFT的運(yùn)算公式：,1) 旋轉(zhuǎn)因子：,2) 系數(shù)：,DITIFFT運(yùn)算流圖,DITIFFT運(yùn)算流圖(防止溢出),如果希望直接調(diào)用FFT子程序計(jì)算IFFT，則可用下面的方法：,對上式兩邊同時(shí)取共軛，得：,例1、如果通

11、用計(jì)算機(jī)的速度為平均每次復(fù)乘需要 5s，每次復(fù)加需要0.5s，用它來計(jì)算512點(diǎn) 的DFTx(n)，問：,1）直接計(jì)算需要多少時(shí)間？ 2）用FFT需要多少時(shí)間？,解：1）用DFT進(jìn)行運(yùn)算：,復(fù)乘：T1=N2510-6=1.31072秒,復(fù)加：T2=N(N-1)0.510-6=0.130816秒,總共：T=T1+T2=1.441536秒,2）用FFT進(jìn)行運(yùn)算：,復(fù)乘：T1=(N/2)log2N510-6=0.01152秒,復(fù)加：T2= Nlog2N 0.510-6=0.002304秒,總共：T=T1+T2=0.013824秒,例2、對一個(gè)連續(xù)時(shí)間信號(hào)xa(t)采樣1秒得到4096個(gè)采 樣點(diǎn)的序列，求：,1）若采樣后沒有發(fā)生混疊現(xiàn)象，xa(t)的最高頻率是 多少？,解：1秒內(nèi)采樣4096個(gè)點(diǎn)，說明采樣頻率是4096Hz。,2）若計(jì)算采樣信號(hào)的4096點(diǎn)DFT，DFT系數(shù)之間 的頻率間隔是多少？,解：(要求解的是頻譜分辨的間隔F),例3、長度為240點(diǎn)的序列x(n)與長度為N點(diǎn)的h(n)卷 積。當(dāng)N=10和240時(shí)，直接進(jìn)行卷積 x(n)*h(n) 和用 IFFTX(K)H(K) 的方法相比，那種方法求 解y(n)的效率更高？,LN1+N2-1,直接進(jìn)行卷積（N=10）：,乘法