DEFINITION MODULE Fourier;

  (*------------------------------------------------------------------------*)
  (* Module provides some routines for Fourier transformations              *)
  (*                                                                        *)
  (* It might be a good idea to have a look on Peter Moylans Fourier module *)
  (* as well.                                                               *)
  (*------------------------------------------------------------------------*)
  (* Implementation : Michael Riedl                                         *)
  (* Licence        : GNU Lesser General Public License (LGPL)              *)
  (*------------------------------------------------------------------------*)

  (* $Id: Fourier.def,v 1.3 2018/03/05 11:21:42 mriedl Exp mriedl $ *)

PROCEDURE RealDFT(VAR A,B : ARRAY OF LONGREAL; (* Koeff.       *)
                  VAR F   : ARRAY OF LONGREAL; (* Werte        *)
                      N   : CARDINAL;          (* Anzahl Werte *)
                      rev : BOOLEAN);

          (*----------------------------------------------------------------*)
          (* Berechnet die diskreten Fourierkoeffizenten der in F "uber-    *)
          (* gebenen Funktionswerte in den Felder A und B. Wenn rev = wahr  *)
          (* wird aud A,B der Vektor F zurueckberechnet                     *)
          (*                                                                *)
          (* Calculates in A & B the discret fourier coefficients of the    *)
          (* function values given in array F, if rev = true the reverse    *)
          (* calculation is performed.                                      *)
          (*                                                                *)
          (* This routine is mainly intended for test purpose               *)
          (*----------------------------------------------------------------*)

PROCEDURE CmplxFT(VAR X   : ARRAY OF LONGCOMPLEX;
                  VAR Y   : ARRAY OF LONGCOMPLEX;
                      N   : CARDINAL;
                      rev : BOOLEAN);

          (*----------------------------------------------------------------*)
          (* If rev = false perform a forward complex Fourier transform     *)
          (* as given by:                                                   *)
          (*                                                                *)
          (*   Y_i = \sum_{j=1}^N  X_j \exp{{2 \pi I J} \over N}            *)
          (*                                                                *)
          (*   \forall i \in [1..N]                                         *)
          (*                                                                *)
          (* If rev = true perform a backward complex Fourier transform     *)
          (* as given by:                                                   *)
          (*                                                                *)
          (*   X_i = {1 \over N} \sum_{j=1}^N Y_j \exp{{2 \pi I J} \over N} *)
          (*                                                                *)
          (*   \forall i \in [1..N]                                         *)
          (*                                                                *)
          (* This routine is mainly intended for test purpose               *)
          (*----------------------------------------------------------------*)

PROCEDURE CmplxFFT(VAR X   : ARRAY OF LONGCOMPLEX;
                       N   : CARDINAL;
                       rev : BOOLEAN);

          (*----------------------------------------------------------------*)
          (* Komplex nach komplex FFT                                       *)
          (*                                                                *)
          (* X    : Array of data points (real parts) for analyse as input, *)
          (*        overwritten by complex fourier cofficents on output.    *)
          (* X    : Array of complex fourier coefficiens for synthesis as   *)
          (*        input, overwritten by re-formed data points (real part) *)
          (*        on outpt.                                               *)
          (* n    : Number of data points in X, must be a power of two      *)
          (* rev  : true  : Analyse (analysis)                              *)
          (*        false : Synthese (synthesis)                            *)
          (* For more explanaions on analysis and synthesis see remarks for *)
          (* procedure CfstFT below as well.                                *)
          (*                                                                *)
          (* This routine seems to be the most precise one in this library  *)
          (*----------------------------------------------------------------*)

PROCEDURE CfstFT(VAR X   : ARRAY OF LONGCOMPLEX;
                     N   : CARDINAL;
                     rev : BOOLEAN);

          (*----------------------------------------------------------------*)
          (*                                                                *)
          (* Calculates the finite Fourier transform of a complex periodic  *)
          (* sequence X[0:n-1] whose period n must be a power of two.       *)
          (*                                                                *)
          (* Either the direct transform (eq. 1)                            *)
          (*                                                                *)
          (*   $X_j = \sum_{k=0}^{n-1} X_k exp({{-i 2\pi jk} \over n})$     *)
          (*                                                                *)
          (* with                                                           *)
          (*                                                                *)
          (*   $ j=0,1,\dots,n-1 $                                          *)
          (*                                                                *)
          (* or the unscaled inverse transform (eq. 2)                      *)
          (*                                                                *)
          (*   $\alpha_k = \sum_{j=0}^{n-1} X_j exp({{i2\pi jk} \over n})$  *)
          (*                                                                *)
          (* with                                                           *)
          (*                                                                *)
          (*   $ k=0,1,\dots,n-1 $                                          *)
          (*                                                                *)
          (* where X_k, \alpha_k and X_j are complex numbers, may be        *)
          (* calculated. To ensure optimum use of storage, the same array   *)
          (* is used for input and output and all intermediate calculations *)
          (* are carried out in this array.                                 *)
          (*                                                                *)
          (*   N   : the number of values in array X, must be a power of 2  *)
          (*   X   : One dimensional array of dimension (0:d),              *)
          (*         where d >= n-1                                         *)
          (*         For rev = true:                                        *)
          (*           On entry, X(k) = X_k (k=0,n-1)                       *)
          (*           On exit,  X(j) = X_j (j=0,n-1), as defined by (1).   *)
          (*         For rev = false:                                       *)
          (*           On entry, X(j) = X_j (j=0,n-1), as defined by (2).   *)
          (*           On exit , X(k) = X_k (k=0,n-1)                       *)
          (*   rev : true  : The direct  transform is calculated.           *)
          (*         false : The inverse transform is calculated.           *)
          (*                                                                *)
          (* This procedure is a Modula-2 translation of CERN Lib CFSTFT    *)
          (* provided by K.S. K\"olbig and H.-H. Umst\"atter                *)
          (*                                                                *)
          (* This routine is about a factor two faster than CmplxFFT        *)
          (*----------------------------------------------------------------*)

PROCEDURE GlFFT(VAR U     : ARRAY OF LONGCOMPLEX;
                    N     : CARDINAL;
                    Invrs : BOOLEAN);

          (*----------------------------------------------------------------*)
          (* Glassmans General N FFT                                        *)
          (*                                                                *)
          (* Berechne die Fouriertransformation des Feldes U oder die       *)
          (* inverse Fouriertransformierte (Invrs=wahr). Die Anzahl der     *)
          (* Datenpunkte muss keine Potenz von 2 sein.                      *)
          (*                                                                *)
          (* Calculates the Fourier transform of U or the inverse Fourier   *)
          (* transformaiton if Invrs = true. Note that the length N does    *)
          (* not need to be a power of 2 !                                  *)
          (*                                                                *)
          (* Input                                                          *)
          (*                                                                *)
          (*  U     : N-vector to be transformed                            *)
          (*  N     : length of vector to be transformed                    *)
          (*  Invrs : false if the forward transform is desired             *)
          (*          true  if the inverse transform is desired             *)
          (*                                                                *)
          (* Output                                                         *)
          (*                                                                *)
          (*  U     : the DFT of U if Invrs is false or N times             *)
          (*          the inverse dft of U if Invrs is true                 *)
          (*                                                                *)
          (* This routine is nearly as fast as CfstFM but less exact, but   *)
          (* the difference in precison should be absolutly OK for normal   *)
          (* usage. It differs from the "usual" Cooley und Tukey algorithm  *)
          (* based routines.                                                *)
          (*                                                                *)
          (* It might be worth the exercise to improve the code !           *)
          (*                                                                *)
          (* This routine is directly derived form the Fortran code         *)
          (* provided in the reference. A "lineraized" Fortran 90 version   *)
          (* can also be provided by the author of this library             *)
          (*----------------------------------------------------------------*)
          (* Ref: Ferguson, Warren E. jr.; "A simple derivation of          *)
          (*      Glassman's general n fast Fourier transform",             *)
          (*      Comp. & Math. with Applc., 8/6, pp 401-311 (1982)         *)
          (*----------------------------------------------------------------*)

END Fourier.
