学步园

以前从没手动实现过FFT算法,实现了一下.

参考<computational frameworks for the fast fourier transform> 

% 2012.3.10
function test
clear;
clc;
close all;
nFFT=2^10;
N=nFFT;
t=linspace(0,1,N)';
fs=N;
x=sin(2*pi*300*t);


freq=fs*((1:nFFT)-1)/nFFT;  % bug fixed 1:nFFT
figure;
disp('time[Matlab fft]:');
tic
y0=fft(x,nFFT);
toc
plot(freq,abs(y0));
title('fft');

disp('time[recursive]:');
tic
y=my_fft(x,nFFT);
toc
figure;
plot(freq,abs(y));
title('my-fft');

figure(3);
plot(freq, abs(y-y0));
title('my-fft-error');

disp('time[inplace]:');
tic
y_inplace=my_ffy_inplace(x,nFFT);
toc
figure(4);
plot(freq,abs(y_inplace-y0));
title('inplace-error');
figure(5);
plot(freq,abs(y_inplace));

    function x=my_ffy_inplace(x,n)
        J=zeros(n,1);
        for k=1:n
            J(k)=bit_reverse(n,k-1);
        end
        J=J+1;  % 1-base
        x=x(J);     
        t=log2(n);
        for q=1:t
            L=2^q; r=n/L; L_2=L/2;
            for j=0:L_2-1
                w=cos(2*pi*j/L)-1i*sin(2*pi*j/L);
                for k=0:r-1
                    tau=w*x((k)*L+j+1+L_2);
                    x((k)*L+j+1+L_2)=x((k)*L+j+1)-tau;
                    x((k)*L+j+1)=    x((k)*L+j+1)+tau;
                end
            end
        end
    end


    function j=bit_reverse(n,k)
        % 0-base index
        t=log2(n);
        j=0;
        for q=1:t
            s=floor(k/2);
            b_q=k-2*s;
            % horner
            j=2*j+b_q;
            k=s;
        end
    end
        

    function y=my_fft(x,n)
        if(n==1)
            y=x;
        else
            m=n/2;
            w=exp(-2*pi*1i/n);
            Omega=zeros(m,m);
            for j=1:m
                Omega(j,j)=w^(j-1);
            end
            Zr=my_fft(x(1:2:n),m);
            Zb=Omega*my_fft(x(2:2:n),m);
            Im=eye(m);
            y=[Im, Im; Im,-Im]*[Zr;Zb];
        end
    end
end