// - Fast Fourier Transform [64FFT] 
// Univ. of the Ryukyus, Okinawa, Japan

function [y]=myfft64(x);
// internal buffer memory and output memory
x1 = zeros(64,1);
x2 = zeros(64,1);
x3 = zeros(64,1);
y  = zeros(64,1);
// stage 1
for mm = 0:1:15
   //radix-4
   twiddle1=exp(-2*%pi*%i*mm*1/64);
   twiddle2=exp(-2*%pi*%i*mm*2/64);
   twiddle3=exp(-2*%pi*%i*mm*3/64);
   x1(mm+1)  =x(mm+1)  +x(mm+17)+x(mm+33)  +x(mm+49);
   x1(mm+17)=(x(mm+1)-%i*x(mm+17)-x(mm+33)+%i*x(mm+49))*twiddle1;
   x1(mm+33)=(x(mm+1)  -x(mm+17)+x(mm+33)  -x(mm+49))*twiddle2;
   x1(mm+49)=(x(mm+1)+%i*x(mm+17)-x(mm+33)-%i*x(mm+49))*twiddle3;
end;
// stage 2
for nn = 0:16:48
   for mm = 0:1:3
      //radix-4
      twiddle1=exp(-2*%pi*%i*mm*4*1/64);
      twiddle2=exp(-2*%pi*%i*mm*4*2/64);
      twiddle3=exp(-2*%pi*%i*mm*4*3/64);
      x2(mm+nn+1)  =x1(mm+nn+1)  +x1(mm+nn+5)+x1(mm+nn+9)  +x1(mm+nn+13);
      x2(mm+nn+5) =(x1(mm+nn+1)-%i*x1(mm+nn+5)-x1(mm+nn+9)+%i*x1(mm+nn+13))*twiddle1;
      x2(mm+nn+9) =(x1(mm+nn+1)  -x1(mm+nn+5)+x1(mm+nn+9)  -x1(mm+nn+13))*twiddle2;
      x2(mm+nn+13)=(x1(mm+nn+1)+%i*x1(mm+nn+5)-x1(mm+nn+9)-%i*x1(mm+nn+13))*twiddle3;
   end;
end;
// stage 3
for nn = 0:4:60
   x3(nn+1) = x2(nn+1)  +x2(nn+2)+x2(nn+3)  +x2(nn+4);
   x3(nn+2) = x2(nn+1)-%i*x2(nn+2)-x2(nn+3)+%i*x2(nn+4);
   x3(nn+3) = x2(nn+1)  -x2(nn+2)+x2(nn+3)  -x2(nn+4);
   x3(nn+4) = x2(nn+1)+%i*x2(nn+2)-x2(nn+3)-%i*x2(nn+4);
end;
// reorder
for n2 = 0:3
   for n1 = 0:3
      for n0 = 0:3
         y(16*n0+4*n1+n2+1)=x3(16*n2+4*n1+n0+1);
      end;
   end;
end;

endfunction
