自己写的FFT总是得不出正确结果。
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; 本帖最后由 wdhd 于 2016-9-8 14:18 编辑
原帖由 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 ...
用 matlab 自带的 fft 函数不行吗?
回复 #2 eight 的帖子
现在好像流行自己写。 要提取循环平稳频率,需要自己写FFT 上面那段代码是转自别人的,但我试了一下好象与MATLAB带的FFT有出入,不知道源创者是否注意到了。 %***************具体的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
页:
[1]