UNIT U_FTG2;                           {Last mod by JFH on 03/18/96}

(* *****************************************************************************

 UNIT FOR PERFORMING GLASSMAN FFT AND INVERSE

 Note that with the Glassman FFT, the number of points does not have to be a
     power of 2.  Instead, the transform will take advantage of any factoring
     that can be done on the number of points.  Which means that if the number
     of points is prime, the transform will not be "fast" at all.

 Rev. 08/30/00 by JFH for Delphi-5 by going to dynamic arrays required no
     change other than the removal of tComplexArray.  Also fixed error in
     CAScale and removed unused variables in GFTInnerLoop.
 Rev. 03/18/96 by JFH for Delphi-2 by removing def of "PI" and going
     to open arrays.
 Pgm. 03/27/89 by J F Herbster

***************************************************************************** *)

INTERFACE

TYPE
  tComplex = record r,i: single;
  end;
  tFactorArray = array [1..32] of integer;

Procedure CxAdd (var sum: tComplex; a,b: tComplex);

Procedure CxMpy (var prd: tComplex; a,b: tComplex);

Procedure CAScale
   (const Inp: array of tComplex; var aOut: array of tComplex;
    A: single; N: integer);

Procedure GetFactors
    (n: integer; var factors: tFactorArray; var nf: integer);

Procedure CFTG
   (const Cinp: array of tComplex; var Cout, CScr: array of tComplex;
    NC: integer);
{ returns a complex Fourier transform of the Cinp array into the Cout array.
  The CScr (scratch) array may be same as the input array. The transform will
  be from the time to the frequency domain if NC is positive of the number
  of complex points to be converted or from the frequency to the time domain
  if NC is the negative of the number of complex points.  The only difference
  between the two is just a divide by the number of points. }

IMPLEMENTATION {===============================================================}

Procedure CxAdd (var sum: tComplex; a,b: tComplex);
    begin sum.r:=a.r+b.r; sum.i:=a.i+b.i end;

Procedure CxMpy (var prd: tComplex; a,b: tComplex);
    begin prd.r:=a.r*b.r-a.i*b.i; prd.i:=a.i*b.r+a.r*b.i end;

Procedure CAScale
   (const Inp: array of tComplex; var aOut: array of tComplex;
    A: single; N: integer);
Var i: word;
Begin
  For i := N-1 {JFH, 08/30/00} downto 0 do
       begin aOut[i].r:=A*Inp[i].r; aOut[i].i:=A*Inp[i].i; end;
End;

Procedure GetFactors(n: integer; var factors: tFactorArray; var nf: integer);
{ Factor n into its nf primes in array factors[]. }
  Var nr,tf,q: word;
Begin
  nr:=n; q:=nr; tf:=2; nf:=1;
  While tf<q do begin
      q:=nr div tf;
      if q*tf=nr
        then begin factors[nf]:=tf; Inc(nf); nr:=q end
        else Inc(tf);
    end{while};
    factors[nf]:=nr;
  End;

Procedure GFTInnerLoop
   (const Inp: array of tComplex; var aOut: array of tComplex;
    NC: integer; X, S: integer);
Var WDL1,WD,CTemp,CTemp2: tComplex;  N,L1,L,K,D,I,U,H: word; a: single;
Begin
  N:=abs(NC); D:=N div S; I:=0;
  a:=(-(2*Pi)*D)/NC; WD.r:=cos(a); WD.i:=sin(a);
  For L1:=1 to S do
    begin
      If L1=1
        then begin WDL1.r := 1.0; WDL1.i := 0.0 end
        else CxMpy(WDL1,WDL1,WD);
      For L:=1 to D do
        begin
          K:= (L+(L1-1)*D*X) mod N;
          U:=K+(X-1)*D;
          CTemp:=Inp[U-1];
          For H:=2 to X do
            begin
              CxMpy (CTemp,CTemp,WDL1);
              U:=U-D; CTemp2:=Inp[U-1];
              CxAdd (CTemp,CTemp,CTemp2);
            end;
          aOut[i]:=CTemp;
          Inc(I);
        end{L-loop};
    end{L1-loop};
End;

Procedure CFTG
   (const Cinp: array of tComplex; var Cout, CScr: array of tComplex;
    NC: integer);
{ CScr array may be same as the Cinp array.}
Var  a: single; j,n,x,s,nf: integer; r: tFactorArray;
Begin

{ Let n be the number of complex elements in the array.}
  n:=abs(NC);
  If n=0
    then exit;

  GetFactors(n,r,nf);          { Factor n into its primes.}

{ Now do the transformation.}
  s:=1;
  For j:=1 to nf do
    begin
      x:=r[j]; s:=s*x;
      If j=1
        then GFTInnerLoop (CInp, COut, NC, X, S)
        else if Odd(j)
          then GFTInnerLoop (CScr, COut, NC, X, S)
          else GFTInnerLoop (COut, CScr, NC, X, S);
    end{j-loop};

{ Finally, check for the required scaling and that data is in out array.}
  If nc < 0 then a:=1.0 else a:=1.0/n;
  If not Odd(nf)
    then CAScale({in}CScr, {out}COut, a, N)
    else if a<>1
      then CAScale(COut, COut, a, N);

End;

END.