快速傅里葉變換

1、2021/3/231n自己能編寫自己能編寫FFT算法算法n理論上理解理論上理解FFT算法算法2021/3/232n直接計算直接計算DFT的問題及改進(jìn)的途徑的問題及改進(jìn)的途徑 2021/3/233 2021/3/234 設(shè)復(fù)序列設(shè)復(fù)序列x(n) 長度為長度為N點(diǎn)點(diǎn),其其DFT為為10( )( )NnkNnX kx n Wk=0,N-1 （1）計算一個）計算一個X(k) 值的運(yùn)算量值的運(yùn)算量復(fù)數(shù)乘法次數(shù)復(fù)數(shù)乘法次數(shù): N復(fù)數(shù)加法次數(shù)復(fù)數(shù)加法次數(shù): N1（2）計算全部）計算全部N個個X(k) 值的運(yùn)算量值的運(yùn)算量復(fù)數(shù)乘法次數(shù)復(fù)數(shù)乘法次數(shù): N2復(fù)數(shù)加法次數(shù)復(fù)數(shù)加法次數(shù): N(N1)2021/3/23

2、5N點(diǎn)點(diǎn)DFT的復(fù)數(shù)乘法次數(shù)舉例的復(fù)數(shù)乘法次數(shù)舉例NN2NN22464404941612816384864256 65 536 16256512 262 144 3210241024 1 048 576 結(jié)論結(jié)論:當(dāng)當(dāng)N很大時很大時,其運(yùn)算量很大其運(yùn)算量很大,對實時性很強(qiáng)的信號處理對實時性很強(qiáng)的信號處理來說來說,要求計算速度快要求計算速度快,因此需要改進(jìn)因此需要改進(jìn)DFT的計算方法的計算方法,以大以大大減少運(yùn)算次數(shù)。大減少運(yùn)算次數(shù)。 2021/3/236 nkNW主要原理是利用系數(shù)主要原理是利用系數(shù) 的以下特性對的以下特性對DFT進(jìn)行分解進(jìn)行分解: nkNW（2）對稱性）對稱性 ()nkNW(

3、)k N nNW（1）周期性）周期性 ()()n N kn kNnkNNNWWW（3）可約性）可約性 mnknkmNNWW/nknk mNN mWW另外另外,12/NNWkNNkNWW)2/(2021/3/237DIT-FFT(Decimation-In-Time)2021/3/2381(2 )( )xrx r設(shè)設(shè)N2M,將將x(n)按按 n 的奇偶分為兩組的奇偶分為兩組: 2(21)( )xrx rr =0,1, 12N10( ) ( )( )NnkNnX kDFT x nx n W則則1010)()(NnnnkNNnnnkNWnxWnx為奇數(shù)為偶數(shù)2021/3/239120)12(1202

4、) 12()2(NrkrNNrrkNWrxWrx1202212021)()(NrrkNkNNrrkNWrxWWrx)()(21kXWkXkN1100( )( )( )NNnknkNNnnnnX kx n Wx n W為偶數(shù)為奇數(shù)式中式中,X1(k)和和X2(k)分別是分別是x1(n)和和x2(n)的的N/2的的DFT。另外另外,式中式中k的取值范圍是的取值范圍是:0,1, ,N/21 。2021/3/2310因此因此, 只能計算出只能計算出X(k)的前一半值。的前一半值。12( )( )( )kNX kX kW Xk后一半后一半X(k) 值值, N/2 , N/2 1, ,N-1 ?1()2N

5、Xk2 1(2)120( )Nr NkNrx r W2 1120( )NrkNrx r W1( )X k同理可得同理可得22()( )2NXkXk2021/3/2311(2)NkkNNWW 因此可得后半部分因此可得后半部分X(k) )2()2()2(221NkXWNkXNkXNkN12( )( )( )kNX kX kW Xk由前半部分由前半部分X(k) )()(21kXWkXkNk=0,1, ,N/21k=0,1, ,N/212021/3/2312k=k=0,10,1, , , ,N/2N/21 1結(jié)論結(jié)論:12( )( )( )kNX kXkW Xk12()( )( )2kNNX kX k

6、W Xk2 11120( )( )NrkNrXkx r W2 12220( )( )NrkNrXkx r W因此因此, ,只要求只要求出出2 2個個N/2N/2點(diǎn)點(diǎn)的的DFTDFT, ,即即X X1 1(k)(k)和和X X2 2(k)(k), ,再經(jīng)再經(jīng)過這種運(yùn)算過這種運(yùn)算就可求出全就可求出全部部X X(k)(k)的值的值2021/3/231312( )( )( )kNX kX kW Xk12()( )( )2kNNX kX kW Xk蝶形運(yùn)算式蝶形運(yùn)算式蝶形運(yùn)算蝶形運(yùn)算信號流圖符號信號流圖符號 X1(k)X2(k)蝶形運(yùn)算的運(yùn)算量蝶形運(yùn)算的運(yùn)算量:每次蝶形含一次復(fù)數(shù)乘和兩每次蝶形含一次復(fù)數(shù)

7、乘和兩 次復(fù)數(shù)加次復(fù)數(shù)加2021/3/23140NW1NW2NW3NW以以N=8為例為例,分分解為解為2個個4點(diǎn)的點(diǎn)的DFT,然后做然后做8/2=4次蝶形運(yùn)次蝶形運(yùn)算即可求出所算即可求出所有有8點(diǎn)點(diǎn)X(k)的值。的值。2021/3/2315復(fù)數(shù)乘法次數(shù)復(fù)數(shù)乘法次數(shù): N2復(fù)數(shù)加法次數(shù)復(fù)數(shù)加法次數(shù): N(N1)復(fù)數(shù)乘法次數(shù)復(fù)數(shù)乘法次數(shù): 2*(N/2)2+N/2=N2/2+N/2復(fù)數(shù)加法次數(shù)復(fù)數(shù)加法次數(shù): 2*(N/2)(N/21)+2*N/2=N2/2nN點(diǎn) 每次蝶形含一次復(fù)數(shù)每次蝶形含一次復(fù)數(shù)乘和兩次復(fù)數(shù)加乘和兩次復(fù)數(shù)加2021/3/2316 由于由于N2M,因而因而N/2仍是偶數(shù)仍是偶數(shù)

8、,可以進(jìn)一步把每個可以進(jìn)一步把每個N/2點(diǎn)點(diǎn)子序列再按其奇偶部分分解為兩個子序列再按其奇偶部分分解為兩個N/4點(diǎn)的子序列。點(diǎn)的子序列。 以以N/2點(diǎn)序列點(diǎn)序列x1(r)為例為例 1314(2 )( )0,1,1(21)( )4xlx lNlxlx l則有則有 rkNNrWrxkX212011)()(klNNllkNNlWlxWlx)12(21401221401) 12()2(lkNNlkNlkNNlWlxWWlx41404241403)()()()(42/3kXWkXkNk=0,1, 14N2021/3/2317且且13/24( )( )4kNNXkXkWXkk=0,1, 14N由此可見由此可

9、見,一個一個N/2點(diǎn)點(diǎn)DFT可分解成兩個可分解成兩個N/4點(diǎn)點(diǎn)DFT。 同理同理,也可對也可對x2(n)進(jìn)行同樣的分解進(jìn)行同樣的分解,求出求出X2(k)。 2021/3/2318概念概念:信號流圖信號流圖2021/3/2319100(0)( )(0)(1)(0)(1)NnkNNnXx n WxxxW x10( )( )0,1NnkNnX kx n Wk 對此例對此例N=8,最后剩下的是最后剩下的是4個個N/4= 2點(diǎn)的點(diǎn)的DFT,2點(diǎn)點(diǎn)DFT也可以由蝶形運(yùn)算來完成。也可以由蝶形運(yùn)算來完成。N=2這說明這說明,N=2,N=2M M的的DFTDFT可全部由蝶形運(yùn)算來完成。可全部由蝶形運(yùn)算來完成。因

10、此因此:N:N必須是必須是2 2的整數(shù)次冪的整數(shù)次冪. . “基基2 2”100(1)( )(0)(1)(0)(1)NnkNNnXx n WxxxW x蝶形運(yùn)算蝶形運(yùn)算2021/3/2320N=8按時間抽取法按時間抽取法FFT信號流圖信號流圖 x(n)抽取抽取X(k)順序順序DIT-FFT(Decimation-In-Time)FFT算法算法2021/3/2321由按時間抽取法FFT的信號流圖可知,當(dāng)N=2M時,共有 級蝶形運(yùn)算;每級都由 個蝶形運(yùn)算組成,而每個蝶形有 次復(fù)乘、 次復(fù)加,因此每級運(yùn)算都需 次復(fù)乘和 次復(fù)加。 MN/2 N/2 12N2021/3/2322FFT算法中算法中,這樣

11、這樣 級運(yùn)算總共需要級運(yùn)算總共需要: M復(fù)數(shù)乘法復(fù)數(shù)乘法: 2log22NNMN復(fù)數(shù)加法復(fù)數(shù)加法: 2logN MNN直接直接DFT算法運(yùn)算量算法運(yùn)算量 復(fù)數(shù)乘法復(fù)數(shù)乘法: 復(fù)數(shù)加法復(fù)數(shù)加法: N2N(N1)直接計算直接計算DFT與與FFT算法的計算量之比為算法的計算量之比為MNNNNNM222log2log2N=2M2021/3/2323NN2計算量之比M NN2計算量之比M 2414.012816 38444836.641644.025665 5361 02464.0864125.4512262 1442 304113.816256328.010241 048 5765 120204.83

12、210288012.820484 194 30411 264372.464404919221.4NN2log2NN2log22021/3/23242021/3/2325N=1024; M=80;x=1:M,zeros(1,N-M);t=cputime;y1=fft(x,N);time_fft=cputime-t;t1=cputime;y2=dft(x,N);time_dft=cputime-t1;time_dft=6.0290time_fft=0.01002021/3/2326L=1L=2L=3倒倒序序2021/3/2327n同址運(yùn)算（原位運(yùn)算）n序列的逆序排列n蝶形運(yùn)算兩節(jié)點(diǎn)間的距離n相同旋

13、轉(zhuǎn)因子節(jié)點(diǎn)間的距離n (旋轉(zhuǎn)因子）的確定pNW2021/3/2328 某一列任何兩個節(jié)點(diǎn)某一列任何兩個節(jié)點(diǎn)k 和和j 的節(jié)點(diǎn)變量進(jìn)行蝶形運(yùn)算的節(jié)點(diǎn)變量進(jìn)行蝶形運(yùn)算后后,得到結(jié)果為下一列得到結(jié)果為下一列k、j兩節(jié)點(diǎn)的節(jié)點(diǎn)變量兩節(jié)點(diǎn)的節(jié)點(diǎn)變量,而和其他節(jié)而和其他節(jié)點(diǎn)變量無關(guān)。這種原位運(yùn)算結(jié)構(gòu)可以節(jié)省存儲單元點(diǎn)變量無關(guān)。這種原位運(yùn)算結(jié)構(gòu)可以節(jié)省存儲單元,降低降低設(shè)備成本。設(shè)備成本。)(kX)2(NkX)(1kX)(2kXkNW運(yùn)算前運(yùn)算前運(yùn)算后運(yùn)算后( )A j)(kA( )A j )(kA例例n同址運(yùn)算（原位運(yùn)算）2021/3/23292021/3/2330)(01221)()(BINMMDECn

14、nnnnn 由于由于 x(n) 被反復(fù)地按奇、偶分組被反復(fù)地按奇、偶分組,所以流圖輸所以流圖輸入端的入端的排列不再是順序的排列不再是順序的,但仍有規(guī)律可循但仍有規(guī)律可循: 因為因為 N=2M , 對于任意對于任意 n（0n N-1),可以用可以用M個個二進(jìn)制碼表示為二進(jìn)制碼表示為:10,01221nnnnnMM n 反復(fù)按奇、偶分解時反復(fù)按奇、偶分解時,即按二進(jìn)制碼的即按二進(jìn)制碼的“0” “1” 分解。分解。n序列的逆序排列2021/3/23312021/3/2332自然順序自然順序 n二進(jìn)制數(shù)二進(jìn)制數(shù)倒位序二進(jìn)制數(shù)倒位序二進(jìn)制數(shù)倒位序順序數(shù)倒位序順序數(shù)00000000100110042010

15、01023011110641000011510110156110011371111117 nIJ規(guī)律規(guī)律:最左端加最左端加1,向右進(jìn)位向右進(jìn)位倒序倒序2021/3/2333I=J 不換不換IJ 不換不換IX= czt(x, M, W, A)n其中其中,x是待變換的時域信號是待變換的時域信號x(n),其長度為其長度為N,M是變換是變換長度長度,W確定變換的步長確定變換的步長,A確定變換的起點(diǎn)。若確定變換的起點(diǎn)。若M= N,A= 1,則則CZT變成變成DFT。2021/3/237610( )( )( )( ) ()Nmy nx nh nx m h nm設(shè)一離散線性移不變系統(tǒng)的沖激響應(yīng)為設(shè)一離散線性

16、移不變系統(tǒng)的沖激響應(yīng)為 ,其輸其輸入信號為入信號為 .其輸出為其輸出為 .并且并且 的長度為的長度為N N點(diǎn)點(diǎn), 的長度為的長度為M M點(diǎn)點(diǎn),計算其線性卷積。計算其線性卷積。)(nx)(nx)(nh)(nh)(ny)(ny)(nx)(nh2021/3/2377用FFT算 的步驟: )(n ny y1.( ), ( )1 ;x n h nLNM將補(bǔ)零點(diǎn)，至少為點(diǎn)2.( )( );H kFFT h nLFFT求, 點(diǎn)3.( )( )X kFFT x n求,L點(diǎn)FFT；求)()()(.4kHkXkY。求)()(.5kYIFFTny2021/3/2378FFTFFTIFFTxx(n)h(n)y(n)X

17、(k)H(k)Y(k)流程圖流程圖程序參考第三章幻燈片程序參考第三章幻燈片2021/3/2379例例4.8.2 設(shè)模擬信號設(shè)模擬信號 ,以以 t= 0.01n (n=0: N-1) 進(jìn)行取樣進(jìn)行取樣,試用試用fft函數(shù)對其做頻譜分析。函數(shù)對其做頻譜分析。N分別為分別為:(1) N=10;(2) N=50;(3) N=100;(2) N=1024。 ( )2sin(20)5cos(40)x ttt程序清單如下程序清單如下 %計算計算N=10的的FFT并繪出其幅頻曲線并繪出其幅頻曲線N=10;n=0:N-1;t=0.01*n;q=n*2*pi/N;x=2*sin(20*pi*t)+5*cos(40*pi*t);y=fft(x,N);figure(1)subplot(2,2,1)plot(q,abs(y)title(FFT N=10)2021/3/2380%計算計算N=50的的FFT并繪出其幅頻曲線并繪出其幅頻曲線N=50;n=0:N-1;t=0.01*n;q=n*2*pi/N;x=2*sin(20*pi*t)+5*cos(40*pi*t);y=fft(x,N);figure(1)subplot(2,2,2)plot(q,abs(y)title(FFT N=50)2021/3/2381%計算計算N=100的的FFT并繪出其幅頻曲線并繪出其幅頻曲線N=100;n=0:N-1;t=0.0