9.7 Unit FFTU

{-------------------------------------------------------------------}

{ Eindimensionale Fast Fourier Transformation ohne Coprozessor }
unit FFTU;

{$N+,E+}

{ Achtung : die Unit FFTU benoetigt waehrend der Kompilation die
            Datei BITREV.OBJ ! }

{-------------------------------------------------------------------}

interface

uses
  Glob;     { Globale Deklarationen }

procedure Tabelle(n : integer);  { Muss vor dem ersten Aufruf von
                                   FFT aufgerufen werden ! }
procedure FFT(var XReal, XImag : Realarray; n, nu : integer;
              hin : boolean);

{-------------------------------------------------------------------}

implementation

type
  RealTypearray = array [0..nmax - 1] of single;

var
  i, q1, q2, q3, q4 : integer;
  Tabarray          : ^RealTypearray;

{-------------------------------------------------------------------}

{$L BITREV}
{$F-}

function BitRev(m, nu : integer): integer; external;

{$F+}

{-------------------------------------------------------------------}

procedure Tabelle(n : integer);

var
  i, nmax : integer;
  s       : single;

begin
  nmax := n div 4;
  s := 2.0 * Pi / n;
  for i := 0 to nmax do Tabarray^[i] := sin(s * i);
  q1 := n div 4;
  q2:= n div 2;
  q3 := q1 + q2;
  q4 := n;
end; { Tabelle }

{-------------------------------------------------------------------}

function Sinus(x : integer): single;

begin
  if x <= q1 then Sinus :=  Tabarray^[x]
  else if x <= q2 then Sinus :=  Tabarray^[q2 - x]
       else if x <= q3 then Sinus :=  - Tabarray^[x - q2]
            else Sinus :=  - Tabarray^[q4 - x]
end; { Sinus }

{-------------------------------------------------------------------}

function Cosinus(x : integer): single;

begin
  if x <= q1 then Cosinus :=  Tabarray^[q1 - x]
  else if x <= q2 then Cosinus := - Tabarray^[x - q1]
       else if x <= q3 then Cosinus := - Tabarray^[q3 - x]
            else Cosinus := Tabarray^[x - q3]
end; { Cosinus }

{-------------------------------------------------------------------}

procedure FFT(var XReal, XImag : Realarray; n, nu : integer;
              hin : boolean);

type
  PointRec = record
    Low, High : word;
  end;

var
  n2, kn2, nu1, i, l, k, k0, k0step, p, pstep :	integer;
  FFTfilename                                 : string [25];
  f                                           : Text;
  TReal, TImag, c, s, rn2, FFTNorm            : single;
  RealPtr1, ImagPtr1, RealPtr2, ImagPtr2      :	^single;

  RRealPtr1 : PointRec absolute RealPtr1;
  RRealPtr2 : PointRec absolute RealPtr2;
  RImagPtr1 : PointRec absolute ImagPtr1;
  RImagPtr2 : PointRec absolute ImagPtr2;

begin
  k := 0;

  while (k < n) do
  begin
    i := BitRev(k, nu);
    if (i > k) then
    begin
      TReal := XReal[k];
      TImag := XImag[k];
      XReal[k] := XReal[i];
      XImag[k] := XImag[i];
      XReal[i] := TReal;
      XImag[i] := TImag
    end;
    inc(k)
  end;

  n2 := 1;
  nu1 := 1;
  k0step := 2;
  pstep := n shr 1;

  for l := 1 to nu do
  begin
    k0 := 0;

    repeat
      k := k0;
      kn2 := k0 + n2;
      p := 0;
      RealPtr1 := @XReal[kn2];
      ImagPtr1 := @XImag[kn2];
      RealPtr2 := @XReal[k];
      ImagPtr2 := @XImag[k];
      for i := 1 to n2 do
      begin
        c := Cosinus(p);
        if hin then s := - Sinus(p) else s := Sinus(p);

        TReal := RealPtr1^ * c + ImagPtr1^ * s;
        TImag := ImagPtr1^ * c - RealPtr1^ * s;

        RealPtr1^ := RealPtr2^ - TReal;
        ImagPtr1^ := ImagPtr2^ - TImag;
        RealPtr2^ := RealPtr2^ + TReal;
        ImagPtr2^ := ImagPtr2^ + TImag;

        inc(RRealPtr1.Low, sizeof(Single));
        inc(RRealPtr2.Low, sizeof(Single));
        inc(RImagPtr1.Low, sizeof(Single));
        inc(RImagPtr2.Low, sizeof(Single));

        inc(p, pstep)
      end;
      inc(k0, k0step)
    until (k0 = n);

    k0step := k0step + k0step;
    inc(nu1);
    n2 := n2 + n2;
    pstep := pstep shr 1
  end;

  if hin then
  begin
    rn2 := 2 / n;
    for i := 0 to n - 1 do
    begin
      XReal[i] := rn2 * XReal[i];
      XImag[i] := rn2 * XImag[i];
    end;
  end
  else
  begin
    for i := 0 to n - 1 do
    begin
      XReal[i] := 0.5 * XReal[i];
      XImag[i] := 0.5 * XImag[i];
    end;
  end;
end; { FFT }

{-------------------------------------------------------------------}

end. { FFTU }

{-------------------------------------------------------------------}