快速傅里叶变换_蝶形运算_按频率抽取基2-fft算法_MATLAB代码.pdf

《快速傅里叶变换_蝶形运算_按频率抽取基2-fft算法_MATLAB代码.pdf》由会员分享，可在线阅读，更多相关《快速傅里叶变换_蝶形运算_按频率抽取基2-fft算法_MATLAB代码.pdf（4页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、function y=MyFFT_FB(x,n)%MYFFT_TB:My Fast Fourier Transform Frequency Based%按频率抽取基 2-fft 算法%input:%x-输入的一维样本%n-变换长度，缺省时 n=length(x)当 n 小于 x 数据长度时，x 数据被截断到第 n 个数据%当 n 大于时，x 数据在尾部补 0 直到 x 含 n 个数据%output:%y-1*n 的向量，快速傅里叶变换结果%variable define:%N-一维数据 x 的长度%xtem-临时储存 x 数据用%m,M-对 N 进行分解 N=2m*M,M 为不能被 2 整除的

2、整数%two_m-2m%adr-变址，1*N 的向量%l-当前蝶形运算的级数%W-长为 N/2 的向量，记录 W(0,N),W(1,N),.W(N/2-1,N)%d-蝶形运算两点间距离%t-第 l 级蝶形运算含有的奇偶数组的个数%mul-标量，乘数%ind1,ind2-标量，下标%tem-标量，用于临时储存%参考文献：%81c 输入参数个数检查 msg=nargchk(1,2,nargin);error(msg);%输入数据截断或加 0 N=length(x);if nargin=2 if Nn%加 0 xtem=x;x=zeros(1,n);x(1:N)=xtem;N=n;else%截断 x

3、tem=x;x=xtem(1:n);N=n;end end%对 N 进行分解 N=2m*M m,M=factorize(N);two_m=N/M;%变换 if m=0%如果 N 可以被 2 整除 adr=address(m,M,two_m);y=x;%蝶形运算级数 l=m 时%计算 W 向量 W=exp(-2*pi*i*(0:N/2-1)/N);%蝶形运算 d=N/2;t=1;for l=1:m%加 for ii=0:t-1 ind1=ii*2*d+1;ind2=ind1+d;for r=0:d-1 tem=y(ind1)+y(ind2);y(ind2)=y(ind1)-y(ind2);y(i

4、nd1)=tem;ind1=ind1+1;ind2=ind2+1;end end%乘 for r=0:d-1 mul=W(r*t+1);for ii=0:t-1 y(ii*2*d+d+1+r)=y(ii*2*d+d+1+r)*mul;end end d=d/2;t=t*2;end%直接傅里叶变换 if M=1%N 分解含有非 2 因数 M 时，对 y 中每 M 个数据做直接傅里叶变换 for ii=1:two_m y(ii-1)*M+1:ii*M)=DDFT(y(ii-1)*M+1:ii*M);end end%变址输出 y=y(adr+1);else%如果 N 不能被 2 整除 y=DDFT(

5、x);end end%内嵌函数=function y=DDFT(x)%直接离散傅里叶变换%input:%x-样本数据，N 维向量%output:%y-N 维向量%参考文献：%结构动力学，克拉夫，P82%variable define%s-sum，用于求和 N=length(x);y=zeros(size(x);for n=1:N s=0;for m=1:N s=s+x(m)*exp(-i*2*pi*(m-1)*(n-1)/N);end y(n)=s;end end function m,M=factorize(N)%对 N 分解 m=0;while true if mod(N,2)=0 m=m

6、+1;N=N/2;else M=N;break;end end end function adr=address(m,M,two_m)%变址%b-2m*m 的矩阵，用来存储二进制数据%ds-数，公差 adr=zeros(two_m,M);b=de2bi(0:two_m-1,m);%转换为 2 进制 注：matlab 中二进制0 1 1=6 b=b(:,end:-1:1);%逆序 adr(:,1)=bi2de(b);%2 进制转换为 10 进制 if M=1 ds=two_m;adr=adr(:,1)*ones(1,M);adr=adr+ds*ones(size(adr,1),1)*(0:M-1);adr=reshape(adr,1,);end end