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.




 
PROGRAM FFT
 
CONST CNT=8
 
!$DYNAMIC
DIM REL[0],IMG[0],CMP[0],V[0]
 
BEGIN
   SIG=INT(LOG(CNT)/LOG(2)+0.9999)
   REAL1=2^SIG
 
   REAL=REAL1-1
   REAL2=INT(REAL1/2)
   REAL4=INT(REAL1/4)
   REAL3=REAL4+REAL2
 
!$DIM REL[REAL1],IMG[REAL1],CMP[REAL3]
 
FOR I=0 TO CNT-1 DO
   READ(REL[I],IMG[I])
END FOR
 
DATA(1,0,1,0,1,0,1,0,0,0,0,0,0,0,0,0)
 
SIG2=INT(SIG/2)
SIG1=SIG-SIG2
CNT1=2^SIG1
CNT2=2^SIG2
 
!$DIM V[CNT1-1]
V[0]=0
DV=1
PTR=CNT1
 
FOR J=1 TO SIG1 DO
  HLFPTR=INT(PTR/2)
  PT=CNT1-HLFPTR
  FOR I=HLFPTR TO PT STEP PTR DO
    V[I]=V[I-HLFPTR]+DV
  END FOR
  DV=2*DV
  PTR=HLFPTR
END FOR
 
K=2*π/REAL1
 
FOR X=0 TO REAL4 DO
   CMP[X]=COS(K*X)
   CMP[REAL2-X]=-CMP[X]
   CMP[REAL2+X]=-CMP[X]
END FOR
 
PRINT("FFT: BIT REVERSAL")
 
FOR I=0 TO CNT1-1 DO
  IP=I*CNT2
  FOR J=0 TO CNT2-1 DO
    H=IP+J
    G=V[J]*CNT2+V[I]
    IF G>H THEN
       SWAP(REL[G],REL[H])
       SWAP(IMG[G],IMG[H])
    END IF
  END FOR
END FOR
 
T=1
FOR STAGE=1 TO SIG DO
  PRINT("STAGE:";STAGE)
  D=INT(REAL2/T)
  FOR II=0 TO T-1 DO
     L=D*II
     LS=L+REAL4
     FOR I=0 TO D-1 DO
       A=2*I*T+II
       B=A+T
       F1=REL[A]
       F2=IMG[A]
       CNT1=CMP[L]*REL[B]
       CNT2=CMP[LS]*IMG[B]
       CNT3=CMP[LS]*REL[B]
       CNT4=CMP[L]*IMG[B]
       REL[A]=F1+CNT1-CNT2
       IMG[A]=F2+CNT3+CNT4
       REL[B]=F1-CNT1+CNT2
       IMG[B]=F2-CNT3-CNT4
     END FOR
  END FOR
  T=2*T
END FOR
 
PRINT("NUM REAL     IMAG")
FOR I=0 TO REAL DO
    IF ABS(REL[I])<1E-5 THEN REL[I]=0 END IF
    IF ABS(IMG[I])<1E-5 THEN IMG[I]=0 END IF
    PRINT(I;"";)
    WRITE("##.###### ##.######";REL[I];IMG[I])
END FOR
END PROGRAM
 