--
-- demo\rosetta\FastFourierTransform.exw
-- =====================================
--
--  Originally written by Robert Craig and posted to EuForum Dec 13, 2001
--

constant REAL = 1, IMAG = 2

type complex(sequence x)
    return length(x)=2 and atom(x[REAL]) and atom(x[IMAG])
end type

function p2round(integer x)
-- rounds x up to a power of two
integer p = 1
    while p<x do
        p += p
    end while
    return p
end function

function log2(atom x)
-- return log2 of x, or -1 if x is not a power of 2
    if x>0 then
        integer p = -1
        while floor(x)=x do
            x /= 2
            p += 1
        end while
        if x=0.5 then
            return p
        end if
    end if
    return -1
end function

function bitrev(sequence a)
-- bitrev an array of complex numbers
integer j=1, n = length(a)
    for i=1 to n-1 do
        if i<j then
            {a[i],a[j]} = {a[j],a[i]}
        end if
        integer k = n/2
        while k<j do
            j -= k
            k /= 2
        end while
        j = j+k
    end for
    return a
end function

function cmult(complex arg1, complex arg2)
-- complex multiply
    return {arg1[REAL]*arg2[REAL]-arg1[IMAG]*arg2[IMAG],
            arg1[REAL]*arg2[IMAG]+arg1[IMAG]*arg2[REAL]}
end function

function ip_fft(sequence a)
-- perform an in-place fft on an array of complex numbers
-- that has already been bit reversed
integer n = length(a)
integer ip, le, le1
complex u, w, t

    for l=1 to log2(n) do
        le = power(2, l)
        le1 = le/2
        u = {1, 0}
        w = {cos(PI/le1), sin(PI/le1)}
        for j=1 to le1 do
            for i=j to n by le do
                ip = i+le1
                t = cmult(a[ip], u)
                a[ip] = sq_sub(a[i],t)
                a[i] = sq_add(a[i],t)
            end for
            u = cmult(u, w)
        end for
    end for
    return a
end function

function fft(sequence a)
integer n = length(a)
    if log2(n)=-1 then
        puts(1, "input vector length is not a power of two, padded with 0's\n\n")
        n = p2round(n)
         -- pad with 0's
        for j=length(a)+1 to n do
            a = append(a,{0, 0})
        end for
    end if
    a = ip_fft(bitrev(a))
    -- reverse output from fft to switch +ve and -ve frequencies
    for i=2 to n/2 do
        integer j = n+2-i
        {a[i],a[j]} = {a[j],a[i]}
    end for
    return a
end function

function ifft(sequence a)
integer n = length(a)
    if log2(n)=-1 then ?9/0 end if -- (or as above?)
    a = ip_fft(bitrev(a))
    -- modifies results to get inverse fft
    for i=1 to n do
        a[i] = sq_div(a[i],n)
    end for
    return a
end function

constant a = {{1, 0},
              {1, 0},
              {1, 0},
              {1, 0},
              {0, 0},
              {0, 0},
              {0, 0},
              {0, 0}}

printf(1, "Results of %d-point fft:\n\n", length(a))
ppOpt({pp_Nest,1,pp_IntFmt,"%10.6f",pp_FltFmt,"%10.6f"})
pp(fft(a))
printf(1, "\nResults of %d-point inverse fft (rounded to 6 d.p.):\n\n", length(a))
pp(ifft(fft(a)))