Digitale bildgebende Verfahren: Transformationen/ FFT
Modul FFT
[Bearbeiten]Das folgende Programm in der Programmiersprache Component Pascal zur Berechnung von zweidimensionalen Fast-Fourier-Transformationen (FFT) kann mit dem Kommando FFT.ExampleCallCommand aufgerufen werden.
MODULE FFT;
(*
Title: Fast-Fourier-Transformation, FFT
Last edits: 4th January 2012 / layout: 12th January 2022
Author: Dr. Markus Bautsch, Berlin
Programming Language: Component Pascal
Reference: Elbert Oran, FFT-Anwendungen, Oldenbourg, München, Wien, 1997
*)
IMPORT
Math, Out;
TYPE
COMPLEX* = RECORD re*, im*: REAL END;
FFTArray* = POINTER TO ARRAY OF COMPLEX;
FFTArray2D* = POINTER TO ARRAY OF FFTArray;
CosSinArray = POINTER TO ARRAY OF REAL;
BitrevArray = POINTER TO ARRAY OF INTEGER;
VAR
Bitrev: BitrevArray;
Cos, Sin: CosSinArray;
N, NU: INTEGER;
(*
Amplitude computes and returns the absolute real value of the complex number c.
*)
PROCEDURE Amplitude* (c: COMPLEX): REAL;
VAR
cRe, cIm: REAL;
result: REAL;
BEGIN
cRe := c.re;
cIm := c.im;
result := Math.Sqrt (cRe * cRe + cIm * cIm);
RETURN result;
END Amplitude;
(*
InitArrays initilialises the global variables.
*)
PROCEDURE InitArrays ();
VAR
bitset: SET;
i, j, nu1: INTEGER;
pi2N, arg: REAL;
BEGIN
pi2N := 2 * Math.Pi () / N;
i := N;
REPEAT
DEC (i);
bitset := {};
nu1 := 0;
j := NU;
REPEAT
DEC (j);
IF j IN BITS (i) THEN INCL (bitset, nu1) END;
INC (nu1);
UNTIL j = 0;
Bitrev [i] := ORD (bitset);
arg := pi2N * i;
Cos [i] := Math.Cos (arg);
Sin [i] := Math.Sin (arg);
UNTIL i = 0;
END InitArrays;
(*
InitFFT prepares and checks the global variables. Finally it calls InitArrays in order to intilalise all global variables.
The parameter n gives the dimension of the FFT-Arrays measured as the power of 2.
*)
PROCEDURE InitFFT (n: INTEGER);
VAR
a: INTEGER;
BEGIN
ASSERT (n >= 2);
NU := 0;
a := n;
REPEAT
ASSERT ((a MOD 2) = 0);
a := a DIV 2;
INC (NU);
UNTIL a <= 1;
N := n;
NEW (Cos, N);
ASSERT (Cos # NIL);
NEW (Sin, N);
ASSERT (Sin # NIL);
NEW (Bitrev, N);
ASSERT (Bitrev # NIL);
InitArrays ();
END InitFFT;
(*
The procedure Sort carries out the bit reversal within an FFTArray H.
*)
PROCEDURE Sort (H: FFTArray);
VAR
i, l, k: INTEGER;
temp: COMPLEX;
BEGIN
(* bit reversal: *)
l := N - 1;
k := 0;
REPEAT
INC (k);
i := Bitrev [k];
IF i > k THEN
temp := H [k];
H [k] := H [i];
H [i] := temp;
END;
UNTIL k = l;
END Sort;
(*
FFT calculates the one-dimensional discrete fourier transformation H of a complex array h.
The result is stored in the same array h because of efficiency reasons.
The two halfs of the result are not exchanged, i.e. the constant offset level for the whole fourier transform is contained in h [0] !
*)
PROCEDURE FFT* (h: FFTArray);
VAR
i, l, k, kn2, n, n2, p: INTEGER;
hc: COMPLEX;
bitset: SET;
cos, sin, tempRe, tempIm, hcre, hcim: REAL;
BEGIN
ASSERT (h # NIL);
n := LEN (h);
IF N # n THEN InitFFT (n) END;
h [0].re := 0.5 * (h [0].re + h [N-1].re);
h [0].im := 0.5 * (h [0].im + h [N-1].im);
n2 := N DIV 2;
FOR l := NU - 1 TO 0 BY -1 DO
bitset := {l};
k := 0;
WHILE k < N DO
kn2 := k + n2;
FOR i := 0 TO n2 - 1 DO
p := Bitrev [k DIV ORD (bitset)];
cos := Cos [p];
sin := Sin [p];
hc := h [kn2];
hcre := hc.re;
hcim := hc.im;
tempRe := hcre * cos + hcim * sin;
tempIm := hcim * cos - hcre * sin;
hc := h [k];
hcre := hc.re;
hcim := hc.im;
h [kn2].re := hcre - tempRe;
h [kn2].im := hcim - tempIm;
h [k].re := hcre + tempRe;
h [k].im := hcim + tempIm;
INC (k);
INC (kn2);
END;
INC (k, n2);
END;
n2 := n2 DIV 2;
END;
Sort (h);
END FFT;
(*
FFT2D calculates the two-dimensional discrete fourier transformation H of an FFTArray2D h2.
The result is stored in the same array h2 because of efficiency reasons.
*)
PROCEDURE FFT2D* (h2: FFTArray2D);
VAR
nx, ny: INTEGER;
i, j: INTEGER;
ht: FFTArray;
BEGIN
ASSERT (h2 # NIL);
ny := LEN (h2);
ASSERT (ny >= 2);
ASSERT (h2 [0] # NIL);
nx := LEN (h2 [0]);
ASSERT (nx >= 2);
FOR j := 1 TO ny - 1 DO
ASSERT (h2 [j] # NIL);
ASSERT (LEN (h2 [j]) = nx);
END;
(* calc columns: *)
FOR j := 0 TO ny - 1 DO FFT (h2 [j]) END;
(* calc rows: *)
NEW (ht, ny);
ASSERT (ht # NIL);
FOR i := 0 TO nx - 1 DO
FOR j := 0 TO ny - 1 DO ht [j] := h2 [j, i] END;
FFT (ht);
FOR j := 0 TO ny - 1 DO h2 [j, i] := ht [j] END;
END;
END FFT2D;
(*
The command ExampleCallCommand is public and can be called by the run time system.
The two-dimensional FFT-Array is stored in Example2DArray.
*)
PROCEDURE ExampleCallCommand*;
VAR
n, i, j: INTEGER;
Example2DArray: FFTArray2D;
BEGIN
n:= 8; (* 2^8 = 256 elements in one dimension *)
(* Build two-dimensional FFT-Array *)
NEW (Example2DArray, n);
FOR j := 0 TO n - 1 DO
NEW (Example2DArray [j], n);
FOR i := 0 TO n-1 DO
Example2DArray [j, i].re := Math.Sin (2 * Math.Pi () * i / n) * Math.Sin (4 * Math.Pi () * j / n) - 0.5;
Example2DArray [j, i].im := 0;
END;
END;
(* Transform two-dimensional FFT-Array *)
FFT2D (Example2DArray);
(* Output transformed two-dimensional FFT-Array *)
FOR j := 0 TO n - 1 DO
FOR i := 0 TO n - 1 DO
Out.Real (Amplitude (Example2DArray [j, i]), 20);
Out.String (" ");
END;
Out.Ln;
END;
END ExampleCallCommand;
BEGIN
(*
Initialisation of the global variables when this module is loaded by the run time system.
*)
N := 0;
NU := 0;
Sin := NIL;
Cos := NIL;
Bitrev := NIL;
END FFT.
(*
Result table of fftArray2D (values rounded to one decimal place) after the call of FFT.ExampleCallCommand
(the four sub blocks have to be exchanged diagonally row 0 and column 0 have been removed):
0.2 0.2 2.0 0.2 2.0 0.2 0.2
1.4 1.4 16.2 1.4 16.2 1.4 1.4
0.2 0.2 2.0 0.2 2.0 0.2 0.2
0.2 0.2 2.0 31.8 2.0 0.2 0.2
0.2 0.2 2.0 0.2 2.0 0.2 0.2
1.4 1.4 16.2 1.4 16.2 1.4 1.4
0.2 0.2 2.0 0.2 2.0 0.2 0.2
*)
Ausgabe
[Bearbeiten]Ausgabe der vom Programm berechneten Werte (die gedoppelten Werte sind durchgestrichen und werden nach der Umsortierung nicht weiter verwendet):
31.82322330470336 2.007797300526125 0.176776695296637 0.17677669529663710.17677669529663590.1767766952966373 0.176776695296637 2.007797300526125 0.176776695296637 2.007797300526125 0.1767766952966368 0.17677669529663660.17677669529663680.1767766952966369 0.1767766952966368 2.007797300526125 1.425219281373923 16.18737934317967 1.425219281373922 1.4252192813739221.4252192813739221.425219281373922 1.425219281373922 16.18737934317967 0.176776695296637 2.007797300526125 0.1767766952966368 0.17677669529663690.17677669529663680.1767766952966367 0.1767766952966368 2.0077973005261250.1767766952966365 2.007797300526125 0.176776695296637 0.1767766952966366 0.1767766952966376 0.1767766952966366 0.176776695296637 2.0077973005261250.176776695296637 2.007797300526125 0.1767766952966368 0.17677669529663660.17677669529663680.1767766952966369 0.1767766952966368 2.007797300526125 1.425219281373923 16.18737934317967 1.425219281373922 1.4252192813739221.4252192813739221.425219281373923 1.425219281373922 16.18737934317967 0.176776695296637 2.007797300526125 0.1767766952966368 0.17677669529663690.17677669529663680.1767766952966367 0.1767766952966368 2.007797300526125