% - Fast Fourier Transform [1024FFT] - (Tue, 7 Nov 2000)
% Univ. of the Ryukyus, Okinawa, Japan

function [y]=myfft1024(x);
x1 = zeros(1024,1);
x2 = zeros(1024,1);
x3 = zeros(1024,1);
x4 = zeros(1024,1);
x5 = zeros(1024,1);
y  = zeros(1024,1);
% stage 1
for mm = 0:1:255
   %radix-4
  %twiddle0=exp(-2*pi*j*mm*0/1024);
   twiddle1=exp(-2*pi*j*mm*1/1024);
   twiddle2=exp(-2*pi*j*mm*2/1024);
   twiddle3=exp(-2*pi*j*mm*3/1024);
   x1(mm+1)  = x(mm+1)  +x(mm+257)+x(mm+513)  +x(mm+769);%twiddle0;
   x1(mm+257)=(x(mm+1)-j*x(mm+257)-x(mm+513)+j*x(mm+769))*twiddle1;
   x1(mm+513)=(x(mm+1)  -x(mm+257)+x(mm+513)  -x(mm+769))*twiddle2;
   x1(mm+769)=(x(mm+1)+j*x(mm+257)-x(mm+513)-j*x(mm+769))*twiddle3;
end;
% stage 2
for nn = 0:256:768
   for mm = 0:1:63
      %radix-4
     %twiddle0=exp(-2*pi*j*mm*4*0/1024);
      twiddle1=exp(-2*pi*j*mm*4*1/1024);
      twiddle2=exp(-2*pi*j*mm*4*2/1024);
      twiddle3=exp(-2*pi*j*mm*4*3/1024);
      x2(mm+nn+1)  = x1(mm+nn+1)  +x1(mm+nn+65)+x1(mm+nn+129)  +x1(mm+nn+193);%twiddle0;
      x2(mm+nn+65) =(x1(mm+nn+1)-j*x1(mm+nn+65)-x1(mm+nn+129)+j*x1(mm+nn+193))*twiddle1;
      x2(mm+nn+129)=(x1(mm+nn+1)  -x1(mm+nn+65)+x1(mm+nn+129)  -x1(mm+nn+193))*twiddle2;
      x2(mm+nn+193)=(x1(mm+nn+1)+j*x1(mm+nn+65)-x1(mm+nn+129)-j*x1(mm+nn+193))*twiddle3;
   end;
end;
% stage 3
for nn = 0:64:960
   for mm = 0:1:15
      %radix-4
     %twiddle0=exp(-2*pi*j*mm*16*0/1024);
      twiddle1=exp(-2*pi*j*mm*16*1/1024);
      twiddle2=exp(-2*pi*j*mm*16*2/1024);
      twiddle3=exp(-2*pi*j*mm*16*3/1024);
      x3(mm+nn+1) = x2(mm+nn+1)  +x2(mm+nn+17)+x2(mm+nn+33)  +x2(mm+nn+49);%twiddle0;
      x3(mm+nn+17)=(x2(mm+nn+1)-j*x2(mm+nn+17)-x2(mm+nn+33)+j*x2(mm+nn+49))*twiddle1;
      x3(mm+nn+33)=(x2(mm+nn+1)  -x2(mm+nn+17)+x2(mm+nn+33)  -x2(mm+nn+49))*twiddle2;
      x3(mm+nn+49)=(x2(mm+nn+1)+j*x2(mm+nn+17)-x2(mm+nn+33)-j*x2(mm+nn+49))*twiddle3;
   end;
end;
% stage 4
for nn = 0:16:1008
   for mm = 0:1:3
      %radix-4
     %twiddle0=exp(-2*pi*j*mm*64*0/1024);
      twiddle1=exp(-2*pi*j*mm*64*1/1024);
      twiddle2=exp(-2*pi*j*mm*64*2/1024);
      twiddle3=exp(-2*pi*j*mm*64*3/1024);
      x4(mm+nn+1) = x3(mm+nn+1)  +x3(mm+nn+5)+x3(mm+nn+9)  +x3(mm+nn+13);%twiddle0;
      x4(mm+nn+5) =(x3(mm+nn+1)-j*x3(mm+nn+5)-x3(mm+nn+9)+j*x3(mm+nn+13))*twiddle1;
      x4(mm+nn+9) =(x3(mm+nn+1)  -x3(mm+nn+5)+x3(mm+nn+9)  -x3(mm+nn+13))*twiddle2;
      x4(mm+nn+13)=(x3(mm+nn+1)+j*x3(mm+nn+5)-x3(mm+nn+9)-j*x3(mm+nn+13))*twiddle3;
   end;
end;
% stage 5
for nn = 0:4:1020
   x5(nn+1) = x4(nn+1)  +x4(nn+2)+x4(nn+3)  +x4(nn+4);
   x5(nn+2) = x4(nn+1)-j*x4(nn+2)-x4(nn+3)+j*x4(nn+4);
   x5(nn+3) = x4(nn+1)  -x4(nn+2)+x4(nn+3)  -x4(nn+4);
   x5(nn+4) = x4(nn+1)+j*x4(nn+2)-x4(nn+3)-j*x4(nn+4);
end;
% reorder
for n4 = 0:3
   for n3 = 0:3
      for n2 = 0:3
         for n1 = 0:3
            for n0 = 0:3
               y(256*n0+64*n1+16*n2+4*n3+n4+1)=x5(256*n4+64*n3+16*n2+4*n1+n0+1);
            end;
         end;
      end;
   end;
end;