第3章.離散傅里葉變換及其快速算法與應(yīng)用

1、n離散傅里葉變換（DFT)的定義nDFT與DTFT和Z變換的關(guān)系nDFT的數(shù)字計算nDFT的性質(zhì)nDFT 的快速算法FFTnDFT(FFT)的應(yīng)用第第3 3章章 離散傅里葉變換及快速算法與應(yīng)用離散傅里葉變換及快速算法與應(yīng)用 傅里葉變換對的數(shù)據(jù)如果連續(xù)或無限，則不適合在計算機上運算。從數(shù)值計算的角度出發(fā)，我們感興趣的是時域及頻域都離散且有限的情況，這就是我們本章要研究的離散傅里葉變換。離散傅里葉變換(Discrete Fourier Transform，簡稱，簡稱DFT)是數(shù)字信號處理中非常重要的一種數(shù)學(xué)變換，除了在理論上相當(dāng)重要之外，因其存在有效的快速算法，故在各種數(shù)字信號處理的算法中起核心作

2、用。引言引言3.1.1 離散傅里葉變換的定義離散傅里葉變換的定義10210(0 -1)(:( ) ( )( )0,1,11( )( )( )0,1,1NknNnjNnNknNkx nM nMx nNDFTX kDFT x nx n WkNWeNMx nIDFT X kX k WnkNN）是長為 （的有限長序列，則定義）的 點其中為旋轉(zhuǎn)因子，離散傅里葉逆變換IDFT：3.1.2 DFT與與DTFT和和Z變換的關(guān)系變換的關(guān)系210101202( ),0,1,-1( ) ( )( )() ( )( )( ) ( )( )()|( )|20,1,1,jkNMnnMjj nnMjknNNnjkz eNx

3、 nM nMX zZT x nx n zX eDTFT x nx n eX kDFT x nx n eX eX zkNkN長表示頻率：結(jié)論：結(jié)論： DFT是有限長序列x(n)的傅里葉變換在區(qū)間0，2上的N點等間隔采樣；DFT是對有限長序列x(n) 進行Z變換后在單位圓上進行等間隔采樣的采樣值 。3.1.2 DFT與與DTFT和和Z變換的關(guān)系變換的關(guān)系kN223.1.4 DFT的數(shù)字計算的數(shù)字計算（直接定義）（直接定義）functionX=dft(x,N)dw=2*pi/N;n=0:N-1;k=0:N-1;M=length(x);if NM for l=M+1:N x(l)=0; endendX

4、=x*(exp(-j*dw).(n*k);DFT計算子函數(shù)：IDFT計算子函數(shù)：3.1.4 DFT的數(shù)字計算的數(shù)字計算（直接計算）（直接計算）functionx=idft(X,N)dw=2*pi/N;n=0:N-1;k=0:N-1;x=X*(exp(j*dw).(n*k)/N;3.1.4 DFT的數(shù)字計算的數(shù)字計算（直接計算）（直接計算）41 ( )( ),( )48.x nR nx nDFT例 求的 點和 點x=1 1 1 1 ;n=0:3;k=-1500:1500;dw=(pi/500);X,w=dift(x,n,dw,k);magX=abs(X);figure(1)plot(w/pi,m

5、agX);N=4;X=dft(x,N);MX=abs(X);figure(2)stem(MX);x1=idft(X,N);figure(3)stem(real(x1)n1）隱含的周期性n 2）循環(huán)移位性質(zhì)3.2 DFT的的性質(zhì)性質(zhì)( )( ).knNWX kx nN因的周期性，和隱含周期性，且周期均為11( ),( )()( )( ),011,( ):mod( ,)NNNNx nMy nx nmR nnNnnMNnnNMnnmatlabn N長為則其循環(huán)移位定義：表示模 對 求余：為整數(shù)循環(huán)移位的實現(xiàn)子程序：function y=cirshift(x,m,N)x=x zeros(1,N-len

6、gth(x);n=0:N-1;n=mod(n-m,N);y=x(n+1);152( )10 (0.8) ,010,(6)nx nnx n例 已知求n=0:10;x=10*(0.8).n;y=cirshift(x,6,15);n=0:14;x=x,zeros(1,4);subplot(2,1,1)stem(n,x);subplot(2,1,2)stem(n,y);024681012140510 xn024681012140246810 x(n-6)3.2 DFT的的性質(zhì)性質(zhì)n2）循環(huán)移位性質(zhì)時域循環(huán)移位定理頻率循環(huán)移位定理( )()( )( )( )NNkmNy nx nmR nY kWX k則

7、( )()( )( )( )NNnlNY kXklR ky nW x n則n3）循環(huán)卷積定理1.兩個有限長序列的循環(huán)卷積3.2 DFT的的性質(zhì)性質(zhì)10( )( ). ( )( )( )( ) () ( )= ( )( )max,LcLLmh nx nNM h nx nLy nh m x nmR nh nx nLN M設(shè)與的長度分別為 和與的的循環(huán)卷積：循環(huán)卷積的計算：function y=circonv(x1,x2,L)x1=x1 zeros(1,L-length(x1);x2=x2 zeros(1,L-length(x2);n=0:L-1;ix2=x2(mod(-n,L)+1);for k=

8、0:L-1 y1=cirshift(ix2,k,L); yc(k+1)=x1*y1;endy=yc;123( )1,2,2 ,( )1,2,3,4 ,04.x nx nn例 已知求其4點、5點和6點循環(huán)卷積。x1=1 2 2;x2=1 2 3 4;y4=circonv(x1,x2,4);y5=circonv(x1,x2,5);y6=circonv(x1,x2,6);y4 = 15 12 9 14y5 = 9 4 9 14 14y6 = 1 4 9 14 14 8n時域循環(huán)卷積定理n頻域循環(huán)卷積定理3.2 DFT的的性質(zhì)性質(zhì)1212122121( )( )max,( )( )( ),( )( )

9、( )( )x nx nNNNN Nx nx nx nx nNDFTX kXkX k和的長度分別為和，則的 點為：1212( )( )( ),1( )( )( )x nx n x nX kX kXkN如果則循環(huán)卷積定理驗證12124( )1,2,2,0 ,( )1,2,3,4 ,04.4( )( )( )( ),x nx nny nX kXkY k例 已知求其 點循環(huán)卷積、并驗證循環(huán)卷積定理。x1=1 2 2 0;x2=1 2 3 4;N=4;y=circonv(x1,x2,N); X1=dft(x1,N);X2=dft(x2,N);Y=dft(y,N); yi=idft(X1.*X2,N);

10、%yi=y?nDFT的對稱性0 -1,/ 2NN因DFT定義區(qū)間為這里的對稱是指關(guān)于的對稱性。( )( ) ( ) ,( )(- )-|( )| |(- )|-arg( )-arg(- )-Nx nX kDFT x nXkX N kX kX N kX kX N k若是實序列，則共軛對稱模復(fù)角15( )1,2,3,4 ,010.x nn例 其DFT:X1 = 10 , -2 + 2i , -2 , -2 - 2i3.3 頻率域采樣頻率域采樣( ),(),0 2(),)( )( )( )()( ),( )( )jjNNiNx nMDTFTX eX eNX kX kx nNxnIDFT X kx n

11、iN R nNM xnx n如果序列的長度為其為在對等間隔抽樣 個點得到抽樣序列（能由（重建嗎？ 應(yīng)滿足何條件？110( )( )11( )( )1NNkkNNMX kX zzX zX kNW z內(nèi)插公式：時，可由從建例 3.3.1M=26;n=0:M;xa=0:M/2;xb=ceil(M/2)-1:-1:0;xn=xa,xb;N=27;k=0:N-1;dw=2*pi/N;X,w=dift(xn,n,dw,k);figure(1)subplot(1,2,1)stem(n,xn);xlabel(xn);x1=idft(X,N);subplot(1,2,2)stem(k,real(x1)xlabe

12、l(x1);010203002468101214xn010203002468101214x1N=27010203002468101214xn05101502468101214x1N=16010203002468101214xn01020304002468101214x1N=32 直接直接DFT存在的問題存在的問題結(jié)論：結(jié)論：直接DFT運算共需N N2 2次復(fù)數(shù)乘法次復(fù)數(shù)乘法, N(N-1)N(N-1)次復(fù)數(shù)加法次復(fù)數(shù)加法。由于一次復(fù)乘需4次實乘和2次實加，一次復(fù)加需2次實加， 因此直接DFT運算共需4N2次實乘及2N(2N-1)次實加。當(dāng)N很大時運算量非常大，如果要進行實時信號處理就必須設(shè)法減

13、少運算量。 1,.,1 , 0 )()(DFT)(10NkWnxnxkXNnknN4.DFT 的快速算法的快速算法FFT 減少減少DFT運算量的方法運算量的方法思路：思路：將長序列分解為短序列，利用DFT系數(shù)的對稱性和周期性，分解計算再合并處理，可大大減少運算量。 方法：方法：根據(jù)分解與合并方法的不同，有多種快速算法?；舅惴煞譃閮深悾?按時間抽取法（Decimation in Time，縮寫為DIT） 按頻率抽取法（Decimation in Frequency，縮寫為DIF）時間抽取法：時間抽取法：將序列按時間順序的奇偶分解為越來越短的子序列進行計算的快速算法，也稱Cooley-Tuke

14、y算法。頻率抽取法：頻率抽取法：將序列按時間順序的前后分解為越來越短的子序列進行計算的快速算法， 也稱Sande-Tukey算法。 FFT( Fast Fourier Transform )并不是一種新的傅里葉變換。FFT是計算離散傅里葉變換（DFT）的快速算法。重點以按時間抽取的基2FFT(DIT-FFT)算法為例，介紹FFT算法的基本原理，歸納FFT算法的運算規(guī)律，總結(jié)FFT算法的編程框圖，列出FFT算法的實現(xiàn)程序，從而建立起一種從算法理論到程序?qū)崿F(xiàn)全過程的完整概念。 另外，簡單介紹進一步減少運算量的措施和其它快速算法。 快速傅里葉變換快速傅里葉變換將x(n)按 n 的奇偶分為兩組：)()

15、()()() 12()2()()()()(211202212021120)12(12021210kXWkXWrxWWrxWrxWrxWnxWnxWnxkXkNNrrkNkNNrrkNNrkrNNrrkNNnnkNNnnkNNnnkN為奇為偶 基基2DIT-FFT2DIT-FFT算法的算法的算法原理算法原理 DIT蝶形運算流圖符號)()()(21kXWkXkXkN)()()2()2()2(212)2(1kXWkXNkXWNkXNkXkNNkN說明：式中得到的X(k)只有N/2個點，實際上X(k)有N個點，要用X1(k)和X2(k)表示全部的X(k)值，還要應(yīng)用系數(shù)的周期性和對稱性。 基基2DIT

16、-FFT2DIT-FFT算法的算法的運算流程運算流程 按時間抽取將一個N點DFT分解為兩個N/2點DFT（N=8）8點DFT二次時域抽取分解運算流圖N=8 點DIT-FFT運算流圖例：計算x(n)=0,1,2,3,4,5,6,7的DFT值X(k) 基基2DIT-FFT2DIT-FFT算法的算法的運算規(guī)律運算規(guī)律 蝶形運算：蝶形運算：蝶形運算是分級進行的；每級的蝶形運算可以按旋轉(zhuǎn)因子的指數(shù)大小排序進行；如果指數(shù)大小一樣則可從上往下依次蝶算。蝶形運算的規(guī)律是編程的基礎(chǔ) 。原位計算：原位計算：當(dāng)數(shù)據(jù)輸入到存儲器以后，每一級運算的結(jié)果仍可儲存在同一組存儲器中，直到最后輸出，中間無需其它存儲器，這叫原位

17、計算。 序列倒序序列倒序 ：運算前要先對輸入的序列進行位序顛倒 。順序數(shù)（I） 順序二進制（IB） 倒序二進制（JB） 倒序數(shù) （J） 0123456700000101001110010111011100010001011000110101111104261537規(guī)律：規(guī)律：倒序二進制數(shù)是順序二進制數(shù)的倒置。 基基2DIT-FFT2DIT-FFT算法的算法的程序框圖程序框圖 基基2DIT-FFT2DIT-FFT算法的實現(xiàn)算法的實現(xiàn)程序程序 clc;close all;clear;format compact;%倒序xn=0,1,2,3,4,5,6,7;%可取任意序列M=nextpow2(len

18、gth(xn), N=2M,%將數(shù)據(jù)輸入到存儲地址A=xn,zeros(1,N-length(xn); disp(輸入到各存儲單元的數(shù)據(jù):),disp(A);%數(shù)據(jù)倒序操作J=0;%給倒序數(shù)賦初值for I=0:N-1;%按順序交換數(shù)據(jù)和算倒序數(shù) if I=K; J=J-K;K=K/2; end J=J+K;enddisp(倒序后各存儲單元的數(shù)據(jù):),disp(A);%分級按序依次進行蝶形運算WN=exp(-j*2*pi/N);%計算常量for L=1:M;%分級計算 disp(運算級次:),disp(L); B=2(L-1); for R=0:B-1;%各級按序計算 P=2(M-L)*R;

19、for K=R:2L:N-2;%每序依次計算 T=A(K+1)+A(K+B+1)*WNP; A(K+B+1)=A(K+1)-A(K+B+1)*WNP; A(K+1)=T; end end disp(本級運算后各存儲單元的數(shù)據(jù):),disp(A); enddisp(輸出各存儲單元的數(shù)據(jù):),Xk=A,disp(調(diào)用fft函數(shù)運算的結(jié)果:),fftxn=fft(xn,N), 基基2DIT-FFT2DIT-FFT算法的運算量分析算法的運算量分析 對于任何一個2的整數(shù)次冪 N=2M 點的DFT運算，總可以通過M 次分解，直到DFT 就是其本身的1點序列。這樣的M次分解，構(gòu)成從x(n)到X(k)的M 級

20、運算。從流程圖可看到，每級運算都由N/2個蝶形運算構(gòu)成。因此每級運算都需要 (N/2) 次復(fù)乘和N次復(fù)加(每個結(jié)作加、減各一次)。故總共需要的運算： 復(fù)乘 (N/2)M=(N/2)log2N 復(fù)加 NM=N log2N 實際運算量與這個數(shù)字稍有出入，但按時間抽取法所需的計算量，大概與N log2N 成正比，而直接運算的運算量與N2成正比，當(dāng)N較大時，F(xiàn)FT算法優(yōu)勢明顯。 IDFT IDFT的快速算法的快速算法方法方法1 1：只要把DFT運算中的蝶形系數(shù)W nkN 改為W -nkN 。結(jié)果再除N，則FFT運算程序可直接進行IDFT運算。FFT算法可用于IDFT運算，簡稱為IFFT。方法方法2 2

21、：不改FFT程序，直接用它作IFFT(1)將X（k）取共軛（虛部乘以-1） (2)然后直接作FFT運算(3)最后對FFT結(jié)果取共軛并除N得x(n) ( )X k 求共軛( )XkFFT 求( )FFT Xk( )FFT XkN 除以( )x n 求共軛1010)(1)()()(NknkNNnnkNWkXNnxWnxkX)(1)(1)()(1)(1010kXDFTNWkXNnxWkXNnxNknkNNknkNn3.4 DFT 的應(yīng)用舉例DFT的廣泛應(yīng)用是在其快速算法FFT出現(xiàn)之后，主要應(yīng)用在卷積和相關(guān)的具體計算，及用DFT近似計算信號的傅里葉變換。本節(jié)介紹DFT的兩個重要應(yīng)用：一是用DFT計算卷積；二是用DFT對連續(xù)信號或序列進行譜分析。 3.4.1 用DFT計算線性卷積補零至 L求FFT)()()(nhnxny補零至 L求FFT)()()(kHkXkY