OpenEuphoria

Re: FFT algorithm

• Posted by Robert Craig <rds at RapidEuphoria.com> Dec 13, 2001
• 491 views

Michael J. Sabal writes:
> Does anybody have a pre-built FFT (Fast Fourier Transform) algorithm
> already built in Euphoria?

I wrote the following code many years ago,
back in the pre-1.0 days of Euphoria.
It was one of the first Euphoria programs I wrote.
I used to be heavily involved in writing a compiler for
a machine that was designed to do FFTs and litte else,
but I've forgotten all that stuff now. I hope the code below 
does approximately what you want.

    ----------------------------------
    --    FAST FOURIER TRANSFORM    --
    ----------------------------------

constant SCREEN = 1,
      pi = 3.141592653589793,
    REAL = 1, IMAG = 2

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

global function p2round(integer x)
-- rounds x up to a power of two
    integer p

    p = 1
    while p < x do
         p = p + p
    end while
    return p
end function


global function log2(atom x)
-- return log2 of x, or -1 if x is not a power of 2

    integer p

    if x < 1 then
         return -1 
    end if
    p = -1
    while floor(x) = x do
         x = x / 2
         p = p + 1
    end while
    if x = 0.5 then
         return p
    else
         return -1
    end if 
end function

function input()
-- input all the complex numbers

    sequence a
    integer n

    -- read in all the numbers, for now...
    a = {{1.0, 1.1},
          {1.0, 1.9},
          {1.0, 1.0},
          {4.1, 1.0},
          {1.0, 1.0},
          {1.0, 6.5},
          {1.5, 1.0},
          {1.0, 8.0}}
    
    n = length(a)
    if log2(n) = -1 then
         puts(SCREEN, "input vector length is not a power of two\n") 
         puts(SCREEN, "--> padded with 0's\n\n")
         n = p2round(n)
         -- pad with 0's 
         for j = length(a)+1 to n do
             a[j] = {0, 0}
         end for
    end if
    return a
end function

function bitrev(sequence a)
-- bitrev an array of complex numbers

    atom nv2
    integer nm1, j, k, n
    complex temp

    n = length(a)
    nv2 = n / 2
    nm1 = n - 1
    j = 1
    for i = 1 to nm1 do
         if i < j then
             temp = a[i]
             a[i] = a[j]
             a[j] = temp
         end if
         k = nv2
         while k < j do
             j = j - k
             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 fft(sequence a)
-- perform an in-place fft on an array of complex numbers
-- that has already been bit reversed

    integer m, ip, le
    integer n
    integer le1
    complex u, w, t

    n = length(a)
    m = log2(n)
    for l = 1 to m 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] = a[i] - t
                  a[i]  = a[i] + t
             end for
             u = cmult(u, w)
         end for
    end for
    return a
end function


function frequency_reverse(sequence a)
-- reverse output from fft to switch +ve and -ve frequencies

    integer n
    sequence temp

    n = length(a)
    for i = 2 to n / 2 do
         temp = a[i]
         a[i] = a[n + 2 - i]
         a[n + 2 - i] = temp
    end for
    return a
end function


procedure output(sequence a)
-- output an array of complex numbers

    for i = 1 to length(a) do
         print(SCREEN, a[i])
    end for
end procedure


function make_inverse(sequence a)
-- modifies results to get inverse fft

    for i = 1 to length(a) do
         a[i] = a[i] / length(a) 
    end for
    return a
end function


global procedure main(integer argv)
-- FFT / IFFT 
    sequence a

    a = fft(bitrev(input()))

    if argv = 1 then
         a = make_inverse(a)
         printf(SCREEN, "Results of %d-point inverse fft:\n\n", length(a))
    else 
         a = frequency_reverse(a)
         printf(SCREEN, "Results of %d-point fft:\n\n", length(a))
    end if
    output(a)
end procedure

main(0)

---------------------------------------------------
Regards,
   Rob Craig
   Rapid Deployment Software
   http://www.RapidEuphoria.com