function Y=fft(x)
n=size(x);
for l=0:12,
if 2^l>=n
N=2^l;
L=l;
break
end
end
for t=n(2)+1:N, %不足2^L的数在末尾补0
x(t)=0;
end
x %x为补0后的序列
nv2=N/2; %倒位序
nm1=N-1;
I=0;
J=0;
while I<nm1
if I<J
t =x(J+1);
x(J+1)=x(I+1);
x(I+1)=t;
end
K=nv2;
while K<=J
J=J-K;
K=K/2;
end
J=J+K;
I=I+1;
end
X=x
M=1;
while M<=L
LE=2^M;
LE1=LE/2;
U=1;
W=exp(-j*pi/LE1);
J=0;
while J<=LE1-1
I=J;
while I<=N-1
IP=I+LE1;
T=x(IP+1)*U;
x(IP+1)=x(I+1)-T;
x(I+1)=x(I+1)+T;
I=I+LE;
end
U=U*W;
J=J+1;
end
M=M+1;
end
Y=x;
原帖由 tianyoume 于 2007-9-11 15:52 发表
function Y=fft(x)
n=size(x);
for l=0:12,
if 2^l>=n
N=2^l;
L=l;
break
end
end
for t=n(2)+1:N, %不足2^L的数在末尾补0
x(t)=0;
end
x %x ...
%***************具体的FFT的蝶形运算*************************************
for L=1:M
B=2^(L-1);
for J=0:B-1
P=J*2^(M-L);
for k=J:2^L:N-1
xn(k+1)=xn(k+1)+xn(k+1+B)*(exp((-1)*i*2*pi/N*P)); % DIT:先复乘后加减
xn(k+1+B)=xn(k+1)-xn(k+1+B)*(exp((-1)*i*2*pi/N*P));
end
end
end