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.