% Name:FFT_DIT
clear,close all;
N=input('Type in the DFT length N=');
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;c=0;
    while k~=M+1&c~=1;
            if J>=N/(2.^k);
                J=J-N/(2.^k);c=0;
            else J=J+N/(2.^k);c=1;                
            end
            k=k+1;            
        end        
     x1(I+1)=x(J+1);
end
%=============================================
% Compute 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)
stem(k,x,'.'),title('The sequence x(n)')
xlabel('Index n'),axis([0,max(k),1.1*min(x),1.1*max(x)])
subplot(212)
stem(k,abs(x1),'r.'),title('The DFT of the sequence x(n)')
xlabel('Index k'),axis([0,max(k),1.1*min(abs(x1)),1.1*max(abs(x1))])