IMPLEMENTATION MODULE Fourier;
(*------------------------------------------------------------------------*)
(* Das Modul enthaelt einige Routinen fuer Fourtransforamtionen *)
(* This Module provides some routines for fourier transformations *)
(*------------------------------------------------------------------------*)
(* Letzte Aenderung: *)
(* *)
(* 30.05.98, MRi: Durchsicht *)
(* 13.10.15, MRi: Ersetzen von LFLOAT(#) durch VAL(LONGREAL,#) *)
(* 02.10.16, MRi: Umsellen auf ISO Modules Complex Typ *)
(* Prozedur four1 geloescht *)
(* 13.10.16, MRi: Umbenennen der Prozeduren, einfuegen von CfstFT *)
(* 03.03.18, MRi: Anpassen der Schnittstelle von CfstFT an CmplxFFT, *)
(* Vereinfachung in RealDFT. *)
(* 04.03.18, MRi: In RealDFT den Parameter "rev" eingefuegt, Prozedur *)
(* CmplxFT eingefuegt *)
(* 05.03.18, MRi: Norming the result implemented in all procedures exept *)
(* CmplxFFT where is was already present *)
(* 12.03.18, MRi: Erstellen der ersten Fortran-Version von GlFFT nach *)
(* Literaturvorgabe *)
(* 13.03.18, MRi: Umstellen von GlFFT auf linearen Index und Uebersetzung *)
(* nach M2, Normierung eingefuegt *)
(*------------------------------------------------------------------------*)
(* Procedure CmplxFFT is based on a Pascal version provided by *)
(* Hans-J\"urgen Hochkamp. He stated that his on is in parts based on a *)
(* routine he found in the magazin C'T (Heise Verlag Hannover). I checked *)
(* the listings provided at on FTP server of Heise but could find a *)
(* similar routine in that repository - so I suppose we cannot figure out *)
(* anymore where it originally comes from. *)
(* *)
(* Procedure CfstFT is a Modula-2 translation of CERN Lib CFSTFT *)
(* provided by K.S. K\"olbig and H.-H. Umst\"atter *)
(*------------------------------------------------------------------------*)
(* Offene Punkte *)
(* *)
(* - Eine Inter zu Interger FFT einfuegen ... z.B. *)
(* Monro, Donald M.; ���A portable integer FFT in FORTRAN���, Computer *)
(* Programs in Biomedicine, 7-4, pp. 267-272 (1977) *)
(*------------------------------------------------------------------------*)
(* Testroutinen in TstFFT.mod *)
(*------------------------------------------------------------------------*)
(* Implementation : Michael Riedl *)
(* Licence : GNU Lesser General Public License (LGPL) *)
(*------------------------------------------------------------------------*)
(* $Id: Fourier.mod,v 1.5 2018/03/05 11:21:46 mriedl Exp mriedl $ *)
FROM SYSTEM IMPORT TSIZE;
FROM Storage IMPORT ALLOCATE,DEALLOCATE;
FROM LongMath IMPORT pi,sqrt,sin,cos;
FROM LongComplexMath IMPORT zero,one,scalarMult,conj,abs;
FROM LMathLib IMPORT Log2;
PROCEDURE RealDFT(VAR A,B : ARRAY OF LONGREAL; (* Koeff. *)
VAR F : ARRAY OF LONGREAL; (* Werte *)
N : CARDINAL; (* Anzahl Werte *)
rev : BOOLEAN);
VAR k,t : CARDINAL;
rezN,tmp,theta : LONGREAL;
Fk,Ak,Bk : LONGREAL;
BEGIN
rezN := 1.0 / LFLOAT(N);
FOR k:=0 TO N-1 DO
tmp := 2.0*pi*LFLOAT(k)*rezN;
IF NOT rev THEN
Ak:=0.0; Bk:=0.0;
FOR t:=0 TO N-1 DO
theta := tmp*LFLOAT(t);
Ak := Ak + F[t]*cos(theta);
Bk := Bk - F[t]*sin(theta);
END;
A[k] := Ak;
B[k] := Bk;
ELSE
Fk:=0.0;
FOR t:=0 TO N-1 DO
theta := tmp*LFLOAT(t);
Fk:=Fk + abs( CMPLX(A[t],B[t]) * CMPLX(cos(theta), sin(theta)) );
END;
F[k] := Fk;
END;
END;
IF NOT rev THEN (* Norm *)
FOR k:=0 TO N-1 DO A[k]:=rezN*A[k]; B[k]:=rezN*B[k]; END;
END;
END RealDFT;
PROCEDURE CmplxFT(VAR X : ARRAY OF LONGCOMPLEX;
VAR Y : ARRAY OF LONGCOMPLEX;
N : CARDINAL;
rev : BOOLEAN);
VAR i,j : CARDINAL;
rezN,tmp,theta : LONGREAL;
val : LONGCOMPLEX;
BEGIN
rezN := 1.0 / LFLOAT(N);
FOR i:=0 TO N-1 DO
val := zero;
tmp := - 2.0*pi*LFLOAT(i)*rezN;
IF NOT rev THEN
FOR j:=0 TO N-1 DO
theta := LFLOAT(j)*tmp;
val:=val + X[j]*CMPLX(cos(theta), -sin(theta));
END;
Y[i] := val;
ELSE
FOR j:=0 TO N-1 DO
theta := LFLOAT(j)*tmp;
val := val + Y[j]*CMPLX(cos(theta), sin(theta));
END;
X[i] := scalarMult(rezN,val);
END;
END;
IF NOT rev THEN (* Norm *)
FOR i:=0 TO N-1 DO Y[i]:=scalarMult(rezN,Y[i]); END;
ELSE
FOR i:=0 TO N-1 DO X[i]:=scalarMult(LFLOAT(N),X[i]); END;
END;
END CmplxFT;
PROCEDURE CmplxFFT(VAR X : ARRAY OF LONGCOMPLEX;
N : CARDINAL;
rev : BOOLEAN);
VAR a,norm : LONGREAL;
t,w : LONGCOMPLEX;
i,j,l,m,mr,n1,nh,istep : CARDINAL;
BEGIN
mr := 0; (* Initialisierung der Variablen *)
n1 := N - 1;
nh := N DIV 2;
FOR m:=1 TO n1 DO (* Bit-Reverse *)
l := nh;
WHILE ((mr + l) > n1) DO l := l DIV 2; END;
mr := (mr MOD l) + l;
IF (mr > m) THEN (* Vertauschen *)
t := X[m]; X[m] := X[mr]; X[mr] := t;
END;
END;
(*
j:=0;
FOR i:=0 TO N-2 DO
IF (i < j) THEN t := X[j]; X[j] := X[i]; X[i] := t; END;
l := nh;
WHILE (l <= j) DO j:=j - l; l:=l DIV 2; END;
INC(j,l);
END;
*)
l:=1;
WHILE (l < N) DO (* fft *)
istep := 2*l;
FOR m:=1 TO l DO
a := - pi*LFLOAT(m-1) / LFLOAT(l);
w := CMPLX(cos(a),sin(a));
IF rev THEN w := conj(w); END;
i := m - 1;
REPEAT
j := i + l;
t := w*X[j];
X[j] := X[i] - t;
X[i] := X[i] + t;
INC(i,istep);
UNTIL (i >= N);
END;
l := istep;
END;
IF NOT rev THEN (* Normierung *)
norm := 1.0 / LFLOAT(N);
FOR i:=0 TO N-1 DO X[i]:=scalarMult(norm,X[i]); END;
END;
END CmplxFFT;
PROCEDURE CfstFT(VAR X : ARRAY OF LONGCOMPLEX;
N : CARDINAL;
rev : BOOLEAN);
VAR c,s,norm : LONGREAL;
i,j,k,l,le,le1,m : CARDINAL;
u,w,t : LONGCOMPLEX;
BEGIN
IF (N # 0) THEN
m := Log2(N);
j := 0;
FOR i:=0 TO N-2 DO
IF (i < j) THEN
t := X[j];
X[j] := X[i];
X[i] := t;
END;
k := N DIV 2;
WHILE (k <= j) DO DEC(j,k); k:=k DIV 2; END;
INC(j,k);
END;
FOR i:=0 TO N-1 BY 2 DO
t := X[i+1];
X[i+1] := X[i] - t;
X[i+0] := X[i] + t;
END;
c := 0.0;
IF rev THEN s:=-1.0; ELSE s:=1.0; END;
le := 2;
FOR l:=2 TO m DO
w := CMPLX(c,s);
u := w;
c := sqrt(0.5*c + 0.5);
s := IM(w) / (c + c);
le1 := le;
le := le1 + le1;
i:=0;
WHILE (i < N) DO
t := X[i+le1];
X[i+le1] := X[i] - t;
X[i] := X[i] + t;
INC(i,le);
END;
FOR j:=2 TO le1 DO
i:=j-1;
WHILE (i <= N-1) DO
t := X[i+le1]*u;
X[i+le1] := X[i] - t;
X[i ] := X[i] + t;
INC(i,le);
END;
u:=u * w;
END;
END; (* FOR l *)
END; (* IF *)
IF NOT rev THEN (* Norm *)
norm := 1.0 / LFLOAT(N);
FOR i:=0 TO N-1 DO X[i]:=scalarMult(norm,X[i]); END;
END;
END CfstFT;
PROCEDURE GlFFT(VAR U : ARRAY OF LONGCOMPLEX;
N : CARDINAL;
Invrs : BOOLEAN);
PROCEDURE Glassman( nA,nB,nC : CARDINAL;
VAR Uin : ARRAY OF LONGCOMPLEX; (* CONST *)
VAR Uout : ARRAY OF LONGCOMPLEX;
Invrs : BOOLEAN);
(*------------------------------------------------------*)
(* In this code I "linearized" the three dimensional *)
(* array found in the original code *)
(*------------------------------------------------------*)
VAR angel : LONGREAL;
ia,ib,ic,jc,jcr : CARDINAL;
delta,omega,sum : LONGCOMPLEX;
BEGIN
angel := 2.0*pi / LFLOAT(nA*nC);
delta := CMPLX(cos(angel),-sin(angel));
IF Invrs THEN
delta := conj(delta);
END;
omega:=one;
FOR ic:=0 TO nC-1 DO
FOR ia:=0 TO nA-1 DO
FOR ib:=0 TO nB-1 DO
(* Uin[ib,nC,ia] *)
sum := Uin[ib + ((nC-1) + ia*nC)*nB];
FOR jcr:=2 TO nC DO
jc := nC - jcr;
(* Uin[ib,jc,ia] *)
sum:=sum*omega + Uin[ib + (jc + ia*nC)*nB];
END;
(* Uout[ib,ia,ic] *)
Uout[ib + (ia + ic*nA)*nB] := sum;
END; (* FOR ib *)
omega:=omega*delta;
END; (* FOR ia *)
END; (* FOR ic *)
END Glassman;
VAR nA,nB,nC,i : CARDINAL;
inu : BOOLEAN;
Work : POINTER TO ARRAY [0..MAX(INTEGER)-1] OF LONGCOMPLEX;
Norm : LONGREAL;
BEGIN
ALLOCATE(Work,N*TSIZE(LONGCOMPLEX));
nA:=1; nB:=N; nC:=1; inu:=TRUE;
WHILE (nB > 1) DO
nA:=nC*nA;
nC := 1;
REPEAT INC(nC); UNTIL ((nB MOD nC) = 0);
nB:=nB DIV nC;
IF (inu) THEN
Glassman(nA,nB,nC,U,Work^,Invrs);
ELSE
Glassman(nA,nB,nC,Work^,U,Invrs);
END;
inu := NOT inu;
END; (* WHILE *)
Norm := 1.0 / LFLOAT(N);
IF inu AND NOT Invrs THEN (* Normiere *)
FOR i:=0 TO N-1 DO U[i]:=scalarMult(Norm,U[i]); END;
ELSE
IF Invrs THEN
FOR i:=0 TO N-1 DO U[i]:=Work^[i]; END;
ELSE (* Normiere *)
FOR i:=0 TO N-1 DO U[i]:=scalarMult(Norm,Work^[i]); END;
END;
END;
DEALLOCATE(Work,N*TSIZE(LONGCOMPLEX));
END GlFFT;
END Fourier.