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

1、第四章 快速傅里葉變換（FFT） 通信1402班 肖太龍 Contents 從DFT到FFT FFT的基本概念 FFT算法的實現(xiàn)原理 FFT的物理意義 DFT與FFT時間復(fù)雜度比較 從DFT到FFT的緣由 依次類推，N/2個 點的時間復(fù)雜度為 多少呢？N/4呢？ 從哪個方面減少運算量呢？ 周期性表現(xiàn)為： 對稱性表現(xiàn)為： 1965年庫利和圖基提出一 種DFT的快速算法，以此來 解決DFT變換的效率問題。 算法的提出的初衷主要是為 了減少乘法次數(shù)，又因為旋 轉(zhuǎn)因子的對稱性和周期性能 達(dá)到這一目標(biāo)，因此在理論 上是可以實現(xiàn)的。 FFT算法就是不斷地把長序 列的DFT分解為短序列的 DFT的算法。最常

2、用的是基 2FFT算法。 22 ()jm lNjm m lNm NN NN WeeW 2 mN mN mm NNNN N m m NN WWWW WW 或者 淺談FFT實現(xiàn)原理 lFT算法基本上分為兩大類： l時域抽取法FFT(Decimation In Time FFT,簡稱DIT-FFT)。 l頻域抽取法FFT(Decimation In Frequency FFT,簡稱 DIFFFT)。 l我們主要來分析DIF-FFT算法。 A 算法理論推導(dǎo) B 一個簡單的算例 8點DFT一次時域抽取分解 C Matlab繪制圖像 DIF-FFT與DIT-FFT算法有什么異同？ 算法理論推導(dǎo)： 設(shè)序列

3、的長度為N，且滿足N=2M ，M為自然數(shù)。 按n奇偶把 分解為兩個N/2點的子序列 則x(n)的DFT為 ( )x n ( )x n 1 2 ( )(2 ),0,1,1 2 ( )(21),0,1,1 2 N x rxrr N x rxrr /2 1/2 1 2(21) 00 /2 1/2 1 2 12 00 ( )( )( ) (2 )(21) ( )( ) knkn NN nn NN krkr NN rr NN kkr NN rr X kx n Wx n W xr WxrW x rWx r W 奇數(shù)偶數(shù) 因為 所以 2 2 2 2 2 /2 jkr N jkr krkr N NN WeeW

4、 /2 1/2 1 1/22/2 00 12 ( )( )( ) ( )( ) NN krkkr NNN rr k N X kx r WWx r W X kW Xk 0,1,2,1kN 其中X1(k)和X2(k)分別為x1(r)和x2(r)的N/2點DFT，即 /2 1 11/21 0 /2 1 22/22 0 ( )( )( ) ( )( )( ) N kr N r N kr N r X kx r WDFT x r Xkx r WDFT x r 由于X1(k)和X2(k)均以N/2為周期，且 2 N k k NN WW 所以X(k)又可表示為 12 12 ( )( )( )0,1,1 2 (

5、)( )( )0,1,1 22 k N k N N X kX kW Xkk NN X kX kW Xkk 至此，一次時域抽取法的理論推導(dǎo)就完成了。從上面的兩個式子 中，我們定義一種新的運算符蝶形運算符。 C A B A BC A BC N/2點 DFT WN 0 N/2點 DFT WN 1 WN 2 WN 3 x(0) X1(0) x(2) x(4) x(6) x(1) x(3) x(5) x(7) 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) 8點DFT一次時域抽取分解運算流

6、圖 完成一次蝶形運算需要一次復(fù)數(shù)乘法和兩次復(fù)數(shù)加 法運算。因此在8點DFT一次時域抽取分解中，共需要計 算兩個N/2點DFT運算和N/2個蝶形運算。 所以按照上圖的計算DFT時，總的復(fù)數(shù)乘法次數(shù)為 復(fù)數(shù)加法次數(shù)為 2 2 1 (1) 2() 2222 N NNN NN 2 (1)2 222 NNN N 類似地，我們將 按奇偶分解成兩個N/4點子 序列 ，其表達(dá)式分別如下： 1( ) x r 34 ( )( )x lx l和 31 41 ( )(2 ) ,0,1,1 ( )(21)4 x lxl N l x lxl 那么， 又可表示為 1( ) X k /4 1/4 1 2(21) 11/21/

7、2 00 /4 1/4 1 3/4/24/4 00 3/24 ( )(2 )(21) ( )( ) ( )( ),0,1,/ 2 1 NN klkl NN ii NN klkkl NNN ii k N X kxl WxlW x l WWx l W XkWXkkN 式中 同理，由 的周期性和 的對稱性 最后得到： /4 1 33/43 0 /4 1 44/44 0 ( )( )( ) ( )( )( ) N kl N i N kl N i x kx l WDFT x l x kx l WDFT x l 13/24 13/24 ( )( )( ) ,0,1,1 4 (/ 4)( )( ) k N

8、k N X kXkWXk N k X kNXkWXk 34 ( )( )XkXk和 /2 m N W /4 /2/2 k Nk NN WW 用同樣的方法可以計算出 其中 25/26 25/26 ( )( )( ) ,0,1,/41 (/4)( ) k N k N XkXkWXk kN XkNX kWXk /4 1 55/45 0 /4 1 66/46 0 52 62 ( )( )( ) ( )( )( ) ( )(2 ) ,0,1,/41 ( )(21) N kl N i N kl N i Xkx l WDFT x l Xkx l WDFT x l x lxl lN x lxl N/4點 DF

9、T WN 1 2 WN 1 2 WN 0 WN 1 WN 2 WN 3 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) x(0) X3(0) X3(1) X4(0) X4(1) x(4) x(2) x(6) x(1) x(5) x(3) x(7) N/4點 DFT N/4點 DFT N/4點 DFT WN 0 2 WN 0 2 8點DFT二次時域抽取分解運算 WN 0 WN 1 WN 2 WN 3 WN 0 WN 2 WN 0 WN 2 WN 0 WN 0 WN 0

10、WN 0 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) A(0) 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) 8點DIT-FFT運算流圖 從上面的蝶形算法，當(dāng) 時，其運算應(yīng)該有 M級蝶形，每一級都由N/2個蝶形運算構(gòu)成。因此每一級運算 都需要N/2次復(fù)數(shù)乘法和N次復(fù)數(shù)加法，所以M

11、級蝶形運算的乘 法次數(shù)為： 加法次數(shù)： 2MN 2 =log 22 N M NN CM 2 log N A CN MN 一個簡單的算例 計算序列 的8點DFT。分別 用基本的DFT算法和FFT算法實現(xiàn)，體會計算過程中的時 間復(fù)雜度。 當(dāng)然針對計算機來說，計算乘法所需要的時間遠(yuǎn)大于加 法。一般的計算機，差不多相差十倍左右。 用計算機產(chǎn)生隨機的1024個點構(gòu)成序列，然后取 N=1024.此時計算時間差距就會加大。N=2048時，時間 差距會更大。 ( )1 1 1 1x n ， ， ， 為了更進(jìn)一步的體會FFT的物理意義，引入一個算例： 假設(shè)對某信號經(jīng)過ADC之后，得到一個序列 ，此時我們 不知道

12、其具體的函數(shù)表達(dá)式。但是我們可以對其做FFT運算。 現(xiàn)在我們做一個驗證： 假設(shè)我們有一個信號，它含有2V的直流分量，頻率為50Hz、相位為- 30度、幅度為3V的交流信號，以及一個頻率為75Hz、相位為90度、 幅度為1.5V的交流信號。用數(shù)學(xué)表達(dá)式就是如下： S=2+3*cos(2*pi*50*t-pi*30/180)+1.5*cos(2*pi*75*t+pi*90/180) 現(xiàn)在對其進(jìn)行變換，取樣點為1024，采樣頻率為1024Hz 注意這里的取樣頻率只要滿足原始頻率的2倍即可，且取樣 點和取樣頻率根據(jù)頻率分辨率來選取。 ( )x n 020040060080010001200 -3 -2

13、 -1 0 1 2 3 4 5 6 7 原 始 信 號 t/s 幅度大小 0100200300400500600 0 0.5 1 1.5 2 2.5 3 X: 0 Y: 2 幅 度 -頻 率 曲 線 圖 f/Hz 幅度 X: 50 Y: 3 X: 75 Y: 1.5 0100200300400500600 -200 -150 -100 -50 0 50 100 150 200 X: 50 Y: -30 相 位 -頻 率 曲 線 圖 f/Hz 相位 X: 75 Y: 90 X: 0 Y: 0 從圖中的結(jié)果可以得出，當(dāng)頻率為0、50、75時，對 應(yīng)的幅度值依次為2、3、1.5，相位依次為0、90、

14、-30。當(dāng) 然，這看上去幾乎是沒有誤差的，因為我們?nèi)宇l率和所取 點數(shù)比較大，當(dāng)N=Fs=256時，會存在誤差，但是這個誤差 完全不影響我們對信號函數(shù)的分析與判斷。從而進(jìn)一步的驗 證了DFT與FFT算法的正確性。 上圖是從模擬信號到頻域離散信號的完整的過程。這其中對 應(yīng)有幾個概念容易混淆。因此對此做出區(qū)分： 數(shù)字頻率、模擬頻率與采樣頻率：模擬頻率通常用 表示，數(shù)字頻率用 表示。此時的數(shù)字頻率主要是針對序列 而言，因此沒有采樣頻率 就不會有數(shù)字頻率的概念，所以 數(shù)字頻率與采樣頻率和模擬頻率一定滿足某種關(guān)系。我們知道 有如下過程： 因此 其中 為取樣間隔。對應(yīng)的采樣頻率 f s F ( )sin(

15、2) a x tft ( )( )sin(2)sin() s as t nT x nx tfnTn 2 s fT s T 1 s s F T DFT取樣點N與采樣頻率和頻率分辨率（步進(jìn)值）的關(guān)系 首先我們根據(jù)奈奎斯特定律得到： ， 為連續(xù)信號的截至（最高）頻率（因為它可能有很多頻率 成分）?，F(xiàn)在問題就是取樣點N的選??？現(xiàn)在我們從序列的 FT說起！ FT完成的任務(wù)就是對時域連續(xù)信號抽樣之后通過序列的傅 里葉變換得到的頻譜關(guān)系，我們知道這個譜是連續(xù)的。此 時的譜是一個周期譜，那么周期是多少呢？是不是與時域 采樣頻率有關(guān)系呢？答案是肯定的，此時的周期就是 （注意此時是頻域，頻譜圖中的橫坐標(biāo)表示的頻率

16、）。 2 sc Ff c f s F 對于周期連續(xù)的譜，計算機還是無法操作，因此我們還 得對FT變換之后的譜進(jìn)行頻域采樣（DFT的實質(zhì)）。類似于 時域采樣，頻域FT譜線的采樣頻率為多少合適呢？顯然此時 只需一個周期即可，即在一個周期中取N個點，滿足 ： 式中 為頻率分辨率， 為采樣頻率。 因此關(guān)于頻率分辨率和取樣點N的關(guān)系就得出來了。 s FN f f s F 首先必須得畫出序列對應(yīng)的FT譜線圖，橫坐標(biāo)為 ，縱坐 標(biāo)為頻率。此時的FT譜線應(yīng)該是對上圖的周期延拓。這就是 我們一般只取其主值-pi,pi得原因。 得到FT譜線圖之后，接下來就是DFT對其進(jìn)行取樣！ 由于DFT的物理意義就是對0,2*

17、pi的等間隔取樣。即 22 0,1,239 40 kkk N 從題目給定的頻譜圖得出模擬頻率應(yīng)該是連續(xù)分 布的。但是離散化之后只能得到離散的模擬頻率， 由 確定。但是不能超過200Hz。 k fk f DFT算法與FFT算法耗時對比基于matlab編程 DFTDFT函數(shù)函數(shù) function Xk=dft(xn,N) n=0:1:N-1; k=0:1:N-1; WN=exp(-j*2*pi/N); nk=n*k; WNnk=WN.nk; Xk=xn*WNnk; end DFTDFT測試代碼：測試代碼： clc clear N=4096; xn=rand(1,4096); tic Xk=dft(

18、xn,4096); toc Elapsed Elapsed time is time is 10.993700 seconds.10.993700 seconds. 由于Matlab工具箱自帶FFT算法，因此調(diào)用庫函數(shù)即可。 但是可能存在一定的問題。庫函數(shù)里面的FFT算法性能可能會 更加優(yōu)化，耗時會更少。 N=4096; xn=rand(1,4096); tic Xk=fft(xn,4096); toc Elapsed time is 0.002649 seconds. 從耗時來看，同樣是4096個點， fft算法所需時間遠(yuǎn)小于DFT算法 時間。（做一次實驗的結(jié)果）由 于FFT具體的算法不得而知，因 此該實驗只能說明，F(xiàn)FT快速算 法