快速傅里叶变换MATLAB.docx
- 文档编号:290121
- 上传时间:2022-10-08
- 格式:DOCX
- 页数:8
- 大小:78.09KB
快速傅里叶变换MATLAB.docx
《快速傅里叶变换MATLAB.docx》由会员分享,可在线阅读,更多相关《快速傅里叶变换MATLAB.docx(8页珍藏版)》请在冰豆网上搜索。
快速傅里叶变换MATLAB
快速傅里叶变换(MATLAB版)
DFT是信号分析与处理中的一种重要变换。
但直接计算DFT的计算量与变换区间长度N的平方成正比,当N较大时,计算量太大,直接用DFT算法进行谱分析和信号的实时处理是不切实际的。
1965年库利和图基发现了DFT的一种快速算法,使DFT的运算效率提高1-2个数量级,为数字信号处理技术应用于各种信号的实时处理创造了条件,推动了数字处理技术的发展。
它的思想是不断地把长序列的DFT分解成几个短序列的DFT,并利用旋转因子的周期性和对称性来减少DFT的运算次数。
假设2M
N=,全部计算分解为M级,利用旋转因子WNk的周期性、对称性进行合并、归类处理,以减少DFT的运算次数。
周期性:
mmlNNN
WW+=对称性:
[]*mmNNWW-=,/2mNmNNWW+=-FFT算法基本上分为两大类:
时域抽取法FFT(简称DIT-FFT)和频域抽取法FFT(简称DIF―FFT)。
下面略去推导和证明,仅以长度为8的序列为例子说明这两种算法。
1.DIT-FFT
N=8点DIT-FFT运算流图x(0)x(4)
x
(2)
x(6)x
(1)x(5)
x(3)
(0)
(1)
(2)
(3)(4)(5)(6)x(7)NL=1级L=2
L=3X(7)
相应的MATLAB程序:
%自编FFT
%基2时域抽取DIT-FFT
%输入x,输出X均为行向量
functionX=myfft1(x)
iflength(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);%旋转因子
fork=1:
K
fori=0:
2^(K-k)-1
forj=0:
2^(k-1)-1
ifmod(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
x(0)x
(1)x
(2)x(3)
x(4)x(5)x(6)x(7)x
(0)
(0)
(4)
(2)
(6)
X
(1)
X(5)
X(3)
X(7)
x(0)
N
相应的MATLAB程序:
%自编FFT
%基2频域抽取DIF-FFT
%输入x,输出X均为行向量
functionX=myfft2(x)
iflength(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);%旋转因子
fork=1:
K
fori=0:
2^(k-1)-1
forj=0:
2^(K-k)-1
ifmod(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^(K-k+1)+2^(K-k)+1)=(X2(j+i*2^(K-k+1)+1)-X2(j+i*2^(K-k+1)+2^(K-k)+1))*W_n^(j*2^(k-1));
end
end
end
end
if(mod(K,2)==1)
X=X2;
else
X=X1;
end
%频域序列倒序
X=X(invertorder([0:
length(X)-1]));
3.倒序处理
%输入为自然序列,从0开始,2^M。
%输出为它的倒序数,从1开始。
functions=invertorder(t)
iflength(t)~=2^fix(log2(length(t)))%
end
N=length(t);
M=log2(N);
s=zeros(1,N);
J=t
(1);%从0开始
s
(1)=1;
fori=2:
N
forj=1:
M
ifJ>=N/2^j
J=J-N/2^j;
else
s(i)=J+N/2^j+1;
J=s(i)-1;
break;
end
end
end
4.二维傅里叶变换
思想:
先对矩阵每行作变换,再对每列变换,但要注意MATLAB里转置是共轭转置,单独的转置是.’。
%自编二维傅里叶变换
functionF=myfftdimention2(X)
%先对每一行向量变换,再对每一列向量变换。
[m,n]=size(X);
if(m~=2^nextpow2(m))
M=2^nextpow2(m);
else
M=m;
end
if(n~=2^nextpow2(n))
N=2^nextpow2(n);
else
N=n;
end
F=zeros(M,N);
x=zeros(1,n);
y=zeros(1,m);
fx=zeros(1,N);
fy=zeros(1,M);
fori=1:
M
x=X(i,:
);
fx=myfft1(x);
F(i,:
)=fx;
end
forj=1:
N
y=F(:
j).';%转置但不共轭
fy=myfft1(y);
F(:
j)=fy;
end
温馨推荐
您可前往XX文库小程序
享受更优阅读体验
不去了
立即体验
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 快速 傅里叶变换 MATLAB