Program:      FFTN

Purpose:      Apply fast Fourier transform in N dimensions, i.e. for a profile (N=1),
              a map or a series of maps (N=2) or a data cube (N=3).

Category:     MANIPULATION

File:         fftn.src

Author:       M.G.R. Vogelaar


Keywords:

   INVERSE=   FFT or Inverse FFT: ... [F]/I

              Set modus of operation: FFT or the inverse transformation.


   PHASEIN=   Input data in amplitude and phase? ... Y/[N]

              If True, then the phases must be entered in radians.


   RINSET=    Set, subset(s) real part of input:

              Note that if your input is not contiguous (e.g.:
              Aurora FREQ 10 14:16 20) then the output will be from
              FREQ 10 to 20 and the subsets that are not included in the input
              are set to blank in the output. The transform is only applied to
              the given subsets.
              

   IINSET=    Set, subset imaginary part of input:


   BOX=       Enter box in ....


   PHASEOUT=  Output data in amplitude and phase? ... Y/[N]

              The output phases are in radians.
              The set names get the extension _amp and _phase
              Else they get extensions _real and _imag


              If INVERSE=n

   SHIFTIN=   Set first element in sample to the center of the input data? ... Y/[N]

   SHIFTOUT=  Shift the zero-freq. component of OUTPUT to the center of the spectrum? ... Y/[N]


              If INVERSE=y
              
   SHIFTIN=   Shift the zero-freq. component of INPUT to the left of the spectrum? ... Y/[N]

   SHIFTOUT=  Set first element of transform to be the center of the output data? ... Y/[N]"


   OUTSET=    Name of output set:


Description:
              The program calculates an N dimensional Fourier transform of the
              input data.
              The number of axes in the selected subset(s) sets the value of N.
              With this program it is possible to do a 3-dimensional FFT on a
              data cube or one or more 1-dimensional FFT's in an image or
              2-dimensional FFT's on channel maps in a data cube etc.


              Gaussians
              ---------
              As a testcase, one often uses a Gaussian function as a real input.
              The FFT of this symmetrical real function should be also gaussian
              and the imaginary part is zero. An inverse FFT should return the
              original image.

              With numpy's fft() function (and therefore in this task), we need
              to do something extra to achieve the same result.
              First of all, the values in the result of an FFT follow so-called
              standard order: If A = fft(a, n),
              then A[0] contains the zero-frequency term (the mean of the signal),
              which is always purely real for real inputs. Then A[1:n/2] contains
              the positive-frequency terms, and A[n/2+1:] contains the negative-
              frequency terms, in order of decreasingly negative frequency.
              The routine fftshift(A) shifts transforms and their
              frequencies to put the zero-frequency components in the middle.
              This action is set with keyword SHIFTOUT=Y

              But if your input is a Gaussian with the peak in the center, your
              output wil not be Gaussian and the imaginary part is non-zero.
              The reason is that the input is not a real-valued, symmetrical time
              series because the peak is not at t=0.
              We have to shift the peak of the Gaussian to time = 0 and wrap
              around the shifted data for t < 0. The function that does this is
              called ifftshift(A) and is executed if we set SHIFTIN=Y

              From Gaussian with the peak in the center to Gaussian with the
              peak in the center, use:

              BOX=
              IINSET=
              INVERSE= n
              OUTSET= gaussout
              PHASEIN=
              PHASEOUT=
              RINSET= gaussin
              SHIFTIN= y
              SHIFTOUT= y

              To get the original Gaussian back:

              BOX=
              IINSET= gaussout_imag
              INVERSE= y
              OUTSET= igauss
              PHASEIN=
              PHASEOUT=
              RINSET= gaussout_real
              SHIFTIN= y
              SHIFTOUT= y

             
Notes:        If you shifted a transform and want to apply the inverse transform,
              you should also set SHIFTIN=Y because the input is not the FFT but
              the shifted FFT.
                          

Updates:      Jul 30, 2014: VOG, Document created.