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.




 
    P =8
    S  =int( log( P) /log( 2) +0.9999)
 
    Pi =3.14159265
    R1 =2^S
 
    R =R1 -1
    R2 =div( R1,  2)
    R4 =div( R1,  4)
    R3 =R4 +R2
 
    Dim Re( R1), Im( R1), Co( R3)
 
    for N =0 to P -1
        read dummy: Re( N) =dummy
        read dummy: Im( N) =dummy
    next N
 
    data    1, 0,      1, 0,      1, 0,      1, 0,      0, 0,     0, 0,      0, 0,       0, 0
 
    S2 =div( S, 2)
    S1 =S -S2
    P1 =2^S1
    P2 =2^S2
 
    dim V( P1 -1)
    V( 0) =0
    DV =1
    DP =P1
 
    for J =1 to S1
        HA =div( DP, 2)
        PT =P1 -HA
        for I =HA to PT step DP
            V( I) =V( I -HA) +DV
        next I
        DV =DV +DV
        DP =HA
    next J
 
    K =2 *Pi /R1
 
    for X =0 to R4
        COX =cos( K *X)
        Co( X) =COX
        Co( R2 -X) =0 -COX
        Co( R2 +X) =0 -COX
    next X
 
    print "FFT: bit reversal"
 
    for I =0 to P1 -1
        IP =I *P2
        for J =0 to P2 -1
            H =IP +J
            G =V( J) *P2 +V( I)
            if G >H then temp =Re( G): Re( G) =Re( H): Re( H) =temp
            if G >H then temp =Im( G): Im( G) =Im( H): Im( H) =temp
        next J
    next I
 
    T =1
 
    for stage =0 to S -1
        print "  Stage:- "; stage
        D =div( R2, T)
        for Z =0 to T -1
            L   =D *Z
            LS  =L +R4
            for I =0 to D -1
                A      =2 *I *T +Z
                B      =A +T
                F1     =Re( A)
                F2     =Im( A)
                P1     =Co( L)  *Re( B)
                P2     =Co( LS) *Im( B)
                P3     =Co( LS) *Re( B)
                P4     =Co( L)  *Im( B)
                Re( A) =F1 +P1 -P2
                Im( A) =F2 +P3 +P4
                Re( B) =F1 -P1 +P2
                Im( B) =F2 -P3 -P4
            next I
        next Z
        T =T +T
    next stage
 
    print "   M          Re( M)       Im( M)"
 
    for M =0 to R
        if abs( Re( M)) <10^-5 then Re( M) =0
        if abs( Im( M)) <10^-5 then Im( M) =0
        print "   "; M, Re( M), Im( M)
    next M
 
    end
 
 
    wait
 
    function div( a, b)
        div =int( a /b)
    end function
 
    end
 