快速傅里叶变换_蝶形运算_按时间抽取基2-fft算法_MATLAB代码

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=2^m*M,M为不能被2整除的整数
% two_m -- 2^m
% 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 -- 标量，用于临时储存
%参考文献：
% /view/fea1e985b9d528ea81c779ee.html
%% 输入参数个数检查
msg=nargchk(1,2,nargin);
error(msg);
%% 输入数据截断或加0
N=length(x);
if nargin==2
if N<n % 加0
xtem=x;
x=zeros(1,n);
x(1:N)=xtem;
N=n;
else % 截断
xtem=x;
x=xtem(1:n);
N=n;
end
end
%% 对N进行分解N=2^m*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_m
y((ii-1)*M+1 : ii*M ) = DDFT( y((ii-1)*M+1 : ii*M ) );
end
end
%% 计算W向量
W=exp(-2*pi*i* ( 0:N/2-1 ) /N);
%% 蝶形运算
d=M;
t=N/2/d;
for l=m:-1:1
% 乘
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
% 加
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(ind1)=tem;
ind1=ind1+1;
ind2=ind2+1;
end
end
d=2*d;t=t/2;
end
else
%% 如果N 不能被2整除
y=DDFT(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+1;
N=N/2;
else
M=N;
break;
end
end
end
function adr=address(m,M,two_m)
%% 变址
% b -- 2^m * 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
联系方式：matrixsuper@。