自己写的FFT总是得不出正确结果。

tianyoume

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;

eight

原帖由 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 函数不行吗？

zhangnan3509

现在好像流行自己写。

tianyoume

要提取循环平稳频率，需要自己写FFT

tianyoume

上面那段代码是转自别人的，但我试了一下好象与MATLAB带的FFT有出入，不知道源创者是否注意到了。

graduate

%***************具体的ＦＦＴ的蝶形运算*************************************
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