------------------------->  GNU Sather - sourcefile  <-------------------------
-- Copyright (C) 1994 by International Computer Science Institute            --
-- This file is part of the GNU Sather library. It is free software; you may --
-- redistribute  and/or modify it under the terms of the GNU Library General --
-- Public  License (LGPL)  as published  by the  Free  Software  Foundation; --
-- either version 2 of the license, or (at your option) any later version.   --
-- This  library  is distributed  in the  hope that it will  be  useful, but --
-- WITHOUT ANY WARRANTY without even the implied warranty of MERCHANTABILITY --
-- or FITNESS FOR A PARTICULAR PURPOSE. See Doc/LGPL for more details.       --
-- The license text is also available from:  Free Software Foundation, Inc., --
-- 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA                     --
-------------->  Please email comments to <bug-sather@gnu.org>  <--------------

-- 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 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 FFT_TEST 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;
