第五章-快速傅立葉變換(FFT)課件

第五章快速傅立葉變換(FFT)5.1引言5.2直接計算DFT的特點及減少運算量的基本途徑5.3基2FFT算法第五章快速傅立葉變換(FFT)5.1引言15.1引言 影響數(shù)字信號處理發(fā)展的最主要因素之一是處理速度。DFT使計算機在頻域處理信號成為可能，但是當(dāng)N很大時，直接計算N點DFT的計算量非常大。快速傅立葉變換(FFT：FastFourierTransform)可使實現(xiàn)DFT的運算量下降幾個數(shù)量級，從而使數(shù)字信號處理的速度大大提高。自從1965年第一篇DFT快速算法的論文發(fā)表以來，人們已經(jīng)研究出多種FFT算法，它們的復(fù)雜度和運算效率各不相同。本章主要介紹最基本的基2FFT算法及其編程方法。5.1引言 影響數(shù)字信號處理發(fā)展的最主要因素之一是處理25.2直接計算DFT的特點及減少運算量的基本途徑 長度為N的序列x(n)的N點DFT為 的周期性：k=0,1,…,N-1(5.2.1)5.2直接計算DFT的特點及減少運算量的基本途徑 長度3 的對稱性： (5.2.2)

(5.2.3) 的對稱性：45.3基2FFT算法 5.3.1DIT-FFT算法 序列x(n)的N(N=2M)點DFT為k=0,1,…,N-15.3基2FFT算法 5.3.1DIT-FFT算法5 將上面的和式按n的奇偶性分解為 令x1(l)=x(2l),x2(l)=x(2l+1)。 因為W2klN=WklN/2,所以上式可寫成(5.3.1) 將上面的和式按n的奇偶性分解為(5.3.1)6 (5.3.1)式說明，按n的奇偶性將x(n)分解為兩個N/2長的序列x1(l)和x2(l)，則N點DFT可分解為兩個N/2點DFT來計算。用X1(k)和X2(k)分別表示

將以上兩式代入(5.3.1)式，并利用和X1(k)、X2(k)的隱含周期性可得到： (5.3.1)式說明，按n的奇偶性將x(n)分解為兩個N7

這樣，就將N點DFT的計算分解為計算兩個N/2點離散傅立葉變換X1(k)和X2(k)。為了將如上分解過程用運算流圖表示，以便估計其運算量，觀察運算規(guī)律，總結(jié)編程方法，先介紹一種表示上式的蝶形圖。 蝶形圖及其運算功能如圖5.3.1所示。

8圖5.3.1蝶形運算圖圖5.3.1蝶形運算圖9

根據(jù)圖5.3.2可以求得第一次分解后的運算量。圖5.3.2包括兩個N/2點DFT和N/2個蝶形，每個點DFT需要次復(fù)數(shù)乘法和次復(fù)數(shù)加法運算，每個蝶形只有一次復(fù)數(shù)乘法運算和兩次復(fù)數(shù)加法運算。所以，總的復(fù)數(shù)乘法次數(shù)為 總的復(fù)數(shù)加法次數(shù)為

根據(jù)圖5.3.2可以求得第一次分解后的運算量。圖5.310圖5.3.28點DFT一次時域抽取分解運算流圖圖5.3.28點DFT一次時域抽取分解運算流圖11 N=8點DIT-FFT的運算流圖如圖5.3.3(a)所示。根據(jù)WkN/m=WkmN，將圖5.3.3(a)轉(zhuǎn)換成如圖5.3.3(b)所示的標(biāo)準(zhǔn)形式的運算流圖。 N=8點DIT-FFT的運算流圖如圖5.3.3(a)所示12（a）圖5.3.38點DIT-FFT運算流圖x(0)x(4)x(2)x(6)x(1)x(5)x(3)x(7)A(0)A(1)A(2)A(3)A(4)A(5)A(6)A(7)A(0)A(1)A(2)A(3)A(4)A(5)A(6)A(7)X(0)X(1)X(2)X(3)X(4)X(5)X(6)X(7)A(0)A(1)A(2)A(3)A(4)A(5)A(6)A(7)A(0)A(1)A(2)A(3)A(4)A(5)A(6)A(7)（a）x(0)x(4)x(2)x(6)x(1)x(5)x(313x(0)x(4)x(2)x(6)x(1)x(5)x(3)x(7)A(0)A(1)A(2)A(3)A(4)A(5)A(6)A(7)A(0)A(1)A(2)A(3)A(4)A(5)A(6)A(7)X(0)X(1)X(2)X(3)X(4)X(5)X(6)X(7)A(0)A(1)A(2)A(3)A(4)A(5)A(6)A(7)A(0)A(1)A(2)A(3)A(4)A(5)A(6)A(7)（b）圖5.3.38點DIT-FFT運算流圖x(0)x(4)x(2)x(6)x(1)x(5)x(3)x(14 5.3.2DIT-FFT的運算效率 DIT-FFT的運算效率指直接計算DFT的運算量與DIT-FFT的運算量之比。 由圖5.3.3可見，N=2M時，其DIT-FFT運算流圖由M級蝶形構(gòu)成，每級有N/2個蝶形。因此，每級需要N/2次復(fù)數(shù)乘法運算和N次復(fù)數(shù)加法運算，M級形共需復(fù)數(shù)乘法次數(shù)CM(2)和復(fù)數(shù)加法次數(shù)CA(2)分別為 (5.3.5)

(5.3.6) 5.3.2DIT-FFT的運算效率15 直接計算N點DFT的復(fù)數(shù)乘法次數(shù)為N2,復(fù)數(shù)加法次數(shù)為(N-1)。當(dāng)N>>1時,N2/CM(2)>>1，所以N越大，DIT-FFT運算效率越高。DIT-FFT算法與DFT所需乘法次數(shù)與N的關(guān)系曲線如圖5.3.4所示。例如，N=210=1024時，DIT-FFT的運算效率為 而當(dāng)N=211=2048時， 直接計算N點DFT的復(fù)數(shù)乘法次數(shù)為N2,復(fù)數(shù)加法次數(shù)為(16圖5.3.4DIT-FFT與DFT所需乘法次數(shù)比較曲線圖5.3.4DIT-FFT與DFT所需乘法次數(shù)比較曲線175.3.3DIT-FFT的運算規(guī)律及編程思想1.原位計算由圖5.3.3可以看出，DIT-FFT的運算過程很有規(guī)律。2.旋轉(zhuǎn)因子的變化規(guī)律觀察圖5.3.3(a)不難發(fā)現(xiàn)，第L級共有2L-1個不同的旋轉(zhuǎn)因子。N=23=8時的各級旋轉(zhuǎn)因子表示如下：L=1時，J=0L=2時，J=0,1L=3時，J=0,1,2,35.3.3DIT-FFT的運算規(guī)律及編程思想18對N=2M的一般情況，第L級的旋轉(zhuǎn)因子為 J=0,1,2,…,2L-1-1由于2L=2M×2L-M=N·2L-M所以 J=0,1,2,…,2L-1-1，(5.3.7)(5.3.8)這樣，就可按(5.3.7)式和(5.3.8)式確定第L級運算的旋轉(zhuǎn)因子(實際編程序時，L為最外層循環(huán)變量)。對N=2M的一般情況，第L級的旋轉(zhuǎn)因子為19 3.蝶形運算規(guī)律 設(shè)序列x(n)經(jīng)時域抽選(倒序)后，存入數(shù)組X中。如果蝶形運算的兩個輸入數(shù)據(jù)相距B個點，應(yīng)用原位計算，則蝶形運算可表示成如下形式：XL(J)XL-1(J)+XL-1(J+B)WpN

XL(J+B)XL-1(J)-XL-1(J+B)WpN 式中

p=J·2M-L;J=0,1,…,2L-1-1;L=1,2,…,M 3.蝶形運算規(guī)律20 設(shè)

式中，下標(biāo)R表示取實部，I表示取虛部， 設(shè)21 4.編程思想及程序框圖 仔細觀察圖5.3.3，還可歸納出一些編程序有用的運算規(guī)律：第L級中，每個蝶形的兩個輸入數(shù)據(jù)相距B=2L-1個點；同一旋轉(zhuǎn)因子對應(yīng)著間隔為2L點的2M-L個蝶形。 總結(jié)上述運算規(guī)律，便可采用下述運算方法。先從輸入端(第1級)開始，逐級進行，共進行M級運算。當(dāng)進行第L級運算時，依次求出2L-1個不同的旋轉(zhuǎn)因子，每求出一個旋轉(zhuǎn)因子，就計算完它對應(yīng)的所有2M-L個蝶形。這樣，我們可用三重循環(huán)程序?qū)崿F(xiàn)DIT-FFT運算，程序框圖如圖5.3.5所示。 4.編程思想及程序框圖22圖5.3.5DIT-FFT運算和程序框圖圖5.3.5DIT-FFT運算和程序框圖23 5.序列的倒序 DIT-FFT算法的輸入序列的排序看起來似乎很亂，但仔細分析就會發(fā)現(xiàn)這種倒序是很有規(guī)律的。由于N=2M，因此順序數(shù)可用M位二進制數(shù)(nM-1

nM-2…n1n0)表示。M次偶奇時域抽選過程如圖5.3.6所示。 第一次按最低位n0的0和1將x(n)分解為偶奇兩組，第二次又按次低位n1的0、1值分別對偶奇組分解；依次類推，第M次按nM-1位分解，最后所得二進制倒序數(shù)如圖5.3.6所示。表5.3.1列出了N=8時以二進制數(shù)表示的順序數(shù)和倒序數(shù)。 5.序列的倒序24圖5.3.6形成例序的樹狀圖(N=23)圖5.3.6形成例序的樹狀圖(N=23)25表5.3.1順序和倒序二進制數(shù)對照表

表5.3.1順序和倒序二進制數(shù)對照表26 形成倒序J后，將原存儲器中存放的輸入序列重新按倒序排列。設(shè)原輸入序列x(n)先按自然順序存入數(shù)組A中。例如，對N=8,A(0)，A(1)，A(2)，…，A(7)中依次存放著x(0),x(1),…,x(7)。對x(n)的重新排序(倒序)規(guī)律如圖5.3.7所示。倒序的程序框圖如圖5.3.8所示，圖5.3.8中的虛線框內(nèi)是完成計算倒序值的運算流程圖。 形成倒序J后，將原存儲器中存放的輸入序列重新按倒序排列。27圖5.3.7倒序規(guī)律圖5.3.7倒序規(guī)律28圖5.3.8倒序程序框圖圖5.3.8倒序程序框圖29 5.3.4DIF-FFT 在基2快速算法中，頻域抽取法FFT也是一種常用的快速算法，簡稱DIF-FFT。 設(shè)序列x(n)長度為N=2M,首先將x(n)前后對半分開，得到兩個子序列，其DFT可表示為如下形式： 5.3.4DIF-FFT30 式中

WkN/2N=(-1)k=

將X(k)分解成偶數(shù)組與奇數(shù)組，當(dāng)k取偶數(shù)(k=2r,r=0,1,…,N/2-1)時，1k=偶數(shù)-1k=奇數(shù)(5.3.9) 式中1k=偶數(shù)(5.331 當(dāng)k取奇數(shù)(k=2r+1,r=0,1,…,N/2-1)時 令(5.3.10) 當(dāng)k取奇數(shù)(k=2r+1,r=0,1,…,N/232 將x1(n)和x2(n)分別代入(5.3.9)式和(5.3.10)式，可得

x1(n)、x2(n)和x

(n)的關(guān)系也可用圖5.3.9所示的蝶形運算流圖符號表示。圖5.3.10表示N=8時一次分解的運算流圖。(5.3.11) 將x1(n)和x2(n)分別代入(5.3.9)式和(5.33圖5.3.9DIF-FFT蝶形運算流圖符號圖5.3.9DIF-FFT蝶形運算流圖符號34圖5.3.10DIF-FFT一次分解運算流圖(N=8)圖5.3.10DIF-FFT一次分解運算流圖(N=835 圖5.3.11示出了N=8時二次分解運算流圖。這樣繼續(xù)分解下去，經(jīng)過M-1次分解，最后分解為2M-1個2點DFT，2點DFT就是一個基本蝶形運算流圖。當(dāng)N=8時，經(jīng)兩次分解，便分解為四個2點DFT,如圖5.3.11所示。N=8的完整DIF-FFT運算流圖如圖5.3.12所示。 圖5.3.11示出了N=8時二次分解運算流圖。這樣繼續(xù)分36圖5.3.11DIF-FFT二次分解運算流圖(N=8)圖5.3.11DIF-FFT二次分解運算流圖(N=8)37圖5.3.12DIF-FFT運算流圖(N=8)圖5.3.12DIF-FFT運算流圖(N=8)38 5.3.5IDFT的高效算法 上述FFT算法流圖也可以用于離散傅立葉逆變換(IDFT:InverseDiscreteFourierTransform)。比較DFT和IDFT的運算公式： 由DIF-FFT運算流圖改成的DIT-IFFT運算流圖如圖5.3.13所示。 5.3.5IDFT的高效算法39圖5.3.13DIT-IFFT運算流圖圖5.3.13DIT-IFFT運算流圖40 在實際中，有時為了防止運算過程中發(fā)生溢出，將1/N分配到每一級蝶形運算中。由于1/N=(1/2)M,因此每級的每個蝶形輸出支路均有一相乘因子1/2，這種運算的蝶形流圖如圖5.3.14所示。由圖可知，乘法次數(shù)比圖5.3.13增加了(N/2)(M-1)次。 在實際中，有時為了防止運算過程中發(fā)生溢出，將1/N分配到41圖5.3.14DIT-IFFT運算流圖(防止溢出)圖5.3.14DIT-IFFT運算流圖(防止溢出)42 如果希望直接調(diào)用FFT子程序計算IFFT，則可用下面的方法： 對上式兩邊同時取共軛，得由于因此 如果希望直接調(diào)用FFT子程序計算IFFT，則可用下面的方43第五章快速傅立葉變換(FFT)5.1引言5.2直接計算DFT的特點及減少運算量的基本途徑5.3基2FFT算法第五章快速傅立葉變換(FFT)5.1引言445.1引言 影響數(shù)字信號處理發(fā)展的最主要因素之一是處理速度。DFT使計算機在頻域處理信號成為可能，但是當(dāng)N很大時，直接計算N點DFT的計算量非常大。快速傅立葉變換(FFT：FastFourierTransform)可使實現(xiàn)DFT的運算量下降幾個數(shù)量級，從而使數(shù)字信號處理的速度大大提高。自從1965年第一篇DFT快速算法的論文發(fā)表以來，人們已經(jīng)研究出多種FFT算法，它們的復(fù)雜度和運算效率各不相同。本章主要介紹最基本的基2FFT算法及其編程方法。5.1引言 影響數(shù)字信號處理發(fā)展的最主要因素之一是處理455.2直接計算DFT的特點及減少運算量的基本途徑 長度為N的序列x(n)的N點DFT為 的周期性：k=0,1,…,N-1(5.2.1)5.2直接計算DFT的特點及減少運算量的基本途徑 長度46 的對稱性： (5.2.2)

(5.2.3) 的對稱性：475.3基2FFT算法 5.3.1DIT-FFT算法 序列x(n)的N(N=2M)點DFT為k=0,1,…,N-15.3基2FFT算法 5.3.1DIT-FFT算法48 將上面的和式按n的奇偶性分解為 令x1(l)=x(2l),x2(l)=x(2l+1)。 因為W2klN=WklN/2,所以上式可寫成(5.3.1) 將上面的和式按n的奇偶性分解為(5.3.1)49 (5.3.1)式說明，按n的奇偶性將x(n)分解為兩個N/2長的序列x1(l)和x2(l)，則N點DFT可分解為兩個N/2點DFT來計算。用X1(k)和X2(k)分別表示

將以上兩式代入(5.3.1)式，并利用和X1(k)、X2(k)的隱含周期性可得到： (5.3.1)式說明，按n的奇偶性將x(n)分解為兩個N50

這樣，就將N點DFT的計算分解為計算兩個N/2點離散傅立葉變換X1(k)和X2(k)。為了將如上分解過程用運算流圖表示，以便估計其運算量，觀察運算規(guī)律，總結(jié)編程方法，先介紹一種表示上式的蝶形圖。 蝶形圖及其運算功能如圖5.3.1所示。

51圖5.3.1蝶形運算圖圖5.3.1蝶形運算圖52

根據(jù)圖5.3.2可以求得第一次分解后的運算量。圖5.3.2包括兩個N/2點DFT和N/2個蝶形，每個點DFT需要次復(fù)數(shù)乘法和次復(fù)數(shù)加法運算，每個蝶形只有一次復(fù)數(shù)乘法運算和兩次復(fù)數(shù)加法運算。所以，總的復(fù)數(shù)乘法次數(shù)為 總的復(fù)數(shù)加法次數(shù)為

根據(jù)圖5.3.2可以求得第一次分解后的運算量。圖5.353圖5.3.28點DFT一次時域抽取分解運算流圖圖5.3.28點DFT一次時域抽取分解運算流圖54 N=8點DIT-FFT的運算流圖如圖5.3.3(a)所示。根據(jù)WkN/m=WkmN，將圖5.3.3(a)轉(zhuǎn)換成如圖5.3.3(b)所示的標(biāo)準(zhǔn)形式的運算流圖。 N=8點DIT-FFT的運算流圖如圖5.3.3(a)所示55（a）圖5.3.38點DIT-FFT運算流圖x(0)x(4)x(2)x(6)x(1)x(5)x(3)x(7)A(0)A(1)A(2)A(3)A(4)A(5)A(6)A(7)A(0)A(1)A(2)A(3)A(4)A(5)A(6)A(7)X(0)X(1)X(2)X(3)X(4)X(5)X(6)X(7)A(0)A(1)A(2)A(3)A(4)A(5)A(6)A(7)A(0)A(1)A(2)A(3)A(4)A(5)A(6)A(7)（a）x(0)x(4)x(2)x(6)x(1)x(5)x(356x(0)x(4)x(2)x(6)x(1)x(5)x(3)x(7)A(0)A(1)A(2)A(3)A(4)A(5)A(6)A(7)A(0)A(1)A(2)A(3)A(4)A(5)A(6)A(7)X(0)X(1)X(2)X(3)X(4)X(5)X(6)X(7)A(0)A(1)A(2)A(3)A(4)A(5)A(6)A(7)A(0)A(1)A(2)A(3)A(4)A(5)A(6)A(7)（b）圖5.3.38點DIT-FFT運算流圖x(0)x(4)x(2)x(6)x(1)x(5)x(3)x(57 5.3.2DIT-FFT的運算效率 DIT-FFT的運算效率指直接計算DFT的運算量與DIT-FFT的運算量之比。 由圖5.3.3可見，N=2M時，其DIT-FFT運算流圖由M級蝶形構(gòu)成，每級有N/2個蝶形。因此，每級需要N/2次復(fù)數(shù)乘法運算和N次復(fù)數(shù)加法運算，M級形共需復(fù)數(shù)乘法次數(shù)CM(2)和復(fù)數(shù)加法次數(shù)CA(2)分別為 (5.3.5)

(5.3.6) 5.3.2DIT-FFT的運算效率58 直接計算N點DFT的復(fù)數(shù)乘法次數(shù)為N2,復(fù)數(shù)加法次數(shù)為(N-1)。當(dāng)N>>1時,N2/CM(2)>>1，所以N越大，DIT-FFT運算效率越高。DIT-FFT算法與DFT所需乘法次數(shù)與N的關(guān)系曲線如圖5.3.4所示。例如，N=210=1024時，DIT-FFT的運算效率為 而當(dāng)N=211=2048時， 直接計算N點DFT的復(fù)數(shù)乘法次數(shù)為N2,復(fù)數(shù)加法次數(shù)為(59圖5.3.4DIT-FFT與DFT所需乘法次數(shù)比較曲線圖5.3.4DIT-FFT與DFT所需乘法次數(shù)比較曲線605.3.3DIT-FFT的運算規(guī)律及編程思想1.原位計算由圖5.3.3可以看出，DIT-FFT的運算過程很有規(guī)律。2.旋轉(zhuǎn)因子的變化規(guī)律觀察圖5.3.3(a)不難發(fā)現(xiàn)，第L級共有2L-1個不同的旋轉(zhuǎn)因子。N=23=8時的各級旋轉(zhuǎn)因子表示如下：L=1時，J=0L=2時，J=0,1L=3時，J=0,1,2,35.3.3DIT-FFT的運算規(guī)律及編程思想61對N=2M的一般情況，第L級的旋轉(zhuǎn)因子為 J=0,1,2,…,2L-1-1由于2L=2M×2L-M=N·2L-M所以 J=0,1,2,…,2L-1-1，(5.3.7)(5.3.8)這樣，就可按(5.3.7)式和(5.3.8)式確定第L級運算的旋轉(zhuǎn)因子(實際編程序時，L為最外層循環(huán)變量)。對N=2M的一般情況，第L級的旋轉(zhuǎn)因子為62 3.蝶形運算規(guī)律 設(shè)序列x(n)經(jīng)時域抽選(倒序)后，存入數(shù)組X中。如果蝶形運算的兩個輸入數(shù)據(jù)相距B個點，應(yīng)用原位計算，則蝶形運算可表示成如下形式：XL(J)XL-1(J)+XL-1(J+B)WpN

XL(J+B)XL-1(J)-XL-1(J+B)WpN 式中

p=J·2M-L;J=0,1,…,2L-1-1;L=1,2,…,M 3.蝶形運算規(guī)律63 設(shè)

式中，下標(biāo)R表示取實部，I表示取虛部， 設(shè)64 4.編程思想及程序框圖 仔細觀察圖5.3.3，還可歸納出一些編程序有用的運算規(guī)律：第L級中，每個蝶形的兩個輸入數(shù)據(jù)相距B=2L-1個點；同一旋轉(zhuǎn)因子對應(yīng)著間隔為2L點的2M-L個蝶形。 總結(jié)上述運算規(guī)律，便可采用下述運算方法。先從輸入端(第1級)開始，逐級進行，共進行M級運算。當(dāng)進行第L級運算時，依次求出2L-1個不同的旋轉(zhuǎn)因子，每求出一個旋轉(zhuǎn)因子，就計算完它對應(yīng)的所有2M-L個蝶形。這樣，我們可用三重循環(huán)程序?qū)崿F(xiàn)DIT-FFT運算，程序框圖如圖5.3.5所示。 4.編程思想及程序框圖65圖5.3.5DIT-FFT運算和程序框圖圖5.3.5DIT-FFT運算和程序框圖66 5.序列的倒序 DIT-FFT算法的輸入序列的排序看起來似乎很亂，但仔細分析就會發(fā)現(xiàn)這種倒序是很有規(guī)律的。由于N=2M，因此順序數(shù)可用M位二進制數(shù)(nM-1

nM-2…n1n0)表示。M次偶奇時域抽選過程如圖5.3.6所示。 第一次按最低位n0的0和1將x(n)分解為偶奇兩組，第二次又按次低位n1的0、1值分別對偶奇組分解；依次類推，第M次按nM-1位分解，最后所得二進制倒序數(shù)如圖5.3.6所示。表5.3.1列出了N=8時以二進制數(shù)表示的順序數(shù)和倒序數(shù)。 5.序列的倒序67圖5.3.6形成例序的樹狀圖(N=23)圖5.3.6形成例序的樹狀圖(N=23)68表5.3.1順序和倒序二進制數(shù)對照表

表5.3.1順序和倒序二進制數(shù)對照表69 形成倒序J后，將原存儲器中存放的輸入序列重新按倒序排列。設(shè)原輸入序列x(n)先按自然順序存入數(shù)組A中。例如，對N=8,A(0)，A(1)，A(2)，…，A(7)中依次存放著x(0),x(1),…,x(7)。對x(n)的重新排序(倒序)規(guī)律如圖5.3.7所示。倒序的程序框圖如圖5.3.8所示，圖5.3.8中的虛線框內(nèi)是完成計算倒序值的運算流程圖。 形成倒序J后，將原存儲器中存放的輸入序列重新按倒序排列。70圖5.3.7倒序規(guī)律圖5.3.7倒序規(guī)律71圖5.3.8倒序程序框圖圖5.3.8倒序程序框圖72 5.3.4DIF-FFT 在基2快速算法中，頻域抽取法FFT也是一種常用的快速算法，簡稱DIF-FFT。 設(shè)序列x(n)長度為N=2M,首先將x(n)前后對半分開，得到兩個子序列，其DFT可表示為如下形式： 5.3.4DIF-FFT73 式中

WkN/2N=(-1)k=

將X(k)分解成偶數(shù)組與奇數(shù)組，當(dāng)k取偶數(shù)(k=2r,r=0,1,…,N/2-1)時，1