fft.sa
Generated by gen_html_sa_files from ICSI. Contact gomes@icsi.berkeley.edu for details
---------------------------> Sather 1.1 source file <--------------------------
-- Copyright (C) International Computer Science Institute, 1994. COPYRIGHT --
-- NOTICE: This code is provided "AS IS" WITHOUT ANY WARRANTY and is subject --
-- to the terms of the SATHER LIBRARY GENERAL PUBLIC LICENSE contained in --
-- the file "Doc/License" of the Sather distribution. The license is also --
-- available from ICSI, 1947 Center St., Suite 600, Berkeley CA 94704, USA. --
--------> Please email comments to sather-bugs@icsi.berkeley.edu. <----------
-- fft.sa: Sather 1.0 Radix-2 FFT routines.
--
-- Author: Jeff A. Bilmes <bilmes@icsi.berkeley.edu>
-- The fft routines themselves have been collected over the years
-- from various sources, but a primary source has been:
-- "Discrete-Time Signal Processing" by Oppenheim and Schafer,
-- Prentice-Hall, ISBN 0-13-216292-X
class FFT
class FFT is
-- FFT class. Useful routines include:
-- fft_basic: Radix-2 FFT, in-place and in-order. in-order result.
-- fft_real : Radix-2 FFT, fft_basic, but optimized for real input data.
-- fft_torl : Radix-2 iFFT, but optimized for real result data.
-- dft: direct DFT routine.
-- fft: various flavors of simplified radix-2 fft routines.
-- fft_2dimensional : Image FFT, for real data.
-- ifft_2dimensional: inverse of the above
--
-- For those fft's that take a 'twid' argument (twiddle factors),
-- you must first assign the complex exponential factors using the
-- routine 'assign_basis(size)'. Some routines take a twid argument
-- thereby optimizing for speed.
private shared twiddles:FLIST{ARRAY{CPX}}; -- prestored twiddle factors.
init is
twiddles := #FLIST{ARRAY{CPX}}(1);
end;
-- Basic support routines. No real need to be private since they might
-- be useful.
pure_real(x:ARRAY{CPX},n:INT):BOOL is
-- Query whether the data in x is all real.
loop n.times!;
if x.elt!.im /= 0.0 then
return false; end;
end;
return true; end;
is_power_of_two(n:INT):BOOL is
-- Query whether n is a 2^k for some k
loop while!(n>1);
if n.mod(2) /= 0 then
break! end;
n := n.rshift(1); end;
return (n = 1); end;
find_table(n:INT):ARRAY{CPX} is
-- search our list of existing twiddle arrays
loop twid:ARRAY{CPX} := twiddles.elt!;
if twid.size = n then
return twid; end; end;
return void; end;
assign_basis(n:INT):ARRAY{CPX} is
-- Return up the complex array to be the array
-- of <cos,sin> pairs (i.e., twiddle factors) corresponding
-- to the desired fft size 'n'. This routine allocates
-- memory for you.
expn:ARRAY{CPX};
i:INT;
expn := find_table(n);
if ~void(expn) then
return expn; end;
expn := #ARRAY{CPX}(n);
i := 0; loop while!(i<n);
a::= 2.0* FLT::pi * i.flt/n.flt;
expn[i] := #CPX(a.cos,-(a.sin));
i := i + 1;
end; -- loop
twiddles := twiddles.push(expn);
return expn; end;
reverse_dig(x:ARRAY{CPX},n:INT,stride:INT) is
-- bit reverse the elements in array x
i,j,k,jj,ii:INT;
tmp:CPX;
i := 0; j := 0; loop while!(i<n-1);
if i < j then
jj := j * stride;
ii := i * stride;
tmp := x[jj];
x[jj] := x[ii];
x[ii] := tmp;
end;
k := n/2;
loop while!(k<=j);
j := j - k;
k := k.rshift(1); end;
j := j + k;
i := i + 1; end; end;
int_2_complex(s:ARRAY{INT},c:ARRAY{CPX}):ARRAY{CPX} is
-- create a complex array fttable from an integer array
i:INT;
if void(c) then
c := #ARRAY{CPX}(s.size); end;
i :=0; loop while!(i<s.size);
c[i] := #CPX(s[i].flt,0.0);
i := i + 1; end;
return c; end;
conj_scale(x:ARRAY{CPX},n:INT,scale:FLT) is
-- Conjugate and scale complex data (e.g. prior to IFFT by FFT)
loop n.times!;
x.set!(x.elt!.conjugate * scale);
end; end;
-- miscale:FLT := -scale;
-- loop while!(n > 0);
-- n := n-1;
-- x[n] := #CPX(x[n].re * scale,x[n].im * miscale);
-- end; end;
-- The core fft routines. I.e., ones that a user probably won't use.
reals(x:ARRAY{CPX},stride:INT,sign:INT,twid:ARRAY{CPX}) is
-- This routine separates out the reals from the result
-- when we do a real n point fft using an n/2 point complex fft.
-- 'sign' is 1 for FFT2real, -1 for torl.
n:INT := x.size;
half:INT := n/2;
quarter:INT := n/4;
m,mm:INT;
a,b,s,t:CPX;
half := half * stride;
n := n * stride;
if sign = 1 then x[half] := x[0]; end; -- only for Yc to Xr
mm := 0; loop while!(mm <= quarter);
m := mm * stride;
s := #CPX( (1.0 + twid[mm].im) / 2.0,
-sign.flt * twid[mm].re / 2.0);
t := #CPX( 1.0 - s.re, -s.im);
a := x[m];
b := x[half-m].conjugate;
a := a * s;
b := b * t;
b := b + a;
a := x[m].conjugate;
x[m] := b;
if m /= 0 then
b := b.conjugate;
x[n - m] := b;
end;
b := x[half - m];
a := a * (t.conjugate);
b := b * (s.conjugate);
b := b+a;
x[half - m] := b;
if m /= 0 then
b := b.conjugate;
x[half + m] := b;
end;
mm := mm + 1; end; end;
fft_raw(x:ARRAY{CPX},n:INT,dilate:INT,stride:INT,twid:ARRAY{CPX}) is
-- Data is x.
-- Data size is n.
-- 'dilate' means: library global expn is the (cos, -j sin)
-- array, EXCEPT for effective data size n/dilate,
-- stride is the offset of each successive data term, as
-- in "fft" above.
j,m,p,q,i,k,n2,n1:INT;
tmp,tmp_twid:CPX;
m := 1;
n2 := n;
loop while!(m < n);
n1 := n2;
n2 := n2.rshift(1);
j := 0; q:= 0; loop while!(j < n2);
tmp_twid := twid[q];
q := q + m*dilate;
k := j; loop while!(k<n);
p := (k + n2) * stride;
i := k * stride;
tmp := x[i] - x[p];
x[i] := x[i] + x[p];
x[p] := tmp_twid * tmp;
k := k + n1; end;
j := j + 1; end;
m := m.lshift(1); end; end;
ifft_stride(x:ARRAY{CPX},stride:INT) is
-- optimized for stride indexing.
ex:ARRAY{CPX};
tmp:CPX;
i:INT;
j,k:INT;
n::= x.size;
ex := assign_basis(n);
-- reverse direction of x
i :=0; loop while!(i<n/2);
j := (n - 1 - i) * stride;
k := i * stride;
tmp := x[j];
x[j] := x[k];
x[k] := tmp;
i := i + 1; end;
fft_basic(x,stride);
-- now do final twiddle
i:=0; loop while!(i<n);
k := i * stride;
x[k] := x[k]*ex[i] / n.flt;
i := i + 1; end; end;
-- The basic fft routines.
fft_real(x:ARRAY{CPX},stride:INT,twid:ARRAY{CPX})
pre is_power_of_two(x.size) is
-- Perform FFT on real input data.
-- Data must be arranged as {re, 0, re, 0, re, 0...}.
-- So input data is assumed to be pure real.
-- The full in-order complex result is left in x.
-- This routine is optimized for real input data by
-- doing an n/2-point complex FFT instead of an n-point
-- complex FFT. See "Discrete-Time Signal Processing"
-- by Oppenheim and Schafer, problem 9.31 on page 656
-- for a detailed description of this method.
half:INT := x.size.rshift(1);
m,mm:INT;
if half.rshift(1) = 0 then return; end;
-- pack length n real signal into length n/2 complex signal
mm := 0;
loop while!(mm < half);
m := mm * stride;
x[m] := #CPX(x[m.lshift(1)].re,x[m.lshift(1) + stride].re);
mm := mm + 1;
end; -- loop
-- do n/2 length fft
fft_raw(x,half,2,stride,twid);
reverse_dig(x,half,stride);
-- extract the length n result
reals(x,stride,1,twid);
end; -- fft_real
fft_real(x:ARRAY{CPX},stride:INT) is
fft_real(x,stride,assign_basis(x.size));
end;
fft_torl(x:ARRAY{CPX},stride:INT,scale:FLT,twid:ARRAY{CPX})
pre is_power_of_two(x.size) is
-- Performs FFT on data assumed complex conjugate symetric to give purely
-- real result. This is essentially the inverse of FFT2real.
half:INT := x.size.rshift(1);
m,mm:INT;
if half.rshift(1) = 0 then return; end;
reals(x,stride,-1,twid);
conj_scale(x, half + 1, 2.0 * scale);
fft_raw(x, half, 2, stride, twid);
reverse_dig(x, half, stride);
mm := half-1;
loop while!(mm >= 0);
m := mm * stride;
x[m*2] := x[m*2].re(x[m].re);
-- We need to conjugate as we separate result for a true ifft.
x[m*2+stride] := x[m*2+stride].re(-(x[m].im));
x[m*2] := x[m*2].im(0.0);
x[m*2+stride] := x[m*2+stride].im(0.0);
mm := mm - 1; end; end;
fft_torl(x:ARRAY{CPX},stride:INT,scale:FLT) is
fft_torl(x,stride,scale,assign_basis(x.size));
end; -- fft_torl
fft_basic(x:ARRAY{CPX},stride:INT,twid:ARRAY{CPX})
pre is_power_of_two(x.size) is
-- Perform FFT for n = a power of 2.
-- Operate and modify array x in place. Result
-- goes in x in order.
-- The relevant data are the complex numbers
-- x[0],x[stride],x[2*stride], ...
fft_raw(x, x.size, 1, stride, twid);
reverse_dig(x, x.size, stride);
end;
fft_basic(x:ARRAY{CPX},stride:INT) is
fft_basic(x,stride,assign_basis(x.size));
end;
dft(data:ARRAY{CPX},res:ARRAY{CPX},twid:ARRAY{CPX}):ARRAY{CPX} is
-- Perform direct Discrete Fourier Transform (DFT).
-- 'data' may be any length.
j, k, m:INT;
n:INT := data.size;
if void(res) then
res := #ARRAY{CPX}(data.size);
end;
j := 0; loop while!(j<n);
res[j] := #CPX(0.0,0.0);
k := 0; loop while!(k<n);
m := (j * k) % n;
res[j] := res[j] + data[k]*twid[m];
k := k + 1; end;
j := j + 1; end;
return res;
end; -- dft
-- Simplified FFT routine interface.
-- for arbitrary complex input data.
fft(x:ARRAY{CPX},stride:INT) is
-- complex data in x, complex result in x in order.
ex:ARRAY{CPX};
ex := assign_basis(x.size);
fft_basic(x,stride,ex);
end;
fft(x:ARRAY{CPX}) is fft(x,1); end;
ifft(x:ARRAY{CPX},stride:INT) is
-- inverse fft, complex data in x, complex result in x.
conj_scale(x,x.size,1.0/x.size.flt);
fft(x,stride);
conj_scale(x,x.size,1.0);
-- alternatively, we could do: ifft_stride(x,stride);
end;
ifft(x:ARRAY{CPX}) is ifft(x,1); end;
-- for real data
fft(x:ARRAY{FLT},stride:INT):ARRAY{CPX} is
-- real data in x, in-order result returned as a new array.
res ::= #ARRAY{CPX}(x.size);
loop res.set!(#CPX(x.elt!,0.0)); end;
fft_real(res,stride);
return res;
end;
fft(x:ARRAY{FLT}):ARRAY{CPX} is return fft(x,1); end;
ifft_real(x:ARRAY{CPX},stride:INT) is
-- assumped complex conj. sym. data in x, from real data.
-- Produces a real result (i.e., .im attributes are all zeros).
fft_torl(x,stride,1.0/x.size.flt);
end;
ifft_real(x:ARRAY{CPX}) is ifft_real(x,1); end;
-- simple 2D routines.
fft(x:ARRAY2{CPX}) is
-- Perform 2D FFT on image (both dims must be powers of 2).
-- We assume here that row order is fastest in 'x'.
i:INT;
tmp:ARRAY{CPX};
twid:ARRAY{CPX};
-- do rows
tmp := #(x.nc);
twid := assign_basis(x.nc);
i:=0; loop i := x.row_ind!;
loop tmp.set!(x.row_elt!(i)); end;
fft_basic(tmp,1,twid);
loop x.set_row!(i,tmp.elt!); end;
end;
-- do cols
tmp := #(x.nr);
twid := assign_basis(x.nr);
i:=0; loop i := x.col_ind!;
loop tmp.set!(x.col_elt!(i)); end;
fft_basic(tmp,1,twid);
loop x.set_col!(i,tmp.elt!); end;
end; end;
ifft(x:ARRAY2{CPX}) is
ifft_2dimensional(x); end;
fft_2dimensional(x:ARRAY2{CPX}) is
-- Perform 2D FFT on image (both dims must be powers of 2).
-- IMAGE IS ASSUMED PURE-REAL (i.e., imaginary portions are zeros).
-- If the image is not real, the 'fft_real' call should be change
-- to 'fft_basic'.
-- We assume here that row order is fastest in 'x'.
i:INT;
tmp:ARRAY{CPX};
twid:ARRAY{CPX};
-- do rows
tmp := #(x.nc);
twid := assign_basis(x.nc);
i:=0; loop i := x.row_ind!;
loop tmp.set!(x.row_elt!(i)); end;
fft_real(tmp,1,twid);
loop x.set_row!(i,tmp.elt!); end;
end;
-- do cols
tmp := #(x.nr);
twid := assign_basis(x.nr);
i:=0; loop i := x.col_ind!;
loop tmp.set!(x.col_elt!(i)); end;
fft_basic(tmp,1,twid);
loop x.set_col!(i,tmp.elt!); end;
end; end;
ifft_2dimensional(x:ARRAY2{CPX}) is
-- Perform 2D inverse FFT on image (both dims must be powers of 2).
-- Assume row indexing is the fastest varying index.
i:INT;
tmp:ARRAY{CPX};
twid:ARRAY{CPX};
-- do rows
tmp := #ARRAY{CPX}(x.nc);
twid := assign_basis(x.nc);
i:=0; loop i := x.row_ind!;
loop tmp.set!(x.row_elt!(i)); end;
-- do the ifft over the row.
conj_scale(tmp,tmp.size,1.0/tmp.size.flt);
fft_basic(tmp,1,twid);
conj_scale(tmp,tmp.size,1.0);
loop x.set_row!(i,tmp.elt!); end;
end;
-- do cols
tmp := #(x.nr);
twid := assign_basis(x.nr);
i:=0; loop i := x.col_ind!;
loop tmp.set!(x.col_elt!(i)); end;
-- do the ifft over the col.
conj_scale(tmp,tmp.size,1.0/tmp.size.flt);
fft_basic(tmp,1,twid);
conj_scale(tmp,tmp.size,1.0);
loop x.set_col!(i,tmp.elt!); end;
end; end;
-- simple 3D routines.
fft(x:ARRAY3{CPX}) is
-- Perform 3D FFT on image (both dims must be powers of 2).
-- We assume here that row order is fastest in 'x'.
i,j:INT;
tmp:ARRAY{CPX};
twid:ARRAY{CPX};
-- do rows
tmp := #(x.nc);
twid := assign_basis(x.nc);
i:=0; loop i := x.row_ind!;
j:=0; loop j := x.elem_ind!;
loop tmp.set!(x.row_elt!(i,j)); end;
fft_basic(tmp,1,twid);
loop x.set_row!(i,j,tmp.elt!); end;
end; end;
-- do cols
tmp := #(x.nr);
twid := assign_basis(x.nr);
i:=0; loop i := x.col_ind!;
j:=0; loop j := x.elem_ind!;
loop tmp.set!(x.col_elt!(i,j)); end;
fft_basic(tmp,1,twid);
loop x.set_col!(i,j,tmp.elt!); end;
end; end;
-- do elements
tmp := #(x.ne);
twid := assign_basis(x.ne);
i:=0; loop i := x.row_ind!;
j:=0; loop j := x.col_ind!;
loop tmp.set!(x.elem_elt!(i,j)); end;
fft_basic(tmp,1,twid);
loop x.set_elem!(i,j,tmp.elt!); end;
end; end;
end;
ifft(x:ARRAY3{CPX}) is
ifft_3dimensional(x); end;
fft_3dimensional(x:ARRAY3{CPX}) is
-- Perform 3D FFT on image (both dims must be powers of 2).
-- IMAGE IS ASSUMED PURE-REAL (i.e., imaginary portions are zeros).
-- If the image is not real, the 'fft_real' call should be change
-- to 'fft_basic'.
-- We assume here that row order is fastest in 'x'.
i,j:INT;
tmp:ARRAY{CPX};
twid:ARRAY{CPX};
-- do rows
tmp := #(x.nc);
twid := assign_basis(x.nc);
i:=0; loop i := x.row_ind!;
j:=0; loop j := x.elem_ind!;
loop tmp.set!(x.row_elt!(i,j)); end;
fft_real(tmp,1,twid);
loop x.set_row!(i,j,tmp.elt!); end;
end; end;
-- do cols
tmp := #(x.nr);
twid := assign_basis(x.nr);
i:=0; loop i := x.col_ind!;
j:=0; loop j := x.elem_ind!;
loop tmp.set!(x.col_elt!(i,j)); end;
fft_basic(tmp,1,twid);
loop x.set_col!(i,j,tmp.elt!); end;
end; end;
-- do elements
tmp := #(x.ne);
twid := assign_basis(x.ne);
i:=0; loop i := x.row_ind!;
j:=0; loop j := x.col_ind!;
loop tmp.set!(x.elem_elt!(i,j)); end;
fft_basic(tmp,1,twid);
loop x.set_elem!(i,j,tmp.elt!); end;
end; end;
end;
ifft_3dimensional(x:ARRAY3{CPX}) is
-- Perform 3D inverse FFT on image (both dims must be powers of 2).
-- Assume row indexing is the fastest varying index.
i,j:INT;
tmp:ARRAY{CPX};
twid:ARRAY{CPX};
-- do rows
tmp := #ARRAY{CPX}(x.nc);
twid := assign_basis(x.nc);
i:=0; loop i := x.row_ind!;
j:=0; loop j := x.elem_ind!;
loop tmp.set!(x.row_elt!(i,j)); end;
-- do the ifft over the row.
conj_scale(tmp,tmp.size,1.0/tmp.size.flt);
fft_basic(tmp,1,twid);
conj_scale(tmp,tmp.size,1.0);
loop x.set_row!(i,j,tmp.elt!); end;
end;
end;
-- do cols
tmp := #(x.nr);
twid := assign_basis(x.nr);
i:=0; loop i := x.col_ind!;
j:=0; loop j := x.elem_ind!;
loop tmp.set!(x.col_elt!(i,j)); end;
-- do the ifft over the col.
conj_scale(tmp,tmp.size,1.0/tmp.size.flt);
fft_basic(tmp,1,twid);
conj_scale(tmp,tmp.size,1.0);
loop x.set_col!(i,j,tmp.elt!); end;
end;
end;
-- do elements
tmp := #(x.ne);
twid := assign_basis(x.ne);
i:=0; loop i := x.row_ind!;
j:=0; loop j := x.col_ind!;
loop tmp.set!(x.elem_elt!(i,j)); end;
-- do the ifft over the col.
conj_scale(tmp,tmp.size,1.0/tmp.size.flt);
fft_basic(tmp,1,twid);
conj_scale(tmp,tmp.size,1.0);
loop x.set_elem!(i,j,tmp.elt!); end;
end;
end;
end;
end; -- class FFT
class TEST_FFT
class TEST_FFT is
const sig_length:INT := 16;
main is
real_sig:ARRAY{FLT};
real_sig_fft:ARRAY{CPX};
real_image:ARRAY2{CPX};
tmp:ARRAY{CPX};
cpx_sig:ARRAY{CPX};
i:INT;
-- test real versions.
#OUT+"Real array\n";
real_sig := #ARRAY{FLT}(sig_length);
loop real_sig.set!(-7.0+15.0*(RND::uniform.flt)); end;
loop #OUT+real_sig.elt!+" "; end; #OUT+"\n";
real_sig_fft := FFT::fft(real_sig);
#OUT+"FFT of above\n";
loop #OUT+real_sig_fft.elt!.str+" " end; #OUT+"\n";
FFT::ifft_real(real_sig_fft);
#OUT+"iFFT of above\n";
loop #OUT+real_sig_fft.elt!.str+" " end; #OUT+"\n";
-- test complex versions
#OUT+"\nComplex array\n";
cpx_sig := #ARRAY{CPX}(sig_length);
loop cpx_sig.set!(#CPX(15.0*(RND::uniform.flt)-7.0,15.0*(RND::uniform.flt)-7.0)); end;
loop #OUT+cpx_sig.elt!.str+" " end; #OUT+"\n";
FFT::fft(cpx_sig);
#OUT+"FFT of above\n";
loop #OUT+cpx_sig.elt!.str+" " end; #OUT+"\n";
tmp := #ARRAY{CPX}(sig_length);
loop tmp.set!(cpx_sig.elt!); end;
FFT::ifft(cpx_sig);
#OUT+"iFFT of above\n";
loop #OUT+cpx_sig.elt!.str+" " end; #OUT+"\n";
FFT::ifft_stride(tmp,1);
#OUT+"iFFT_stride of above\n";
loop #OUT+tmp.elt!.str+" " end; #OUT+"\n";
-- test dft
#OUT+"\nDFT\n";
cpx_sig := #ARRAY{CPX}(sig_length+1);
loop cpx_sig.set!(#CPX(15.0*(RND::uniform.flt)-7.0,15.0*(RND::uniform.flt)-7.0)); end;
loop #OUT+cpx_sig.elt!.str+" " end; #OUT+"\n";
cpx_sig := FFT::dft(cpx_sig,void,FFT::assign_basis(cpx_sig.size));
#OUT+"DFT of above\n";
loop #OUT+cpx_sig.elt!.str+" " end; #OUT+"\n";
-- test 2d FFT
real_image := #ARRAY2{CPX}(8,8);
loop real_image.set!(#CPX(15.0*(RND::uniform.flt)-7.0,0.0)); end;
#OUT+"\nImage array\n";
loop i := real_image.row_ind!;
loop #OUT+real_image.row_elt!(i).str+" " end; #OUT+"\n"; end;
#OUT+"Transform of above\n";
FFT::fft_2dimensional(real_image);
loop i := real_image.row_ind!;
loop #OUT+real_image.row_elt!(i).str+" " end; #OUT+"\n"; end;
#OUT+"Inverse Transform of above\n";
FFT::ifft_2dimensional(real_image);
loop i := real_image.row_ind!;
loop #OUT+real_image.row_elt!(i).str+" " end; #OUT+"\n"; end;
end; -- main
end;