第四章-快速傅里葉變換(FFT)

1、第四章第四章快速傅里葉變換快速傅里葉變換（FFT）2021/8/61主要內容主要內容qDIT-FFT算法算法 qDIF-FFT算法算法qIFFT算法算法qChirp-FFT算法算法q線性卷積的線性卷積的FFT算法算法4.1 引言引言q FFT: Fast Fourier Transformq 1965年，年，Cooley-Turky 發(fā)表文章發(fā)表文章機器計算傅機器計算傅里葉級數的一種算法里葉級數的一種算法，提，提出出FFT算法，解決算法，解決DFT運算量太大，在實際使用中受限制的問題。運算量太大，在實際使用中受限制的問題。q FFT的應用。頻譜分析、濾波器實現、實時信的應用。頻譜分析、濾波器實

2、現、實時信號處理等。號處理等。q DSP芯片實現。芯片實現。TI公司的公司的TMS 320c30，10MHz時鐘，基時鐘，基2-FFT1024點點FFT時間時間15ms。q 典型應用：信號頻譜計算、系統分析等典型應用：信號頻譜計算、系統分析等)()(kXnxDFT )()()(nynhnxFFTnhnyIFFTFFTnx)()()( 系統分析系統分析 頻譜分析與功率譜計算頻譜分析與功率譜計算4.2 直接計算直接計算DFT的問題及改進途徑的問題及改進途徑10)()(NnknNWnxkX10)(1)(NkknNWkXNnx1、 DFT與與IDFT( )Nx n點有限長序列2、DFT與與IDFT運算

3、特點運算特點復數乘法復數乘法復數加法復數加法一個一個X(k)NN 1N個個X(k)(N點點DFT)N 2N (N 1)10( )NnkNnx n Wajbcjdacbdj adcb同理：同理：IDFT運算量與運算量與DFT相同。相同。實數乘法實數乘法實數加法實數加法一次復乘一次復乘42一次復加一次復加2一個一個X (k) 4N2N+2 (N 1)=2 (2N 1)N個個X (k)(N點點DFT)4N 22N (2N 1)3、降低、降低DFT運算量的考慮運算量的考慮nkNW 的特性*()() ()nknkNn kn NkNNNNWWWW對稱性()() nkN n kn N kNNNWWW周期性

4、nkmnkNmNWW可約性/nknk mNN mWW0/2(/2) 11Nk NkNNNNWWWW 特殊點：2jnknkNNWeNknkNNWWnNnkNNWW2jmnkmNe221NjjNee FFT算法分類算法分類:q 時間抽選法時間抽選法DIT: Decimation-In-Timeq 頻率抽選法頻率抽選法DIF: Decimation-In-FrequencyFFTDFTDFTDFTDFT算法的基本思想： 利用系數的特性，合并運算中的某些項， 把長序列短序列，從而減少其運算量。4.3 按時間抽?。ò磿r間抽?。―IT）的）的FFT算法算法12/.210) 12()()2()(21Nrrx

5、rxrxrx，（Decimation In Time）1、算法原理、算法原理設序列點數設序列點數 N = 2L，L 為整數。為整數。 若不滿足，則補零若不滿足，則補零將序列將序列x(n)按按n的奇偶分成兩組：的奇偶分成兩組：N為為2的整數冪的的整數冪的FFT算法稱算法稱基基-2FFT算法算法。將將N點點DFT定義式分解為兩個長度為定義式分解為兩個長度為N/2的的DFT10)()()(NnknNWnxnxDFTkXkrNnNrrkNnNrWrxWrx)12(12/0212/0) 12()2( 為奇為偶 )(12/02/2)(2/12/0121)()(kXNrrkNkNkXrkNNrWrxWWrx

6、)()()(21kXWkXkXkN記：記： （1 1）rkNrkNWW2/2（這一步利用：（這一步利用： ）,0,1,./2 1r kN再利用周期性求再利用周期性求X(k)的后半部分的后半部分/22NkNkkNNNNWWWW 又)(2)()()(222112/02/112/0)2/(2/11kXkNXkXWrxWrxkNXNrrkNNrkNrNrkNkNrNWW2/)2/(2/)2()2()2()2(12/,.2 , 1 , 0)()()(2)2/(121kNXWkNXkNXNkkXWkXkXkNNkN，12/,.2 , 1 , 0)()(21NkkXWkXkN，將上式表達的運算用一個專用將上

7、式表達的運算用一個專用“蝶形蝶形”信流圖表示。信流圖表示。)(1kX)(2kX)()(21kXWkXkN)()(21kXWkXkNkNW1212( )( )( )()( )( )2kNkNX kXkW XkNX kXkW Xk0,1,.,/2 1kN注：注：a. 上支路為加法，下支路為減法；上支路為加法，下支路為減法； b. 乘法運算的支路標箭頭和系數。乘法運算的支路標箭頭和系數。用用“蝶形結蝶形結”表示上面運算的分解：表示上面運算的分解： 328N)0(x)1 (x)2(x)3(x)4(x)5(x)6(x)7(x)0(X) 1 (X)2(X)3(X)4(X)5(X)6(X)7(X1NW0NW

8、2NW3NW)0(1X)1 (1X)2(1X)3(1X)0(2X)1 (2X(3)2X)2(2XDFTN點2DFTN點2分解后的運算量：分解后的運算量：復數乘法復數乘法復數加法復數加法一個一個N/2點點DFT(N/2)2N/2 (N/2 1)兩個兩個N/2點點DFTN2/2N (N/2 1)一個蝶形一個蝶形12N/2個蝶形個蝶形N/2N總計總計22/2/2/2NNN2/2 1/2N NNN運算量減少了近一半運算量減少了近一半進一步分解進一步分解MN2122MN2N4N由于由于 ， 仍為偶數，因此，兩個仍為偶數，因此，兩個 點點DFTDFT又可同樣進一步分解為又可同樣進一步分解為4 4個個 點的

9、點的DFTDFT。1314(2 )( )(21)( )xlx lxlx l0,1,.,/4 1lN13/2413/24( )( )( )()( )( )4kNkNX kXkWXkNX kXkWXk0,1,.,14Nk 02/NW12/NW)(3lx)(4lx)2(x)4(x)6(x)0(x)0(1X) 1 (1X)2(1X) 3(1X) 0(3X) 1 (3X)0(4X) 1 (4XDFTN點4DFTN點4“蝶形蝶形”信流圖表示信流圖表示 N點點DFT分解為四個分解為四個N/4點的點的DFTDFTN點4DFTN點4DFTN點4DFTN點4)2(x)4( x)6( x)0( x) 1 ( x)

10、3 ( x)5(x)7( x0NW2NW0NW2NW1NW0NW2NW3NW)0(X) 1 (X) 2(X) 3(X) 4(X) 5(X)6(X)7(X)(.kX)(.nxq 類似的分解一直繼續(xù)下去，直到分解為最類似的分解一直繼續(xù)下去，直到分解為最后的兩類蝶形運算為止后的兩類蝶形運算為止(2點點DFT).q 如上述如上述N=8=23，N/4=2點中：點中： 類似進一步分解類似進一步分解1點點DFTx(0)1點點DFTx(4)X3(0)X3(1)02W進一步簡化為蝶形流圖：進一步簡化為蝶形流圖：0NWX3(0)X3(1)x(0)x(4)4()0()4()0()0(004/3xWxxWxXNN)4

11、()0()4()0() 1 (014/3xWxxWxXNN因此因此8 8點點FFTFFT時間抽取方法的信流圖如下時間抽取方法的信流圖如下)2(x)4(x)6(x)0(x) 1 ( x) 3 ( x)5(x)7(x0NW0NW0NW0NW第一級.0NW2NW0NW2NW 第二級.)(0kX1m)(1kX)(2kX)(3kX2m3m1NW0NW2NW3NW)0(X) 1 (X)2(X)3(X)4(X)5(X)6(X)7(X 第三級.FFT運算量與運算特點運算量與運算特點 1 N=2L時，共有時，共有L=log2N級運算；每一級有級運算；每一級有N/2個蝶形結。個蝶形結。2每一級有每一級有N個數據（

12、中間數據），且每級只用到個數據（中間數據），且每級只用到本級的轉入中間數據，適合于迭代運算。本級的轉入中間數據，適合于迭代運算。3計算量：計算量： 每級每級N/2次復乘法，次復乘法，N次復加。（每蝶形只乘一次，次復加。（每蝶形只乘一次，加減各一次）。共有加減各一次）。共有L*N/2=N/2log2N 次復乘法；次復乘法；復加法復加法L*N=Nlog2N 次。與直接次。與直接DFT定義式運算定義式運算量相比量相比(倍數倍數) N2/(Nlog2N) 。當。當 N大時，此倍數大時，此倍數很大。很大。222()2()loglog2FFmDFTNNNmFFTNN比較比較DFT 參考參考P150 表表4

13、-1 圖圖4-6可以直觀看出，當點數可以直觀看出，當點數N越大時，越大時，FFT的優(yōu)點更突出。的優(yōu)點更突出。按時間抽取按時間抽取FFT蝶形運算特點蝶形運算特點 1、關于、關于FFT運算的混序與順序處理（位倒序處理）運算的混序與順序處理（位倒序處理） 由于輸入序列按時間序位的奇偶抽取，故輸入序由于輸入序列按時間序位的奇偶抽取，故輸入序列是混序的，為此需要先進行混序處理。列是混序的，為此需要先進行混序處理。混序規(guī)律：混序規(guī)律： x(n)按按n位置進行碼位（二進制）倒置位置進行碼位（二進制）倒置規(guī)律輸入，而非自然排序，即得到混序排列。所規(guī)律輸入，而非自然排序，即得到混序排列。所以稱為位倒序處理。以稱

14、為位倒序處理。位倒序實現：位倒序實現：（1）DSP實現采用位倒序尋址實現采用位倒序尋址（2）通用計算機實現可以有兩個方法：一是嚴格按）通用計算機實現可以有兩個方法：一是嚴格按照位倒序含義進行；二是倒進位的加照位倒序含義進行；二是倒進位的加N/2。倒位序倒位序自然序自然序00000000100410010102201011063011001141001015510101136110111771112 102( )()x nnn n n倒位序倒位序例例 計算計算 ， 。計算。計算 點點FFTFFT。用時間抽取輸入倒序算法，。用時間抽取輸入倒序算法，問倒序前寄存器的數問倒序前寄存器的數 和倒序后和倒

15、序后 的的數據值？數據值？2)(nnx31.2 , 1 , 0n32N)13( x)13( x16913)13(2x5232N210)01101()13(102)22()10110(48422)13(2x解：倒序前解：倒序前 倒序倒序 倒序為倒序為 倒序后倒序后 DIT FFT中最主要的蝶形運算實現中最主要的蝶形運算實現（1）參與蝶形運算的兩類結點（信號）間）參與蝶形運算的兩類結點（信號）間“距離距離”（碼地址）與其所處的第幾級蝶形有關；第（碼地址）與其所處的第幾級蝶形有關；第m級的級的“結距離結距離”為為 (即原位計算迭代）即原位計算迭代）（2）每級迭形結構為）每級迭形結構為，12m)2(,

16、.2 , 1LNLmrNmmmmmrNmmmmWkXkXkXWkXkXkX)2()()2()2()()(1111111q 蝶形運算兩節(jié)點的第一個節(jié)點為蝶形運算兩節(jié)點的第一個節(jié)點為k值，表值，表示成示成L位二進制數，左移位二進制數，左移L m位，把右位，把右邊空出的位置補零，結果為邊空出的位置補零，結果為r的二進制數。的二進制數。2( )2L mrk（3） 的確定：的確定： 第第m級的級的r取值：取值：rNWkNWDIT算法的其他形式流圖算法的其他形式流圖q 輸入倒位序輸出自然序輸入倒位序輸出自然序q 輸入自然序輸出倒位序輸入自然序輸出倒位序q 輸入輸出均自然序輸入輸出均自然序q 相同幾何形狀相

17、同幾何形狀q 輸入倒位序輸出自然序輸入倒位序輸出自然序q 輸入自然序輸出倒位序輸入自然序輸出倒位序參考參考P154-155時間抽取、時間抽取、 輸入自然順序、輸入自然順序、 輸出倒位序的輸出倒位序的FFTFFT流圖流圖 0NW0NW0NW2NW1110NW10NW2NW1111X(0)x(0)X(4)x(1)X(2)x(2)X(6)x(3)X(1)x(4)X(5)x(5)X(3)x(6)X(7)x(7)110NW0NW2NW1NW3NW11 例例 用用FFT算法處理一幅算法處理一幅NN點的二維圖像，如用每秒可點的二維圖像，如用每秒可做做10萬次復數乘法的計算機，當萬次復數乘法的計算機，當N=1

18、024時，問需要多少時間時，問需要多少時間（不考慮加法運算時間）？（不考慮加法運算時間）？ 解解 當當N=1024點時，點時，FFT算法處理一幅二維圖像所需復數乘算法處理一幅二維圖像所需復數乘法約為法約為 次，僅為直接計算次，僅為直接計算DFT所需時間的所需時間的10萬萬分之一。分之一。 即原需要即原需要3000小時，現在只需要小時，現在只需要2 分鐘。分鐘。 722210log2NN4.4 按頻率抽?。ò搭l率抽取（DIF）的）的FFT算法算法q 與與DIT-FFT算法類似分解，但是抽取的是算法類似分解，但是抽取的是X(k)。即分解即分解X(k)成奇數與偶數序號的兩個序列。成奇數與偶數序號的兩

19、個序列。q 設：設： N = 2L，L 為整數。將為整數。將X(k)按按k的奇偶分的奇偶分組前，先將輸入組前，先將輸入x(n)按按n的順序分成前后兩半：的順序分成前后兩半：（Decimation In Frequency)一、算法原理一、算法原理12/12/0)()()(NNnnkNNnnkNWnxWnxkX12/0)(212/02)()(NnknNNNnnkNNWnxWnx12/022/)()(NnnkNNkNNWnxWnxkNkNW) 1(2/2 10( )( 1)2NknkNnNx nx nW 0,1,.,1kN下面討論下面討論：的)(12,2kXrrk12/02/212/022) 1

20、()()()()()2(NnnrNNNnrnNNWnxnxWnxnxrX12/02/212/0)12(2)2()()()()() 12(NnnrNnNNNnnrNNWWnxnxWnxnxrX按按k k的奇偶將的奇偶將X(k)X(k)分成兩部分：分成兩部分：顯然：顯然：。點的對應兩個DFTNrXrX2/) 12(),2(nNNNWnxnxnxnxnxnx)()()()()()(2221)()(2NnxnxnNNNWnxnxnxnx)()()()(22nNW令：令：用蝶型結構圖表示為：用蝶型結構圖表示為：x1(0)x1(1)-1x1(2)x1(3)-1x2(0)x2(1)-1x2(2)x2(3)-

21、1N/2點DFTN/2點DFTx(0)x(7)x(1)x(2)x(3)x(4)x(5)x(6)X1(0)=X(0)X2(0)=X(1)X1(1)=X(2)X1(2)=X(4)X1(3)=X(6)X2(1)=X(3)X2(2)=X(5)X2(3)=X(7)1NW0NW2NW3NW311411/2( )( )(/4)( ) ( )(/4)nNx nx nx nNx nx nx nNW0,1,.,14Nn 313414( )(2 ) ( )( )(21)( )X kXkDFT x nX kXkDFT x n0,1,.,14Nk N/2仍為偶數，進一步分解：仍為偶數，進一步分解：N/2 N/4x3(0

22、)x3(1)-1-1x4(0)x4(1)N/4點DFTN/4點DFTx1(0)x1(1)x1(2)x1(3)X3(0)=X1(0)=X(0)X4(0)=X1(1)=X(2)X3(1)=X1(2)=X(4)X4(1)=X1(3)=X(6)0/2NW1/2NWq 按照以上思路繼續(xù)分解，即一個按照以上思路繼續(xù)分解，即一個N/2的的DFT分解成兩個分解成兩個N/4點點DFT，直到只計算，直到只計算2點的點的DFT，這就是，這就是DIF-FFT算法。算法。2個個1點的點的DFT蝶形流圖蝶形流圖 進一步簡化為蝶形流圖：進一步簡化為蝶形流圖：1點點DFTx3(0)1點點DFTx3(1)X(0)X(4)02W

23、02WX(0)X(4)x3(0)x3(1)2(x)4(x)6(x)0(x)1(x)3(x)5(x)7(x)0(X)1(X)2(X)3(X)4(X)5(X)6(X)7(X0NW0NW1NW2NW3NW0NW0NW0NW2NW0NW2NW0NW1m第一級： 2m第二級：3m第三級：二、按頻率抽取二、按頻率抽取FFT蝶形運算特點蝶形運算特點1）原位計算）原位計算1111( )( )( )( )( )( )mmmrmmmNXkXkXjXjXkXj WrNW1( )mXk1( )mXj( )mXk( )mXj-1L級蝶形運算，每級級蝶形運算，每級N/2個蝶形，每個蝶形結構：個蝶形，每個蝶形結構：m表示第

24、表示第m級迭代，級迭代，k，j表示數據所在的行數表示數據所在的行數2）蝶形運算）蝶形運算對對N=2L點點FFT，輸入自然序，輸出倒位序，輸入自然序，輸出倒位序，兩節(jié)點距離：兩節(jié)點距離：2L-m=N / 2m1111( )( )2( )22mmmmrmmmNmmNXkXkXkNNXkXkXkW第第m級運算：級運算：q 蝶形運算兩節(jié)點的第一個節(jié)點為蝶形運算兩節(jié)點的第一個節(jié)點為k值，表示值，表示成成L位二進制數，左移位二進制數，左移m-1位，把右邊空出的位，把右邊空出的位置補零，結果為位置補零，結果為r的二進制數。的二進制數。rNW 的確定12( )2mrk存儲單元存儲單元輸入序列輸入序列x(n)

25、: N個存儲單元個存儲單元rNW系數系數 ：N / 2個存儲單元個存儲單元三、三、DIT與與DIF的異同的異同q 基本蝶形不同基本蝶形不同2log2FNmN2logFaNNq DIT: 先復乘后加減先復乘后加減q DIF: 先減后復乘先減后復乘q 運算量相同運算量相同q 都可原位運算都可原位運算q DIT和和DIF的基本蝶形互為轉置的基本蝶形互為轉置4.5 IDFT的的FFT算法算法（FFT應用一）應用一） 一、從定義比較分析一、從定義比較分析knNNkWkXNkXIDFTnx10)(1)()(10)()()(NnknNWnxnxDFTkX與與DFT的比較：的比較： 1）、旋轉因子）、旋轉因子

26、WN-kn 的不同；的不同； 2）、結果還要乘）、結果還要乘 1/N。 )(10*10*)(1)(1)()(kXDFTknNNkknNNkWkXNWkXNkXIDFTnx二、實現算法二、實現算法直接使用直接使用FFT程序的算法程序的算法*)(1)(kXDFTNnx共軛共軛FFT共軛共軛乘乘1/ N( )X k*( )Xk( )x n直接調用直接調用FFT子程序計算子程序計算IFFT的方法：的方法：4.6 線性調頻線性調頻Z變換（變換（Chirp-Z變換）算變換）算法法 （FFT應用二）應用二）單位圓與非單位圓采樣單位圓與非單位圓采樣（a） 沿單位圓采樣沿單位圓采樣; (b) 沿沿AB弧采樣弧采

27、樣 ooooABX(ej)RezRezjImzjImzAB(a)(b)X(ej) 螺線采樣螺線采樣 0jImzRezo0A01.0zM1A0W01z0(M I)0zk=AW-k k=0, 1, , M-1 0000jjeWWeAAChirp-Z變換的線性系統表示變換的線性系統表示 x(n)2/ 2 )(nWnh2/ 2nnWAg(n)h(n)1 / h(n)X(zn)由于系統的單位脈沖響應由于系統的單位脈沖響應 可以想象為頻率可以想象為頻率隨時間隨時間(n)(n)呈線性增長的復指數序列。在雷達系統中呈線性增長的復指數序列。在雷達系統中, ,這這種信號稱為種信號稱為線性調頻信號（線性調頻信號（C

28、hirp SignalChirp Signal），因此，這，因此，這里的變換稱為里的變換稱為線性調頻線性調頻Z Z變換變換。22)(nWnh一、基本算法思路一、基本算法思路10)()()()()(MmmnxmhnhnxnyLMMd)1()(nMhnh2/LMmd4.7 線性卷積的線性卷積的FFT算法算法（FFT應用三）應用三）若若L點點x(n)，M點點h(n)，則直接計算其線性卷積則直接計算其線性卷積y(n)需運算量：需運算量：若系統滿足線性相位，即：若系統滿足線性相位，即：則需運算量：則需運算量：FFT法：以圓周卷積代替線性卷積法：以圓周卷積代替線性卷積21mNML令 ( )01( )01x

29、 nnLx nLnN( )01( )01h nnMh nMnNN( )( )* ( )( ) ( )y nx nh nx nh n則 NN2log2NN2log2NN2log2)()() 1nhFFTkH)()()2nxFFTkX)()()()3kXkHkY)()()4kYIFFTnyN總運算量：總運算量： 次乘法次乘法NNNmF2log23比較直接計算和比較直接計算和FFT法計算的運算量法計算的運算量22(1 3/2*log)dmFmMLKmNN222413/2*(1log)106logmMMKMMM223logmMKL討論：討論：ML12NMLM 則1）當）當1NMLL 則2）當）當mLK

30、 需采用分段卷積重疊相加法重疊保留法ML 見見P183 表表 x(n)長度很長時，將長度很長時，將x(n)分為分為L長的若干長的若干小的片段，小的片段，L與與M可比擬?？杀葦M。nLiniLnxnxi，其它，01) 1()()(iinxnx)()()()()(nhnxnyiinhnx)()(1 1、重疊相加法、重疊相加法iiny)( 則：則： 輸出：輸出：)()()(nhnxnyii1MLN其中：其中：可以用圓周卷積計算：可以用圓周卷積計算：mN2 選選 ，上面圓周卷積可用，上面圓周卷積可用FFTFFT計算。計算。 )()()(nhnxnyiiN 由于由于yi(n)長度為長度為N，而，而xi(n

31、)長度長度L ，必有，必有M-1點重疊，點重疊， yi(n)應相加才能構成最后應相加才能構成最后y(n)的。的。iinyny)()(h(n)0N 1M 1x(n)0L2L3LnnnnnL 10 x0(n)N 10 x1(n)L2L 1LN 13L 10 x2(n)2L2LN 1重疊相加法圖形重疊相加法圖形nnnN 10y0(n)x0(n) h (n)N2L2L N 100L N 1Ly1(n)x1(n) h (n)Ny2(n)x2(n) h (n)N和上面的討論一樣，和上面的討論一樣， 用用FFT法實現重疊相加法的步驟如下法實現重疊相加法的步驟如下: 計算計算N點點FFT， H(k)=DFTh

32、(n)； 計算計算N點點FFT，Xi(k)=DFTxi(n)； 相乘，相乘，Yi(k)=Xi(k)H(k)； 計算計算N點點IFFT，yi(n)=IDFTYi(k)； 將各段將各段yi(n)（包括重疊部分）相加，（包括重疊部分）相加， 。重疊相加的名稱是由于各輸出段的重疊部分相加而得名的。重疊相加的名稱是由于各輸出段的重疊部分相加而得名的。 )()(0nynyii例例 已知序列已知序列xn=n+2,0 n 12, hn=1,2,1試試利用重疊相加法計算線性卷積利用重疊相加法計算線性卷積, 取取L=5 。yn=2, 7, 12, 16, 20, 24, 28, 32, 36, 40, 44, 48, 52, 41, 14解解: 重疊相加法重疊相加法x1n=2, 3, 4, 5, 6x2n=7, 8, 9, 10, 11x3n