unit Un_FFT;

interface

 const
  TNArraySize = 1023;

 type
  Float       = real;
  TNvector    = array[0..TNArraySize] of Float;
  TNvectorPtr = ^TNvector;
  Analyses = (RF, RN, RC, RA, CF, CN, CC, CA);

 var
  WhichAnalysis : Analyses;         { Indicates which application  }
                                    { will be run                  }
  NumPoints : integer;              { Number of points  }
  XReal, XImag : TNvectorPtr;       { One set of complex data points  }
  HReal, HImag : TNvectorPtr;       { Another set of complex data points  }
  Inverse : boolean;                { False ==> forward Fourier transform  }
                                    { True ===> inverse Fourier transform  }
  Auto : boolean;                   { False ==> crosscorrelation  }
                                    { True ===> autocorrelation  }
  Error : byte;                     { Flags if something went wrong  }
  OutFile: text;

procedure TestInput(NumPoints    : integer;
                var NumberOfBits : byte;
                var Error        : byte);

{-------------------------------------------------------------}
{- Input: NumPoints                                          -}
{- Output: NumberOfBits, Error                               -}
{-                                                           -}
{- This procedure checks the input.  If the number of points -}
{- (NumPoints) is less than two or is not a multiple of two  -}
{- then an error is returned.  NumberOfBits is the number of -}
{- bits necessary to represent NumPoints in binary (e.g. if  -}
{- NumPoints = 16, NumberOfBits = 4).                        -}
{-------------------------------------------------------------}

procedure MakeSinCosTable(NumPoints : integer;
                      var SinTable  : TNvectorPtr;
                      var CosTable  : TNvectorPtr);

{--------------------------------------------------------}
{- Input: NumPoints                                     -}
{- Output: SinTable, CosTable                           -}
{-                                                      -}
{- This procedure fills in a table with sin and cosine  -}
{- values.  It is faster to pull data out of this       -}
{- table than it is to calculate the sines and cosines. -}
{--------------------------------------------------------}

procedure FFT(NumberOfBits : byte;
              NumPoints    : integer;
              Inverse      : boolean;
          var XReal        : TNvectorPtr;
          var XImag        : TNvectorPtr;
          var SinTable     : TNvectorPtr;
          var CosTable     : TNvectorPtr);

{-----------------------------------------------------}
{- Input: NumberOfBits, NumPoints, Inverse, XReal,   -}
{-        XImag, SinTable, CosTable                  -}
{- Output: XReal, XImag                              -}
{-                                                   -}
{- This procedure implements the actual fast Fourier -}
{- transform routine.  The vector X, which must be   -}
{- entered in bit-inverted order, is transformed in  -}
{- place.  The transformation uses the Cooley-Tukey  -}
{- algorithm.                                        -}
{-----------------------------------------------------}

procedure RealFFT(NumPoints : integer;
                  Inverse   : boolean;
              var XReal     : TNvectorPtr;
              var XImag     : TNvectorPtr;
              var Error     : byte);

{---------------------------------------------------------------------------}
{-                                                                         -}
{-    Input: NumPoints, Inverse, XReal, XImag,                             -}
{-    Output: XReal, XImag, Error                                          -}
{-                                                                         -}
{-    Purpose:  This procedure uses the complex Fourier transform          -}
{-              routine (FFT) to transform real data.  The real data       -}
{-              is in the vector XReal.  Appropriate shuffling of indices  -}
{-              changes the real vector into two vectors (representing     -}
{-              complex data) which are only half the size of the original -}
{-              vector.  Appropriate unshuffling at the end produces the   -}
{-              transform of the real data.                                -}
{-                                                                         -}
{-  User Defined Types:                                                    -}
{-         TNvector = array[0..TNArraySize] of real                        -}
{-      TNvectorPtr = ^TNvector                                            -}
{-                                                                         -}
{- Global Variables:  NumPoints   : integer     Number of data             -}
{-                                              points in X                -}
{-                    Inverse     : boolean     False => forward transform -}
{-                                              True ==> inverse transform -}
{-                    XReal,XImag : TNvectorPtr Data points                -}
{-                    Error       : byte        Indicates an error         -}
{-                                                                         -}
{-             Errors:  0: No Errors                                       -}
{-                      1: NumPoints < 2                                   -}
{-                      2: NumPoints not a power of two                    -}
{-                         (or 4 for radix-4 transforms)                   -}
{-                                                                         -}
{---------------------------------------------------------------------------}

procedure Start;

implementation

procedure TestInput{(NumPoints    : integer;
                 var NumberOfBits : byte;
                 var Error        : byte)};

type
  ShortArray = array[1..13] of integer;

var
  Term : integer;

const
  PowersOfTwo : ShortArray = (2, 4, 8, 16, 32, 64, 128, 256,
                              512, 1024, 2048, 4096, 8192);

begin
  Error := 2; { Assume NumPoints not a power of two  }
  if NumPoints < 2 then
    Error := 1;     { NumPoints < 2  }
  Term := 1;
  while (Term <= 13) and (Error = 2) do
  begin
    if NumPoints = PowersOfTwo[Term] then
    begin
      NumberOfBits := Term;
      Error := 0;  { NumPoints is a power of two }
    end;
    Term := Succ(Term);
  end;
end; { procedure TestInput }

procedure MakeSinCosTable{(NumPoints : integer;
                       var SinTable  : TNvectorPtr;
                       var CosTable  : TNvectorPtr)};

var
  RealFactor, ImagFactor : Float;
  Term : integer;
  TermMinus1 : integer;
  UpperLimit : integer;

begin
  RealFactor :=  Cos(2 * Pi / NumPoints);
  ImagFactor := -Sqrt(1 - Sqr(RealFactor));
  CosTable^[0] := 1;
  SinTable^[0] := 0;
  CosTable^[1] := RealFactor;
  SinTable^[1] := ImagFactor;
  UpperLimit := NumPoints shr 1 - 1;
  for Term := 2 to UpperLimit do
  begin
    TermMinus1 := Term - 1;
    CosTable^[Term] := CosTable^[TermMinus1] * RealFactor -
                       SinTable^[TermMinus1] * ImagFactor;
    SinTable^[Term] := CosTable^[TermMinus1] * ImagFactor +
                       SinTable^[TermMinus1] * RealFactor;
  end;
end; { procedure MakeSinCosTable }

procedure FFT{(NumberOfBits : byte;
               NumPoints    : integer;
               Inverse      : boolean;
           var XReal        : TNvectorPtr;
           var XImag        : TNvectorPtr;
           var SinTable     : TNvectorPtr;
           var CosTable     : TNvectorPtr)};

const
  RootTwoOverTwo = 0.707106781186548;

var
  Term : byte;
  CellSeparation : integer;
  NumberOfCells : integer;
  NumElementsInCell : integer;
  NumElInCellLess1 : integer;
  NumElInCellSHR1 : integer;
  NumElInCellSHR2 : integer;
  CosTerm, SumTerm, DifTerm : Float;
  Element : integer;
  CellElements : integer;
  ElementInNextCell : integer;
  Index : integer;
  RealDummy, ImagDummy : Float;
  Dummy1, Dummy2 : Float;

procedure BitInvert(NumberOfBits : byte;
                    NumPoints    : integer;
                var XReal        : TNvectorPtr;
                var XImag        : TNvectorPtr);

{-----------------------------------------------------------}
{- Input: NumberOfBits, NumPoints                          -}
{- Output: XReal, XImag                                    -}
{-                                                         -}
{- This procedure bit inverts the order of data in the     -}
{- vector X.  Bit inversion reverses the order of the      -}
{- binary representation of the indices; thus 2 indices    -}
{- will be switched.  For example, if there are 16 points, -}
{- Index 7 (binary 0111) would be switched with Index 14   -}
{- (binary 1110).  It is necessary to bit invert the order -}
{- of the data so that the transformation comes out in the -}
{- correct order.                                          -}
{-----------------------------------------------------------}

var
  Term : integer;
  Invert : integer;
  Hold : Float;
  NumPointsDiv2, K : integer;

begin
  NumPointsDiv2 := NumPoints shr 1;
  Invert := 0;
  for Term := 0 to NumPoints - 2 do
  begin
    if Term < Invert then   { Switch these two indices  }
    begin
      Hold := XReal^[Invert];
      XReal^[Invert] := XReal^[Term];
      XReal^[Term] := Hold;
      Hold := XImag^[Invert];
      XImag^[Invert] := XImag^[Term];
      XImag^[Term] := Hold;
    end;
    K := NumPointsDiv2;
    while K <= Invert do
    begin
      Invert := Invert - K;
      K := K shr 1;
    end;
    Invert := Invert + K;
  end;
end; { procedure BitInvert }

begin { procedure FFT }
  { The data must be entered in bit inverted order }
  { for the transform to come out in proper order  }
  BitInvert(NumberOfBits, NumPoints, XReal, XImag);

  if Inverse then
    { Conjugate the input  }
    for Element := 0 to NumPoints - 1 do
      XImag^[Element] := -XImag^[Element];

  NumberOfCells := NumPoints;
  CellSeparation := 1;
  for Term := 1 to NumberOfBits do
  begin
    { NumberOfCells halves; equals 2^(NumberOfBits - Term)  }
    NumberOfCells := NumberOfCells shr 1;
    { NumElementsInCell doubles; equals 2^(Term-1)  }
    NumElementsInCell := CellSeparation;
    { CellSeparation doubles; equals 2^Term  }
    CellSeparation := CellSeparation SHL 1;
    NumElInCellLess1 := NumElementsInCell - 1;
    NumElInCellSHR1 := NumElementsInCell shr 1;
    NumElInCellSHR2 := NumElInCellSHR1 shr 1;

    { Special case: RootOfUnity = EXP(-i 0)  }
    Element := 0;
    while Element < NumPoints do
    begin
      { Combine the X[Element] with the element in  }
      { the identical location in the next cell     }
      ElementInNextCell := Element + NumElementsInCell;
      RealDummy := XReal^[ElementInNextCell];
      ImagDummy := XImag^[ElementInNextCell];
      XReal^[ElementInNextCell] := XReal^[Element] - RealDummy;
      XImag^[ElementInNextCell] := XImag^[Element] - ImagDummy;
      XReal^[Element] := XReal^[Element] + RealDummy;
      XImag^[Element] := XImag^[Element] + ImagDummy;
      Element := Element + CellSeparation;
    end;

    for CellElements := 1 to NumElInCellSHR2 - 1 do
    begin
      Index := CellElements * NumberOfCells;
      CosTerm := CosTable^[Index];
      SumTerm := CosTerm + SinTable^[Index];
      DifTerm := SinTable^[Index] - CosTerm;
      Element := CellElements;

      while Element < NumPoints do
      begin
        { Combine the X[Element] with the element in  }
        { the identical location in the next cell     }
        ElementInNextCell := Element + NumElementsInCell;
        Dummy1 := (XReal^[ElementInNextCell] +
        XImag^[ElementInNextCell]) * CosTerm;
        Dummy2 := XReal^[ElementInNextCell] * DifTerm;
        RealDummy := Dummy1 - XImag^[ElementInNextCell] * SumTerm;
        ImagDummy := Dummy1 + Dummy2;;
        XReal^[ElementInNextCell] := XReal^[Element] - RealDummy;
        XImag^[ElementInNextCell] := XImag^[Element] - ImagDummy;
        XReal^[Element] := XReal^[Element] + RealDummy;
        XImag^[Element] := XImag^[Element] + ImagDummy;
        Element := Element + CellSeparation;
      end;
    end; { for }

    { Special case: RootOfUnity = EXP(-i PI/4)  }
    if Term > 2 then
    begin
      Element := NumElInCellSHR2;
      while Element < NumPoints do
      begin
        { Combine the X[Element] with the element in  }
        { the identical location in the next cell     }
        ElementInNextCell := Element + NumElementsInCell;
        RealDummy := RootTwoOverTwo * (XReal^[ElementInNextCell] +
                                       XImag^[ElementInNextCell]);
        ImagDummy := RootTwoOverTwo * (XImag^[ElementInNextCell] -
                                       XReal^[ElementInNextCell]);
        XReal^[ElementInNextCell] := XReal^[Element] - RealDummy;
        XImag^[ElementInNextCell] := XImag^[Element] - ImagDummy;
        XReal^[Element] := XReal^[Element] + RealDummy;
        XImag^[Element] := XImag^[Element] + ImagDummy;
        Element := Element + CellSeparation;
      end;
    end;

    for CellElements := NumElInCellSHR2 + 1 to NumElInCellSHR1 - 1 do
    begin
      Index := CellElements * NumberOfCells;
      CosTerm := CosTable^[Index];
      SumTerm := SinTable^[Index] + CosTerm;
      DifTerm := SinTable^[Index] - CosTerm;
      Element := CellElements;
      while Element < NumPoints do
      begin
        { Combine the X[Element] with the element in  }
        { the identical location in the next cell     }
        ElementInNextCell := Element + NumElementsInCell;
        Dummy1 := (XReal^[ElementInNextCell] +
                   XImag^[ElementInNextCell]) * CosTerm;
        Dummy2 := XReal^[ElementInNextCell] * DifTerm;
        RealDummy := Dummy1 - XImag^[ElementInNextCell] * SumTerm;
        ImagDummy := Dummy1 + Dummy2;;
        XReal^[ElementInNextCell] := XReal^[Element] - RealDummy;
        XImag^[ElementInNextCell] := XImag^[Element] - ImagDummy;
        XReal^[Element] := XReal^[Element] + RealDummy;
        XImag^[Element] := XImag^[Element] + ImagDummy;
        Element := Element + CellSeparation;
      end;
    end; { for }

    { Special case: RootOfUnity = EXP(-i PI/2)  }
    if Term > 1 then
    begin
      Element := NumElInCellSHR1;
      while Element < NumPoints do
      begin
        { Combine the X[Element] with the element in  }
        { the identical location in the next cell     }
        ElementInNextCell := Element + NumElementsInCell;
        RealDummy :=  XImag^[ElementInNextCell];
        ImagDummy := -XReal^[ElementInNextCell];
        XReal^[ElementInNextCell] := XReal^[Element] - RealDummy;
        XImag^[ElementInNextCell] := XImag^[Element] - ImagDummy;
        XReal^[Element] := XReal^[Element] + RealDummy;
        XImag^[Element] := XImag^[Element] + ImagDummy;
        Element := Element + CellSeparation;
      end;
    end;

    for CellElements := NumElInCellSHR1 + 1 to
                        NumElementsInCell - NumElInCellSHR2 - 1 do
    begin
      Index := CellElements * NumberOfCells;
      CosTerm := CosTable^[Index];
      SumTerm := SinTable^[Index] + CosTerm;
      DifTerm := SinTable^[Index] - CosTerm;
      Element := CellElements;
      while Element < NumPoints do
      begin
        { Combine the X[Element] with the element in  }
        { the identical location in the next cell     }
        ElementInNextCell := Element + NumElementsInCell;
        Dummy1 := (XReal^[ElementInNextCell] +
                   XImag^[ElementInNextCell]) * CosTerm;
        Dummy2 := XReal^[ElementInNextCell] * DifTerm;
        RealDummy := Dummy1 - XImag^[ElementInNextCell] * SumTerm;
        ImagDummy := Dummy1 + Dummy2;;
        XReal^[ElementInNextCell] := XReal^[Element] - RealDummy;
        XImag^[ElementInNextCell] := XImag^[Element] - ImagDummy;
        XReal^[Element] := XReal^[Element] + RealDummy;
        XImag^[Element] := XImag^[Element] + ImagDummy;
        Element := Element + CellSeparation;
      end;
    end; { for }

    { Special case: RootOfUnity = EXP(-i 3PI/4)  }
    if Term > 2 then
    begin
      Element := NumElementsInCell - NumElInCellSHR2;
      while Element < NumPoints do
      begin
        { Combine the X[Element] with the element in  }
        { the identical location in the next cell     }
        ElementInNextCell := Element + NumElementsInCell;
        RealDummy := -RootTwoOverTwo * (XReal^[ElementInNextCell] -
                                        XImag^[ElementInNextCell]);
        ImagDummy := -RootTwoOverTwo * (XReal^[ElementInNextCell] +
                                        XImag^[ElementInNextCell]);
        XReal^[ElementInNextCell] := XReal^[Element] - RealDummy;
        XImag^[ElementInNextCell] := XImag^[Element] - ImagDummy;
        XReal^[Element] := XReal^[Element] + RealDummy;
        XImag^[Element] := XImag^[Element] + ImagDummy;
        Element := Element + CellSeparation;
      end;
    end;

    for CellElements := NumElementsInCell - NumElInCellSHR2 + 1 to
                                             NumElInCellLess1 do
    begin
      Index := CellElements * NumberOfCells;
      CosTerm := CosTable^[Index];
      SumTerm := SinTable^[Index] + CosTerm;
      DifTerm := SinTable^[Index] - CosTerm;
      Element := CellElements;
      while Element < NumPoints do
      begin
        { Combine the X[Element] with the element in  }
        { the identical location in the next cell     }
        ElementInNextCell := Element + NumElementsInCell;
        Dummy1 := (XReal^[ElementInNextCell] +
                   XImag^[ElementInNextCell]) * CosTerm;
        Dummy2 := XReal^[ElementInNextCell] * DifTerm;
        RealDummy := Dummy1 - XImag^[ElementInNextCell] * SumTerm;
        ImagDummy := Dummy1 + Dummy2;;
        XReal^[ElementInNextCell] := XReal^[Element] - RealDummy;
        XImag^[ElementInNextCell] := XImag^[Element] - ImagDummy;
        XReal^[Element] := XReal^[Element] + RealDummy;
        XImag^[Element] := XImag^[Element] + ImagDummy;
        Element := Element + CellSeparation;
      end;
    end; { for }
  end;

  {----------------------------------------------------}
  {-  Divide all the values of the transformation     -}
  {-  by the square root of NumPoints. If taking the  -}
  {-  inverse, conjugate the output.                  -}
  {----------------------------------------------------}

  if Inverse then
    ImagDummy := -1 / Sqrt(NumPoints)
  else
    ImagDummy :=  1 / Sqrt(NumPoints);
  RealDummy := ABS(ImagDummy);
  for Element := 0 to NumPoints - 1 do
  begin
    XReal^[Element] := XReal^[Element] * RealDummy;
    XImag^[Element] := XImag^[Element] * ImagDummy;
  end;
end; { procedure FFT }

procedure RealFFT{(NumPoints : integer;
                   Inverse   : boolean;
               var XReal     : TNvectorPtr;
               var XImag     : TNvectorPtr;
               var Error     : byte)};
var
  SinTable, CosTable : TNvectorPtr;   { Tables of sine and cosine values  }
  NumberOfBits : byte;                { Number of bits necessary to     }
                                      { represent the number of points  }

procedure MakeRealDataComplex(NumPoints : integer;
                          var XReal     : TNvectorPtr;
                          var XImag     : TNvectorPtr);

{-----------------------------------------------------------}
{- Input: NumPoints, XReal                                 -}
{- Output: XReal, XImag                                    -}
{-                                                         -}
{- This procedure shuffles the real data.  There are       -}
{- 2*NumPoints real data points in the vector XReal.  The  -}
{- data is shuffled so that there are NumPoints complex    -}
{- data points.  The real part of the complex data is      -}
{- made up of those points whose original array Index was  -}
{- even; the imaginary part of the complex data is made    -}
{- up of those points whose original array Index was odd.  -}
{-----------------------------------------------------------}

var
  Index, NewIndex : integer;
  DummyReal, DummyImag : TNvectorPtr;

begin
  New(DummyReal);
  New(DummyImag);
  for Index := 0 to NumPoints - 1 do
  begin
    NewIndex := Index shl 1;
    DummyReal^[Index] := XReal^[NewIndex];
    DummyImag^[Index] := XReal^[NewIndex + 1];
  end;
  XReal^ := DummyReal^;
  XImag^ := DummyImag^;
  Dispose(DummyReal);
  Dispose(DummyImag);
end; { procedure MakeRealDataComplex }

procedure UnscrambleComplexOutput(NumPoints : integer;
                              var SinTable  : TNvectorPtr;
                              var CosTable  : TNvectorPtr;
                              var XReal     : TNvectorPtr;
                              var XImag     : TNvectorPtr);

{----------------------------------------------------------}
{- Input: NumPoints, SinTable, CosTable, XReal, XImag     -}
{- Output: XReal, XImag                                   -}
{-                                                        -}
{- This procedure unshuffles the complex transform.       -}
{- The transform has NumPoints elements.  This procedure  -}
{- unshuffles the transform so that it is 2*NumPoints     -}
{- elements long.  The resulting vector is symmetric      -}
{- about the element NumPoints.                           -}
{- Both the forward and inverse transforms are defined    -}
{- with a 1/Sqrt(NumPoints) factor.  Since the real FFT   -}
{- algorithm operates on vectors of length NumPoints/2,   -}
{- the unscrambled vectors must be divided by Sqrt(2).    -}
{----------------------------------------------------------}

var
  PiOverNumPoints : Float;
  Index : integer;
  indexSHR1 : integer;
  NumPointsMinusIndex : integer;
  SymmetricIndex : integer;
  Multiplier : Float;
  Factor : Float;
  CosFactor, SinFactor : Float;
  RealSum, ImagSum, RealDif, ImagDif : Float;
  RealDummy, ImagDummy : TNvectorPtr;
  NumPointsSHL1 : integer;

begin
  New(RealDummy);
  New(ImagDummy);
  RealDummy^ := XReal^;
  ImagDummy^ := XImag^;
  PiOverNumPoints := Pi / NumPoints;
  NumPointsSHL1 := NumPoints shl 1;
  RealDummy^[0] := (XReal^[0] + XImag^[0]) / Sqrt(2);
  ImagDummy^[0] := 0;
  RealDummy^[NumPoints] := (XReal^[0] - XImag^[0]) / Sqrt(2);
  ImagDummy^[NumPoints] := 0;
  for Index := 1 to NumPoints - 1 do
  begin
    Multiplier := 0.5 / Sqrt(2);
    Factor := PiOverNumPoints * Index;
    NumPointsMinusIndex := NumPoints - Index;
    SymmetricIndex := NumPointsSHL1 - Index;
    if Odd(Index) then
      begin
        CosFactor :=  Cos(Factor);
        SinFactor := -Sin(Factor);
      end
    else
      begin
        indexSHR1 := Index shr 1;
        CosFactor := CosTable^[indexSHR1];
        SinFactor := SinTable^[indexSHR1];
      end;

    RealSum := XReal^[Index] + XReal^[NumPointsMinusIndex];
    ImagSum := XImag^[Index] + XImag^[NumPointsMinusIndex];
    RealDif := XReal^[Index] - XReal^[NumPointsMinusIndex];
    ImagDif := XImag^[Index] - XImag^[NumPointsMinusIndex];

    RealDummy^[Index] := Multiplier * (RealSum + CosFactor * ImagSum
                         + SinFactor * RealDif);
    ImagDummy^[Index] := Multiplier * (ImagDif + SinFactor * ImagSum
                         - CosFactor * RealDif);
    RealDummy^[SymmetricIndex] :=  RealDummy^[Index];
    ImagDummy^[SymmetricIndex] := -ImagDummy^[Index];
  end;  { for }

  XReal^ := RealDummy^;
  XImag^ := ImagDummy^;
  Dispose(RealDummy);
  Dispose(ImagDummy);
end; { procedure UnscrambleComplexOutput }

begin { procedure RealFFT }

  { The number of complex data points will  }
  { be half the number of real data points  }
  NumPoints := NumPoints shr 1;

  TestInput(NumPoints, NumberOfBits, Error);

  if Error = 0 then
  begin
    New(SinTable);
    New(CosTable);
    MakeRealDataComplex(NumPoints, XReal, XImag);
    MakeSinCosTable(NumPoints, SinTable, CosTable);
    FFT(NumberOfBits, NumPoints, Inverse, XReal, XImag, SinTable, CosTable);
    UnscrambleComplexOutput(NumPoints, SinTable, CosTable, XReal, XImag);
    NumPoints := NumPoints shl 1;    { The number of complex points  }
                                     { in the transform will be the  }
                                     { same as the number of real    }
                                     { points in input data.         }
    Dispose(SinTable);
    Dispose(CosTable);
  end;
end; { procedure RealFFT }

procedure Initialize(var NumPoints : integer;
                     var XReal     : TNvectorPtr;
                     var XImag     : TNvectorPtr;
                     var HReal     : TNvectorPtr;
                     var HImag     : TNvectorPtr;
                     var Error     : byte);

{----------------------------------------------------------}
{- Output: NumPoints, XReal, XImag, HReal, HImag, Error   -}
{-                                                        -}
{- This procedure initializes the above variables to zero -}
{----------------------------------------------------------}

begin
  NumPoints := 0;
  New(XReal);
  New(XImag);
  New(HReal);
  New(HImag);
  FillChar(XReal^, SizeOf(XReal^), 0);
  FillChar(XImag^, SizeOf(XImag^), 0);
  FillChar(HReal^, SizeOf(HReal^), 0);
  FillChar(HImag^, SizeOf(HImag^), 0);
  Error := 0;
end; { procedure Initialize }

procedure GetData(var NumPoints     : integer;
                  var WhichAnalysis : Analyses;
                  var Auto          : boolean;
                  var XReal         : TNvectorPtr;
                  var XImag         : TNvectorPtr;
                  var HReal         : TNvectorPtr;
                  var HImag         : TNvectorPtr);

{--------------------------------------------------------------}
{- Output: NumPoints, WhichAnalysis, Auto, XReal, XImag,      -}
{-         HReal, HImag                                       -}
{-                                                            -}
{- This procedure reads in data from either the keyboard      -}
{- or a data file.                                            -}
{--------------------------------------------------------------}

var
  Ch : char;
  NumPoints1 : integer;
  NumPoints2 : integer;

function TestForPowersOfTwo(NumPoints : integer) : boolean;

{---------------------------------------------------------}
{- Input: NumPoints                                      -}
{- Output: TestForPowersOfTwo                            -}
{-                                                       -}
{- This procedure checks if NumPoints is a power of two. -}
{- It returns True if it is, False if it isn't.          -}
{---------------------------------------------------------}

type
  ShortArray = array[1..13] of integer;

const
  PowersOfTwo : ShortArray = (2, 4, 8, 16, 32, 64, 128, 256,
                              512, 1024, 2048, 4096, 8192);

var
  Term : integer;
  Test : boolean;

begin
  Test := false; { Assume NumPoints not a power of two  }
  Term := 1;
  while (Term <= 13) and (not Test) do
  begin
    if NumPoints = PowersOfTwo[Term] then
      Test := true; { NumPoints is a power of two  }
    Term := Succ(Term);
  end;
  TestForPowersOfTwo := Test;
end; { function TestForPowersOfTwo }

procedure GetRealVectorFromFile(var NumPoints : integer;
                                var XReal     : TNvectorPtr);

{-----------------------------------------------}
{- Output: NumPoints, X                        -}
{-                                             -}
{- This procedure reads in a real vector of    -}
{- data points from a data file.               -}
{-----------------------------------------------}

var
  FileName : string[255];
  InFile : text;

begin
 FileName:='fft.dat';
 AssignFile(InFile, FileName);
 Reset(InFile);
 NumPoints := 0;
 while not EOF(InFile)  do
  begin
   Readln(InFile, XReal^[NumPoints]);
   NumPoints := Succ(NumPoints);
  end;
 Close(InFile);
end; { procedure GetRealVectorFromFile }

procedure GetComplexVectorFromFile(var NumPoints : integer;
                                   var XReal     : TNvectorPtr;
                                   var XImag     : TNvectorPtr);

{-----------------------------------------------}
{- Output: NumPoints, XReal, XImag             -}
{-                                             -}
{- This procedure reads in a complex vector    -}
{- of data points from a data file.            -}
{-----------------------------------------------}

var
  FileName : string[255];
  InFile : text;

begin
 FileName:='fftcmplx.dat';
 Assign(InFile, FileName);
 Reset(InFile);
 NumPoints := 0;
 while not EOF(InFile) do
  begin
   Readln(InFile, XReal^[NumPoints], XImag^[NumPoints]);
   NumPoints := Succ(NumPoints);
  end;
 Close(InFile);
end; { procedure GetComplexVectorFromFile }

procedure GetInverse(var Inverse : boolean);

var
  Ch : char;

begin
//  if Ch = 'F' then
    Inverse := false  { Forward transform  }
//  else
//    Inverse := true;  { Inverse transform  }
end; { procedure GetInverse }

begin { procedure GetData }
{  Writeln('   1) Real Fast Fourier Transform');
  Writeln('   2) Real Convolution');
  Writeln('   3) Real Autocorrelation');
  Writeln('   4) Real Crosscorrelation');
  Writeln('   5) Complex Fast Fourier Transform');
  Writeln('   6) Complex Convolution');
  Writeln('   7) Complex Autocorrelation');
  Writeln('   8) Complex Crosscorrelation');
  Writeln;
  Write('   Select a number (1-8): ');
  repeat
    Ch := ReadKey;
  until Ch in ['1'..'8'];

  Writeln;
  case Ch of
    '1' : begin
            Writeln('********* Real fast Fourier transform *********');}
            GetInverse(Inverse);
            WhichAnalysis := RF;
{          end;

    '2' : begin
            Writeln('********* Real Convolution *********');
            WhichAnalysis := RN;
          end;

    '3' : begin
            Writeln('********* Real Autocorrelation *********');
            Auto := true;            { Autocorrelation  }
 {           WhichAnalysis := RA;
          end;

    '4' : begin
            Writeln('********* Real Crosscorrelation *********');
            Auto := false;           { Crosscorrelation  }
  {          WhichAnalysis := RC;
          end;

    '5' : begin
            Writeln('********* Complex fast Fourier transform *********');
            GetInverse(Inverse);
            WhichAnalysis := CF;
          end;

    '6' : begin
            Writeln('********* Complex Convolution *********');
            WhichAnalysis := CN;
          end;

    '7' : begin
            Writeln('********* Complex Autocorrelation *********');
            Auto := true;            { Autocorrelation  }
   {         WhichAnalysis := CA;
          end;

    '8' : begin
            Writeln('********* Complex Crosscorrelation *********');
            Auto := false;           { Crosscorrelation  }
    {        WhichAnalysis := CC;
          end;
  end; { case }

{  Writeln;
  Ch := InputChannel('Input Data From');
  if Ch = 'K' then
    repeat  { Until NumPoints1 := NumPoints2  }
 {     NumPoints1 := NumPoints2;
      case WhichAnalysis of
        RF, RA : GetRealVectorFromKeyboard(NumPoints, XReal);

        RN, RC : begin
                   Writeln;
                   Writeln('The first function:');
                   GetRealVectorFromKeyboard(NumPoints1, XReal);
                   Writeln;
                   Writeln('The 2nd function:');
                   GetRealVectorFromKeyboard(NumPoints2, HReal);
                   if NumPoints1 <> NumPoints2 then
                   begin
                     Writeln;
                     Write('The number of points in these two files');
                     Writeln(' are different.');
                   end;
                   NumPoints := NumPoints1;
                 end;

        CF, CA : GetComplexVectorFromKeyboard(NumPoints, XReal, XImag);

        CN, CC : begin
                   Writeln;
                   Writeln('The first function:');
                   GetComplexVectorFromKeyboard(NumPoints1, XReal, XImag);
                   Writeln;
                   Writeln('The 2nd function:');
                   GetComplexVectorFromKeyboard(NumPoints2, HReal, HImag);
                   if NumPoints1 <> NumPoints2 then
                   begin
                     Writeln;
                     Write('The number of points in these two files');
                     Writeln(' are different.');
                   end;
                   NumPoints := NumPoints1;
                 end;
      end { case }
 {   until NumPoints1 = NumPoints2
  else { Ch = 'F' }
    repeat
      NumPoints2 := NumPoints1;
      case WhichAnalysis of
        RF, RA : GetRealVectorFromFile(NumPoints, XReal);

       { RN, RC : begin
                   Writeln;
                   Writeln('The first function:');
                   GetRealVectorFromFile(NumPoints1, XReal);
                   Writeln;
                   Writeln('The 2nd function:');
                   GetRealVectorFromFile(NumPoints2, HReal);
                   if NumPoints1 <> NumPoints2 then
                   begin
                     Writeln;
                     Write('The number of points in these two files');
                     Writeln(' are different.');
                   end;
                   NumPoints := NumPoints1;
                 end;

        CF, CA : GetComplexVectorFromFile(NumPoints, XReal, XImag);

        CN, CC : begin
                   Writeln;
                   Writeln('The first function:');
                   GetComplexVectorFromFile(NumPoints1, XReal, XImag);
                   Writeln;
                   Writeln('The 2nd function:');
                   GetComplexVectorFromFile(NumPoints2, HReal, HImag);
                   if NumPoints1 <> NumPoints2 then
                   begin
                     Writeln;
                     Write('The number of points in these two files');
                     Writeln(' are different.');
                   end;
                   NumPoints := NumPoints1;
                 end; }
      end; { case }
    until NumPoints1 = NumPoints2;
//  GetOutputFile(OutFile);
end; { procedure GetData }

procedure Results(NumPoints     : integer;
                  WhichAnalysis : Analyses;
              var XReal         : TNvectorPtr;
              var XImag         : TNvectorPtr;
                  Error         : byte);

{------------------------------------------------------------}
{- This procedure outputs the results to the device OutFile -}
{------------------------------------------------------------}

var
  Index : integer;

begin
 AssignFile(OutFile, 'fftout.dat');
 Rewrite(OutFile);
 Writeln(OutFile);
 Writeln(OutFile);
// if Error >= 1 then
//    DisplayError;

  case Error of
    0 : begin
          case WhichAnalysis of
            RF : Writeln(OutFile, 'Results of real Fourier transform:');
            RN : Writeln(OutFile, 'Results of real convolution:');
            RC : Writeln(OutFile, 'Results of real crosscorrelation:');
            RA : Writeln(OutFile, 'Results of real autocorrelation:');
            CF : Writeln(OutFile, 'Results of complex Fourier transform:');
            CN : Writeln(OutFile, 'Results of complex convolution:');
            CC : Writeln(OutFile, 'Results of complex crosscorrelation:');
            CA : Writeln(OutFile, 'Results of complex autocorrelation:');
          end; { case }

          for Index := 0 to NumPoints - 1 do
            Writeln(OutFile, XReal^[Index]:0:5, '    ',  XImag^[Index]:0:5);
          Writeln(OutFile);
        end;

{    1 : begin
          Writeln(OutFile, 'The number of data points: ', NumPoints);
          Writeln(OutFile, 'There are too few data points.');
        end;

    2 : begin
          Writeln(OutFile, 'The number of data points: ', NumPoints);
          Writeln(OutFile,
                  'The number of data points must be a power of two for a');
          Writeln(OutFile,
             'radix-2 transform or a power of four for a radix-4 transform.');
        end;}
  end; { case }
end; { procedure Results }

procedure Start;
 begin
  Initialize(NumPoints, XReal, XImag, HReal, HImag, Error);
  GetData(NumPoints, WhichAnalysis, Auto, XReal, XImag, HReal, HImag);

  {case WhichAnalysis of
    RF :} RealFFT(NumPoints, Inverse, XReal, XImag, Error);
 {   RN : RealConvolution(NumPoints, XReal, XImag, HReal, Error);
    RC : RealCorrelation(NumPoints, Auto, XReal, XImag, HReal, Error);
    RA : RealCorrelation(NumPoints, Auto, XReal, XImag, XReal, Error);
    CF : ComplexFFT(NumPoints, Inverse, XReal, XImag, Error);
    CN : ComplexConvolution(NumPoints, XReal, XImag, HReal, HImag, Error);
    CC : ComplexCorrelation(NumPoints, Auto, XReal, XImag,
                            HReal, HImag, Error);
    CA : ComplexCorrelation(NumPoints, Auto, XReal, XImag,
                            XReal, XImag, Error);
  end; { case }

  Results(NumPoints, WhichAnalysis, XReal, XImag, Error);
  Close(OutFile);
  Dispose(XReal);
  Dispose(XImag);
  Dispose(HReal);
  Dispose(HImag);
end;

end.
