快速傅里叶变换MATLAB

快速傅里叶变换（MATLAB版）DFT是信号分析与处理中的一种重要变换。但直接计算DFT的计算量与变换区间长度N的平方成正比，当N较大时，计算量太大，直接用DFT算法进行谱分析和信号的实时处理是不切实际的。1965年库利和图基发现了DFT的一种快速算法，使DFT的运算效率提高1-2个数量级，为数字信号处理技术应用于各种信号的实时处理创造了条件，推动了数字处理技术的发展。它的思想是不断地把长序列的DFT分解成几个短序列的DFT，并利用旋转因子的周期性和对称性来减少DFT的运算次数。假设N=2M，全部计算分解为M级，利用旋转因子WNk的周期性、对称性进行合并、归类处理，以减少DFT的运算次数。周期性：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); %旋转因子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%偶数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); %旋转因子fork=1:Kfori=0:2A(k-1)-1TOC\o"1-5"\h\zforj=0:2A(K-k)-1ifmod(k,2)==1%奇数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%偶数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]));倒序处理

%%%实现十进制序列的二进制倒序输出输入为自然序列，从0开始，2AM0输出为它的倒序数，从%%%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里转置是共轭转置，单独的转置是.’。%自编二维傅里叶变换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