home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 9 / 09.iso / l / l042 / 2.ddi / CHAP10.ARC / FFTB2.PAS < prev    next >
Encoding:
Pascal/Delphi Source File  |  1987-12-30  |  29.6 KB  |  647 lines
  1. unit FFTB2;
  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 non 8087 version                            -}
  10. {-                                                                          -}
  11. {----------------------------------------------------------------------------}
  12.  
  13. {$N-} { Doesn't use an 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.  
  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.  
  291. type
  292.   ShortArray = array[1..13] of integer;
  293.  
  294. var
  295.   Term : integer;
  296.  
  297. const
  298.   PowersOfTwo : ShortArray = (2, 4, 8, 16, 32, 64, 128, 256,
  299.                               512, 1024, 2048, 4096, 8192);
  300.  
  301. begin
  302.   Error := 2; { Assume NumPoints not a power of two  }
  303.   if NumPoints < 2 then
  304.     Error := 1;     { NumPoints < 2  }
  305.   Term := 1;
  306.   while (Term <= 13) and (Error = 2) do
  307.   begin
  308.     if NumPoints = PowersOfTwo[Term] then
  309.     begin
  310.       NumberOfBits := Term;
  311.       Error := 0;  { NumPoints is a power of two }
  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(2 * Pi / NumPoints);
  329.   ImagFactor := -Sqrt(1 - Sqr(RealFactor));
  330.   CosTable^[0] := 1;
  331.   SinTable^[0] := 0;
  332.   CosTable^[1] := RealFactor;
  333.   SinTable^[1] := ImagFactor;
  334.   UpperLimit := NumPoints shr 1 - 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{(NumberOfBits : 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.   NumElInCellSHR1 : integer;
  363.   NumElInCellSHR2 : integer;
  364.   CosTerm, SumTerm, DifTerm : Float;
  365.   Element : integer;
  366.   CellElements : integer;
  367.   ElementInNextCell : integer;
  368.   Index : integer;
  369.   RealDummy, ImagDummy : Float;
  370.   Dummy1, Dummy2 : Float;
  371.  
  372. procedure BitInvert(NumberOfBits : byte;
  373.                     NumPoints    : integer;
  374.                 var XReal        : TNvectorPtr;
  375.                 var XImag        : TNvectorPtr);
  376.  
  377. {-----------------------------------------------------------}
  378. {- Input: NumberOfBits, NumPoints                          -}
  379. {- Output: XReal, XImag                                    -}
  380. {-                                                         -}
  381. {- This procedure bit inverts the order of data in the     -}
  382. {- vector X.  Bit inversion reverses the order of the      -}
  383. {- binary representation of the indices; thus 2 indices    -}
  384. {- will be switched.  For example, if there are 16 points, -}
  385. {- Index 7 (binary 0111) would be switched with Index 14   -}
  386. {- (binary 1110).  It is necessary to bit invert the order -}
  387. {- of the data so that the transformation comes out in the -}
  388. {- correct order.                                          -}
  389. {-----------------------------------------------------------}
  390.  
  391. var
  392.   Term : integer;
  393.   Invert : integer;
  394.   Hold : Float;
  395.   NumPointsDiv2, K : integer;
  396.  
  397. begin
  398.   NumPointsDiv2 := NumPoints shr 1;
  399.   Invert := 0;
  400.   for Term := 0 to NumPoints - 2 do
  401.   begin
  402.     if Term < Invert then   { Switch these two indices  }
  403.     begin
  404.       Hold := XReal^[Invert];
  405.       XReal^[Invert] := XReal^[Term];
  406.       XReal^[Term] := Hold;
  407.       Hold := XImag^[Invert];
  408.       XImag^[Invert] := XImag^[Term];
  409.       XImag^[Term] := Hold;
  410.     end;
  411.     K := NumPointsDiv2;
  412.     while K <= Invert do
  413.     begin
  414.       Invert := Invert - K;
  415.       K := K shr 1;
  416.     end;
  417.     Invert := Invert + K;
  418.   end;
  419. end; { procedure BitInvert }
  420.  
  421. begin { procedure FFT }
  422.   { The data must be entered in bit inverted order }
  423.   { for the transform to come out in proper order  }
  424.   BitInvert(NumberOfBits, NumPoints, XReal, XImag);
  425.  
  426.   if Inverse then
  427.     { Conjugate the input  }
  428.     for Element := 0 to NumPoints - 1 do
  429.       XImag^[Element] := -XImag^[Element];
  430.  
  431.   NumberOfCells := NumPoints;
  432.   CellSeparation := 1;
  433.   for Term := 1 to NumberOfBits do
  434.   begin
  435.     { NumberOfCells halves; equals 2^(NumberOfBits - Term)  }
  436.     NumberOfCells := NumberOfCells shr 1;
  437.     { NumElementsInCell doubles; equals 2^(Term-1)  }
  438.     NumElementsInCell := CellSeparation;
  439.     { CellSeparation doubles; equals 2^Term  }
  440.     CellSeparation := CellSeparation SHL 1;
  441.     NumElInCellLess1 := NumElementsInCell - 1;
  442.     NumElInCellSHR1 := NumElementsInCell shr 1;
  443.     NumElInCellSHR2 := NumElInCellSHR1 shr 1;
  444.  
  445.     { Special case: RootOfUnity = EXP(-i 0)  }
  446.     Element := 0;
  447.     while Element < NumPoints do
  448.     begin
  449.       { Combine the X[Element] with the element in  }
  450.       { the identical location in the next cell     }
  451.       ElementInNextCell := Element + NumElementsInCell;
  452.       RealDummy := XReal^[ElementInNextCell];
  453.       ImagDummy := XImag^[ElementInNextCell];
  454.       XReal^[ElementInNextCell] := XReal^[Element] - RealDummy;
  455.       XImag^[ElementInNextCell] := XImag^[Element] - ImagDummy;
  456.       XReal^[Element] := XReal^[Element] + RealDummy;
  457.       XImag^[Element] := XImag^[Element] + ImagDummy;
  458.       Element := Element + CellSeparation;
  459.     end;
  460.  
  461.     for CellElements := 1 to NumElInCellSHR2 - 1 do
  462.     begin
  463.       Index := CellElements * NumberOfCells;
  464.       CosTerm := CosTable^[Index];
  465.       SumTerm := CosTerm + SinTable^[Index];
  466.       DifTerm := SinTable^[Index] - CosTerm;
  467.       Element := CellElements;
  468.  
  469.       while Element < NumPoints do
  470.       begin
  471.         { Combine the X[Element] with the element in  }
  472.         { the identical location in the next cell     }
  473.         ElementInNextCell := Element + NumElementsInCell;
  474.         Dummy1 := (XReal^[ElementInNextCell] +
  475.         XImag^[ElementInNextCell]) * CosTerm;
  476.         Dummy2 := XReal^[ElementInNextCell] * DifTerm;
  477.         RealDummy := Dummy1 - XImag^[ElementInNextCell] * SumTerm;
  478.         ImagDummy := Dummy1 + Dummy2;;
  479.         XReal^[ElementInNextCell] := XReal^[Element] - RealDummy;
  480.         XImag^[ElementInNextCell] := XImag^[Element] - ImagDummy;
  481.         XReal^[Element] := XReal^[Element] + RealDummy;
  482.         XImag^[Element] := XImag^[Element] + ImagDummy;
  483.         Element := Element + CellSeparation;
  484.       end;
  485.     end; { for }
  486.  
  487.     { Special case: RootOfUnity = EXP(-i PI/4)  }
  488.     if Term > 2 then
  489.     begin
  490.       Element := NumElInCellSHR2;
  491.       while Element < NumPoints do
  492.       begin
  493.         { Combine the X[Element] with the element in  }
  494.         { the identical location in the next cell     }
  495.         ElementInNextCell := Element + NumElementsInCell;
  496.         RealDummy := RootTwoOverTwo * (XReal^[ElementInNextCell] +
  497.                                        XImag^[ElementInNextCell]);
  498.         ImagDummy := RootTwoOverTwo * (XImag^[ElementInNextCell] -
  499.                                        XReal^[ElementInNextCell]);
  500.         XReal^[ElementInNextCell] := XReal^[Element] - RealDummy;
  501.         XImag^[ElementInNextCell] := XImag^[Element] - ImagDummy;
  502.         XReal^[Element] := XReal^[Element] + RealDummy;
  503.         XImag^[Element] := XImag^[Element] + ImagDummy;
  504.         Element := Element + CellSeparation;
  505.       end;
  506.     end;
  507.  
  508.     for CellElements := NumElInCellSHR2 + 1 to NumElInCellSHR1 - 1 do
  509.     begin
  510.       Index := CellElements * NumberOfCells;
  511.       CosTerm := CosTable^[Index];
  512.       SumTerm := SinTable^[Index] + CosTerm;
  513.       DifTerm := SinTable^[Index] - CosTerm;
  514.       Element := CellElements;
  515.       while Element < NumPoints do
  516.       begin
  517.         { Combine the X[Element] with the element in  }
  518.         { the identical location in the next cell     }
  519.         ElementInNextCell := Element + NumElementsInCell;
  520.         Dummy1 := (XReal^[ElementInNextCell] +
  521.                    XImag^[ElementInNextCell]) * CosTerm;
  522.         Dummy2 := XReal^[ElementInNextCell] * DifTerm;
  523.         RealDummy := Dummy1 - XImag^[ElementInNextCell] * SumTerm;
  524.         ImagDummy := Dummy1 + Dummy2;;
  525.         XReal^[ElementInNextCell] := XReal^[Element] - RealDummy;
  526.         XImag^[ElementInNextCell] := XImag^[Element] - ImagDummy;
  527.         XReal^[Element] := XReal^[Element] + RealDummy;
  528.         XImag^[Element] := XImag^[Element] + ImagDummy;
  529.         Element := Element + CellSeparation;
  530.       end;
  531.     end; { for }
  532.  
  533.     { Special case: RootOfUnity = EXP(-i PI/2)  }
  534.     if Term > 1 then
  535.     begin
  536.       Element := NumElInCellSHR1;
  537.       while Element < NumPoints do
  538.       begin
  539.         { Combine the X[Element] with the element in  }
  540.         { the identical location in the next cell     }
  541.         ElementInNextCell := Element + NumElementsInCell;
  542.         RealDummy :=  XImag^[ElementInNextCell];
  543.         ImagDummy := -XReal^[ElementInNextCell];
  544.         XReal^[ElementInNextCell] := XReal^[Element] - RealDummy;
  545.         XImag^[ElementInNextCell] := XImag^[Element] - ImagDummy;
  546.         XReal^[Element] := XReal^[Element] + RealDummy;
  547.         XImag^[Element] := XImag^[Element] + ImagDummy;
  548.         Element := Element + CellSeparation;
  549.       end;
  550.     end;
  551.  
  552.     for CellElements := NumElInCellSHR1 + 1 to
  553.                         NumElementsInCell - NumElInCellSHR2 - 1 do
  554.     begin
  555.       Index := CellElements * NumberOfCells;
  556.       CosTerm := CosTable^[Index];
  557.       SumTerm := SinTable^[Index] + CosTerm;
  558.       DifTerm := SinTable^[Index] - CosTerm;
  559.       Element := CellElements;
  560.       while Element < NumPoints do
  561.       begin
  562.         { Combine the X[Element] with the element in  }
  563.         { the identical location in the next cell     }
  564.         ElementInNextCell := Element + NumElementsInCell;
  565.         Dummy1 := (XReal^[ElementInNextCell] +
  566.                    XImag^[ElementInNextCell]) * CosTerm;
  567.         Dummy2 := XReal^[ElementInNextCell] * DifTerm;
  568.         RealDummy := Dummy1 - XImag^[ElementInNextCell] * SumTerm;
  569.         ImagDummy := Dummy1 + Dummy2;;
  570.         XReal^[ElementInNextCell] := XReal^[Element] - RealDummy;
  571.         XImag^[ElementInNextCell] := XImag^[Element] - ImagDummy;
  572.         XReal^[Element] := XReal^[Element] + RealDummy;
  573.         XImag^[Element] := XImag^[Element] + ImagDummy;
  574.         Element := Element + CellSeparation;
  575.       end;
  576.     end; { for }
  577.  
  578.     { Special case: RootOfUnity = EXP(-i 3PI/4)  }
  579.     if Term > 2 then
  580.     begin
  581.       Element := NumElementsInCell - NumElInCellSHR2;
  582.       while Element < NumPoints do
  583.       begin
  584.         { Combine the X[Element] with the element in  }
  585.         { the identical location in the next cell     }
  586.         ElementInNextCell := Element + NumElementsInCell;
  587.         RealDummy := -RootTwoOverTwo * (XReal^[ElementInNextCell] -
  588.                                         XImag^[ElementInNextCell]);
  589.         ImagDummy := -RootTwoOverTwo * (XReal^[ElementInNextCell] +
  590.                                         XImag^[ElementInNextCell]);
  591.         XReal^[ElementInNextCell] := XReal^[Element] - RealDummy;
  592.         XImag^[ElementInNextCell] := XImag^[Element] - ImagDummy;
  593.         XReal^[Element] := XReal^[Element] + RealDummy;
  594.         XImag^[Element] := XImag^[Element] + ImagDummy;
  595.         Element := Element + CellSeparation;
  596.       end;
  597.     end;
  598.  
  599.     for CellElements := NumElementsInCell - NumElInCellSHR2 + 1 to
  600.                                              NumElInCellLess1 do
  601.     begin
  602.       Index := CellElements * NumberOfCells;
  603.       CosTerm := CosTable^[Index];
  604.       SumTerm := SinTable^[Index] + CosTerm;
  605.       DifTerm := SinTable^[Index] - CosTerm;
  606.       Element := CellElements;
  607.       while Element < NumPoints do
  608.       begin
  609.         { Combine the X[Element] with the element in  }
  610.         { the identical location in the next cell     }
  611.         ElementInNextCell := Element + NumElementsInCell;
  612.         Dummy1 := (XReal^[ElementInNextCell] +
  613.                    XImag^[ElementInNextCell]) * CosTerm;
  614.         Dummy2 := XReal^[ElementInNextCell] * DifTerm;
  615.         RealDummy := Dummy1 - XImag^[ElementInNextCell] * SumTerm;
  616.         ImagDummy := Dummy1 + Dummy2;;
  617.         XReal^[ElementInNextCell] := XReal^[Element] - RealDummy;
  618.         XImag^[ElementInNextCell] := XImag^[Element] - ImagDummy;
  619.         XReal^[Element] := XReal^[Element] + RealDummy;
  620.         XImag^[Element] := XImag^[Element] + ImagDummy;
  621.         Element := Element + CellSeparation;
  622.       end;
  623.     end; { for }
  624.   end;
  625.  
  626.   {----------------------------------------------------}
  627.   {-  Divide all the values of the transformation     -}
  628.   {-  by the square root of NumPoints. If taking the  -}
  629.   {-  inverse, conjugate the output.                  -}
  630.   {----------------------------------------------------}
  631.  
  632.   if Inverse then
  633.     ImagDummy := -1 / Sqrt(NumPoints)
  634.   else
  635.     ImagDummy :=  1 / Sqrt(NumPoints);
  636.   RealDummy := ABS(ImagDummy);
  637.   for Element := 0 to NumPoints - 1 do
  638.   begin
  639.     XReal^[Element] := XReal^[Element] * RealDummy;
  640.     XImag^[Element] := XImag^[Element] * ImagDummy;
  641.   end;
  642. end; { procedure FFT }
  643.  
  644. {$I FFT.inc}     { Include procedure code }
  645.  
  646. end. { FFTB2 }
  647.