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