1、2第第4章章 快速傅里葉變換快速傅里葉變換4.1 引言引言4.2 直接計算直接計算DFT的問題及改進的途徑的問題及改進的途徑4.3 按時間抽取（按時間抽取（DIT）的基）的基-2FFT算法算法4.4 按頻率抽取（按頻率抽取（DIF）的基）的基-2FFT算法算法4.5 IDFT的高效算法的高效算法4.6 實序列的實序列的FFT算法算法 4.7 線性調頻線性調頻Z變換（變換（CZT）34.1 引言引言 有限長序列通過離散傅里葉變換（有限長序列通過離散傅里葉變換（DFT）將其頻域離）將其頻域離散化成有限長序列，但其計算量太大，很難實時處理，因此散化成有限長序列，但其計算量太大，很難實時處理，因此引出

2、了快速傅里葉變換（引出了快速傅里葉變換（FFT）。）。 FFT并不是一種新的變換形式，它只是并不是一種新的變換形式，它只是DFT的一種的一種快快速算法速算法，并且根據對序列分解與選取方法的不同產生了多種，并且根據對序列分解與選取方法的不同產生了多種算法。算法。 FFT在離散傅里葉反變換、線性卷積和線性相關等方在離散傅里葉反變換、線性卷積和線性相關等方面也有重要應用。面也有重要應用。4 44.2.1 直接計算直接計算DFT的運算量問題的運算量問題 若若x(n)是是N點序列，實現點序列，實現x(n)的的DFT，即求出，即求出N個個X(k)，需要需要N2次復數乘法，次復數乘法，N(N-1)次復數加法

3、。次復數加法。 10)()()(NnknNDFTWnxkXnx1, 1 , 0 Nk 10)(1)()(NkknNIDFTWkXNnxkX1, 1 , 0 Nn5 5一個一個復數乘法復數乘法包括包括4個實數乘法和個實數乘法和2個實數加法。個實數加法。一次復數一次復數加法需二次實數加法。加法需二次實數加法。例如：例如： x(n) N=1024，N2=1048576)Re)(ImIm)(ReIm)(ImRe)(ReIm)Re(Im)(Re)()(101010nkNnkNnkNNnnkNNnnkNnkNNnnkNWnxWnxjWnxWnxWjWnxjnxWnxkX 因而每運算一個因而每運算一個X(k

4、)需需4N次實數乘法和次實數乘法和2N+2(N-1)=2(2N-1)次實數加法。次實數加法。 所以，整個所以，整個DFT運算總共需要運算總共需要4N2次實數乘法和次實數乘法和2N(2N-1)次實數加法次實數加法。 (a+jb)(c+jd)=(ac-bd)+j(bc+ad)6 解決耗時的乘法問題是將數字信號處理理論用于實際的解決耗時的乘法問題是將數字信號處理理論用于實際的關鍵問題。關鍵問題。特別是特別是40年前，計算機的速度相當慢。因此，很年前，計算機的速度相當慢。因此，很多學者對解決多學者對解決DFT的快速計算問題產生了極大的興趣。的快速計算問題產生了極大的興趣。 Cooley J W, Tu

5、key J W. An algorithm for the machine computation of complex Fourier series. Mathematics of Computation, 1965, pp29730174.2.2 改善途徑改善途徑 FFT的思路：的思路： 10)()(NnknNWnxkX1, 1 , 0 Nk 如何充如何充分利用這些分利用這些關系關系1. 10 W1, 1. 2 mNNWWrrNWW . 31. 42/ NWrrNWW 2/. 5對稱性對稱性周期性周期性8以四點以四點DFT為例為例11111111(0)11(1)(2)1111(3)11xW

6、WxxxWW0000012302460369(0)(0)(1)(1)(2)(2)(3)(3)WWWWXxXWWWWxXxWWWWXxWWWW 304)()(nknWnxkX3 , 2 , 1 , 0 k11111111(0)11(2)(1)1111(3)11xWWxxxWW911(0) (0)(2) (1)(3)(1) (0)(2) (1)(3)(2) (0)(2) (1)(3)(3) (0)(2) (1)(3)XxxxxXxxxxWXxxxxXxxxxW11(0)(2)xx(1)(3)xx(0)(1)XX(2)(3)XX111W(0)(2)xx(0)(2)xx(1)(3)xx(1)(3)xx

7、10 利用利用WNnk的對稱性和周期性，將大點數的的對稱性和周期性，將大點數的DFT分解成若分解成若干個小點數的干個小點數的DFT，FFT正是基于這個基本思路發展起來的。正是基于這個基本思路發展起來的。 分類：按時間抽取（分類：按時間抽取（DIT）算法和按頻率抽取（）算法和按頻率抽取（DIF）算法。算法。 問題是問題是如何分最有效？如何分最有效？ N點DFTN/2點DFTN/4點 DFT2點 DFT 1個 2個 4個 N/2個Decimation in Time Decimation in Frequency 11knNW的周期性周期性knNW的對稱性對稱性nkNnkNWW *)(nkNNnk

8、NWW )2/()()(NknNkNnNnkNWWW 124.3 按時間抽取（按時間抽取（DIT）的基）的基-2FFT算法算法 1.算法原理算法原理 設序列設序列x(n)長度為長度為N2M，M為整數，將為整數，將x(n)（n=0,1N-1）按）按n的奇偶分解成兩組的奇偶分解成兩組N/2點的子序列點的子序列 x1(r)= x(2r) x2(r)= x(2r+1) r =0,1N/2-1（長度為（長度為N/2）則則 10)()(NnnkNWnxkX 為為奇奇數數為為偶偶數數nnkNnnkNWnxWnx)()( 12/0)12(12/02)12()2(NrkrNNrrkNWrxWrx 12/02/1

9、2/02/)12()2(NrrkNkNNrrkNWrxWWrx13 12/02/12/02/)12()2()(NrrkNkNNrrkNWrxWWrxkX)(1kX)(2kX12/,.,1 , 0)()()2/(21 NkkXWkXNkXkN12/,.,1 , 0)()()(21 NkkXWkXkXkN X1(k)、X2(k)都是都是 N/2 點的點的 DFT，它們各自又可分成，它們各自又可分成 N/4 點的點的DFT，如此繼續分下去，直至兩點，如此繼續分下去，直至兩點DFT。兩點。兩點DFT不需要乘法運算：不需要乘法運算：)1()0()1()1()0()0(xxXxxX 2點點DFT的運算流圖

10、的運算流圖15 由于每一步分解都是按輸入序列在時間上的奇偶次序，由于每一步分解都是按輸入序列在時間上的奇偶次序，分解成兩個半長的子序列，所以稱分解成兩個半長的子序列，所以稱“按時間抽取算法按時間抽取算法”。162.運算量比較運算量比較 對對N=2M，共可分，共可分M次次，即，即m =0,1M-1，每一級有，每一級有N/2個如下的蝶形單元個如下的蝶形單元 一個蝶形單元只需一個蝶形單元只需一一次復數乘法，次復數乘法，兩兩次復數加法。次復數加法。 總共復數乘法：總共復數乘法： 復數加法：復數加法： )()()()()()(1111jXWkXjXjXWkXkXmrNmmmrNmm NNMN2log22

11、 NNMN2log 可以共享可以共享17直接計算直接計算DFT與與FFT算法的計算量之比為算法的計算量之比為 NNNMNN222log22 183.算法特點算法特點1）原位運算原位運算 兩點構成一個蝶形單元，并且這兩點不再參與別的蝶形兩點構成一個蝶形單元，并且這兩點不再參與別的蝶形單元的運算。單元的運算。 所謂原位運算，就是當數據輸入到存儲器后，每一級運所謂原位運算，就是當數據輸入到存儲器后，每一級運算的結果仍然貯存在這同一組存儲器中，直到最后輸出，中算的結果仍然貯存在這同一組存儲器中，直到最后輸出，中間無需其它存儲器。間無需其它存儲器。192）旋轉因子的變化規律）旋轉因子的變化規律 如上所述

12、，如上所述，N點點DITFFT運算流圖中，每級都有運算流圖中，每級都有N/2個蝶個蝶形。每個蝶形都要乘以因子形。每個蝶形都要乘以因子 ，稱其為旋轉因子，稱其為旋轉因子，p稱為旋稱為旋轉因子的指數。轉因子的指數。 不難發現，不難發現，第第L級級共有共有2L-1個不同的旋轉因子。個不同的旋轉因子。N=23=8時的各級旋轉因子表示如下：時的各級旋轉因子表示如下： pNW3,2,1,0,31,0,20,1222/24/ JWWWLJWWWLJWWWLJJNpNJJNpNJJNpNLLL一般情況，第一般情況，第L級的旋轉因子指數為級的旋轉因子指數為L為從左到右的運算級數，編程時為從左到右的運算級數，編程

13、時L為循環變量。為循環變量。LMJp 2MLJL,.,2 , 1, 12,.,2 , 1 , 01 20 如果蝶形運算的兩個輸入數據相距如果蝶形運算的兩個輸入數據相距B個點，應用原位運個點，應用原位運算，則算，則pNLLLpNLLLWBJXJXBJXWBJXJXJX)()()()()()(1111 式中式中LMJp2MLJL,.,2 , 1, 12,.,2 , 1 , 01B 兩蝶形節點的兩蝶形節點的“距離距離”（蝶距）（蝶距）12LB 3）蝶形運算）蝶形運算21 4）倒位序）倒位序 輸入序列輸入序列x(n)的排列次序不符合自然順序。此現象是由的排列次序不符合自然順序。此現象是由于按于按n的奇

14、偶分組進行的奇偶分組進行DFT運算而造成的，這種排列方式稱運算而造成的，這種排列方式稱為為“倒位序倒位序”。 所謂所謂“倒位序倒位序”，是指按二進制表示的數字首尾位置顛，是指按二進制表示的數字首尾位置顛倒，重新按十進制讀數。倒，重新按十進制讀數。自然順序二進制順序碼位倒置碼位倒讀順序0 000000001 100110042 201001023 301111064 410000115 510110156 611001137 7111111722倒位序的實現倒位序的實現 在實際運算時，先按自然順序將輸入序列存入存儲單元，在實際運算時，先按自然順序將輸入序列存入存儲單元，再通過再通過變址運算變址運

15、算將自然順序變換成按將自然順序變換成按DIT-FFT算法要求的順算法要求的順序。序。 變址的過程可以用程序安排加以實現變址的過程可以用程序安排加以實現- 稱為稱為“整序整序”或或“重排重排”（采用碼位倒讀）（采用碼位倒讀）注意注意:(1)當當I=J時，數據不必調換；時，數據不必調換；(2)當當IJ時，必須將原來存放數據時，必須將原來存放數據x(I)送入暫存器送入暫存器R，再將，再將x(J)送入送入x(I)，R中內容送中內容送x(J)，進行數據對調。進行數據對調。(3)為了避免再次調換已調換過的數據，保證調換只進行一次，為了避免再次調換已調換過的數據，保證調換只進行一次，否則又變回原狀否則又變回

16、原狀。JI（或（或I1, 向內盤旋（內縮）；向內盤旋（內縮）； W01, 向外盤旋（外伸）向外盤旋（外伸） 表示相鄰分析點間的夾角表示相鄰分析點間的夾角。0 0)(00000000 kjkjkkjkkeWAeWeAAWz 。變變換換求求出出該該序序列列即即由由，為為均均勻勻分分布布在在單單位位園園上上此此時時等等分分）時時，（。即即即即（當當滿滿足足下下面面特特殊殊條條件件：DFTCZTzNNMWeeWWcAeAAbNMakNjjj 221,)(01,1)()00200000045 10)()(NnnknkWAnxzX1,.,1 , 0 Mk， )(00000000 kjkjkkjkkeWAe

17、WeAAWz 10)()(NnnznxzX為了提高運算速度為了提高運算速度 )(21222nkknnk 1022)(2222)()(NnknknnkWWWAnxzX 102)(22222)(NnnknnkWWAnxW46 102)(22222)()(NnnknnkkWWAnxWzX22)()(nnWAnxng 1,.,1 , 0Nn22)(nWnh 102)()()(2NnkknkhngWzX)(*)(22khkgWk 1,.,1 , 0 Mk， 47472.CZT計算框計算框圖圖)(kzX)(nx22nnWA 22kW22)(nWnh )(ng)()(kVnV 102)()()(2Nnkkn

18、khngWzX其中：卷積可以通過圓周卷積來實現，這樣可其中：卷積可以通過圓周卷積來實現，這樣可借用借用FFT快速算法。快速算法。 計算步驟計算步驟 自學自學22)()(nnWAnxng 22)(nWnh )(nh48說明說明變變換換）（線線性性調調頻頻此此變變換換稱稱為為線線性性調調頻頻信信號號（號號為為在在雷雷達達系系統統中中稱稱這這種種信信列列成成線線性性增增長長的的復復指指數數序序隨隨時時間間可可以以想想象象為為頻頻率率序序列列ZCZTSignalChirpnWeeWnhnWnhnjnjnn)21, 1()()()(00022202022 49用函數用函數chirp產生線性調頻信號，并畫出波形。產生線性調頻信號，并畫出波形。t=0:0.001:2; % 2 secs 1kHz sample ratey=chirp(t,0,1,150); % Start DC, cross 150Hz at t=1s plot(t,y),axis(0 0.5 -1 1) 00.050.10.1