fft64.m

64-point radix-2 fft module

% 256-points fft  (DIF radix-2)  only support 256-fft
clc;clear all;
N = 256;
x =[randint(1,N)] ;
reg1 = [];
reg2 = [];
reg3 = [];
reg4 = [];
reg5 = [];
reg6 = [];
reg7 = [];
reg8 = [];
output = [];

%----stage1-----
for i = 1:N/2
    reg1(i)    = x(i) + x(i+N/2) ;
    reg1(i+N/2) = (x(i) - x(i+N/2))*exp(-j*2*pi*(i-1)/N);
end;
%----stage2----- 
for i = 1:N/4
    for k = 0:N/2:N-1
        reg2(i+k)    = reg1(i+k) + reg1(i+N/4+k);
        reg2(i+N/4+k)= (reg1(i+k) - reg1(i+N/4+k))*exp(-j*2*pi*2*(i-1)/N);
    end;
end;

%----stage3-----
for i = 1:N/8
    for k = 0:N/4:N-1;
        reg3(i+k)    = reg2(i+k) + reg2(i+N/8+k); 
        reg3(i+N/8+k)= (reg2(i+k) - reg2(i+N/8+k))*exp(-j*2*pi*4*(i-1)/N);
    end
end;
%----stage4-----
for i = 1:N/16
    for k = 0:N/8:N-1;
        reg4(i+k)     = reg3(i+k) + reg3(i+N/16+k); 
        reg4(i+N/16+k)= (reg3(i+k) - reg3(i+N/16+k))*exp(-j*2*pi*8*(i-1)/N);
    end
end;    
%----stage5-----
for i = 1:N/32
    for k = 0:N/16:N-1;
        reg5(i+k)    = reg4(i+k) + reg4(i+N/32+k); 
        reg5(i+N/32+k)  = (reg4(i+k) - reg4(i+N/32+k))*exp(-j*2*pi*N/16*(i-1)/N);
    end
end;  
%----stage6-----
for i = 1:N/64
    for k = 0:N/32:N-1;
        reg6(i+k)    = reg5(i+k) + reg5(i+N/64+k); 
        reg6(i+N/64+k)  = (reg5(i+k) - reg5(i+N/64+k))*exp(-j*2*pi*N/8*(i-1)/N);
    end
end;   
%----stage7-----
for i = 1:N/128
    for k = 0:N/64:N-1;
        reg7(i+k)    = reg6(i+k) + reg6(i+N/128+k); 
        reg7(i+N/128+k)  = (reg6(i+k) - reg6(i+N/128+k))*exp(-j*2*pi*N/4*(i-1)/N);
    end
end;   
%----stage8-----
for k=1:2:N,
reg8(k)=reg7(k)+reg7(k+1);
end;
for k=2:2:N,
reg8(k)=(reg7(k-1)-reg7(k));
end;
output = reg8;
%----bit reversal---
num=0:N-1;
temp=dec2bin(0:N-1);
z=bin2dec(temp(:,end:-1:1))';
output=output([z+1]);
%-----check-----
figure;
stem(real(fft(x,N)),N);
hold on;
stem(real(output),'rx',N);
legend('func','matlab')