《快速傅里叶变换蝶形运算按时间抽取基fft算法MATLAB代码.docx》由会员分享，可在线阅读，更多相关《快速傅里叶变换蝶形运算按时间抽取基fft算法MATLAB代码.docx（4页珍藏版）》请在三一文库上搜索。
1、function y=MyFFT_TB(x,n)%MYFFT_TB:My Fast Fourier Transform Time 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=2Am*M,M 为不能被2整除的整数% two_m - 2Am%
2、 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 -标量，用于临时储存%参考文献：% %输入参数个数检查msg=nargchk(1,2,nargin);error(msg);%输入数据截断或加0N=length(x);if nargin=2if N<n % 加 0xtem=x;x=zeros(1,n);x(1:N)=xtem;N=n;else %截断x
3、tem=x;x=xtem(1:n);N=n;% 对N进行分解 N=2Am*M m,M=factorize(N);two_m=N/M;%变换if m=0% 如果N可以被2整除adr=address(m,M,two_m);y=x(adr+1); %蝶形运算级数 l=m时if M=1 % N 分解含有非2因数M时，对y中每M个数据做直接傅里叶变换 for ii=1:two_my(ii-1)*M+1 : ii*M ) = DDFT( y(ii-1)*M+1 : ii*M );endend%计算W向量W=exp(-2*pi*i* ( 0:N/2-1 ) /N);%蝶形运算d=M;t=N/2/d;for
4、l=m:-1:1%乘for r=0:d-1mul=W(r*t+1);for ii=0:t-1y(ii*2*d+d+1+r) = y(ii*2*d+d+1+r)*mul;endend%加for ii=0:t-1ind1=ii*2*d+1;ind2=ind1+d;for r=0:d-1tem=y(ind1)+y(ind2);y(ind2)=y(ind1)-y(ind2);y(ind1)=tem;ind1=ind1+1;ind2=ind2+1;endendd=2*d;t=t/2;endelse%如果N不能被2整除y=DDFT(x);endend% 内嵌函数=function y=DDFT(x)%直接
5、离散傅里叶变换%input:% x -样本数据，N维向量%output:% y - N维向量%参考文献：% 结构动力学，克拉夫，P82% variable define% s - sum,用于求和N=length(x);y=zeros(size(x);for n=1:Ns=0;for m=1:Ns=s+x(m)*exp( -i*2*pi*(m-1)*(n-1)/N );endy(n)=s;endendfunction m,M=factorize(N)%对N分解m=0;while trueif mod(N,2)=0m=m+1;N=N/2;elseM=N;break;end endendfunct
6、ion adr=address(m,M,two_m)%变址% b - 2Am * m 的矩阵，用来存储二进制数据% ds -数，公差adr=zeros(two_m,M);b=de2bi(0:two_m-1,m);% 转换为 2 进制 注：matlab 中二进制0 1 1=6b=b(:,end:-1:1);% 逆序adr(:,1)=bi2de(b);%2 进制转换为 10 进制if M=1ds=two_m;adr=adr(:,1)*ones(1,M);adr=adr+ds*ones(size(adr,1),1)*(0:M-1);adr=reshape(adr',1,);endend联系方式:matrixsuper