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.




cnt  = 8
sig  = int(log(cnt) /log(2) +0.9999)
 
pi    = 3.14159265
real1 = 2^sig
 
real  = real1 -1
real2 = int(real1 /  2)
real4 = int(real1 /  4)
real3 = real4 +real2
 
dim rel(real1)
dim img(real1)
dim cmp(real3)
 
for i = 0 to cnt -1
    read rel(i)
    read img(i)
next i
 
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
    hlfPtr = int(ptr / 2)
    pt     = cnt1 - hlfPtr
    for i = hlfPtr to pt step ptr
        v(i) = v(i -hlfPtr) + dv
    next i
    dv = dv + dv
    ptr = hlfPtr
next j
 
k = 2 *pi /real1
 
for x = 0 to real4
    cmp(x)         = cos(k *x)
    cmp(real2 - x) = 0 - cmp(x)
    cmp(real2 + x) = 0 - cmp(x)
next x
 
print "fft: bit reversal"
 
for i = 0 to cnt1 -1
    ip = i *cnt2
    for j = 0 to cnt2 -1
        h = ip +j
        g = v(j) *cnt2 +v(i)
        if g >h then 
                temp   = rel(g)
                rel(g) = rel(h)
                rel(h) = temp
                temp   = img(g)
                img(g) = img(h)
                img(h) = temp
         end if
    next j
next i
 
t = 1
for stage = 1 to sig
    print "  stage:- "; stage
    d = int(real2 / t)
    for ii = 0 to t -1
        l   = d *ii
        ls  = l +real4
        for i = 0 to d -1
            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
        next i
    next ii
    t = t +t
next stage
 
print "  Num   real   imag"
for i = 0 to real
    if abs(rel(i)) <10^-5 then rel(i) = 0
    if abs(img(i)) <10^-5 then img(i) = 0
    print "   "; i;"   ";using("##.#",rel(i));"    ";img(i)
next i
end