129 lines
3 KiB
Text
129 lines
3 KiB
Text
--
|
|
-- 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)))
|