fft_dit_yeshan.m

经典《信号与系统》教程的matlab例程,对深入理解信号与系统相关概念有很大帮助

% Name:FFT_DIT
clear,close all;
fs=12800;
N=256;
M=log2(N);
n=0:N-1;
x=input('Type in the sequence x(n)=');%sin(pi*n/128)+sin(2*pi*n/128)+sin(3*pi*n/128)+sin(4*pi*n/128)
%==========================================
% Bit-Reversed Order
%==========================================
k=1:M-1;
x1=zeros(1,N);
x1(1)=x(1);
x1(N)=x(N);
J=0;
for I=1:N-2;
    k=1;s=0;
    while k~=M+1&s~=1;
            if J>=N/(2.^k);
                J=J-N/(2.^k);
            else J=J+N/(2.^k);s=1;                
            end
            k=k+1;            
        end        
     x1(I+1)=x(J+1);
end
%=============================================
% Computation of the complex number WN^p
%=============================================
m=0:N/2-1;
W=exp(-i*2*pi*m/N);
%=============================================
% Butterfly computation
%=============================================
for L=1:M
    B=2^(L-1);
    for J=0:B-1
        p=2^(M-L)*J;
        for k=J+1:2^L:N-1;
            q=x1(k);
            x1(k)=q+x1(k+B)*W(p+1);
            x1(k+B)=q-x1(k+B)*W(p+1);
        end              
    end
end
x1=2*x1/N;
k=0:N-1;
subplot(211)
plot(k,x),title('The original sequence x(n)')
xlabel('Index n'),axis([0,max(k),min(x),1.1*max(x)])
subplot(212)
stem(k,abs(x1),'r.'),title('The DFT of the original sequence x(n)')
xlabel('Index k'),axis([0,N/2,min(abs(x1)),1.1*max(abs(x1))])