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.




      @% = &60A
 
      DIM Complex{r#, i#}
      DIM in{(7)} = Complex{}, out{(7)} = Complex{}
      DATA 1, 1, 1, 1, 0, 0, 0, 0
 
      PRINT "Input (real, imag):"
      FOR I% = 0 TO 7
        READ in{(I%)}.r#
        PRINT in{(I%)}.r# "," in{(I%)}.i#
      NEXT
 
      PROCfft(out{()}, in{()}, 0, 1, DIM(in{()},1)+1)
 
      PRINT "Output (real, imag):"
      FOR I% = 0 TO 7
        PRINT out{(I%)}.r# "," out{(I%)}.i#
      NEXT
      END
 
      DEF PROCfft(b{()}, o{()}, B%, S%, N%)
      LOCAL I%, t{} : DIM t{} = Complex{}
      IF S% < N% THEN
        PROCfft(o{()}, b{()}, B%, S%*2, N%)
        PROCfft(o{()}, b{()}, B%+S%, S%*2, N%)
        FOR I% = 0 TO N%-1 STEP 2*S%
          t.r# = COS(-PI*I%/N%)
          t.i# = SIN(-PI*I%/N%)
          PROCcmul(t{}, o{(B%+I%+S%)})
          b{(B%+I% DIV 2)}.r# = o{(B%+I%)}.r# + t.r#
          b{(B%+I% DIV 2)}.i# = o{(B%+I%)}.i# + t.i#
          b{(B%+(I%+N%) DIV 2)}.r# = o{(B%+I%)}.r# - t.r#
          b{(B%+(I%+N%) DIV 2)}.i# = o{(B%+I%)}.i# - t.i#
        NEXT
      ENDIF
      ENDPROC
 
      DEF PROCcmul(c{},d{})
      LOCAL r#, i#
      r# = c.r#*d.r# - c.i#*d.i#
      i# = c.r#*d.i# + c.i#*d.r#
      c.r# = r#
      c.i# = i#
      ENDPROC
 