functions in fft_utils.i -

             __fft(x);  
 
       -or- __fft(x, dir);  
     Replacement for Yorick's fft to speed up fast Fourier transforms (FFT)  
     in,  e.g.,  iteratives  algorithms  that  necessitate  computation  of  
     several  FFT's  of  arrays  with  same dimension  list.   The  FFT  is  
     performed on  all dimensions of  X and DIR  must be a  scalar (default  
     +1). FFT workspace must be initialized by __fft_init.  

             __fft_init, dimlist;  
 
     Initializes  FFT  workspace  for  further  calls to  __fft  (to  see).  
     DIMLIST is the dimension list of the arrays to transform.  The routine  
     defines 2 external symbols:       
       __fft_setup  - used to store the FFT workspace)  
       __fft_number - used to keep track of the number of calls to __fft  
     In order to avoid namespace pollution/clash, a routine that uses __fft  
     should declare these symbols as local before calling __fft_init, e.g.:  
       local __fft_setup, __fft_number;  
       __fft_init, dimlist;  

             __fft_pl2d_limits, z, scale;  
 
     Private routine used by fft_pli, fft_plc and fft_plfc.  

             _fft_centroid(a1)  
 
       -or- _fft_centroid(a1, repeat)  
     Working routine for fft_centroid: return position of centroid of 1D array  
     A1 assuming FFT geometry for the coordinates.  

             abs2(x)  
 
     Returns abs(X)^2  

             fft_best_dim(len);  
 
     Return the smallest integer which is greater or equal LEN and which is  
     a multiple of powers of 2, 3 and/or 5  

             fft_centroid(a)  
 
       -or- fft_centroid(a, repeat)  
     Return  the  position  of  centroid  of N-dimensional  array  A  assuming  
     coordinates  along  dimensions  of  A  are  wrapped  as  in  a  FFT  (see  
     fft_indgen).  The  algorithm proceeds by computing the  center of gravity  
     of A around its  central element which is the maximum of  A for the first  
     iteration  and  the  closest  to  the previously  computed  centroid  for  
     subsequent  iterations.   The  maximum  number  of  iteration  is  REPEAT  
     (default: 3;  in any  cases, at least  one iteration is  performed).  The  
     Nyquist frequency along each even dimension is omitted to avoid a bias.  

             fft_convolve(orig, psf);  
 
       -or- fft_convolve(orig, psf, do_not_roll);  
     Return discrete convolution (computed by FFT) of array ORIG by point  
     spread function PSF.  Unless argument DO_NOT_ROLL is true, PSF is  
     rolled before.  Note: ORIG and PSF must have same dimension list.  
       

             fft_dist(dimlist);  
 
       -or- fft_dist(dim1, dim2, ...);  
     Returns Euclidian lenght of spatial frequencies in frequel units for a  
     FFT of dimensions DIMLIST.  
     If keyword  NYQUIST is true,  the frequel coordinates get  rescaled so  
     that the Nyquist frequency is  equal to NYQUIST along every dimension.  
     This is obtained by using coordinates:  
        (2.0*NYQUIST/DIM(i))*fft_indgen(DIM(i))  
     along i-th dimension of lenght DIM(i).  
     If  keyword  SQUARE is  true,  the square  of  the  Euclidian norm  is  
     returned instead.  
       

 fft_fine_shift  
 

             ptr = fft_freqlist(dimlist)  
 
     returns a vector of DIMLIST(1) pointers with normalized FFT  
     frequencies along all dimensions of DIMLIST (must be  
     dimsof(SOME_ARRAY)) with "adequate" geometry:  
       *ptr(1) =   (2*pi/dimlist(2))*fft_indgen(dimlist(2));  
       *ptr(2) =  [(2*pi/dimlist(3))*fft_indgen(dimlist(3))];  
       *ptr(3) = [[(2*pi/dimlist(4))*fft_indgen(dimlist(4))]];  
        ...  
       

             fft_gaussian_psf(dimlist, fwhm)  
 
       -or- fft_gaussian_mtf(dimlist, fwhm)  
     Returns   normalized   Gaussian  point   spread   function  (PSF)   or  
     corresponding modulation  transfer function (MTF)  with dimension list  
     DIMLIST  and full  width at  half maximum  equals to  FWHM  along each  
     dimensions (in  the PSF space).  Up  to errors due  to limited support  
     and/or numerical precision, the PSF and the MTF obey:  
       
        sum(PSF) = MTF(1) = 1                       (normalization)  
        MTF = fft(PSF, +1)  
        PSF = fft(MTF, -1)/numberof(MTF)  
          
     where MTF(1) is the 0-th frequency in the MTF.  The standard deviation  
     SIGMA and the FWHM are related by:  
        FWHM = sqrt(8*log(2))*SIGMA  
             ~ 2.354820045031*SIGMA  
          
     Note that, owing  to the limited size of  the support and/or numerical  
     precision, these properties may not be perfectly met; for that reason,  
     _always_ compute directly  what you need, e.g. do not  take the FFT of  
     the PSF if what  you need is the MTF.  Also note  that the geometry is  
     that of  the FFT and that  for unequal dimension lengths,  the PSF has  
     the same width (in "pixels") along every dimension but not the MTF.  

             fft_indgen(len)  
 
     Return FFT frequencies along a dimension of length LEN.  

 fft_interp  
 

 fft_interp_complex  
 

             fft_interp(a, off, setup=)  
 
       -or- fft_interp_real(z, phasor)  
       -or- fft_interp_complex(z, phasor)  
     returns value obtained by Fourier interpolation of A at offset OFF.  
     The function fft_interp computes the forward FFT of A and can make use  
     of pre-computed FFT workspace specified by keyword SETUP.  The two  
     other functions (fft_interp_complex, if A is complex; fft_interp_real  
     otherwise) are usefull when the forward FFT of A and/or the complex  
     phasor corresponding to the shift are already computed, their  
     arguments are:  
       Z = fft(A, +1);  
       PHASOR = fft_shift_phasor(OFF, dimsof(A));  
       

             fft_of_two_real_arrays, a, b, ft_a, ft_b, direction;  
 
       -or- fft_of_two_real_arrays, a, b, ft_a, ft_b, ljdir, rjdir;  
     Computes the FFT  of arrays A and  B and stores them in  TF_A and FT_B  
     respectively.  A and B must have same dimension list.  A single FFT is  
     needed.  Agrguments  DIRECTION, LJDIR,  RJDIR, and keyword  SETUP have  
     the same meaning as for the fft function (which see).  
       

             fft_plc, a;  
 
     Plot contour levels  of a 2-D FFT array A, taking  care of "rolling" A  
     and setting  correct world boundaries.   Keyword SCALE can be  used to  
     indicate the  "frequel" scale along both  axis (SCALE is  a scalar) or  
     along each axis (SCALE  is a 2-element vector: SCALE=[XSCALE,YSCALE]);  
     by default, SCALE=[1.0, 1.0].  Other  keywords have same meaning as in  
     plc routine.  
   KEYWORDS scale, levs, type, width, color, smooth,  
            legend, hide, marks, marker, mspace, mphase.  

             fft_plfc, a;  
 
     Plot  filled contour  levels of  a  2-D FFT  array A,  taking care  of  
     "rolling" A  and setting correct world boundaries.   Keyword SCALE can  
     be used  to indicate the "frequel"  scale along both axis  (SCALE is a  
     scalar)   or  along   each  axis   (SCALE  is   a   2-element  vector:  
     SCALE=[XSCALE,YSCALE]); by default,  SCALE=[1.0, 1.0].  Other keywords  
     have  same meaning  as in  plfc routine.   As with  plfc  routine, the  
     actual level values get saved in external symbol plfc_levs.  
   KEYWORDS scale, levs, colors.  

 fft_plg  
 

             fft_plg, y;  
 
        -or fft_plh, y;  
     Plot 1-D FFT array Y as a curve, taking care of "rolling" Y and setting  
     correct coordinates.  Keyword SCALE can be used to indicate the  
     "frequel" scale along X-axis (SCALE is a scalar); by default,  
     SCALE=1.0.  
   KEYWORDS legend, hide, type, width, color, closed, smooth  
            marks, marker, mspace, mphase.  

             fft_pli, a;  
 
     Plot 2-D FFT array A as an image, taking care of "rolling" A and setting  
     correct world boundaries.  Keyword SCALE can be used to indicate the  
     "frequel" scale along both axis (SCALE is a scalar) or along each axis  
     (SCALE is a 2-element vector: SCALE=[XSCALE,YSCALE]); by default,  
     SCALE=[1.0, 1.0].  
   KEYWORDS legend, hide, top, cmin, cmax.  

             fft_recenter(x, template)  
 
       -or- fft_recenter(x, template, reverse)  
     Returns array X rolled so that it matches most closely array TEMPLATE.  
     X and TEMPLATE may be arrays  of real or complex numbers but must have  
     the same dimension lists.  The  returned value is roll(X, S) where the  
     offsets S minimize:  
         sum(abs(roll(X, S) - TEMPLATE)^2)  
     If optional  argument REVERSE is true,  X is also allowed  to have all  
     its dimensions reversed in order to math TEMPLATE.  
       

             fft_recenter_at_max(z)  
 
     
     Return Z  rolled so  that its element  with maximum value  (or maximum  
     absolute value if  Z is complex) is at the  origin.  If keyword MIDDLE  
     is true  (non-zero and non-nil) the  center is at the  middle of every  
     dimension  otherwise the  center is  the first  element of  the output  
     array (as assumed by the FFT).  

             fft_roll_1d(v, off)  
 
       -or- fft_roll_2d(m, off1, off2)  
     "rolls" dimensions of the vector V (1D array) or matrix M (2D array)  
     and return a result with same data type than original array.  
       

 fft_roll_2d  
 

             fft_shift_phasor(off, dimlist)  
 
       -or- fft_shift_phasor(off, fft_freqlist(dimlist))  
     returns complex phasor to apply in FFT space for a shift by OFF cells  
     in the real space.  DIMLIST is a list of dimensions -- the second  
     callinf sequence is to allow for computing the normalized FFT  
     frequencies only once.  The offset OFF must have as many elements as  
     PTR or as many as dimensions in DIMLIST (i.e. a shift for each  
     dimension) and may be fractionnal.  
     This function is intended for Fourier interpolation.  For instance,  
     assuming A is 2-D real array:  
         z = fft(a, +1);  
         u = fft_freqlist(dimsof(a));  
         q = fft_unphasor([0.33, -0.47], u);  
     Then:  
         fft_interp(z, q, real=1);  
     yields the value of A interpolated at coordinate (0.33, -0.47) in FFT  
     frame, i.e. center of lower left cell is at (0,0).  The shfited version  
     of A by (0.33, -0.47) can be obtained by:  
       
         fft_shift(z, q, real=1);  

             fft_smooth(a, fwhm)  
 
       -or- fft_smooth(a, fwhm, setup=workspace)  
     Returns array A smoothed along  all its dimensions by convolution by a  
     gaussian  with  full  width  at  half maximum  equals  to  FWHM.   See  
     fft_setup for the meaning of keyword SETUP.  

             fft_symmetric_index(dimlist)  
 
       -or- fft_symmetric_index(dim1, dim2, ...);    
     Returns  indices  of  hermitian-symmetry  transform  for  a  FFT  with  
     dimension list DIMLIST.  For instance,  if A is a N-dimensional array,  
     then:  
        AP= A(fft_symmetric_index(dimsof(A)))     
     is  equal to array  A with  its coordinates  negated according  to FFT  
     convention:  
        AP(X1, X2, ..., XN) = A(-X1, -X2, ..., -XN)  
     consequently if A is hermitian then:  
        AP= conj(A).  
       

             fft_fine_shift(a, off, setup=)  
 
       -or- fft_unphasor(z, phasor, setup=, real=)  
     The function fft_fine_shift returns array A shifted by offset OFF  
     which can be fractionnal.  Alternatively, the function fft_unphasor  
     can be used when the forward FFT of A and/or the complex phasor  
     corresponding to the shift are already computed:  
       fft_unphasor(fft(A, +1), fft_shift_phasor(OFF, dimsof(A)),  
                    real=(structof(A) != complex))  
     yields the same result as fft_fine_shift(A, OFF).  
     These functions can make use of pre-computed FFT workspace specified  
     by keyword SETUP (see fft_setup).  
   
     TO-DO: Improve code by not Fourier transforming along direction  
            where OFF is zero (or equal to an integer times the length  
            of the dimension).  
       

            : FFT utility routines in "fft_utils.i"  
 
     This package  is mainly written  to deal with the  particular indexing  
     rules in FFT transformed arrays.  The following routines are provided:  
       abs2         - squared absolute value.  
       fft_best_dim - get best dimension to compute FFT.  
       fft_centroid - get centroid in FFT arrays.  
       fft_convolve - compute discrete convolution thanks to FFT.  
       fft_dist     - compute length of FFT frequencies/coordinates.  
       fft_fine_shift, fft_unphasor - shift/roll an array by non-integer  
                      offset by means of Fourier interpolation.  
       fft_gaussian_mtf - compute Gaussian modulation transfer function.  
       fft_gaussian_psf - compute Gaussian point spread function.  
       fft_indgen   - generate index of FFT frequencies/coordinates.  
       fft_interp, fft_interp_complex, fft_interp_real - interpolate array  
                      at non-integer offsets.  
       fft_plc      - plot contours of 2D FFT.  
       fft_plfc     - plot filled contours of 2D FFT.  
       fft_plg      - plot 1D FFT as curve.  
       fft_plh      - plot 1D FFT with stairs.  
       fft_pli      - plot 2D FFT as image.  
       fft_recenter - recenter array with respect to a template.  
       fft_recenter_at_max - recenter FFT arrays at their maximum.  
       fft_roll_1d  - roll dimension of 1D arrays.  
       fft_roll_2d  - roll dimension of 2D arrays.  
       fft_shift_phasor - get complex phasor for arbitrary shift.  
       fft_smooth   - smooth an array by convolution with a gaussian.  
       fft_symmetric_index - get hermitian-symmetry index for FFT arrays.  
       reverse_dims - reverse all dimensions of an array.  
       __fft        - expert driver for repeated FFT's with same dimensions.  
       __fft_init   - initialization for __fft.  

             reverse_dims(a)  
 
     Returns array A with all its dimensions reversed.