1. FFT algorithm
Does anybody have a pre-built FFT (Fast Fourier Transform) algorithm
already built in Euphoria?
Michael J. Sabal
2. Re: FFT algorithm
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
3. Re: FFT algorithm
Thank you. I need to convert the amplitude wave form from the sound
card into a frequency chart, which I can cross-reference into
alphanumerics. I'll also need to do the reverse. The Linux soundmodem
source files I have use the FFT algorithm to accomplish this, but it
uses so much bit shifting, I had trouble following how to do the same
more efficiently in Euphoria. Your code should help greatly.
>>> rds at RapidEuphoria.com 12/13/01 05:45PM >>>
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.
--SNIP--
4. Re: FFT algorithm
Where'd the coed go?
Alvin
Homepage
http://ka9qlq.tripod.com/home/
Where I live.
http://ka9qlq.tripod.com/CCC/
----- Original Message -----
From: <Sabal.Mike at notations.com>
To: "EUforum" <EUforum at topica.com>
Sent: Friday, December 14, 2001 9:16 AM
Subject: Re: FFT algorithm
>
> Thank you. I need to convert the amplitude wave form from the sound
> card into a frequency chart, which I can cross-reference into
> alphanumerics. I'll also need to do the reverse. The Linux soundmodem
> source files I have use the FFT algorithm to accomplish this, but it
> uses so much bit shifting, I had trouble following how to do the same
> more efficiently in Euphoria. Your code should help greatly.
>
> >>> rds at RapidEuphoria.com 12/13/01 05:45PM >>>
>
> 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.
> --SNIP--
>
>
>
