快速傅里葉變換MATLAB

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

end%時(shí)域序列倒序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⑺相應(yīng)的MATLAB程序：%自編FFT%基2頻域抽取DIF-FFT%輸入乂，輸出X均為行向量functionX=myfft2(x)iflength(x)~=2Afix(log2(length(x)))%如果長(zhǎng)度超出，補(bǔ)足下一個(gè)幕的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]));倒序處理

%%%實(shí)現(xiàn)十進(jìn)制序列的二進(jì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思想：先對(duì)矩陣每行作變換，再對(duì)每列變換，但要注意MATLAB里轉(zhuǎn)置是共軛轉(zhuǎn)置，單獨(dú)的轉(zhuǎn)置是.’。%自編二維傅里葉變換functionF=myfftdimention2(X)%先對(duì)每一行向量變換，再對(duì)每一列向量變換。[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