home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 9 / 09.iso / l / l042 / 2.ddi / CHAP10.ARC / FFT87B2.PAS < prev    next >
Encoding:
Pascal/Delphi Source File  |  1987-12-30  |  29.4 KB  |  636 lines
  1. unit FFT87B2;
  2.  
  3. {----------------------------------------------------------------------------}
  4. {-                                                                          -}
  5. {-     Turbo Pascal Numerical Methods Toolbox                               -}
  6. {-     Copyright (c) 1986, 87 by Borland International, Inc.                -}
  7. {-                                                                          -}
  8. {-  This unit provides procedures for performing real and complex fast      -}
  9. {-  fourier transforms. Radix-2 8087 version                                -}
  10. {-                                                                          -}
  11. {----------------------------------------------------------------------------}
  12.  
  13. {$N+} { Requires the 8087 math co-processor }
  14.  
  15. interface
  16.  
  17. const
  18.   TNArraySize = 1023;
  19.  
  20. type
  21.   Float       = double;
  22.   TNvector    = array[0..TNArraySize] of Float;
  23.   TNvectorPtr = ^TNvector;
  24.  
  25. procedure TestInput(NumPoints    : integer;
  26.                 var NumberOfBits : byte;
  27.                 var Error        : byte);
  28.  
  29. {-------------------------------------------------------------}
  30. {- Input: NumPoints                                          -}
  31. {- Output: NumberOfBits, Error                               -}
  32. {-                                                           -}
  33. {- This procedure checks the input.  If the number of points -}
  34. {- (NumPoints) is less than two or is not a multiple of two  -}
  35. {- then an error is returned.  NumberOfBits is the number of -}
  36. {- bits necessary to represent NumPoints in binary (e.g. if  -}
  37. {- NumPoints = 16, NumberOfBits = 4).                        -}
  38. {-------------------------------------------------------------}
  39.  
  40. procedure MakeSinCosTable(NumPoints : integer;
  41.                       var SinTable  : TNvectorPtr;
  42.                       var CosTable  : TNvectorPtr);
  43.  
  44. {--------------------------------------------------------}
  45. {- Input: NumPoints                                     -}
  46. {- Output: SinTable, CosTable                           -}
  47. {-                                                      -}
  48. {- This procedure fills in a table with sin and cosine  -}
  49. {- values.  It is faster to pull data out of this       -}
  50. {- table than it is to calculate the sines and cosines. -}
  51. {--------------------------------------------------------}
  52.  
  53. procedure FFT(NumberOfBits : byte;
  54.               NumPoints    : integer;
  55.               Inverse      : boolean;
  56.           var XReal        : TNvectorPtr;
  57.           var XImag        : TNvectorPtr;
  58.           var SinTable     : TNvectorPtr;
  59.           var CosTable     : TNvectorPtr);
  60.  
  61. {-----------------------------------------------------}
  62. {- Input: NumberOfBits, NumPoints, Inverse, XReal,   -}
  63. {-        XImag, SinTable, CosTable                  -}
  64. {- Output: XReal, XImag                              -}
  65. {-                                                   -}
  66. {- This procedure implements the actual fast Fourier -}
  67. {- transform routine.  The vector X, which must be   -}
  68. {- entered in bit-inverted order, is transformed in  -}
  69. {- place.  The transformation uses the Cooley-Tukey  -}
  70. {- algorithm.                                        -}
  71. {-----------------------------------------------------}
  72.  
  73. procedure RealFFT(NumPoints : integer;
  74.                   Inverse   : boolean;
  75.               var XReal     : TNvectorPtr;
  76.               var XImag     : TNvectorPtr;
  77.               var Error     : byte);
  78.  
  79.  
  80. {---------------------------------------------------------------------------}
  81. {-                                                                         -}
  82. {-    Input: NumPoints, Inverse, XReal, XImag,                             -}
  83. {-    Output: XReal, XImag, Error                                          -}
  84. {-                                                                         -}
  85. {-    Purpose:  This procedure uses the complex Fourier transform          -}
  86. {-              routine (FFT) to transform real data.  The real data       -}
  87. {-              is in the vector XReal.  Appropriate shuffling of indices  -}
  88. {-              changes the real vector into two vectors (representing     -}
  89. {-              complex data) which are only half the size of the original -}
  90. {-              vector.  Appropriate unshuffling at the end produces the   -}
  91. {-              transform of the real data.                                -}
  92. {-                                                                         -}
  93. {-  User Defined Types:                                                    -}
  94. {-         TNvector = array[0..TNArraySize] of real                        -}
  95. {-      TNvectorPtr = ^TNvector                                            -}
  96. {-                                                                         -}
  97. {- Global Variables:  NumPoints   : integer     Number of data             -}
  98. {-                                              points in X                -}
  99. {-                    Inverse     : boolean     False => forward transform -}
  100. {-                                              True ==> inverse transform -}
  101. {-                    XReal,XImag : TNvectorPtr Data points                -}
  102. {-                    Error       : byte        Indicates an error         -}
  103. {-                                                                         -}
  104. {-             Errors:  0: No Errors                                       -}
  105. {-                      1: NumPoints < 2                                   -}
  106. {-                      2: NumPoints not a power of two                    -}
  107. {-                         (or 4 for radix-4 transforms)                   -}
  108. {-                                                                         -}
  109. {---------------------------------------------------------------------------}
  110.  
  111. procedure RealConvolution(NumPoints : integer;
  112.                       var XReal     : TNvectorPtr;
  113.                       var XImag     : TNvectorPtr;
  114.                       var HReal     : TNvectorPtr;
  115.                       var Error     : byte);
  116.  
  117. {-------------------------------------------------------------------}
  118. {-                                                                 -}
  119. {-   Input: NumPoints, XReal, XImag, HReal                         -}
  120. {-   Output: XReal, XImag, Error                                   -}
  121. {-                                                                 -}
  122. {-   Purpose: This procedure performs a convolution of the         -}
  123. {-            real data XReal and HReal.  The result is returned   -}
  124. {-            in the complex vector XReal, XImag.                  -}
  125. {-                                                                 -}
  126. {- User Defined Types:                                             -}
  127. {-         TNvector = array[0..TNArraySize] of real                -}
  128. {-      TNvectorPtr = ^TNvector                                    -}
  129. {-                                                                 -}
  130. {- Global Variables:  NumPoints : integer     Number of data       -}
  131. {-                                            points in X          -}
  132. {-                    XReal     : TNvectorPtr Data points          -}
  133. {-                    HReal     : TNvectorPtr Data points          -}
  134. {-                    Error     : byte        Indicates an error   -}
  135. {-                                                                 -}
  136. {-             Errors:  0: No Errors                               -}
  137. {-                      1: NumPoints < 2                           -}
  138. {-                      2: NumPoints not a power of two            -}
  139. {-                                                                 -}
  140. {-------------------------------------------------------------------}
  141.  
  142. procedure RealCorrelation(NumPoints : integer;
  143.                       var Auto      : boolean;
  144.                       var XReal     : TNvectorPtr;
  145.                       var XImag     : TNvectorPtr;
  146.                       var HReal     : TNvectorPtr;
  147.                       var Error     : byte);
  148.  
  149. {-------------------------------------------------------------------}
  150. {-                                                                 -}
  151. {-   Input: NumPoints, Auto, XReal, XImag, HReal                   -}
  152. {-   Output: XReal, XImag, Error                                   -}
  153. {-                                                                 -}
  154. {-   Purpose: This procedure performs a correlation (auto or       -}
  155. {-            cross) of the real data XReal and HReal. The         -}
  156. {-            correlation is returned in the complex vector        -}
  157. {-            XReal, XImag.                                        -}
  158. {-                                                                 -}
  159. {- User Defined Types:                                             -}
  160. {-         TNvector = array[0..TNArraySize] OF real                -}
  161. {-      TNvectorPtr = ^TNvector                                    -}
  162. {-                                                                 -}
  163. {- Global Variables:  NumPoints : integer     Number of data       -}
  164. {-                                            points in X          -}
  165. {-                    Auto      : boolean     True => auto-        -}
  166. {-                                                    correlation  -}
  167. {-                                            False=> cross-       -}
  168. {-                                                    correlation  -}
  169. {-                    XReal     : TNvectorPtr First sample         -}
  170. {-                    HReal     : TNvectorPtr Second sample        -}
  171. {-                    Error     : byte        Indicates an error   -}
  172. {-                                                                 -}
  173. {-             Errors:  0: No Errors                               -}
  174. {-                      1: NumPoints < 2                           -}
  175. {-                      2: NumPoints not a power of two            -}
  176. {-                                                                 -}
  177. {-------------------------------------------------------------------}
  178.  
  179. procedure ComplexFFT(NumPoints : integer;
  180.                      Inverse   : boolean;
  181.                  var XReal     : TNvectorPtr;
  182.                  var XImag     : TNvectorPtr;
  183.                  var Error     : byte);
  184.  
  185. {-------------------------------------------------------------------}
  186. {-                                                                 -}
  187. {-   Input: NumPoints, Inverse, XReal, XImag                       -}
  188. {-   Output: XReal, XImag, Error                                   -}
  189. {-                                                                 -}
  190. {-   Purpose: This procedure performs a fast Fourier transform     -}
  191. {-            of the complex data XReal, XImag. The vectors XReal  -}
  192. {-            and XImag are transformed in place.                  -}
  193. {-                                                                 -}
  194. {- User Defined Types:                                             -}
  195. {-         TNvector = array[0..TNArraySize] of real                -}
  196. {-         TNvectorPtr = ^TNvector                                 -}
  197. {-                                                                 -}
  198. {- Global Variables:  NumPoints : integer      Number of data      -}
  199. {-                                             points in X         -}
  200. {-                    Inverse   : BOOLEAN      FALSE => Forward    -}
  201. {-                                                      Transform  -}
  202. {-                                             TRUE => Inverse     -}
  203. {-                                                     Transform   -}
  204. {-                    XReal,                                       -}
  205. {-                    XImag     : TNvectorPtr  Data points         -}
  206. {-                    Error     : byte         Indicates an error  -}
  207. {-                                                                 -}
  208. {-             Errors:  0: No Errors                               -}
  209. {-                      1: NumPoints < 2                           -}
  210. {-                      2: NumPoints not a power of two            -}
  211. {-                                                                 -}
  212. {-------------------------------------------------------------------}
  213.  
  214. procedure ComplexConvolution(NumPoints : integer;
  215.                          var XReal     : TNvectorPtr;
  216.                          var XImag     : TNvectorPtr;
  217.                          var HReal     : TNvectorPtr;
  218.                          var HImag     : TNvectorPtr;
  219.                          var Error     : byte);
  220.  
  221. {-------------------------------------------------------------------}
  222. {-                                                                 -}
  223. {-   Input: NumPoints, XReal, XImag, HReal, HImag                  -}
  224. {-   Output: XReal, XImag, Error                                   -}
  225. {-                                                                 -}
  226. {-   Purpose: This procedure performs a convolution of the         -}
  227. {-            data XReal, XImag and the data HReal and HImag. The  -}
  228. {-            vectors XReal, XImag, HReal and HImag are            -}
  229. {-            transformed in place.                                -}
  230. {-                                                                 -}
  231. {- User Defined Types:                                             -}
  232. {-         TNvector = array[0..TNArraySize] of real                -}
  233. {-      TNvectorPtr = ^TNvector                                    -}
  234. {-                                                                 -}
  235. {- Global Variables:  NumPoints   : integer     Number of data     -}
  236. {-                                              points in X        -}
  237. {-                    XReal,XImag : TNvectorPtr Data points        -}
  238. {-                    HReal,HImag : TNvectorPtr Data points        -}
  239. {-                    Error       : byte        Indicates an error -}
  240. {-                                                                 -}
  241. {-             Errors:  0: No Errors                               -}
  242. {-                      1: NumPoints < 2                           -}
  243. {-                      2: NumPoints not a power of two            -}
  244. {-                                                                 -}
  245. {-------------------------------------------------------------------}
  246.  
  247. procedure ComplexCorrelation(NumPoints : integer;
  248.                          var Auto      : boolean;
  249.                          var XReal     : TNvectorPtr;
  250.                          var XImag     : TNvectorPtr;
  251.                          var HReal     : TNvectorPtr;
  252.                          var HImag     : TNvectorPtr;
  253.                          var Error     : byte);
  254.  
  255. {-------------------------------------------------------------------}
  256. {-                                                                 -}
  257. {-   Input: NumPoints, Auto, XReal, XImag, HReal, HImag            -}
  258. {-   Output: XReal, XImag, Error                                   -}
  259. {-                                                                 -}
  260. {-   Purpose: This procedure performs a correlation (auto or       -}
  261. {-            cross) of the complex data XReal, XImag and the      -}
  262. {-            complex data HReal, HImag. The vectors XReal, XImag, -}
  263. {-            HReal, and HImag are transformed in place.           -}
  264. {-                                                                 -}
  265. {- User Defined Types:                                             -}
  266. {-         TNvector = array[0..TNArraySize] of real                -}
  267. {-      TNvectorPtr = ^TNvector                                    -}
  268. {-                                                                 -}
  269. {- Global Variables:  NumPoints   : integer   Number of data       -}
  270. {-                                            points in X          -}
  271. {-                    Auto        : boolean   True => auto-        -}
  272. {-                                                    correlation  -}
  273. {-                                            False=> cross-       -}
  274. {-                                                    correlation  -}
  275. {-                    XReal,XImag : TNvectorPtr First sample       -}
  276. {-                    HReal,HImag : TNvectorPtr Second sample      -}
  277. {-                    Error       : byte        Indicates an error -}
  278. {-                                                                 -}
  279. {-             Errors:  0: No Errors                               -}
  280. {-                      1: NumPoints < 2                           -}
  281. {-                      2: NumPoints not a power of two            -}
  282. {-                                                                 -}
  283. {-------------------------------------------------------------------}
  284.  
  285. implementation
  286.  
  287. procedure TestInput{(NumPoints    : integer;
  288.                  var NumberOfBits : byte;
  289.                  var Error        : byte)};
  290. type
  291.   ShortArray = array[1..13] of integer;
  292.  
  293. var
  294.   Term : integer;
  295.  
  296. const
  297.   PowersOfTwo : ShortArray = (2, 4, 8, 16, 32, 64, 128, 256,
  298.                               512, 1024, 2048, 4096, 8192);
  299.  
  300. begin
  301.   Error := 2;            { Assume NumPoints not a power of two  }
  302.   if NumPoints < 2 then
  303.     Error := 1;     { NumPoints < 2  }
  304.   Term := 1;
  305.   while (Term <= 13) and (Error = 2) do
  306.   begin
  307.     if NumPoints = PowersOfTwo[Term] then
  308.     begin
  309.       NumberOfBits := Term;
  310.       Error := 0;  { NumPoints is a power of two  }
  311.     end;
  312.     Term := Succ(Term);
  313.   end;
  314. end; { procedure TestInput }
  315.  
  316. procedure MakeSinCosTable{(NumPoints : integer;
  317.                        var SinTable  : TNvectorPtr;
  318.                        var CosTable  : TNvectorPtr)};
  319. var
  320.   RealFactor, ImagFactor : Float;
  321.   Term : integer;
  322.   TermMinus1 : integer;
  323.   UpperLimit : integer;
  324.  
  325. begin
  326.   RealFactor :=  Cos(2 * Pi / NumPoints);
  327.   ImagFactor := -Sqrt(1 - Sqr(RealFactor));
  328.   CosTable^[0] := 1;
  329.   SinTable^[0] := 0;
  330.   CosTable^[1] := RealFactor;
  331.   SinTable^[1] := ImagFactor;
  332.   UpperLimit := NumPoints shr 1 - 1;
  333.   for Term := 2 to UpperLimit do
  334.   begin
  335.     TermMinus1 := Term - 1;
  336.     CosTable^[Term] :=  CosTable^[TermMinus1] * RealFactor -
  337.                         SinTable^[TermMinus1] * ImagFactor;
  338.     SinTable^[Term] :=  CosTable^[TermMinus1] * ImagFactor +
  339.                         SinTable^[TermMinus1] * RealFactor;
  340.   end;
  341. end; { procedure MakeSinCosTable }
  342.  
  343. procedure FFT{(NumberOfBits : byte;
  344.                NumPoints    : integer;
  345.                Inverse      : boolean;
  346.            var XReal        : TNvectorPtr;
  347.            var XImag        : TNvectorPtr;
  348.            var SinTable     : TNvectorPtr;
  349.            var CosTable     : TNvectorPtr)};
  350.  
  351. const
  352.   RootTwoOverTwo = 0.707106781186548;
  353.  
  354. var
  355.   Term : byte;
  356.   CellSeparation : integer;
  357.   NumberOfCells : integer;
  358.   NumElementsInCell : integer;
  359.   NumElInCellLess1 : integer;
  360.   NumElInCellSHR1 : integer;
  361.   NumElInCellSHR2 : integer;
  362.   RealRootOfUnity, ImagRootOfUnity : Float;
  363.   Element : integer;
  364.   CellElements : integer;
  365.   ElementInNextCell : integer;
  366.   Index : integer;
  367.   RealDummy, ImagDummy : Float;
  368.  
  369. procedure BitInvert(NumberOfBits : byte;
  370.                     NumPoints    : integer;
  371.                 var XReal        : TNvectorPtr;
  372.                 var XImag        : TNvectorPtr);
  373.  
  374. {-----------------------------------------------------------}
  375. {- Input: NumberOfBits, NumPoints                          -}
  376. {- Output: XReal, XImag                                    -}
  377. {-                                                         -}
  378. {- This procedure bit inverts the order of data in the     -}
  379. {- vector X.  Bit inversion reverses the order of the      -}
  380. {- binary representation of the indices; thus 2 indices    -}
  381. {- will be switched.  For example, if there are 16 points, -}
  382. {- Index 7 (binary 0111) would be switched with Index 14   -}
  383. {- (binary 1110).  It is necessary to bit invert the order -}
  384. {- of the data so that the transformation comes out in the -}
  385. {- correct order.                                          -}
  386. {-----------------------------------------------------------}
  387.  
  388. var
  389.   Term : integer;
  390.   Invert : integer;
  391.   Hold : Float;
  392.   NumPointsDiv2, K : integer;
  393.  
  394. begin
  395.   NumPointsDiv2 := NumPoints shr 1;
  396.   Invert := 0;
  397.   for Term := 0 to NumPoints - 2 do
  398.   begin
  399.     if Term < Invert then   { Switch these two indices  }
  400.     begin
  401.       Hold := XReal^[Invert];
  402.       XReal^[Invert] := XReal^[Term];
  403.       XReal^[Term] := Hold;
  404.       Hold := XImag^[Invert];
  405.       XImag^[Invert] := XImag^[Term];
  406.       XImag^[Term] := Hold;
  407.     end;
  408.     K := NumPointsDiv2;
  409.     while K <= Invert do
  410.     begin
  411.       Invert := Invert - K;
  412.       K := K shr 1;
  413.     end;
  414.     Invert := Invert + K;
  415.   end;
  416. end; { procedure BitInvert }
  417.  
  418. begin { procedure FFT }
  419.   { The data must be entered in bit inverted order }
  420.   { for the transform to come out in proper order  }
  421.   BitInvert(NumberOfBits, NumPoints, XReal, XImag);
  422.  
  423.   if Inverse then
  424.     { Conjugate the input  }
  425.     for Element := 0 to NumPoints - 1 do
  426.       XImag^[Element] := -XImag^[Element];
  427.  
  428.   NumberOfCells := NumPoints;
  429.   CellSeparation := 1;
  430.   for Term := 1 to NumberOfBits do
  431.   begin
  432.     { NumberOfCells halves; equals 2^(NumberOfBits - Term)  }
  433.     NumberOfCells := NumberOfCells shr 1;
  434.     { NumElementsInCell doubles; equals 2^(Term-1)  }
  435.     NumElementsInCell := CellSeparation;
  436.     { CellSeparation doubles; equals 2^Term  }
  437.     CellSeparation := CellSeparation SHL 1;
  438.     NumElInCellLess1 := NumElementsInCell - 1;
  439.     NumElInCellSHR1 := NumElementsInCell shr 1;
  440.     NumElInCellSHR2 := NumElInCellSHR1 shr 1;
  441.  
  442.     { Special case: RootOfUnity = EXP(-i 0)  }
  443.     Element := 0;
  444.     while Element < NumPoints do
  445.     begin
  446.       { Combine the X[Element] with the element in  }
  447.       { the identical location in the next cell     }
  448.       ElementInNextCell := Element + NumElementsInCell;
  449.       RealDummy := XReal^[ElementInNextCell];
  450.       ImagDummy := XImag^[ElementInNextCell];
  451.       XReal^[ElementInNextCell] := XReal^[Element] - RealDummy;
  452.       XImag^[ElementInNextCell] := XImag^[Element] - ImagDummy;
  453.       XReal^[Element] := XReal^[Element] + RealDummy;
  454.       XImag^[Element] := XImag^[Element] + ImagDummy;
  455.       Element := Element + CellSeparation;
  456.     end;
  457.  
  458.     for CellElements := 1 to NumElInCellSHR2 - 1 do
  459.     begin
  460.       Index := CellElements * NumberOfCells;
  461.       RealRootOfUnity := CosTable^[Index];
  462.       ImagRootOfUnity := SinTable^[Index];
  463.       Element := CellElements;
  464.  
  465.       while Element < NumPoints do
  466.       begin
  467.         { Combine the X[Element] with the element in  }
  468.         { the identical location in the next cell     }
  469.         ElementInNextCell := Element + NumElementsInCell;
  470.         RealDummy := XReal^[ElementInNextCell] * RealRootOfUnity -
  471.                      XImag^[ElementInNextCell] * ImagRootOfUnity;
  472.         ImagDummy := XReal^[ElementInNextCell] * ImagRootOfUnity +
  473.                      XImag^[ElementInNextCell] * RealRootOfUnity;
  474.         XReal^[ElementInNextCell] := XReal^[Element] - RealDummy;
  475.         XImag^[ElementInNextCell] := XImag^[Element] - ImagDummy;
  476.         XReal^[Element] := XReal^[Element] + RealDummy;
  477.         XImag^[Element] := XImag^[Element] + ImagDummy;
  478.         Element := Element + CellSeparation;
  479.       end;
  480.     end;
  481.  
  482.     { Special case: RootOfUnity = EXP(-i PI/4)  }
  483.     if Term > 2 then
  484.     begin
  485.       Element := NumElInCellSHR2;
  486.       while Element < NumPoints do
  487.       begin
  488.         { Combine the X[Element] with the element in  }
  489.         { the identical location in the next cell     }
  490.         ElementInNextCell := Element + NumElementsInCell;
  491.         RealDummy := RootTwoOverTwo * (XReal^[ElementInNextCell] +
  492.                      XImag^[ElementInNextCell]);
  493.         ImagDummy := RootTwoOverTwo * (XImag^[ElementInNextCell] -
  494.                      XReal^[ElementInNextCell]);
  495.         XReal^[ElementInNextCell] := XReal^[Element] - RealDummy;
  496.         XImag^[ElementInNextCell] := XImag^[Element] - ImagDummy;
  497.         XReal^[Element] := XReal^[Element] + RealDummy;
  498.         XImag^[Element] := XImag^[Element] + ImagDummy;
  499.         Element := Element + CellSeparation;
  500.       end;
  501.     end;
  502.  
  503.     for CellElements := NumElInCellSHR2 + 1 to NumElInCellSHR1 - 1 do
  504.     begin
  505.       Index := CellElements * NumberOfCells;
  506.       RealRootOfUnity := CosTable^[Index];
  507.       ImagRootOfUnity := SinTable^[Index];
  508.       Element := CellElements;
  509.       while Element < NumPoints do
  510.       begin
  511.         { Combine the X[Element] with the element in  }
  512.         { the identical location in the next cell     }
  513.         ElementInNextCell := Element + NumElementsInCell;
  514.         RealDummy := XReal^[ElementInNextCell] * RealRootOfUnity -
  515.                      XImag^[ElementInNextCell] * ImagRootOfUnity;
  516.         ImagDummy := XReal^[ElementInNextCell] * ImagRootOfUnity +
  517.                      XImag^[ElementInNextCell] * RealRootOfUnity;
  518.         XReal^[ElementInNextCell] := XReal^[Element] - RealDummy;
  519.         XImag^[ElementInNextCell] := XImag^[Element] - ImagDummy;
  520.         XReal^[Element] := XReal^[Element] + RealDummy;
  521.         XImag^[Element] := XImag^[Element] + ImagDummy;
  522.         Element := Element + CellSeparation;
  523.       end;
  524.     end;
  525.  
  526.     { Special case: RootOfUnity = EXP(-i PI/2)  }
  527.     if Term > 1 then
  528.     begin
  529.       Element := NumElInCellSHR1;
  530.       while Element < NumPoints do
  531.       begin
  532.         { Combine the X[Element] with the element in  }
  533.         { the identical location in the next cell     }
  534.         ElementInNextCell := Element + NumElementsInCell;
  535.         RealDummy :=  XImag^[ElementInNextCell];
  536.         ImagDummy := -XReal^[ElementInNextCell];
  537.         XReal^[ElementInNextCell] := XReal^[Element] - RealDummy;
  538.         XImag^[ElementInNextCell] := XImag^[Element] - ImagDummy;
  539.         XReal^[Element] := XReal^[Element] + RealDummy;
  540.         XImag^[Element] := XImag^[Element] + ImagDummy;
  541.         Element := Element + CellSeparation;
  542.       end;
  543.     end;
  544.  
  545.     for CellElements := NumElInCellSHR1 + 1 to
  546.                         NumElementsInCell - NumElInCellSHR2 - 1 do
  547.     begin
  548.       Index := CellElements * NumberOfCells;
  549.       RealRootOfUnity := CosTable^[Index];
  550.       ImagRootOfUnity := SinTable^[Index];
  551.       Element := CellElements;
  552.       while Element < NumPoints do
  553.       begin
  554.         { Combine the X[Element] with the element in  }
  555.         { the identical location in the next cell     }
  556.         ElementInNextCell := Element + NumElementsInCell;
  557.         RealDummy := XReal^[ElementInNextCell] * RealRootOfUnity -
  558.                      XImag^[ElementInNextCell] * ImagRootOfUnity;
  559.         ImagDummy := XReal^[ElementInNextCell] * ImagRootOfUnity +
  560.                      XImag^[ElementInNextCell] * RealRootOfUnity;
  561.         XReal^[ElementInNextCell] := XReal^[Element] - RealDummy;
  562.         XImag^[ElementInNextCell] := XImag^[Element] - ImagDummy;
  563.         XReal^[Element] := XReal^[Element] + RealDummy;
  564.         XImag^[Element] := XImag^[Element] + ImagDummy;
  565.         Element := Element + CellSeparation;
  566.       end;
  567.     end;
  568.  
  569.     { Special case: RootOfUnity = EXP(-i 3PI/4)  }
  570.     if Term > 2 then
  571.     begin
  572.       Element := NumElementsInCell - NumElInCellSHR2;
  573.       while Element < NumPoints do
  574.       begin
  575.         { Combine the X[Element] with the element in  }
  576.         { the identical location in the next cell     }
  577.         ElementInNextCell := Element + NumElementsInCell;
  578.         RealDummy := -RootTwoOverTwo * (XReal^[ElementInNextCell] -
  579.                                         XImag^[ElementInNextCell]);
  580.         ImagDummy := -RootTwoOverTwo * (XReal^[ElementInNextCell] +
  581.                                         XImag^[ElementInNextCell]);
  582.         XReal^[ElementInNextCell] := XReal^[Element] - RealDummy;
  583.         XImag^[ElementInNextCell] := XImag^[Element] - ImagDummy;
  584.         XReal^[Element] := XReal^[Element] + RealDummy;
  585.         XImag^[Element] := XImag^[Element] + ImagDummy;
  586.         Element := Element + CellSeparation;
  587.       end;
  588.     end;
  589.  
  590.     for CellElements := NumElementsInCell - NumElInCellSHR2 + 1 to
  591.                                             NumElInCellLess1 do
  592.     begin
  593.       Index := CellElements * NumberOfCells;
  594.       RealRootOfUnity := CosTable^[Index];
  595.       ImagRootOfUnity := SinTable^[Index];
  596.       Element := CellElements;
  597.       while Element < NumPoints do
  598.       begin
  599.         { Combine the X[Element] with the element in  }
  600.         { the identical location in the next cell     }
  601.         ElementInNextCell := Element + NumElementsInCell;
  602.         RealDummy := XReal^[ElementInNextCell] * RealRootOfUnity -
  603.                      XImag^[ElementInNextCell] * ImagRootOfUnity;
  604.         ImagDummy := XReal^[ElementInNextCell] * ImagRootOfUnity +
  605.                      XImag^[ElementInNextCell] * RealRootOfUnity;
  606.         XReal^[ElementInNextCell] := XReal^[Element] - RealDummy;
  607.         XImag^[ElementInNextCell] := XImag^[Element] - ImagDummy;
  608.         XReal^[Element] := XReal^[Element] + RealDummy;
  609.         XImag^[Element] := XImag^[Element] + ImagDummy;
  610.         Element := Element + CellSeparation;
  611.       end;
  612.     end;
  613.   end;
  614.  
  615.   {----------------------------------------------------}
  616.   {-  Divide all the values of the transformation     -}
  617.   {-  by the square root of NumPoints. If taking the  -}
  618.   {-  inverse, conjugate the output.                  -}
  619.   {----------------------------------------------------}
  620.  
  621.   if Inverse then
  622.     ImagDummy := -1/Sqrt(NumPoints)
  623.   else
  624.     ImagDummy :=  1/Sqrt(NumPoints);
  625.   RealDummy := ABS(ImagDummy);
  626.   for Element := 0 to NumPoints - 1 do
  627.   begin
  628.     XReal^[Element] := XReal^[Element] * RealDummy;
  629.     XImag^[Element] := XImag^[Element] * ImagDummy;
  630.   end;
  631. end; { procedure FFT }
  632.  
  633. {$I FFT.inc} { Include procedure code }
  634.  
  635. end. { FFT87B2 }
  636.