home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 9 / 09.iso / l / l042 / 2.ddi / CHAP10.ARC / FFT.INC next >
Encoding:
Text File  |  1987-12-30  |  23.4 KB  |  662 lines
  1. {----------------------------------------------------------------------------}
  2. {-                                                                          -}
  3. {-     Turbo Pascal Numerical Methods Toolbox                               -}
  4. {-     Copyright (c) 1986, 87 by Borland International, Inc.                -}
  5. {-                                                                          -}
  6. {----------------------------------------------------------------------------}
  7.  
  8. procedure RealFFT{(NumPoints : integer;
  9.                    Inverse   : boolean;
  10.                var XReal     : TNvectorPtr;
  11.                var XImag     : TNvectorPtr;
  12.                var Error     : byte)};
  13. var
  14.   SinTable, CosTable : TNvectorPtr;   { Tables of sine and cosine values  }
  15.   NumberOfBits : byte;                { Number of bits necessary to     }
  16.                                       { represent the number of points  }
  17.  
  18. procedure MakeRealDataComplex(NumPoints : integer;
  19.                           var XReal     : TNvectorPtr;
  20.                           var XImag     : TNvectorPtr);
  21.  
  22. {-----------------------------------------------------------}
  23. {- Input: NumPoints, XReal                                 -}
  24. {- Output: XReal, XImag                                    -}
  25. {-                                                         -}
  26. {- This procedure shuffles the real data.  There are       -}
  27. {- 2*NumPoints real data points in the vector XReal.  The  -}
  28. {- data is shuffled so that there are NumPoints complex    -}
  29. {- data points.  The real part of the complex data is      -}
  30. {- made up of those points whose original array Index was  -}
  31. {- even; the imaginary part of the complex data is made    -}
  32. {- up of those points whose original array Index was odd.  -}
  33. {-----------------------------------------------------------}
  34.  
  35. var
  36.   Index, NewIndex : integer;
  37.   DummyReal, DummyImag : TNvectorPtr;
  38.  
  39. begin
  40.   New(DummyReal);
  41.   New(DummyImag);
  42.   for Index := 0 to NumPoints - 1 do
  43.   begin
  44.     NewIndex := Index shl 1;
  45.     DummyReal^[Index] := XReal^[NewIndex];
  46.     DummyImag^[Index] := XReal^[NewIndex + 1];
  47.   end;
  48.   XReal^ := DummyReal^;
  49.   XImag^ := DummyImag^;
  50.   Dispose(DummyReal);
  51.   Dispose(DummyImag);
  52. end; { procedure MakeRealDataComplex }
  53.  
  54. procedure UnscrambleComplexOutput(NumPoints : integer;
  55.                               var SinTable  : TNvectorPtr;
  56.                               var CosTable  : TNvectorPtr;
  57.                               var XReal     : TNvectorPtr;
  58.                               var XImag     : TNvectorPtr);
  59.  
  60. {----------------------------------------------------------}
  61. {- Input: NumPoints, SinTable, CosTable, XReal, XImag     -}
  62. {- Output: XReal, XImag                                   -}
  63. {-                                                        -}
  64. {- This procedure unshuffles the complex transform.       -}
  65. {- The transform has NumPoints elements.  This procedure  -}
  66. {- unshuffles the transform so that it is 2*NumPoints     -}
  67. {- elements long.  The resulting vector is symmetric      -}
  68. {- about the element NumPoints.                           -}
  69. {- Both the forward and inverse transforms are defined    -}
  70. {- with a 1/Sqrt(NumPoints) factor.  Since the real FFT   -}
  71. {- algorithm operates on vectors of length NumPoints/2,   -}
  72. {- the unscrambled vectors must be divided by Sqrt(2).    -}
  73. {----------------------------------------------------------}
  74.  
  75. var
  76.   PiOverNumPoints : Float;
  77.   Index : integer;
  78.   indexSHR1 : integer;
  79.   NumPointsMinusIndex : integer;
  80.   SymmetricIndex : integer;
  81.   Multiplier : Float;
  82.   Factor : Float;
  83.   CosFactor, SinFactor : Float;
  84.   RealSum, ImagSum, RealDif, ImagDif : Float;
  85.   RealDummy, ImagDummy : TNvectorPtr;
  86.   NumPointsSHL1 : integer;
  87.  
  88. begin
  89.   New(RealDummy);
  90.   New(ImagDummy);
  91.   RealDummy^ := XReal^;
  92.   ImagDummy^ := XImag^;
  93.   PiOverNumPoints := Pi / NumPoints;
  94.   NumPointsSHL1 := NumPoints shl 1;
  95.   RealDummy^[0] := (XReal^[0] + XImag^[0]) / Sqrt(2);
  96.   ImagDummy^[0] := 0;
  97.   RealDummy^[NumPoints] := (XReal^[0] - XImag^[0]) / Sqrt(2);
  98.   ImagDummy^[NumPoints] := 0;
  99.   for Index := 1 to NumPoints - 1 do
  100.   begin
  101.     Multiplier := 0.5 / Sqrt(2);
  102.     Factor := PiOverNumPoints * Index;
  103.     NumPointsMinusIndex := NumPoints - Index;
  104.     SymmetricIndex := NumPointsSHL1 - Index;
  105.     if Odd(Index) then
  106.       begin
  107.         CosFactor :=  Cos(Factor);
  108.         SinFactor := -Sin(Factor);
  109.       end
  110.     else
  111.       begin
  112.         indexSHR1 := Index shr 1;
  113.         CosFactor := CosTable^[indexSHR1];
  114.         SinFactor := SinTable^[indexSHR1];
  115.       end;
  116.  
  117.     RealSum := XReal^[Index] + XReal^[NumPointsMinusIndex];
  118.     ImagSum := XImag^[Index] + XImag^[NumPointsMinusIndex];
  119.     RealDif := XReal^[Index] - XReal^[NumPointsMinusIndex];
  120.     ImagDif := XImag^[Index] - XImag^[NumPointsMinusIndex];
  121.  
  122.     RealDummy^[Index] := Multiplier * (RealSum + CosFactor * ImagSum
  123.                          + SinFactor * RealDif);
  124.     ImagDummy^[Index] := Multiplier * (ImagDif + SinFactor * ImagSum
  125.                          - CosFactor * RealDif);
  126.     RealDummy^[SymmetricIndex] :=  RealDummy^[Index];
  127.     ImagDummy^[SymmetricIndex] := -ImagDummy^[Index];
  128.   end;  { for }
  129.  
  130.   XReal^ := RealDummy^;
  131.   XImag^ := ImagDummy^;
  132.   Dispose(RealDummy);
  133.   Dispose(ImagDummy);
  134. end; { procedure UnscrambleComplexOutput }
  135.  
  136. begin { procedure RealFFT }
  137.  
  138.   { The number of complex data points will  }
  139.   { be half the number of real data points  }
  140.   NumPoints := NumPoints shr 1;
  141.  
  142.   TestInput(NumPoints, NumberOfBits, Error);
  143.  
  144.   if Error = 0 then
  145.   begin
  146.     New(SinTable);
  147.     New(CosTable);
  148.     MakeRealDataComplex(NumPoints, XReal, XImag);
  149.     MakeSinCosTable(NumPoints, SinTable, CosTable);
  150.     FFT(NumberOfBits, NumPoints, Inverse, XReal, XImag, SinTable, CosTable);
  151.     UnscrambleComplexOutput(NumPoints, SinTable, CosTable, XReal, XImag);
  152.     NumPoints := NumPoints shl 1;    { The number of complex points  }
  153.                                      { in the transform will be the  }
  154.                                      { same as the number of real    }
  155.                                      { points in input data.         }
  156.     Dispose(SinTable);
  157.     Dispose(CosTable);
  158.   end;
  159. end; { procedure RealFFT }
  160.  
  161. procedure RealConvolution{(NumPoints : integer;
  162.                        var XReal     : TNvectorPtr;
  163.                        var XImag     : TNvectorPtr;
  164.                        var HReal     : TNvectorPtr;
  165.                        var Error     : byte)};
  166.  
  167. var
  168.   SinTable, CosTable : TNvectorPtr;    { Tables of sine and cos values  }
  169.   Inverse : boolean;                   { Indicates inverse transform  }
  170.   NumberOfBits : byte;                 { Number of bits required to  }
  171.                                        { represent NumPoints    }
  172.   HImag : TNvectorPtr;                 { Imaginary part of the   }
  173.                                        { FFT transform of HReal  }
  174.  
  175. procedure Multiply(NumPoints : integer;
  176.                var XReal     : TNvectorPtr;
  177.                var XImag     : TNvectorPtr;
  178.                var HReal     : TNvectorPtr;
  179.                var HImag     : TNvectorPtr);
  180.  
  181. {----------------------------------------------}
  182. {- Input: NumPoints, XReal, XImag, HReal,     -}
  183. {-        HImag                               -}
  184. {- Output: XReal, XImag                       -}
  185. {-                                            -}
  186. {- This procedure multiplies each element in  -}
  187. {- X by the corresponding element in H.  The  -}
  188. {- product is returned in X.                  -}
  189. {----------------------------------------------}
  190.  
  191. var
  192.   Term : integer;
  193.   Dummy : Float;
  194.  
  195. begin
  196.   for Term := 0 to NumPoints - 1 do
  197.   begin
  198.     Dummy := XReal^[Term] * HReal^[Term] - XImag^[Term] * HImag^[Term];
  199.     XImag^[Term] := XReal^[Term] * HImag^[Term] +
  200.                     XImag^[Term] * HReal^[Term];
  201.     XReal^[Term] := Dummy;
  202.   end;
  203. end; { procedure Multiply }
  204.  
  205. procedure SeparateTransforms(NumPoints : integer;
  206.                          var XReal     : TNvectorPtr;
  207.                          var XImag     : TNvectorPtr;
  208.                          var HReal     : TNvectorPtr;
  209.                          var HImag     : TNvectorPtr);
  210.  
  211. {-------------------------------------------------------}
  212. {- Input: NumPoints, XReal, XImag                      -}
  213. {- Output: HReal, HImag                                -}
  214. {-                                                     -}
  215. {- The transforms of the real data XReal and HReal are -}
  216. {- stored in the vector XReal, XImag.  By using the    -}
  217. {- symmetries of the transforms of real data, this     -}
  218. {- procedure extracts the complex transforms.  The     -}
  219. {- transform of XReal is returned in XReal, XImag. The -}
  220. {- transform of HReal is returned in HReal, HImag.     -}
  221. {-------------------------------------------------------}
  222.  
  223. var
  224.   Term : integer;
  225.   EndTerm : integer;
  226.   DummyReal : Float;
  227.   DummyImag : Float;
  228.   NumPointsSHR1 : integer;
  229.   NumPointsMinusTerm : integer;
  230.  
  231. begin
  232.   HReal^[0] := XImag^[0];
  233.   HImag^[0] := 0;
  234.   XImag^[0] := 0;
  235.   NumPointsSHR1 := NumPoints SHR 1;
  236.   HReal^[NumPointsSHR1] := XImag^[NumPointsSHR1];
  237.   HImag^[NumPointsSHR1] := 0;
  238.   XImag^[NumPointsSHR1] := 0;
  239.   EndTerm := NumPointsSHR1 - 1;
  240.   for Term := 1 to EndTerm do
  241.   begin
  242.     NumPointsMinusTerm := NumPoints - Term;
  243.     HReal^[Term] := 0.5 * (XImag^[Term] + XImag^[NumPointsMinusTerm]);
  244.     HImag^[Term] := 0.5 * (XReal^[NumPointsMinusTerm] - XReal^[Term]);
  245.     DummyReal := 0.5 * (XReal^[Term] + XReal^[NumPointsMinusTerm]);
  246.     DummyImag := 0.5 * (XImag^[Term] - XImag^[NumPointsMinusTerm]);
  247.     XReal^[Term] := DummyReal;
  248.     XImag^[Term] := DummyImag;
  249.   end;
  250.  
  251.   for Term := 1 to EndTerm do
  252.   begin  { Make use of symmetries to calculate  }
  253.          { the rest of the transform            }
  254.     NumPointsMinusTerm := NumPoints - Term;
  255.     XReal^[NumPointsMinusTerm] :=  XReal^[Term];
  256.     XImag^[NumPointsMinusTerm] := -XImag^[Term];
  257.     HReal^[NumPointsMinusTerm] :=  HReal^[Term];
  258.     HImag^[NumPointsMinusTerm] := -HImag^[Term];
  259.   end;
  260. end; { procedure SeparateTransforms }
  261.  
  262. begin { procedure RealConvolution }
  263.   TestInput(NumPoints, NumberOfBits, Error);
  264.   if Error = 0 then
  265.   begin
  266.     New(HImag);
  267.     New(SinTable);
  268.     New(CosTable);
  269.  
  270.     { Combine the two real transforms  }
  271.     { into one complex transform       }
  272.     XImag^ := HReal^;
  273.  
  274.     MakeSinCosTable(NumPoints, SinTable, CosTable);
  275.  
  276.     { Take the transform of XReal, XImag  }
  277.     Inverse := false;
  278.  
  279.     FFT(NumberOfBits, NumPoints, Inverse, XReal, XImag, SinTable, CosTable);
  280.  
  281.     { Separate out the X transform and H transform  }
  282.     SeparateTransforms(NumPoints, XReal, XImag, HReal, HImag);
  283.  
  284.     { Multiply the two transforms  }
  285.     Multiply(NumPoints, XReal, XImag, HReal, HImag);
  286.  
  287.     { Take the inverse transform of the product  }
  288.     Inverse := true;
  289.     FFT(NumberOfBits, NumPoints, Inverse, XReal, XImag, SinTable, CosTable);
  290.     Dispose(SinTable);
  291.     Dispose(CosTable);
  292.     Dispose(HImag);
  293.   end;
  294. end; { procedure RealConvolution }
  295.  
  296. procedure RealCorrelation{(NumPoints : integer;
  297.                       var Auto       : boolean;
  298.                       var XReal      : TNvectorPtr;
  299.                       var XImag      : TNvectorPtr;
  300.                       var HReal      : TNvectorPtr;
  301.                       var Error      : byte)};
  302.  
  303. var
  304.   SinTable, CosTable : TNvectorPtr;  { Tables of sine and cosine values  }
  305.   Inverse : boolean;                 { Indicates inverse transform  }
  306.   NumberOfBits : byte;               { Number of bits necessary to     }
  307.                                      { represent the number of points  }
  308.   HImag : TNvectorPtr;               { Imaginary part of the   }
  309.                                      { FFT transform of HReal  }
  310.  
  311. procedure SeparateTransforms(NumPoints : integer;
  312.                          var XReal     : TNvectorPtr;
  313.                          var XImag     : TNvectorPtr;
  314.                          var HReal     : TNvectorPtr;
  315.                          var HImag     : TNvectorPtr);
  316.  
  317. {-------------------------------------------------------}
  318. {- Input: NumPoints, XReal, XImag                      -}
  319. {- Output: HReal, HImag                                -}
  320. {-                                                     -}
  321. {- The transforms of the real data XReal and HReal are -}
  322. {- stored in the vector XReal, XImag.  By using the    -}
  323. {- symmetries of the transforms of real data, this     -}
  324. {- procedure extracts the complex transforms.  The     -}
  325. {- transform of XReal is returned in XReal, XImag. The -}
  326. {- transform of HReal is returned in HReal, HImag.     -}
  327. {-------------------------------------------------------}
  328.  
  329. var
  330.   Term : integer;
  331.   EndTerm : integer;
  332.   DummyReal : Float;
  333.   DummyImag : Float;
  334.   NumPointsSHR1 : integer;
  335.   NumPointsMinusTerm : integer;
  336.  
  337. begin
  338.   HReal^[0] := XImag^[0];
  339.   HImag^[0] := 0;
  340.   XImag^[0] := 0;
  341.   NumPointsSHR1 := NumPoints SHR 1;
  342.   HReal^[NumPointsSHR1] := XImag^[NumPointsSHR1];
  343.   HImag^[NumPointsSHR1] := 0;
  344.   XImag^[NumPointsSHR1] := 0;
  345.   EndTerm := NumPointsSHR1 - 1;
  346.   for Term := 1 to EndTerm do
  347.   begin
  348.     NumPointsMinusTerm := NumPoints - Term;
  349.     HReal^[Term] := 0.5 * (XImag^[Term] + XImag^[NumPointsMinusTerm]);
  350.     HImag^[Term] := 0.5 * (XReal^[NumPointsMinusTerm] - XReal^[Term]);
  351.     DummyReal := 0.5 * (XReal^[Term] + XReal^[NumPointsMinusTerm]);
  352.     DummyImag := 0.5 * (XImag^[Term] - XImag^[NumPointsMinusTerm]);
  353.     XReal^[Term] := DummyReal;
  354.     XImag^[Term] := DummyImag;
  355.   end;
  356.  
  357.   for Term := 1 to EndTerm do
  358.   begin  { Make use of symmetries to calculate  }
  359.          { the rest of the transform            }
  360.     NumPointsMinusTerm := NumPoints - Term;
  361.     XReal^[NumPointsMinusTerm] :=  XReal^[Term];
  362.     XImag^[NumPointsMinusTerm] := -XImag^[Term];
  363.     HReal^[NumPointsMinusTerm] :=  HReal^[Term];
  364.     HImag^[NumPointsMinusTerm] := -HImag^[Term];
  365.   end;
  366. end; { procedure SeparateTransforms }
  367.  
  368. procedure Multiply(NumPoints : integer;
  369.                var XReal     : TNvectorPtr;
  370.                var XImag     : TNvectorPtr;
  371.                var HReal     : TNvectorPtr;
  372.                var HImag     : TNvectorPtr);
  373.  
  374. {----------------------------------------------------}
  375. {- Input: NumPoints, XReal, XImag, HReal, HImag     -}
  376. {- Output: XReal, XImag                             -}
  377. {-                                                  -}
  378. {- This procedure performs the following operation: -}
  379. {-                                                  -}
  380. {-     Dummy[f] := X[f] * H[-f]                     -}
  381. {-                                                  -}
  382. {- where -f is represented by NumPoints - f         -}
  383. {- (circular correlation).  Because x and h were    -}
  384. {- real functions, this operation is identical      -}
  385. {- to:                                              -}
  386. {-                           *                      -}
  387. {-    Dummy[f] := X[f] - H[f]                       -}
  388. {-                                                  -}
  389. {- The product is returned in X.                    -}
  390. {----------------------------------------------------}
  391.  
  392. var
  393.   Term : integer;
  394.   NumPointsMinusTerm : integer;
  395.   DummyReal : TNvectorPtr;
  396.   DummyImag : TNvectorPtr;
  397.  
  398. begin
  399.   New(DummyReal);
  400.   New(DummyImag);
  401.   DummyReal^[0] := XReal^[0] * HReal^[0] - XImag^[0] * HImag^[0];
  402.   DummyImag^[0] := XReal^[0] * HImag^[0] + XImag^[0] * HReal^[0];
  403.   for Term := 1 to NumPoints - 1 do
  404.   begin
  405.     NumPointsMinusTerm := NumPoints - Term;
  406.     DummyReal^[Term] := XReal^[Term] * HReal^[NumPointsMinusTerm] -
  407.                         XImag^[Term] * HImag^[NumPointsMinusTerm];
  408.     DummyImag^[Term] := XImag^[Term] * HReal^[NumPointsMinusTerm] +
  409.                         XReal^[Term] * HImag^[NumPointsMinusTerm];
  410.   end;
  411.   XReal^ := DummyReal^;
  412.   XImag^ := DummyImag^;
  413.   Dispose(DummyReal);
  414.   Dispose(DummyImag);
  415. end; { procedure Multiply }
  416.  
  417. begin { procedure RealCorrelation }
  418.   TestInput(NumPoints, NumberOfBits, Error);
  419.   if Error = 0 then
  420.   begin
  421.     if Auto then
  422.       FillChar(XImag^, SizeOf(XImag^), 0)
  423.     else    { Combine the two real transforms  }
  424.             { into one complex transform       }
  425.       XImag^ := HReal^;
  426.     New(HImag);
  427.     New(SinTable);
  428.     New(CosTable);
  429.  
  430.     MakeSinCosTable(NumPoints, SinTable, CosTable);
  431.  
  432.     { Take the transform of XReal, XImag  }
  433.     Inverse := false;
  434.  
  435.     FFT(NumberOfBits, NumPoints, Inverse, XReal, XImag, SinTable, CosTable);
  436.  
  437.     if Auto then
  438.       begin
  439.         HReal^ := XReal^;
  440.         HImag^ := XImag^;
  441.       end
  442.     else
  443.       SeparateTransforms(NumPoints, XReal, XImag, HReal, HImag);
  444.  
  445.     { Multiply the two transforms together  }
  446.     Multiply(NumPoints, XReal, XImag, HReal, HImag);
  447.  
  448.     { Take the inverse transform of the product  }
  449.     Inverse := true;
  450.  
  451.     FFT(NumberOfBits, NumPoints, Inverse, XReal, XImag, SinTable, CosTable);
  452.  
  453.     Dispose(SinTable);
  454.     Dispose(CosTable);
  455.     Dispose(HImag);
  456.   end;
  457. end; { procedure RealCorrelation }
  458.  
  459. procedure ComplexFFT{(NumPoints : integer;
  460.                       Inverse   : boolean;
  461.                   var XReal     : TNvectorPtr;
  462.                   var XImag     : TNvectorPtr;
  463.                   var Error     : byte)};
  464.  
  465. var
  466.   SinTable, CosTable : TNvectorPtr;      { Tables of sin and cosine values  }
  467.   NumberOfBits : byte;                   { Number of bits to represent the  }
  468.                                          { number of data points.           }
  469.  
  470. begin { procedure ComplexFFT }
  471.  
  472.   TestInput(NumPoints, NumberOfBits, Error);
  473.  
  474.   if Error = 0 then
  475.   begin
  476.     New(SinTable);
  477.     New(CosTable);
  478.  
  479.     MakeSinCosTable(NumPoints, SinTable, CosTable);
  480.     FFT(NumberOfBits, NumPoints, Inverse, XReal, XImag, SinTable, CosTable);
  481.  
  482.     Dispose(SinTable);
  483.     Dispose(CosTable);
  484.   end;
  485. end; { procedure ComplexFFT }
  486.  
  487. procedure ComplexConvolution{(NumPoints : integer;
  488.                           var XReal     : TNvectorPtr;
  489.                           var XImag     : TNvectorPtr;
  490.                           var HReal     : TNvectorPtr;
  491.                           var HImag     : TNvectorPtr;
  492.                           var Error     : byte)};
  493.  
  494. var
  495.   SinTable, CosTable : TNvectorPtr;    { Tables of sine and cosine values  }
  496.   Inverse : boolean;                   { Indicates inverse transform  }
  497.   NumberOfBits : byte;                 { # of bits required to  }
  498.                                        { represent NumPoints    }
  499.  
  500. procedure Multiply(NumPoints : integer;
  501.                var XReal     : TNvectorPtr;
  502.                var XImag     : TNvectorPtr;
  503.                var HReal     : TNvectorPtr;
  504.                var HImag     : TNvectorPtr);
  505.  
  506. {----------------------------------------------}
  507. {- Input: NumPoints, XReal, XImag, HReal,     -}
  508. {-        HImag                               -}
  509. {- Output: XReal, XImag                       -}
  510. {-                                            -}
  511. {- This procedure multiplies each element in  -}
  512. {- X by the corresponding element in H.  The  -}
  513. {- product is returned in X.                  -}
  514. {----------------------------------------------}
  515.  
  516. var
  517.   Term : integer;
  518.   Dummy : Float;
  519.  
  520. begin
  521.   for Term := 0 to NumPoints - 1 do
  522.   begin
  523.     Dummy := XReal^[Term] * HReal^[Term] - XImag^[Term] * HImag^[Term];
  524.     XImag^[Term] := XReal^[Term] * HImag^[Term] +
  525.                     XImag^[Term] * HReal^[Term];
  526.     XReal^[Term] := Dummy;
  527.   end;
  528. end; { procedure Multiply }
  529.  
  530. begin { procedure ComplexConvolution }
  531.  
  532.   TestInput(NumPoints, NumberOfBits, Error);
  533.  
  534.   if Error = 0 then
  535.   begin
  536.     New(SinTable);
  537.     New(CosTable);
  538.  
  539.     MakeSinCosTable(NumPoints, SinTable, CosTable);
  540.  
  541.     { Take the transform of XReal, XImag  }
  542.     Inverse := false;
  543.  
  544.     FFT(NumberOfBits, NumPoints, Inverse, XReal, XImag, SinTable, CosTable);
  545.  
  546.     { Take the transform of HReal, HImag  }
  547.  
  548.     FFT(NumberOfBits, NumPoints, Inverse, HReal, HImag, SinTable, CosTable);
  549.  
  550.     { Multiply the two transforms  }
  551.     Multiply(NumPoints, XReal, XImag, HReal, HImag);
  552.  
  553.     { Take the inverse transform of the product  }
  554.     Inverse := true;
  555.  
  556.     FFT(NumberOfBits, NumPoints, Inverse, XReal, XImag, SinTable, CosTable);
  557.  
  558.     Dispose(SinTable);
  559.     Dispose(CosTable);
  560.   end;
  561. end; { procedure ComplexConvolution }
  562.  
  563. procedure ComplexCorrelation{(NumPoints : integer;
  564.                           var Auto      : boolean;
  565.                           var XReal     : TNvectorPtr;
  566.                           var XImag     : TNvectorPtr;
  567.                           var HReal     : TNvectorPtr;
  568.                           var HImag     : TNvectorPtr;
  569.                           var Error     : byte)};
  570.  
  571. var
  572.   SinTable, CosTable : TNvectorPtr;  { Tables of sine and cosine values  }
  573.   Inverse : boolean;                 { Indicates inverse transform  }
  574.   NumberOfBits : byte;               { Number of bits necessary to     }
  575.                                      { represent the number of points  }
  576.  
  577. procedure Multiply(NumPoints : integer;
  578.                var XReal     : TNvectorPtr;
  579.                var XImag     : TNvectorPtr;
  580.                var HReal     : TNvectorPtr;
  581.                var HImag     : TNvectorPtr);
  582.  
  583. {----------------------------------------------------}
  584. {- Input: NumPoints, XReal, XImag, HReal, HImag     -}
  585. {- Output: XReal, XImag                             -}
  586. {-                                                  -}
  587. {- This procedure performs the following operation: -}
  588. {-                                                  -}
  589. {-     Dummy[f] := X[f] * H[-f]                     -}
  590. {-                                                  -}
  591. {- where -f is represented by NumPoints - f         -}
  592. {- (circular correlation).                          -}
  593. {-                                                  -}
  594. {- The product is returned in X.                    -}
  595. {----------------------------------------------------}
  596.  
  597. var
  598.   Term : integer;
  599.   NumPointsMinusTerm : integer;
  600.   DummyReal : TNvectorPtr;
  601.   DummyImag : TNvectorPtr;
  602.  
  603. begin
  604.   New(DummyReal);
  605.   New(DummyImag);
  606.   DummyReal^[0] := XReal^[0] * HReal^[0] - XImag^[0] * HImag^[0];
  607.   DummyImag^[0] := XReal^[0] * HImag^[0] + XImag^[0] * HReal^[0];
  608.   for Term := 1 to NumPoints - 1 do
  609.   begin
  610.     NumPointsMinusTerm := NumPoints - Term;
  611.     DummyReal^[Term] := XReal^[Term] * HReal^[NumPointsMinusTerm] -
  612.                         XImag^[Term] * HImag^[NumPointsMinusTerm];
  613.     DummyImag^[Term] := XImag^[Term] * HReal^[NumPointsMinusTerm] +
  614.                         XReal^[Term] * HImag^[NumPointsMinusTerm];
  615.   end;
  616.   XReal^ := DummyReal^;
  617.   XImag^ := DummyImag^;
  618.   Dispose(DummyReal);
  619.   Dispose(DummyImag);
  620. end; { procedure Multiply }
  621.  
  622. begin { procedure ComplexCorrelation }
  623.  
  624.   TestInput(NumPoints, NumberOfBits, Error);
  625.  
  626.   if Error = 0 then
  627.   begin
  628.     New(SinTable);
  629.     New(CosTable);
  630.  
  631.     MakeSinCosTable(NumPoints, SinTable, CosTable);
  632.  
  633.     { Take the transform of XReal, XImag  }
  634.     Inverse := false;
  635.  
  636.     FFT(NumberOfBits, NumPoints, Inverse, XReal, XImag, SinTable, CosTable);
  637.  
  638.     if not Auto then
  639.       { Take the transform of HReal, HImag  }
  640.       begin
  641.         FFT(NumberOfBits, NumPoints, Inverse, HReal, HImag, SinTable, CosTable);
  642.       end
  643.     else
  644.       { Autocorrelation; set H equal to X  }
  645.       begin
  646.         HReal := XReal;
  647.         HImag := XImag;
  648.       end;
  649.  
  650.     { Multiply the two transforms together  }
  651.     Multiply(NumPoints, XReal, XImag, HReal, HImag);
  652.  
  653.     { Take the inverse transform of the product  }
  654.     Inverse := true;
  655.  
  656.     FFT(NumberOfBits, NumPoints, Inverse, XReal, XImag, SinTable, CosTable);
  657.  
  658.     Dispose(SinTable);
  659.     Dispose(CosTable);
  660.   end;
  661. end; { procedure ComplexCorrelation }
  662.