快速傅里葉變換MATLAB

快速傅里葉變換（MATLAB版）DFT是信號分析與處理中的一種重要變換。但直接計算DFT的計算量與變換區(qū)間長度N的平方成正比，當N較大時，計算量太大，直接用DFT算法進行譜分析和信號的實時處理是不切實際的。1965年庫利和圖基發(fā)現(xiàn)了DFT的一種快速算法，使DFT的運算效率提高1-2個數(shù)量級，為數(shù)字信號處理技術(shù)應用于各種信號的實時處理創(chuàng)造了條件，推動了數(shù)字處理技術(shù)的發(fā)展。它的思想是不斷地把長序列的DFT分解成幾個短序列的DFT，并利用旋轉(zhuǎn)因子的周期性和對稱性來減少DFT的運算次數(shù)。假設N=2M，全部計算分解為M級，利用旋轉(zhuǎn)因子WNk的周期性、對稱性進行合并、歸類處理，以減少DFT的運算次數(shù)。周期性：Wm=Wm+lN對稱性：[聽^]*=Wm，Wm+N/2=-WmFFT算法基本上分為兩大類：時域抽取法FFT（簡稱DIT-FFT）和頻域抽取法FFT（簡稱DIF一FFT）。下面略去推導和證明，僅以長度為8的序列為例子說明這兩種算法。1.DIT-FFTx(0)x(4)x⑵x⑹x⑴xx(0)x(4)x⑵x⑹x⑴x⑸x⑶x⑺WN°Wn2X■WN0L=1級L=2L=3X(0)X⑴X⑵X(3)X(4)XU)X⑹XV)N=8,5、DIT-FFT運算流圖相應的MATLAB程序：%自編FFT%基2時域抽取DIT-FFT%輸入乂，輸出X均為行向量functionX=myfftl(x)iflength(x)?=2Afix(log2(length(x)))%如果長度超出，補足下一個幕的0。x=[x,zeros(1,2Aceil(log2(length(x)))-length(x))];

end%時域序列倒序x=x(invertorder([0:length(x)-1]));N=length(x);K=log2(N);X2=zeros(1,N);X1=x;W_n=exp(-1j*2*pi/N); %旋轉(zhuǎn)因子fork=1:Kfori=0:2A(K-k)-1forj=0:2A(k-1)-1X1(j+i*2Ak+1) +X1(j+i*2Ak+1) -X2(j+i*2Ak+1) X1(j+i*2Ak+1) +X1(j+i*2Ak+1) -X2(j+i*2Ak+1) +X2(j+i*2Ak+1) -X2(j+i*2Ak+2A(k-1)+1)W_nA(j*2A(K-k))*X1(j+i*2Ak+2A(k-1)+1);else%偶數(shù)X1(j+i*2Ak+1)W_nA(j*2A(K-k))*X2(j+i*2Ak+2A(k-1)+1);X1(j+i*2Ak+2A(k-1)+1)W_nA(j*2A(K-k))*X2(j+i*2Ak+2A(k-1)+1);endendendendif(mod(K,2)==1)X=X2;elseX=X1;endDIF-FFTx(0)x⑴x⑵x⑶x(4)x(5)x⑥x⑺滬xi(0)x?1)x3x(0)x⑴x⑵x⑶x(4)x(5)x⑥x⑺滬xi(0)x?1)x3(0)x5(0)0x2(0)x2(3)Wn2x2⑴x2⑵如："w 二x6(0) WN0%⑴二"二X(0)X(4)X⑵X⑥X⑴X(5)X6)X⑺相應的MATLAB程序：%自編FFT%基2頻域抽取DIF-FFT%輸入乂，輸出X均為行向量functionX=myfft2(x)iflength(x)~=2Afix(log2(length(x)))%如果長度超出，補足下一個幕的0。x=[x,zeros(1,2Aceil(log2(length(x)))-length(x))];endN=length(x);K=log2(N);X2=zeros(1,N);X1=x;W_n=exp(-1j*2*pi/N); %旋轉(zhuǎn)因子fork=1:Kfori=0:2A(k-1)-1TOC\o"1-5"\h\zforj=0:2A(K-k)-1ifmod(k,2)==1%奇數(shù)X2(j+i*2A(K-k+1)+1) = X1(j+i*2A(K-k+1)+1) +X1(j+i*2A(K-k+1)+2A(K-k)+1);X2(j+i*2A(K-k+1)+2A(K-k)+1) = (X1(j+i*2A(K-k+1)+1) -X1(j+i*2A(K-k+1)+2A(K-k)+1))*W_nA(j*2A(k-1));else%偶數(shù)X1(j+i*2A(K-k+1)+1) = X2(j+i*2A(K-k+1)+1) +X2(j+i*2A(K-k+1)+2A(K-k)+1);X1(j+i*2A(K-k+1)+2A(K-k)+1) = (X2(j+i*2A(K-k+1)+1) -X2(j+i*2A(K-k+1)+2A(K-k)+1))*W_nA(j*2A(k-1));endendendendif(mod(K,2)==1)X=X2;elseX=X1;end%頻域序列倒序X=X(invertorder([0:length(X)-1]));倒序處理

%%%實現(xiàn)十進制序列的二進制倒序輸出輸入為自然序列，從0開始，2AM0輸出為它的倒序數(shù)，從%%%LHuNJuLHeil(lTuX(I)A(LHuNJuLHeil(lTuX(I)A(I)uX(J)A(J)uTN=length(t);M=log2(N);s=zeros(1,N);J=t(1);%從0開始s(1)=1；fori=2:Nforj=1:MifJ>=N/2AjJ=J-N/2Aj;elses(i)=J+N/2Aj+1;J=s(i)-1;break;endendendKuLHJuJ-KJuJ+KuLHJuJ-KJuJ+K思想：先對矩陣每行作變換，再對每列變換，但要注意MATLAB里轉(zhuǎn)置是共軛轉(zhuǎn)置，單獨的轉(zhuǎn)置是.’。%自編二維傅里葉變換functionF=myfftdimention2(X)%先對每一行向量變換，再對每一列向量變換。[m,n]=size(X);if(m~=2Anextpow2(m))M=2Anextpow2(m);elseM=m;endif(n~=2Anextpow2(n))N=2Anextpow2(n);elseN=n;endF=zeros(M,N);x=zeros(1,n