Task

Calculate the   FFT   (Fast Fourier Transform)   of an input sequence.

The most general case allows for complex numbers at the input
and results in a sequence of equal length, again of complex numbers.
If you need to restrict yourself to real numbers, the output should
be the magnitude   (i.e.:   sqrt(re2 + im2))   of the complex result.

The classic version is the recursive Cooley–Tukey FFT. Wikipedia has pseudo-code for that.
Further optimizations are possible but not required.




test: PROCEDURE OPTIONS (MAIN, REORDER); /* Derived from Fortran Rosetta Code */
 
  /* In-place Cooley-Tukey FFT */
FFT: PROCEDURE (x) RECURSIVE;
   DECLARE  x(*) COMPLEX FLOAT (18);
   DECLARE  t    COMPLEX FLOAT (18);
   DECLARE ( N, Half_N ) FIXED BINARY (31);
   DECLARE ( i, j ) FIXED BINARY (31);
   DECLARE (even(*), odd(*)) CONTROLLED COMPLEX FLOAT (18);
   DECLARE pi FLOAT (18) STATIC INITIAL ( 3.14159265358979323E0);
 
   N = HBOUND(x);
 
   if N <= 1 THEN return;
 
   allocate odd((N+1)/2), even(N/2);
 
    /* divide */
   do j = 1 to N by 2; odd((j+1)/2) = x(j); end;
   do j = 2 to N by 2; even(j/2)    = x(j); end;
 
    /* conquer */
   call fft(odd);
   call fft(even);
 
    /* combine */
   half_N = N/2;
   do i=1 TO half_N;
      t = exp(COMPLEX(0, -2*pi*(i-1)/N))*even(i);
      x(i)        = odd(i) + t;
      x(i+half_N) = odd(i) - t;
   end;
 
   FREE odd, even;
 
END fft;
 
 
   DECLARE data(8)  COMPLEX FLOAT (18) STATIC INITIAL (
                    1, 1, 1, 1, 0, 0, 0, 0);
   DECLARE ( i ) FIXED BINARY (31);
 
   call fft(data);
 
   do i=1 TO 8;
      PUT SKIP LIST ( fixed(data(i), 25, 12) );
   end;
 
END test;