% Name:FFTp
clear
N=64;%input('Type in the DFT length N=');
M=log2(N);
n=0:N-1;
x=u(n)-u(n-5);%input('Type in the sequence x(n)=');
%==========================================
% 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
x1;   
%=============================================
% Computation of the complex number WN^p
%=============================================
m=0:N/2-1;
W=exp(-i*2*pi*m/N);
%=============================================
%
%=============================================
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);
            if L~=M|k~=N/2
                x1(k+B)=q-x1(k+B)*W(p+1);
            end
        end
              
    end
end
for k=1:N/2
    x1(N-k+1)=conj(x1(k+1));
end
k=0:N-1;
subplot(211)
stem(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,max(k),min(abs(x1)),1.1*max(abs(x1))])