FAST FOURIER TRANSFORMS

To load the routines do LOADFILE(FFT,FASL,DSK,SHARE);

The basic funtions are:
	FFT --> fast fourier transform
	IFT --> inverse fast fourier transform

These functions perform a (complex) fast fourier transform on either
1 or 2 dimensional FLOATING-POINT arrays, obtained by:
ARRAY(<ary>,FLOAT,<dim1>); or ARRAY(<ary>,FLOAT,<dim1>,<dim2>);
For 1D arrays <dim1> must equal 2^n-1, and for 2D arrays
<dim1>=<dim2>=2^n-1 (i.e. the array is square).  (Recall that
MACSYMA arrays are indexed from a 0 origin so that there will be
2^n and (2^n)^2 arrays elements in the above two cases.)

Calling sequence is:
FFT(<real-array>,<imag-array>); or IFT(<real-array>,<imag-array>);

The real and imaginary arrays must of course be the same size.
The transforms are done in place so that <real-array> and <imag-
array> will contain the real and imaginary parts of the transform.
(If you want to keep the transformed and un-transfromed arrays
separate copy the arrays before calling FFT or IFT using the
FILLARRAY function - see SHARE;ARRAY USAGE.)

The definitions of the Fast Fourier Transform and its inverse
are given here.  Here A is the array to be transformed and AT is
its transform.  Both A and AT are complex arrays, although as noted
above FFT and IFT can only deal with separate real arrays for
the real and imaginary parts of A and AT.  N (or N^2) is the number
of elements in A in the 1D (or 2D) case.  (In fact these definitions
are not of the FFTs but of the discrete Fourier transforms.  The
FFT and IFT functions merely provided efficient algorithms for the
implementation of these definitions.)

1D case:

      N - 1
      ====		       - 1
      \		 2 %I %PI I K N
AT  =  >    A  %E		   
  K   /	     I
      ====
      I = 0

	  N - 1
	  ====			      - 1
      - 1 \	      - 2 %I %PI I K N
A  = N	   >    AT  %E			  
 I	  /	  K
	  ====
	  K = 0

2D case:

	 N - 1 N - 1
	 ====  ====				   - 1
	 \     \	     2 %I %PI (I K + J L) N
AT     =  >     >    A	   %E			       
  K, L	 /     /      I, J
	 ====  ====
	 I = 0 J = 0

	     N - 1 N - 1
	     ====  ====					  - 1
	 - 2 \	   \		  - 2 %I %PI (I K + J L) N
A     = N     >	    >	 AT     %E			      
 I, J	     /	   /	   K, L
	     ====  ====
	     K = 0 L = 0

Other functions included in this file are:

POLARTORECT(<magnitude-array>,<phase-array>);
converts from magnitude/phase form into real/imaginary form
putting the real part in the magnitude array and the imaginary part
into the phase array.  (<real>=<magnitude>*COS(<phase>)
and <imaginary>=<magnitude>*SIN(<phase>).)

RECTTOPOLAR(<real-array>,<imag-array>);
undoes POLARTORECT.  The phase is given in the range
from -%PI to %PI.

Like FFT and IFT these function accept 1 or 2 dimensional
arrays.  However, the array dimensions need not be a power of 2,
nor need the 2D arrays be square.

(The above 4 functions return a list of their arguments)

EXAMPLE - do LOADFILE(FFT,FASL,DSK,SHARE); DEMO(FFT,DEMO,DSK,SHARE);

DISCLAIMER - I didn't write these routines, but copied them from
AI:LIBDOC;FFT BWOOD1.  However you can report bugs to me.

						CFFK@MC