快速傅里叶变换MATLAB..docVIP

快速傅里叶变换（MATLAB版） DFT是信号分析与处理中的一种重要变换。但直接计算DFT的计算量与变换区间长度N的平方成正比，当N较大时，计算量太大，直接用DFT算法进行谱分析和信号的实时处理是不切实际的。 1965年库利和图基发现了DFT的一种快速算法，使DFT的运算效率提高1－2个数量级，为数字信号处理技术应用于各种信号的实时处理创造了条件，推动了数字处理技术的发展。它的思想是不断地把长序列的DFT分解成几个短序列的DFT，并利用旋转因子的周期性和对称性来减少DFT的运算次数。 假设，全部计算分解为级，利用旋转因子WNk的周期性、对称性进行合并、归类处理，以减少DFT的运算次数。 周期性： 对称性：， FFT算法基本上分为两大类：时域抽取法FFT(简称DIT-FFT)和频域抽取法FFT(简称DIF―FFT)。下面略去推导和证明，仅以长度为8的序列为例子说明这两种算法。 1. DIT-FFT 相应的MATLAB程序： % 自编FFT % 基2时域抽取DIT-FFT % 输入x，输出X均为行向量 function X = myfft1(x) if length(x) ~= 2^fix(log2(length(x))) % 如果长度超出，补足下一个幂的0。 x=[x,zeros(1,2^ceil(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); % 旋转因子 for k=1:K for i=0:2^(K-k)-1 for j=0:2^(k-1)-1 if mod(k,2)==1 % 奇数 X2(j+i*2^k+1) = X1(j+i*2^k+1) + W_n^(j*2^(K-k))*X1(j+i*2^k+2^(k-1)+1); X2(j+i*2^k+2^(k-1)+1) = X1(j+i*2^k+1) - W_n^(j*2^(K-k))*X1(j+i*2^k+2^(k-1)+1); else % 偶数 X1(j+i*2^k+1) = X2(j+i*2^k+1) + W_n^(j*2^(K-k))*X2(j+i*2^k+2^(k-1)+1); X1(j+i*2^k+2^(k-1)+1) = X2(j+i*2^k+1) - W_n^(j*2^(K-k))*X2(j+i*2^k+2^(k-1)+1); end end end end if(mod(K,2) == 1) X=X2; else X=X1; end 2. DIF―FFT 相应的MATLAB程序： % 自编FFT % 基2频域抽取DIF-FFT % 输入x，输出X均为行向量 function X = myfft2(x) if length(x) ~= 2^fix(log2(length(x))) % 如果长度超出，补足下一个幂的0。 x=[x,zeros(1,2^ceil(log2(length(x)))-length(x))]; end N=length(x); K=log2(N); X2=zeros(1,N); X1=x; W_n=exp(-1j*2*pi/N); % 旋转因子 for k=1:K for i=0:2^(k-1)-1 for j=0:2^(K-k)-1 if mod(k,2)==1 % 奇数 X2(j+i*2^(K-k+1)+1) = X1(j+i*2^(K-k+1)+1) + X1(j+i*2^(K-k+1)+2^(K-k)+1); X2(j+i*2^(K-k+1)+2^(K-k)+1) = (X1(j+i*2^(K-k+1)+1) - X1(j+i*2^(K-k+1)+2^(K-k)+1)) * W_n^(j*2^(k-1)); else % 偶数 X1(j+i*2^(K-k+1)+1) = X2(j+i*2^(K-k+1)+1) + X2(j+i*2^(K-k+1)+2^(K-k)+1); X1(j+i*2^(