home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 9 / 09.iso / l / l042 / 2.ddi / CHAP10.ARC / FFT87B4.PAS < prev    next >
Encoding:
Pascal/Delphi Source File  |  1987-12-30  |  39.7 KB  |  883 lines
  1. unit FFT87B4;
  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-4 8087 version                                -}
  10. {-                                                                          -}
  11. {----------------------------------------------------------------------------}
  12.  
  13. {$N+} { Requires the 8087 numeric 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 NumberOfTwoBits : byte;
  27.                 var Error           : byte);
  28.  
  29. {-------------------------------------------------------------}
  30. {- Input: NumPoints                                          -}
  31. {- Output: NumberOfTwoBits, Error                            -}
  32. {-                                                           -}
  33. {- This procedure checks the input.  If the number of points -}
  34. {- (NumPoints) is less than four or is not a multiple of     -}
  35. {- four then an error is returned.  NumberOfTwoBits is the   -}
  36. {- number of twobits (i.e. two bits) necessary to represent  -}
  37. {- NumPoints in base 4 (e.g. if NumPoints = 16,              -}
  38. {- NumberOfTwoBits = 3).                                     -}
  39. {-------------------------------------------------------------}
  40.  
  41. procedure MakeSinCosTable(NumPoints : integer;
  42.                       var SinTable  : TNvectorPtr;
  43.                       var CosTable  : TNvectorPtr);
  44.  
  45. {--------------------------------------------------------}
  46. {- Input: NumPoints                                     -}
  47. {- Output: SinTable, CosTable                           -}
  48. {-                                                      -}
  49. {- This procedure fills in a table with sin and cosine  -}
  50. {- values.  It is faster to pull data out of this table -}
  51. {- than it is to recalculate the sines and cosines.     -}
  52. {--------------------------------------------------------}
  53.  
  54. procedure FFT(NumberOfTwoBits : byte;
  55.               NumPoints       : integer;
  56.               Inverse         : boolean;
  57.           var XReal           : TNvectorPtr;
  58.           var XImag           : TNvectorPtr;
  59.           var SinTable        : TNvectorPtr;
  60.           var CosTable        : TNvectorPtr);
  61.  
  62. {-----------------------------------------------------}
  63. {- Input: NumberOfTwoBits, NumPoints, Inverse,       -}
  64. {-        XReal, XImag, SinTable, CosTable           -}
  65. {- Output: XReal, XImag                              -}
  66. {-                                                   -}
  67. {- This procedure implements the actual fast fourier -}
  68. {- transform routine.  The vector X, which must be   -}
  69. {- entered in twobit-inverted order, is transformed  -}
  70. {- in place.  The transformation uses the            -}
  71. {- Cooley-Tukey algorithm.                           -}
  72. {-----------------------------------------------------}
  73.  
  74. procedure RealFFT(NumPoints : integer;
  75.                   Inverse   : boolean;
  76.               var XReal     : TNvectorPtr;
  77.               var XImag     : TNvectorPtr;
  78.               var Error     : byte);
  79.  
  80.  
  81. {---------------------------------------------------------------------------}
  82. {-                                                                         -}
  83. {-    Input: NumPoints, Inverse, XReal, XImag,                             -}
  84. {-    Output: XReal, XImag, Error                                          -}
  85. {-                                                                         -}
  86. {-    Purpose:  This procedure uses the complex Fourier transform          -}
  87. {-              routine (FFT) to transform real data.  The real data       -}
  88. {-              is in the vector XReal.  Appropriate shuffling of indices  -}
  89. {-              changes the real vector into two vectors (representing     -}
  90. {-              complex data) which are only half the size of the original -}
  91. {-              vector.  Appropriate unshuffling at the end produces the   -}
  92. {-              transform of the real data.                                -}
  93. {-                                                                         -}
  94. {-  User Defined Types:                                                    -}
  95. {-         TNvector = array[0..TNArraySize] of real                        -}
  96. {-      TNvectorPtr = ^TNvector                                            -}
  97. {-                                                                         -}
  98. {- Global Variables:  NumPoints   : integer     Number of data             -}
  99. {-                                              points in X                -}
  100. {-                    Inverse     : boolean     False => forward transform -}
  101. {-                                              True ==> inverse transform -}
  102. {-                    XReal,XImag : TNvectorPtr Data points                -}
  103. {-                    Error       : byte        Indicates an error         -}
  104. {-                                                                         -}
  105. {-             Errors:  0: No Errors                                       -}
  106. {-                      1: NumPoints < 2                                   -}
  107. {-                      2: NumPoints not a power of two                    -}
  108. {-                         (or 4 for radix-4 transforms)                   -}
  109. {-                                                                         -}
  110. {---------------------------------------------------------------------------}
  111.  
  112. procedure RealConvolution(NumPoints : integer;
  113.                       var XReal     : TNvectorPtr;
  114.                       var XImag     : TNvectorPtr;
  115.                       var HReal     : TNvectorPtr;
  116.                       var Error     : byte);
  117.  
  118. {-------------------------------------------------------------------}
  119. {-                                                                 -}
  120. {-   Input: NumPoints, XReal, XImag, HReal                         -}
  121. {-   Output: XReal, XImag, Error                                   -}
  122. {-                                                                 -}
  123. {-   Purpose: This procedure performs a convolution of the         -}
  124. {-            real data XReal and HReal.  The result is returned   -}
  125. {-            in the complex vector XReal, XImag.                  -}
  126. {-                                                                 -}
  127. {- User Defined Types:                                             -}
  128. {-         TNvector = array[0..TNArraySize] of real                -}
  129. {-      TNvectorPtr = ^TNvector                                    -}
  130. {-                                                                 -}
  131. {- Global Variables:  NumPoints : integer     Number of data       -}
  132. {-                                            points in X          -}
  133. {-                    XReal     : TNvectorPtr Data points          -}
  134. {-                    HReal     : TNvectorPtr Data points          -}
  135. {-                    Error     : byte        Indicates an error   -}
  136. {-                                                                 -}
  137. {-             Errors:  0: No Errors                               -}
  138. {-                      1: NumPoints < 2                           -}
  139. {-                      2: NumPoints not a power of two            -}
  140. {-                                                                 -}
  141. {-------------------------------------------------------------------}
  142.  
  143. procedure RealCorrelation(NumPoints : integer;
  144.                       var Auto      : boolean;
  145.                       var XReal     : TNvectorPtr;
  146.                       var XImag     : TNvectorPtr;
  147.                       var HReal     : TNvectorPtr;
  148.                       var Error     : byte);
  149.  
  150. {-------------------------------------------------------------------}
  151. {-                                                                 -}
  152. {-   Input: NumPoints, Auto, XReal, XImag, HReal                   -}
  153. {-   Output: XReal, XImag, Error                                   -}
  154. {-                                                                 -}
  155. {-   Purpose: This procedure performs a correlation (auto or       -}
  156. {-            cross) of the real data XReal and HReal. The         -}
  157. {-            correlation is returned in the complex vector        -}
  158. {-            XReal, XImag.                                        -}
  159. {-                                                                 -}
  160. {- User Defined Types:                                             -}
  161. {-         TNvector = array[0..TNArraySize] OF real                -}
  162. {-      TNvectorPtr = ^TNvector                                    -}
  163. {-                                                                 -}
  164. {- Global Variables:  NumPoints : integer     Number of data       -}
  165. {-                                            points in X          -}
  166. {-                    Auto      : boolean     True => auto-        -}
  167. {-                                                    correlation  -}
  168. {-                                            False=> cross-       -}
  169. {-                                                    correlation  -}
  170. {-                    XReal     : TNvectorPtr First sample         -}
  171. {-                    HReal     : TNvectorPtr Second sample        -}
  172. {-                    Error     : byte        Indicates an error   -}
  173. {-                                                                 -}
  174. {-             Errors:  0: No Errors                               -}
  175. {-                      1: NumPoints < 2                           -}
  176. {-                      2: NumPoints not a power of two            -}
  177. {-                                                                 -}
  178. {-------------------------------------------------------------------}
  179.  
  180. procedure ComplexFFT(NumPoints : integer;
  181.                      Inverse   : boolean;
  182.                  var XReal     : TNvectorPtr;
  183.                  var XImag     : TNvectorPtr;
  184.                  var Error     : byte);
  185.  
  186. {-------------------------------------------------------------------}
  187. {-                                                                 -}
  188. {-   Input: NumPoints, Inverse, XReal, XImag                       -}
  189. {-   Output: XReal, XImag, Error                                   -}
  190. {-                                                                 -}
  191. {-   Purpose: This procedure performs a fast Fourier transform     -}
  192. {-            of the complex data XReal, XImag. The vectors XReal  -}
  193. {-            and XImag are transformed in place.                  -}
  194. {-                                                                 -}
  195. {- User Defined Types:                                             -}
  196. {-         TNvector = array[0..TNArraySize] of real                -}
  197. {-         TNvectorPtr = ^TNvector                                 -}
  198. {-                                                                 -}
  199. {- Global Variables:  NumPoints : integer      Number of data      -}
  200. {-                                             points in X         -}
  201. {-                    Inverse   : BOOLEAN      FALSE => Forward    -}
  202. {-                                                      Transform  -}
  203. {-                                             TRUE => Inverse     -}
  204. {-                                                     Transform   -}
  205. {-                    XReal,                                       -}
  206. {-                    XImag     : TNvectorPtr  Data points         -}
  207. {-                    Error     : byte         Indicates an error  -}
  208. {-                                                                 -}
  209. {-             Errors:  0: No Errors                               -}
  210. {-                      1: NumPoints < 2                           -}
  211. {-                      2: NumPoints not a power of two            -}
  212. {-                                                                 -}
  213. {-------------------------------------------------------------------}
  214.  
  215. procedure ComplexConvolution(NumPoints : integer;
  216.                          var XReal     : TNvectorPtr;
  217.                          var XImag     : TNvectorPtr;
  218.                          var HReal     : TNvectorPtr;
  219.                          var HImag     : TNvectorPtr;
  220.                          var Error     : byte);
  221.  
  222. {-------------------------------------------------------------------}
  223. {-                                                                 -}
  224. {-   Input: NumPoints, XReal, XImag, HReal, HImag                  -}
  225. {-   Output: XReal, XImag, Error                                   -}
  226. {-                                                                 -}
  227. {-   Purpose: This procedure performs a convolution of the         -}
  228. {-            data XReal, XImag and the data HReal and HImag. The  -}
  229. {-            vectors XReal, XImag, HReal and HImag are            -}
  230. {-            transformed in place.                                -}
  231. {-                                                                 -}
  232. {- User Defined Types:                                             -}
  233. {-         TNvector = array[0..TNArraySize] of real                -}
  234. {-      TNvectorPtr = ^TNvector                                    -}
  235. {-                                                                 -}
  236. {- Global Variables:  NumPoints   : integer     Number of data     -}
  237. {-                                              points in X        -}
  238. {-                    XReal,XImag : TNvectorPtr Data points        -}
  239. {-                    HReal,HImag : TNvectorPtr Data points        -}
  240. {-                    Error       : byte        Indicates an error -}
  241. {-                                                                 -}
  242. {-             Errors:  0: No Errors                               -}
  243. {-                      1: NumPoints < 2                           -}
  244. {-                      2: NumPoints not a power of two            -}
  245. {-                                                                 -}
  246. {-------------------------------------------------------------------}
  247.  
  248. procedure ComplexCorrelation(NumPoints : integer;
  249.                          var Auto      : boolean;
  250.                          var XReal     : TNvectorPtr;
  251.                          var XImag     : TNvectorPtr;
  252.                          var HReal     : TNvectorPtr;
  253.                          var HImag     : TNvectorPtr;
  254.                          var Error     : byte);
  255.  
  256. {-------------------------------------------------------------------}
  257. {-                                                                 -}
  258. {-   Input: NumPoints, Auto, XReal, XImag, HReal, HImag            -}
  259. {-   Output: XReal, XImag, Error                                   -}
  260. {-                                                                 -}
  261. {-   Purpose: This procedure performs a correlation (auto or       -}
  262. {-            cross) of the complex data XReal, XImag and the      -}
  263. {-            complex data HReal, HImag. The vectors XReal, XImag, -}
  264. {-            HReal, and HImag are transformed in place.           -}
  265. {-                                                                 -}
  266. {- User Defined Types:                                             -}
  267. {-         TNvector = array[0..TNArraySize] of real                -}
  268. {-      TNvectorPtr = ^TNvector                                    -}
  269. {-                                                                 -}
  270. {- Global Variables:  NumPoints   : integer   Number of data       -}
  271. {-                                            points in X          -}
  272. {-                    Auto        : boolean   True => auto-        -}
  273. {-                                                    correlation  -}
  274. {-                                            False=> cross-       -}
  275. {-                                                    correlation  -}
  276. {-                    XReal,XImag : TNvectorPtr First sample       -}
  277. {-                    HReal,HImag : TNvectorPtr Second sample      -}
  278. {-                    Error       : byte        Indicates an error -}
  279. {-                                                                 -}
  280. {-             Errors:  0: No Errors                               -}
  281. {-                      1: NumPoints < 2                           -}
  282. {-                      2: NumPoints not a power of two            -}
  283. {-                                                                 -}
  284. {-------------------------------------------------------------------}
  285.  
  286. implementation
  287.  
  288. procedure TestInput{(NumPoints       : integer;
  289.                  var NumberOfTwoBits : byte;
  290.                  var Error           : byte)};
  291.  
  292. type
  293.   ShortArray = array[1..6] of integer;
  294.  
  295. var
  296.   Term : integer;
  297.  
  298. const
  299.   PowersOfFour : ShortArray = (4, 16, 64, 256, 1024, 4096);
  300.  
  301. begin
  302.   Error := 2;  { Assume NumPoints not a power of four  }
  303.   if NumPoints < 4 then
  304.     Error := 1;     { NumPoints < 4  }
  305.   Term := 1;
  306.   while (Term <= 6) and (Error = 2) do
  307.   begin
  308.     if NumPoints = PowersOfFour[Term] then
  309.     begin
  310.       NumberOfTwoBits := Term;
  311.       Error := 0;  { NumPoints is a power of four  }
  312.     end;
  313.     Term := Succ(Term);
  314.   end;
  315. end; { procedure TestInput }
  316.  
  317. procedure MakeSinCosTable{(NumPoints : integer;
  318.                        var SinTable  : TNvectorPtr;
  319.                        var CosTable  : TNvectorPtr)};
  320.  
  321. var
  322.   RealFactor, ImagFactor : Float;
  323.   Term : integer;
  324.   TermMinus1 : integer;
  325.   UpperLimit : integer;
  326.  
  327. begin
  328.   RealFactor :=  Cos(Pi / (NumPoints SHR 1));
  329.   ImagFactor := -Sqrt(1 - Sqr(RealFactor));
  330.   CosTable^[0] := 1;
  331.   SinTable^[0] := 0;
  332.   CosTable^[1] := RealFactor;
  333.   SinTable^[1] := ImagFactor;
  334.   UpperLimit := 3 * NumPoints SHR 2 - 1;
  335.   for Term := 2 to UpperLimit do
  336.   begin
  337.     TermMinus1 := Term - 1;
  338.     CosTable^[Term] := CosTable^[TermMinus1] * RealFactor -
  339.                        SinTable^[TermMinus1] * ImagFactor;
  340.     SinTable^[Term] := CosTable^[TermMinus1] * ImagFactor +
  341.                        SinTable^[TermMinus1] * RealFactor;
  342.   end;
  343. end; { procedure MakeSinCosTable }
  344.  
  345. procedure FFT{(NumberOfTwoBits : byte;
  346.                NumPoints       : integer;
  347.                Inverse         : boolean;
  348.            var XReal           : TNvectorPtr;
  349.            var XImag           : TNvectorPtr;
  350.            var SinTable        : TNvectorPtr;
  351.            var CosTable        : TNvectorPtr)};
  352.  
  353. const
  354.   RootTwoOverTwo = 0.707106781186548;
  355.  
  356. var
  357.   Term : byte;
  358.   CellSeparation : integer;
  359.   NumberOfCells : integer;
  360.   NumElementsInCell : integer;
  361.   NumElInCellLess1 : integer;
  362.   NumElInCellDIV2 : integer;
  363.   NumElInCellDIV4 : integer;
  364.   CellElements : integer;
  365.   Index : integer;
  366.   RealRootOfUnity1, ImagRootOfUnity1 : Float;
  367.   RealRootOfUnity2, ImagRootOfUnity2 : Float;
  368.   RealRootOfUnity3, ImagRootOfUnity3 : Float;
  369.   Element0 : integer;
  370.   Element1 : integer;
  371.   Element2 : integer;
  372.   Element3 : integer;
  373.   RealDummy0, ImagDummy0 : Float;
  374.   RealDummy1, ImagDummy1 : Float;
  375.   RealDummy2, ImagDummy2 : Float;
  376.   RealDummy3, ImagDummy3 : Float;
  377.   RealSum02, ImagSum02 : Float;
  378.   RealDif02, ImagDif02 : Float;
  379.   RealSum13, ImagSum13 : Float;
  380.   RealDifi13, ImagDifi13 : Float;
  381.  
  382. procedure BitInvert(NumberOfTwoBits : byte;
  383.                     NumPoints       : integer;
  384.                 var XReal           : TNvectorPtr;
  385.                 var XImag           : TNvectorPtr);
  386.  
  387. {-----------------------------------------------------------}
  388. {- Input: NumberOfBits, NumPoints                          -}
  389. {- Output: XReal, XImag                                    -}
  390. {-                                                         -}
  391. {- This procedure twobit inverts the order of data in the  -}
  392. {- vector X.  Twobit inversion reverses the order of the   -}
  393. {- base 4 representation of the indices; thus 2 indices    -}
  394. {- will be switched.  For example, if there are 16 points, -}
  395. {- Index 11 (23 base 4) would be switched with Index 14    -}
  396. {- (32 base 4).  It is necessary to twobit invert the      -}
  397. {- order of the data so that the transformation comes out  -}
  398. {- in the correct order.                                   -}
  399. {-----------------------------------------------------------}
  400.  
  401. var
  402.   Term : integer;
  403.   Invert : integer;
  404.   Hold : Float;
  405.   DummyTerm, TwoBits : integer;
  406.  
  407. begin
  408.   for Term := 1 to NumPoints - 1 do
  409.   begin
  410.     Invert := 0;
  411.     DummyTerm := Term;
  412.     for TwoBits := 1 to NumberOfTwoBits do
  413.     begin
  414.       Invert := Invert SHL 2 + DummyTerm MOD 4;
  415.       DummyTerm := DummyTerm SHR 2;
  416.     end;
  417.     if Invert > Term then     { Switch the two indices  }
  418.     begin
  419.       Hold := XReal^[Invert];
  420.       XReal^[Invert] := XReal^[Term];
  421.       XReal^[Term] := Hold;
  422.       Hold := XImag^[Invert];
  423.       XImag^[Invert] := XImag^[Term];
  424.       XImag^[Term] := Hold;
  425.     end;
  426.   end;
  427. end; { procedure BitInvert }
  428.  
  429. begin { procedure FFT }
  430.   { The data must be entered in bit inverted order }
  431.   { for the transform to come out in proper order  }
  432.   BitInvert(NumberOfTwoBits, NumPoints, XReal, XImag);
  433.  
  434.   if Inverse then
  435.     { Conjugate the input  }
  436.     for Index := 0 to NumPoints - 1 do
  437.       XImag^[Index] := -XImag^[Index];
  438.  
  439.   NumberOfCells := NumPoints;
  440.   CellSeparation := 1;
  441.   for Term := 1 to NumberOfTwoBits do
  442.   begin
  443.     NumberOfCells := NumberOfCells SHR 2;
  444.     NumElementsInCell := CellSeparation;
  445.     CellSeparation := CellSeparation SHL 2;
  446.     NumElInCellLess1 := NumElementsInCell - 1;
  447.     NumElInCellDIV2 := NumElementsInCell SHR 1;
  448.     NumElInCellDIV4 := NumElInCellDIV2 SHR 1;
  449.  
  450.     { Special case: RootOfUnity1 = EXP(-i 0)  }
  451.     Element0 := 0;
  452.     while Element0 < NumPoints do
  453.     begin
  454.       { Combine the X[Element] with the element in  }
  455.       { the identical location in the next cell     }
  456.       Element1 := Element0 + NumElementsInCell;
  457.       Element2 := Element1 + NumElementsInCell;
  458.       Element3 := Element2 + NumElementsInCell;
  459.  
  460.       RealDummy0 := XReal^[Element0];
  461.       ImagDummy0 := XImag^[Element0];
  462.       RealDummy1 := XReal^[Element1];
  463.       ImagDummy1 := XImag^[Element1];
  464.       RealDummy2 := XReal^[Element2];
  465.       ImagDummy2 := XImag^[Element2];
  466.       RealDummy3 := XReal^[Element3];
  467.       ImagDummy3 := XImag^[Element3];
  468.  
  469.       RealSum02 := RealDummy0 + RealDummy2;
  470.       ImagSum02 := ImagDummy0 + ImagDummy2;
  471.       RealSum13 := RealDummy1 + RealDummy3;
  472.       ImagSum13 := ImagDummy1 + ImagDummy3;
  473.       RealDif02 := RealDummy0 - RealDummy2;
  474.       ImagDif02 := ImagDummy0 - ImagDummy2;
  475.       RealDifi13 := ImagDummy3 - ImagDummy1;
  476.       ImagDifi13 := RealDummy1 - RealDummy3;
  477.  
  478.       XReal^[Element0] := RealSum02 + RealSum13;
  479.       XImag^[Element0] := ImagSum02 + ImagSum13;
  480.       XReal^[Element1] := RealDif02 - RealDifi13;
  481.       XImag^[Element1] := ImagDif02 - ImagDifi13;
  482.       XReal^[Element2] := RealSum02 - RealSum13;
  483.       XImag^[Element2] := ImagSum02 - ImagSum13;
  484.       XReal^[Element3] := RealDif02 + RealDifi13;
  485.       XImag^[Element3] := ImagDif02 + ImagDifi13;
  486.  
  487.       Element0 := Element0 + CellSeparation;
  488.     end; { while }
  489.  
  490.     for CellElements := 1 to  NumElInCellDIV4 - 1 do
  491.     begin
  492.       Index := CellElements * NumberOfCells;
  493.       RealRootOfUnity1 := CosTable^[Index];
  494.       ImagRootOfUnity1 := SinTable^[Index];
  495.       RealRootOfUnity2 := CosTable^[2*Index];
  496.       ImagRootOfUnity2 := SinTable^[2*Index];
  497.       RealRootOfUnity3 := CosTable^[3*Index];
  498.       ImagRootOfUnity3 := SinTable^[3*Index];
  499.       Element0 := CellElements;
  500.       while Element0 < NumPoints do
  501.       begin
  502.         { Combine the X[Element] with the element in  }
  503.         { the identical location in the next cell     }
  504.         Element1 := Element0 + NumElementsInCell;
  505.         Element2 := Element1 + NumElementsInCell;
  506.         Element3 := Element2 + NumElementsInCell;
  507.  
  508.         RealDummy0 := XReal^[Element0];
  509.         ImagDummy0 := XImag^[Element0];
  510.         RealDummy1 := XReal^[Element1] * RealRootOfUnity1 -
  511.                       XImag^[Element1] * ImagRootOfUnity1;
  512.         ImagDummy1 := XReal^[Element1] * ImagRootOfUnity1 +
  513.                       XImag^[Element1] * RealRootOfUnity1;
  514.         RealDummy2 := XReal^[Element2] * RealRootOfUnity2 -
  515.                       XImag^[Element2] * ImagRootOfUnity2;
  516.         ImagDummy2 := XReal^[Element2] * ImagRootOfUnity2 +
  517.                       XImag^[Element2] * RealRootOfUnity2;
  518.         RealDummy3 := XReal^[Element3] * RealRootOfUnity3 -
  519.                       XImag^[Element3] * ImagRootOfUnity3;
  520.         ImagDummy3 := XReal^[Element3] * ImagRootOfUnity3 +
  521.                       XImag^[Element3] * RealRootOfUnity3;
  522.  
  523.         RealSum02 := RealDummy0 + RealDummy2;
  524.         ImagSum02 := ImagDummy0 + ImagDummy2;
  525.         RealSum13 := RealDummy1 + RealDummy3;
  526.         ImagSum13 := ImagDummy1 + ImagDummy3;
  527.         RealDif02 := RealDummy0 - RealDummy2;
  528.         ImagDif02 := ImagDummy0 - ImagDummy2;
  529.         RealDifi13 := ImagDummy3 - ImagDummy1;
  530.         ImagDifi13 := RealDummy1 - RealDummy3;
  531.  
  532.         XReal^[Element0] := RealSum02 + RealSum13;
  533.         XImag^[Element0] := ImagSum02 + ImagSum13;
  534.         XReal^[Element1] := RealDif02 - RealDifi13;
  535.         XImag^[Element1] := ImagDif02 - ImagDifi13;
  536.         XReal^[Element2] := RealSum02 - RealSum13;
  537.         XImag^[Element2] := ImagSum02 - ImagSum13;
  538.         XReal^[Element3] := RealDif02 + RealDifi13;
  539.         XImag^[Element3] := ImagDif02 + ImagDifi13;
  540.  
  541.         Element0 := Element0 + CellSeparation;
  542.       end; { while }
  543.     end; { for }
  544.  
  545.     { special case: RootOfUnity = EXP(-i PI/8)  }
  546.     if Term > 1 then
  547.     begin
  548.       Index := NumElInCellDIV4 * NumberOfCells;
  549.       RealRootOfUnity1 := CosTable^[Index];
  550.       ImagRootOfUnity1 := SinTable^[Index];
  551.       RealRootOfUnity3 := -ImagRootOfUnity1;
  552.       ImagRootOfUnity3 := -RealRootOfUnity1;
  553.  
  554.       Element0 := NumElInCellDIV4;
  555.       while Element0 < NumPoints do
  556.       begin
  557.         { Combine the X[Element] with the element in  }
  558.         { the identical location in the next cell     }
  559.         Element1 := Element0 + NumElementsInCell;
  560.         Element2 := Element1 + NumElementsInCell;
  561.         Element3 := Element2 + NumElementsInCell;
  562.  
  563.         RealDummy0 := XReal^[Element0];
  564.         ImagDummy0 := XImag^[Element0];
  565.         RealDummy1 := XReal^[Element1] * RealRootOfUnity1 -
  566.                       XImag^[Element1] * ImagRootOfUnity1;
  567.         ImagDummy1 := XReal^[Element1] * ImagRootOfUnity1 +
  568.                       XImag^[Element1] * RealRootOfUnity1;
  569.         RealDummy2 := RootTwoOverTwo * (XReal^[Element2] + XImag^[Element2]);
  570.         ImagDummy2 := RootTwoOverTwo * (XImag^[Element2] - XReal^[Element2]);
  571.         RealDummy3 := XReal^[Element3] * RealRootOfUnity3 -
  572.                       XImag^[Element3] * ImagRootOfUnity3;
  573.         ImagDummy3 := XReal^[Element3] * ImagRootOfUnity3 +
  574.                       XImag^[Element3] * RealRootOfUnity3;
  575.  
  576.         RealSum02 := RealDummy0 + RealDummy2;
  577.         ImagSum02 := ImagDummy0 + ImagDummy2;
  578.         RealSum13 := RealDummy1 + RealDummy3;
  579.         ImagSum13 := ImagDummy1 + ImagDummy3;
  580.         RealDif02 := RealDummy0 - RealDummy2;
  581.         ImagDif02 := ImagDummy0 - ImagDummy2;
  582.         RealDifi13 := ImagDummy3 - ImagDummy1;
  583.         ImagDifi13 := RealDummy1 - RealDummy3;
  584.  
  585.         XReal^[Element0] := RealSum02 + RealSum13;
  586.         XImag^[Element0] := ImagSum02 + ImagSum13;
  587.         XReal^[Element1] := RealDif02 - RealDifi13;
  588.         XImag^[Element1] := ImagDif02 - ImagDifi13;
  589.         XReal^[Element2] := RealSum02 - RealSum13;
  590.         XImag^[Element2] := ImagSum02 - ImagSum13;
  591.         XReal^[Element3] := RealDif02 + RealDifi13;
  592.         XImag^[Element3] := ImagDif02 + ImagDifi13;
  593.  
  594.         Element0 := Element0 + CellSeparation;
  595.       end; { while }
  596.     end;
  597.  
  598.     for CellElements := NumElInCellDIV4 + 1 to NumElInCellDIV2 - 1 do
  599.     begin
  600.       Index := CellElements * NumberOfCells;
  601.       RealRootOfUnity1 := CosTable^[Index];
  602.       ImagRootOfUnity1 := SinTable^[Index];
  603.       RealRootOfUnity2 := CosTable^[2*Index];
  604.       ImagRootOfUnity2 := SinTable^[2*Index];
  605.       RealRootOfUnity3 := CosTable^[3*Index];
  606.       ImagRootOfUnity3 := SinTable^[3*Index];
  607.       Element0 := CellElements;
  608.       while Element0 < NumPoints do
  609.       begin
  610.         { Combine the X[Element] with the element in  }
  611.         { the identical location in the next cell     }
  612.         Element1 := Element0 + NumElementsInCell;
  613.         Element2 := Element1 + NumElementsInCell;
  614.         Element3 := Element2 + NumElementsInCell;
  615.  
  616.         RealDummy0 := XReal^[Element0];
  617.         ImagDummy0 := XImag^[Element0];
  618.         RealDummy1 := XReal^[Element1] * RealRootOfUnity1 -
  619.                       XImag^[Element1] * ImagRootOfUnity1;
  620.         ImagDummy1 := XReal^[Element1] * ImagRootOfUnity1 +
  621.                       XImag^[Element1] * RealRootOfUnity1;
  622.         RealDummy2 := XReal^[Element2] * RealRootOfUnity2 -
  623.                       XImag^[Element2] * ImagRootOfUnity2;
  624.         ImagDummy2 := XReal^[Element2] * ImagRootOfUnity2 +
  625.                       XImag^[Element2] * RealRootOfUnity2;
  626.         RealDummy3 := XReal^[Element3] * RealRootOfUnity3 -
  627.                       XImag^[Element3] * ImagRootOfUnity3;
  628.         ImagDummy3 := XReal^[Element3] * ImagRootOfUnity3 +
  629.                       XImag^[Element3] * RealRootOfUnity3;
  630.  
  631.         RealSum02 := RealDummy0 + RealDummy2;
  632.         ImagSum02 := ImagDummy0 + ImagDummy2;
  633.         RealSum13 := RealDummy1 + RealDummy3;
  634.         ImagSum13 := ImagDummy1 + ImagDummy3;
  635.         RealDif02 := RealDummy0 - RealDummy2;
  636.         ImagDif02 := ImagDummy0 - ImagDummy2;
  637.         RealDifi13 := ImagDummy3 - ImagDummy1;
  638.         ImagDifi13 := RealDummy1 - RealDummy3;
  639.  
  640.         XReal^[Element0] := RealSum02 + RealSum13;
  641.         XImag^[Element0] := ImagSum02 + ImagSum13;
  642.         XReal^[Element1] := RealDif02 - RealDifi13;
  643.         XImag^[Element1] := ImagDif02 - ImagDifi13;
  644.         XReal^[Element2] := RealSum02 - RealSum13;
  645.         XImag^[Element2] := ImagSum02 - ImagSum13;
  646.         XReal^[Element3] := RealDif02 + RealDifi13;
  647.         XImag^[Element3] := ImagDif02 + ImagDifi13;
  648.  
  649.         Element0 := Element0 + CellSeparation;
  650.       end; { while }
  651.     end; { for }
  652.  
  653.     { Special case: RootOfUnity1 := EXP(-i PI/4)  }
  654.     if Term > 1 then
  655.     begin
  656.       Element0 := NumElInCellDIV2;
  657.       while Element0 < NumPoints do
  658.       begin
  659.         { Combine the X[Element] with the element in  }
  660.         { the identical location in the next cell     }
  661.         Element1 := Element0 + NumElementsInCell;
  662.         Element2 := Element1 + NumElementsInCell;
  663.         Element3 := Element2 + NumElementsInCell;
  664.  
  665.         RealDummy0 := XReal^[Element0];
  666.         ImagDummy0 := XImag^[Element0];
  667.         RealDummy1 := RootTwoOverTwo * (XReal^[Element1] + XImag^[Element1]);
  668.         ImagDummy1 := RootTwoOverTwo * (XImag^[Element1] - XReal^[Element1]);
  669.         RealDummy2 :=  XImag^[Element2];
  670.         ImagDummy2 := -XReal^[Element2];
  671.         RealDummy3 := -RootTwoOverTwo * (XReal^[Element3] - XImag^[Element3]);
  672.         ImagDummy3 := -RootTwoOverTwo * (XReal^[Element3] + XImag^[Element3]);
  673.  
  674.         RealSum02 := RealDummy0 + RealDummy2;
  675.         ImagSum02 := ImagDummy0 + ImagDummy2;
  676.         RealSum13 := RealDummy1 + RealDummy3;
  677.         ImagSum13 := ImagDummy1 + ImagDummy3;
  678.         RealDif02 := RealDummy0 - RealDummy2;
  679.         ImagDif02 := ImagDummy0 - ImagDummy2;
  680.         RealDifi13 := ImagDummy3 - ImagDummy1;
  681.         ImagDifi13 := RealDummy1 - RealDummy3;
  682.  
  683.         XReal^[Element0] := RealSum02 + RealSum13;
  684.         XImag^[Element0] := ImagSum02 + ImagSum13;
  685.         XReal^[Element1] := RealDif02 - RealDifi13;
  686.         XImag^[Element1] := ImagDif02 - ImagDifi13;
  687.         XReal^[Element2] := RealSum02 - RealSum13;
  688.         XImag^[Element2] := ImagSum02 - ImagSum13;
  689.         XReal^[Element3] := RealDif02 + RealDifi13;
  690.         XImag^[Element3] := ImagDif02 + ImagDifi13;
  691.  
  692.         Element0 := Element0 + CellSeparation;
  693.       end; { while }
  694.     end;
  695.  
  696.     for CellElements := NumElInCellDIV2 + 1 to
  697.         (NumElementsInCell - NumElInCellDIV4 - 1) do
  698.     begin
  699.       Index := CellElements * NumberOfCells;
  700.       RealRootOfUnity1 := CosTable^[Index];
  701.       ImagRootOfUnity1 := SinTable^[Index];
  702.       RealRootOfUnity2 := CosTable^[2*Index];
  703.       ImagRootOfUnity2 := SinTable^[2*Index];
  704.       RealRootOfUnity3 := CosTable^[3*Index];
  705.       ImagRootOfUnity3 := SinTable^[3*Index];
  706.       Element0 := CellElements;
  707.       while Element0 < NumPoints do
  708.       begin
  709.         { Combine the X[Element] with the element in  }
  710.         { the identical location in the next cell     }
  711.         Element1 := Element0 + NumElementsInCell;
  712.         Element2 := Element1 + NumElementsInCell;
  713.         Element3 := Element2 + NumElementsInCell;
  714.  
  715.         RealDummy0 := XReal^[Element0];
  716.         ImagDummy0 := XImag^[Element0];
  717.         RealDummy1 := XReal^[Element1] * RealRootOfUnity1 -
  718.                       XImag^[Element1] * ImagRootOfUnity1;
  719.         ImagDummy1 := XReal^[Element1] * ImagRootOfUnity1 +
  720.                       XImag^[Element1] * RealRootOfUnity1;
  721.         RealDummy2 := XReal^[Element2] * RealRootOfUnity2 -
  722.                       XImag^[Element2] * ImagRootOfUnity2;
  723.         ImagDummy2 := XReal^[Element2] * ImagRootOfUnity2 +
  724.                       XImag^[Element2] * RealRootOfUnity2;
  725.         RealDummy3 := XReal^[Element3] * RealRootOfUnity3 -
  726.                       XImag^[Element3] * ImagRootOfUnity3;
  727.         ImagDummy3 := XReal^[Element3] * ImagRootOfUnity3 +
  728.                       XImag^[Element3] * RealRootOfUnity3;
  729.  
  730.         RealSum02 := RealDummy0 + RealDummy2;
  731.         ImagSum02 := ImagDummy0 + ImagDummy2;
  732.         RealSum13 := RealDummy1 + RealDummy3;
  733.         ImagSum13 := ImagDummy1 + ImagDummy3;
  734.         RealDif02 := RealDummy0 - RealDummy2;
  735.         ImagDif02 := ImagDummy0 - ImagDummy2;
  736.         RealDifi13 := ImagDummy3 - ImagDummy1;
  737.         ImagDifi13 := RealDummy1 - RealDummy3;
  738.  
  739.         XReal^[Element0] := RealSum02 + RealSum13;
  740.         XImag^[Element0] := ImagSum02 + ImagSum13;
  741.         XReal^[Element1] := RealDif02 - RealDifi13;
  742.         XImag^[Element1] := ImagDif02 - ImagDifi13;
  743.         XReal^[Element2] := RealSum02 - RealSum13;
  744.         XImag^[Element2] := ImagSum02 - ImagSum13;
  745.         XReal^[Element3] := RealDif02 + RealDifi13;
  746.         XImag^[Element3] := ImagDif02 + ImagDifi13;
  747.  
  748.         Element0 := Element0 + CellSeparation;
  749.       end; { while }
  750.     end; { for }
  751.  
  752.     { Special case: RootOfUnity = EXP(-i 3*PI/8)  }
  753.     if Term > 1 then
  754.     begin
  755.       Element0 := NumElementsInCell - NumElInCellDIV4;
  756.       Index := Element0 * NumberOfCells;
  757.       RealRootOfUnity1 := CosTable^[Index];
  758.       ImagRootOfUnity1 := SinTable^[Index];
  759.       RealRootOfUnity3 := ImagRootOfUnity1;
  760.       ImagRootOfUnity3 := RealRootOfUnity1;
  761.  
  762.       while Element0 < NumPoints do
  763.       begin
  764.         { Combine the X[Element] with the element in  }
  765.         { the identical location in the next cell     }
  766.         Element1 := Element0 + NumElementsInCell;
  767.         Element2 := Element1 + NumElementsInCell;
  768.         Element3 := Element2 + NumElementsInCell;
  769.  
  770.         RealDummy0 := XReal^[Element0];
  771.         ImagDummy0 := XImag^[Element0];
  772.         RealDummy1 := XReal^[Element1] * RealRootOfUnity1 -
  773.                       XImag^[Element1] * ImagRootOfUnity1;
  774.         ImagDummy1 := XReal^[Element1] * ImagRootOfUnity1 +
  775.                       XImag^[Element1] * RealRootOfUnity1;
  776.         RealDummy2 := -RootTwoOverTwo * (XReal^[Element2] - XImag^[Element2]);
  777.         ImagDummy2 := -RootTwoOverTwo * (XReal^[Element2] + XImag^[Element2]);
  778.         RealDummy3 := XReal^[Element3] * RealRootOfUnity3 -
  779.                       XImag^[Element3] * ImagRootOfUnity3;
  780.         ImagDummy3 := XReal^[Element3] * ImagRootOfUnity3 +
  781.                       XImag^[Element3] * RealRootOfUnity3;
  782.  
  783.         RealSum02 := RealDummy0 + RealDummy2;
  784.         ImagSum02 := ImagDummy0 + ImagDummy2;
  785.         RealSum13 := RealDummy1 + RealDummy3;
  786.         ImagSum13 := ImagDummy1 + ImagDummy3;
  787.         RealDif02 := RealDummy0 - RealDummy2;
  788.         ImagDif02 := ImagDummy0 - ImagDummy2;
  789.         RealDifi13 := ImagDummy3 - ImagDummy1;
  790.         ImagDifi13 := RealDummy1 - RealDummy3;
  791.  
  792.         XReal^[Element0] := RealSum02 + RealSum13;
  793.         XImag^[Element0] := ImagSum02 + ImagSum13;
  794.         XReal^[Element1] := RealDif02 - RealDifi13;
  795.         XImag^[Element1] := ImagDif02 - ImagDifi13;
  796.         XReal^[Element2] := RealSum02 - RealSum13;
  797.         XImag^[Element2] := ImagSum02 - ImagSum13;
  798.         XReal^[Element3] := RealDif02 + RealDifi13;
  799.         XImag^[Element3] := ImagDif02 + ImagDifi13;
  800.  
  801.         Element0 := Element0 + CellSeparation;
  802.       end; { while }
  803.     end;
  804.  
  805.     for CellElements := (NumElementsInCell - NumElInCellDIV4 + 1) to
  806.                                                  NumElInCellLess1 do
  807.     begin
  808.       Index := CellElements * NumberOfCells;
  809.       RealRootOfUnity1 := CosTable^[Index];
  810.       ImagRootOfUnity1 := SinTable^[Index];
  811.       RealRootOfUnity2 := CosTable^[2*Index];
  812.       ImagRootOfUnity2 := SinTable^[2*Index];
  813.       RealRootOfUnity3 := CosTable^[3*Index];
  814.       ImagRootOfUnity3 := SinTable^[3*Index];
  815.       Element0 := CellElements;
  816.       while Element0 < NumPoints do
  817.       begin
  818.         { Combine the X[Element] with the element in  }
  819.         { the identical location in the next cell     }
  820.         Element1 := Element0 + NumElementsInCell;
  821.         Element2 := Element1 + NumElementsInCell;
  822.         Element3 := Element2 + NumElementsInCell;
  823.  
  824.         RealDummy0 := XReal^[Element0];
  825.         ImagDummy0 := XImag^[Element0];
  826.         RealDummy1 := XReal^[Element1] * RealRootOfUnity1 -
  827.                       XImag^[Element1] * ImagRootOfUnity1;
  828.         ImagDummy1 := XReal^[Element1] * ImagRootOfUnity1 +
  829.                       XImag^[Element1] * RealRootOfUnity1;
  830.         RealDummy2 := XReal^[Element2] * RealRootOfUnity2 -
  831.                       XImag^[Element2] * ImagRootOfUnity2;
  832.         ImagDummy2 := XReal^[Element2] * ImagRootOfUnity2 +
  833.                       XImag^[Element2] * RealRootOfUnity2;
  834.         RealDummy3 := XReal^[Element3] * RealRootOfUnity3 -
  835.                       XImag^[Element3] * ImagRootOfUnity3;
  836.         ImagDummy3 := XReal^[Element3] * ImagRootOfUnity3 +
  837.                       XImag^[Element3] * RealRootOfUnity3;
  838.  
  839.         RealSum02 := RealDummy0 + RealDummy2;
  840.         ImagSum02 := ImagDummy0 + ImagDummy2;
  841.         RealSum13 := RealDummy1 + RealDummy3;
  842.         ImagSum13 := ImagDummy1 + ImagDummy3;
  843.         RealDif02 := RealDummy0 - RealDummy2;
  844.         ImagDif02 := ImagDummy0 - ImagDummy2;
  845.         RealDifi13 := ImagDummy3 - ImagDummy1;
  846.         ImagDifi13 := RealDummy1 - RealDummy3;
  847.  
  848.         XReal^[Element0] := RealSum02 + RealSum13;
  849.         XImag^[Element0] := ImagSum02 + ImagSum13;
  850.         XReal^[Element1] := RealDif02 - RealDifi13;
  851.         XImag^[Element1] := ImagDif02 - ImagDifi13;
  852.         XReal^[Element2] := RealSum02 - RealSum13;
  853.         XImag^[Element2] := ImagSum02 - ImagSum13;
  854.         XReal^[Element3] := RealDif02 + RealDifi13;
  855.         XImag^[Element3] := ImagDif02 + ImagDifi13;
  856.  
  857.         Element0 := Element0 + CellSeparation;
  858.       end; { while }
  859.     end; { for }
  860.   end; { for }
  861.  
  862.   {----------------------------------------------------}
  863.   {-  Divide all the values of the transformation     -}
  864.   {-  by the square root of NumPoints. If taking the  -}
  865.   {-  inverse, conjugate the output.                  -}
  866.   {----------------------------------------------------}
  867.  
  868.   if Inverse then
  869.     ImagDummy1 := -1 / Sqrt(NumPoints)
  870.   else
  871.     ImagDummy1 :=  1 / Sqrt(NumPoints);
  872.   RealDummy1 := ABS(ImagDummy1);
  873.   for Element0 := 0 to NumPoints - 1 do
  874.   begin
  875.     XReal^[Element0] := XReal^[Element0] * RealDummy1;
  876.     XImag^[Element0] := XImag^[Element0] * ImagDummy1;
  877.   end;
  878. end; { procedure FFT }
  879.  
  880. {$I FFT.inc} { Include procedure code }
  881.  
  882. end. { FFT87B4 }
  883.