home *** CD-ROM | disk | FTP | other *** search
/ Hacker Chronicles 2 / HACKER2.BIN / 156.TRANSFOR.PAS < prev    next >
Pascal/Delphi Source File  |  1988-05-24  |  16KB  |  439 lines
  1. UNIT Transforms;
  2.  
  3. INTERFACE
  4.  
  5. USES
  6.    DOS,
  7. {$IFDEF DOSCrt}
  8.    DOSCrt,
  9. {$ELSE}
  10.    Crt,
  11. {$ENDIF}
  12.    TextOps,
  13.    Global;
  14.  
  15. CONST
  16.    Fourward = TRUE;          {Forward transform}
  17.    Inverse  = FALSE;         {Inverse transform}
  18.    Four     = TRUE;          {Fourier transform}
  19.    Hart     = FALSE;         {Hartley transform}
  20.  
  21.  
  22. PROCEDURE Transform ( kind      : BOOLEAN;
  23.                       direction : BOOLEAN);
  24.  
  25. {----------------------------------------------------------------------------}
  26. {-  NOTE: Transform is currently only set up for FFT's, not for FHT's;      -}
  27. {-  however, routines for both transforms are included in this file.        -}
  28. {-                                                                          -}
  29. {-  kind = Four : perform fast Fourier transform                            -}
  30. {-  kind = Hart : perform fast Hartley transform                            -}
  31. {-                                                                          -}
  32. {-  direction = Fourward : forward transform from time to frequency domain  -}
  33. {-  direction = Inverse  : reverse transform from frequency to time domain  -}
  34. {----------------------------------------------------------------------------}
  35.  
  36. (****************************************************************************)
  37.  
  38. IMPLEMENTATION
  39.  
  40.  
  41. TYPE
  42.    polar = RECORD
  43.               Mag : REAL;
  44.               Ph  : REAL;
  45.            END;   {RECORD}
  46.  
  47. VAR
  48.    n           : INTEGER;                   (* n = 2*NumFreqs-1             *)
  49.    power_index : 1..MaxPower;               (* Power of two of FFT          *)
  50.    delta       : REAL;                      (* delta = 1/beta               *)
  51.  
  52. {----------------------------------------------------------------------------}
  53. {-    Arg = the counterclockwise angle between the positive half of the     -}
  54. {-    real axis, and the ray from the origin to the complex number.         -}
  55. {-                                                                          -}
  56. {----------------------------------------------------------------------------}
  57.  
  58. FUNCTION Arg (v:polar) : REAL;
  59.  
  60.    CONST
  61.       epsilon = 1e-9;
  62.  
  63.    VAR
  64.       shift : REAL;
  65.  
  66.    BEGIN   {Arg}
  67.       IF v.Mag>=0
  68.          THEN shift:=0
  69.       ELSE IF v.Ph>=0
  70.          THEN shift:=PI
  71.          ELSE shift:=-PI;
  72.       IF (abs(v.Mag) <= epsilon)
  73.          THEN IF v.Ph>=0
  74.             THEN Arg:=+PI/2
  75.             ELSE Arg:=-PI/2
  76.          ELSE Arg:=arctan(v.Ph/v.Mag)+shift;
  77.    END;   {Arg}
  78.  
  79.  
  80. {----------------------------------------------------------------------------}
  81. {-                                                                          -}
  82. {-    FFT and FHT results are given as real and imaginary parts. Polarform  -}
  83. {-    converts those values to magnitude and phase (radians) components.    -}
  84. {----------------------------------------------------------------------------}
  85.  
  86. PROCEDURE PolarForm;
  87.  
  88.    CONST
  89.       epsilon = 1e-06;      (* A tolerance. Any angle <= epsilon2 is 0. *)
  90.  
  91.    VAR
  92.       i     : INTEGER;      (* A counter variable.                      *)
  93.       temp  : polar;        (* A temporary variable.                    *)
  94.       konst : REAL;         (* Magnitude scale factor for transform.    *)
  95.  
  96.  
  97.    BEGIN   {PolarForm}
  98.       konst:=1/(NumPoints*delta);
  99.       FOR i:=0 TO n DO BEGIN
  100.          temp.Mag:=TempXPtr^[i];
  101.          temp.Ph :=TempYPtr^[i];
  102.          TempXPtr^[i]:=sqrt(sqr(temp.Mag)+sqr(temp.Ph))*konst;
  103.          TempYPtr^[i]:=Arg (temp);
  104.       END;   {FOR}
  105.   END;   {PolarForm}
  106.  
  107. {----------------------------------------------------------------------------}
  108. {-  Timer reads the DOS time; it is used to time the FFT and FHT routines.  -}
  109. {----------------------------------------------------------------------------}
  110.  
  111. FUNCTION timer : REAL;
  112.  
  113.    VAR
  114.       hour,min,sec,ms : word;
  115.  
  116.    BEGIN   {timer}
  117.       GetTime (hour,min,sec,ms);
  118.       timer:=hour*3600 + min*60 + sec + ms/100;
  119.    END;   {timer}
  120.  
  121.  
  122. {----------------------------------------------------------------------------}
  123. {-                                                                          -}
  124. {-    Cooley_Tukey_FFT is modified from the FFT algorithm of Numerical      -}
  125. {-    Analysis, 3rd ed. by Burden & Faires, p. 404. The algorithm is        -}
  126. {-    based on the Cooley-Tukey FFT algorithm.                              -}
  127. {-                                                                          -}
  128. {----------------------------------------------------------------------------}
  129.  
  130.  
  131. PROCEDURE Cooley_Tukey_FFT;
  132.  
  133.    TYPE
  134.       binary      = 0..1;
  135.       BinaryArray = ARRAY [0..MaxPower] OF binary;
  136.  
  137.    VAR
  138.       big_M    : INTEGER;
  139.       big_K    : INTEGER;
  140.       q        : INTEGER;
  141.       L        : INTEGER;
  142.       k        : BinaryArray;
  143.       K2       : INTEGER;
  144.       eta      : polar;
  145.       temp     : polar;
  146.       i        : INTEGER;
  147.       j        : INTEGER;
  148.  
  149.    BEGIN   {Cooley_Tukey_FFT}
  150.       OutArrayX^[0]:=1;
  151.       OutArrayY^[0]:=0;
  152.       FOR i:=1 TO NumFreqs DO BEGIN
  153.          OutArrayX^[i]:=cos(i*PI/NumFreqs);
  154.          OutArrayY^[i]:=sin(i*PI/NumFreqs);
  155.          OutArrayX^[i+NumFreqs]:=-OutArrayX^[i];
  156.          OutArrayY^[i+NumFreqs]:=-OutArrayY^[i];
  157.       END;   {FOR}
  158.       big_K:=0;
  159.       big_M:=NumFreqs;
  160.       q:=power_index;
  161.       FOR L:=1 TO succ(power_index) DO BEGIN
  162.          WHILE big_K < n DO BEGIN
  163.             FOR j:=1 TO big_M DO BEGIN
  164.                FOR i:=0 TO power_index DO
  165.                   k[i]:=binary((big_K AND (1 shl i)) shr i);
  166.                K2:=0;
  167.                FOR i:=power_index DOWNTO q DO
  168.                   K2:=K2+(k[i] shl (power_index-i+q));
  169.                eta.Mag:=TempXPtr^[big_K+big_M]*OutArrayX^[K2]
  170.                        -TempYPtr^[big_K+big_M]*OutArrayY^[K2];
  171.                eta.Ph :=TempXPtr^[big_K+big_M]*OutArrayY^[K2]
  172.                        +TempYPtr^[big_K+big_M]*OutArrayX^[K2];
  173.                TempXPtr^[big_K+big_M]:=TempXPtr^[big_K]-eta.Mag;
  174.                TempYPtr^[big_K+big_M]:=TempYPtr^[big_K]-eta.Ph;
  175.                TempXPtr^[big_K]:=TempXPtr^[big_K]+eta.Mag;
  176.                TempYPtr^[big_K]:=TempYPtr^[big_K]+eta.Ph;
  177.                big_K:=succ(big_K);
  178.             END;   {FOR}
  179.             big_K:=big_K+big_M;
  180.          END;   {WHILE}
  181.          big_K:=0;
  182.          big_M:=big_M div 2;
  183.          q:=pred(q);
  184.       END;   {FOR}
  185.       WHILE big_K < n DO BEGIN
  186.          FOR i:=0 TO power_index DO
  187.             k[i]:=binary((big_K AND (1 shl i)) shr i);
  188.          j:=0;
  189.          FOR i:=0 TO power_index DO
  190.             j:=j+(k[i] shl (power_index-i));
  191.          IF j>big_K THEN BEGIN
  192.             temp.Mag:=TempXPtr^[j];
  193.             temp.Ph :=TempYPtr^[j];
  194.             TempXPtr^[j]:=TempXPtr^[big_K];
  195.             TempYPtr^[j]:=TempYPtr^[big_k];
  196.             TempXPtr^[big_K]:=temp.Mag;
  197.             TempYPtr^[big_K]:=temp.Ph;
  198.          END;   {IF}
  199.          big_K:=succ(big_K);
  200.       END;   {WHILE}
  201.    END;   {Cooley_Tukey_FFT}
  202.  
  203.  
  204. {----------------------------------------------------------------------------}
  205. {-                                                                          -}
  206. {-      Fast Hartley Transform (FHT) routine, version 1.00,                 -}
  207. {-        dated 22 November 1986. Appeared in BYTE magazine, April          -}
  208. {-        1988. Modified for use in program FFT on 13 April 1988.           -}
  209. {-                                                                          -}
  210. {----------------------------------------------------------------------------}
  211.  
  212. {----------------------------------------------------------------------------}
  213. {-  Start of the Hartley butterfly transform.                               -}
  214. {----------------------------------------------------------------------------}
  215. {----------------------------------------------------------------------------}
  216. {-  FOR k = 1 TO power DO                                                   -}
  217. {-     1. Calculate the address of the retrograde index for the sine term   -}
  218. {-        for the dual place algorithm, if it is required.                  -}
  219. {-     2. Butterfly transform an index pair.                                -}
  220. {-     3. INCrement trg_ind, j.                                             -}
  221. {-  END - FOR                                                               -}
  222. {----------------------------------------------------------------------------}
  223. {----------------------------------------------------------------------------}
  224. {-  End of Hartley butterfly. The results are scaled if necessary, and      -}
  225. {-  then placed back into the data array.                                   -}
  226. {----------------------------------------------------------------------------}
  227.  
  228.  
  229. PROCEDURE fht;
  230.  
  231.    TYPE
  232.       TNDoubleArray    = ARRAY [1..TNArraySize,1..2] OF REAL;
  233.       TNDoubleArrayPtr = ^TNDoubleArray;
  234.  
  235.    VAR
  236.       i        : INTEGER;
  237.       j        : INTEGER;
  238.       k        : INTEGER;
  239.       trg_ind  : INTEGER;
  240.       trg_inc  : INTEGER;
  241.       power    : INTEGER;
  242.       t_a      : INTEGER;
  243.       f_a      : INTEGER;
  244.       i_temp   : INTEGER;
  245.       section  : INTEGER;
  246.       s_start  : INTEGER;
  247.       s_end    : INTEGER;
  248.       cosine   : REAL;
  249.       sine     : REAL;
  250.       accu     : TNDoubleArrayPtr;
  251.       modify   : INTEGER;
  252.  
  253. {----------------------------------------------------------------------------}
  254. {-  Permutation routine. This routine reorders the data before the          -}
  255. {-  butterfly transform routine is called.                                  -}
  256. {----------------------------------------------------------------------------}
  257.  
  258.    FUNCTION permute (index : INTEGER) : INTEGER;
  259.  
  260.       VAR
  261.          i : INTEGER;
  262.          j : INTEGER;
  263.  
  264.       BEGIN   {permute}
  265.          j:=0;
  266.          DEC (index,1);
  267.          FOR i:=1 TO power_index DO BEGIN
  268.             IF Odd (index)
  269.                THEN INC (j,SUCC(j))
  270.                ELSE INC (j,j);
  271.             index:=index div 2;
  272.          END;   {FOR}
  273.          permute:=SUCC(j);
  274.       END;   {permute}
  275.  
  276. {----------------------------------------------------------------------------}
  277. {-  Main program for the fast Hartley transform.                            -}
  278. {----------------------------------------------------------------------------}
  279.  
  280.    BEGIN   {fht}
  281.       DISPOSE (OutArrayX);
  282.       DISPOSE (OutArrayY);
  283.       NEW (accu);
  284.       FillChar (accu^,SizeOf(accu^),0);
  285.       power:=1;
  286.       f_a:=1;
  287.       t_a:=2;
  288.       FOR i:=1 TO NumPoints DO
  289.          accu^[permute(i),f_a]:=TempXPtr^[i-1];
  290.       FOR i:=1 TO power_index DO BEGIN
  291.          j:=1;
  292.          section:=1;
  293.          trg_inc:=NumPoints div (power+power);
  294.          REPEAT
  295.             trg_ind:=1;
  296.             s_start:=SUCC (section*power);
  297.             s_end:=SUCC (section)*power;
  298.             FOR k:=1 TO power DO BEGIN
  299.                IF (s_start = (j+power)) OR (power < 3)
  300.                   THEN modify:=j+power
  301.                   ELSE modify:=s_start+s_end-(j+power)+1;
  302.                cosine:=cos(2*PI*PRED(trg_ind)/NumPoints);
  303.                sine  :=sin(2*PI*PRED(trg_ind)/NumPoints);
  304.                accu^[j,t_a]:=accu^[j,f_a] +
  305.                     accu^[(j+power),f_a]*cosine + accu^[modify,f_a]*sine;
  306.                accu^[(j+power),t_a]:=accu^[j,f_a] -
  307.                     accu^[(j+power),f_a]*cosine - accu^[modify,f_a]*sine;
  308.                INC (trg_ind,trg_inc);
  309.                INC (j,1);
  310.             END;   {FOR}
  311.             INC (j,power);
  312.             INC (section,2);
  313.          UNTIL j > NumPoints;
  314.          INC (power,power);
  315.          i_temp:=t_a;
  316.          t_a:=f_a;
  317.          f_a:=i_temp;
  318.       END;   {FOR}
  319.       FOR i:=1 TO NumPoints DO TempXPtr^[i-1]:=accu^[i,f_a];
  320.       DISPOSE (accu);
  321.       NEW (OutArrayX);
  322.       NEW (OutArrayY);
  323.    END;   {fht}
  324.  
  325.  
  326. {----------------------------------------------------------------------------}
  327. {-  Compute the Fourier transform from the FHT algorithm.                   -}
  328. {----------------------------------------------------------------------------}
  329.  
  330. PROCEDURE Get_FFT;
  331.  
  332.    VAR
  333.       i     : INTEGER;
  334.  
  335.    BEGIN   {Get_FFT}
  336.       i:=1;
  337.       WHILE i < NumPoints DO BEGIN
  338.          TempXPtr^[i]:=TempXPtr^[i]+TempXPtr^[NumPoints-PRED (i)];
  339.          TempYPtr^[i]:=TempXPtr^[i]-TempXPtr^[NumPoints-PRED (i)];
  340.          INC (i,2);
  341.       END;   {FOR}
  342.    END;   {Get_FFT}
  343.  
  344.  
  345. {----------------------------------------------------------------------------}
  346. {-  Compute the inverse Fourier transform from the FHT algorithm.           -}
  347. {----------------------------------------------------------------------------}
  348.  
  349. PROCEDURE Get_IFT;
  350.  
  351.    VAR
  352.       i     : INTEGER;
  353.  
  354.    BEGIN   {Get_IFT}
  355.       FOR i:=1 TO NumFreqs DO BEGIN
  356.          TempXPtr^[i-1]        :=0.5*(TempXPtr^[i]+TempYPtr^[i]);
  357.          TempXPtr^[NumPoints-i]:=0.5*(TempXPtr^[i]-TempYPtr^[i]);
  358.       END;   {FOR}
  359.    END;   {Get_IFT}
  360.  
  361.  
  362. {----------------------------------------------------------------------------}
  363.  
  364. PROCEDURE Transform ( kind      : BOOLEAN;
  365.                       direction : BOOLEAN);
  366.  
  367.    VAR
  368.       beta        : REAL;                   (* beta = t,max - t,min         *)
  369.       i           : INTEGER;                (* counter variable             *)
  370.       j           : INTEGER;                (* counter variable             *)
  371.       ans         : char;                   (* Apply window?                *)
  372.       StartTime   : REAL;                   (* Time at beginning of FFT     *)
  373.       EndTime     : REAL;                   (* Time when finished FFT       *)
  374.  
  375.  
  376.    BEGIN   {Transform}
  377.       power_index:=round(ln(NumFreqs)/ln(2));
  378.       n:=PRED (NumPoints);
  379.  
  380.       {---  Assign data to temporary variables, apply window if desired.  ---}
  381.       IF direction = Fourward
  382.          THEN BEGIN
  383.             WriteXY ('Apply window [H/R/N]? ',StartColumn,24);
  384.             REPEAT
  385.                ans:=ReadKey;
  386.             UNTIL UpCase(ans) IN ['N',ENTER,'H','R'];
  387.             IF ans = ENTER
  388.                THEN Write ('N')
  389.                ELSE Write (ans);
  390.             FillChar (TempYPtr^,SizeOf(TempYPtr^),0);
  391.             TempXPtr^:=ampl^;
  392.             CASE UpCase(ans) OF
  393.                'H': FOR i:=0 TO n DO
  394.                        TempXPtr^[i]:=ampl^[i]*0.5*(1-cos(PI*i/NumFreqs));
  395.                'R': FOR i:=0 TO n DO
  396.                        TempXPtr^[i]:=ampl^[i]*0.5*(1-cos(PI*(1+i/NumPoints)));
  397.             END;   {CASE}
  398.             delta:=1/(time^[n]-time^[0]);
  399.          END   {THEN}
  400.          ELSE BEGIN
  401.             FOR i:=0 TO n DO BEGIN
  402.                TempXPtr^[i]:=mag^[i]*cos(phase^[i]);
  403.                TempYPtr^[i]:=mag^[i]*sin(phase^[i]);
  404.             END;   {FOR}
  405.             beta:=(freq^[n]-freq^[0])/n;
  406.             delta:=1/(NumPoints*beta);
  407.          END;   {ELSE}
  408.          WriteXY ('Transform in progress ...',StartColumn,24); ClrEOL;
  409.  
  410.       {---  Perform transform on data.  ---}
  411.       StartTime:=timer;
  412.       Cooley_Tukey_FFT;
  413.       EndTime:=timer;
  414.       GotoXY (StartColumn,22);
  415.       Write ('Transform performed in ',(EndTime-StartTime):5:1,' seconds.');
  416.  
  417.       {---  Prepare output arrays.  ---}
  418.       IF direction = Fourward
  419.          THEN BEGIN
  420.             PolarForm;
  421.             FOR i:=0 TO n DO
  422.                freq^[i] :=i*delta;
  423.             mag^  :=TempXPtr^;
  424.             phase^:=TempYPtr^;
  425.          END   {THEN}
  426.          ELSE BEGIN
  427.             time^[0]:=0;
  428.             ampl^[0]:=TempXPtr^[0]*beta;
  429.             FOR i:=1 TO n DO BEGIN
  430.                time^[i]:=i*delta;
  431.                ampl^[i]:=TempXPtr^[n-i+1]*beta;
  432.             END;   {FOR}
  433.          END;   {ELSE}
  434.    END;   {Transform}
  435.  
  436. (****************************************************************************)
  437.  
  438. BEGIN   {Initialization}
  439. END.   {UNIT Transforms}