home *** CD-ROM | disk | FTP | other *** search
/ Hacker Chronicles 2 / HACKER2.BIN / 151.SIGNALPR.PAS < prev    next >
Pascal/Delphi Source File  |  1988-06-15  |  16KB  |  373 lines
  1. UNIT SignalProcessing;
  2.  
  3. INTERFACE
  4.  
  5. USES
  6.    DOS,
  7. {$IFDEF DOSCrt}
  8.    DOSCrt,
  9. {$ELSE}
  10.    Crt,
  11. {$ENDIF}
  12.    Extended_Reals,
  13.    TextOps,
  14.    Global,
  15.    Transforms;
  16.  
  17. {----------------------------------------------------------------------------}
  18. {-                                                                          -}
  19. {-    SPS prints the SPS menu, and gives the user the choice of one of      -}
  20. {-    several signal processing functins:                                   -}
  21. {-                                                                          -}
  22. {-       1. Integration of time domain data.                                -}
  23. {-       2. Differentiation of time domain data.                            -}
  24. {-       3. Fast Fourier Transform of time domain data.                     -}
  25. {-       4. Adjust (+,-,*,/) time domain amplitude by a constant.           -}
  26. {-       5. Multiply time base by a factor.                                 -}
  27. {-       6. Subtract the mean of the time domain waveform.                  -}
  28. {-                                                                          -}
  29. {-    SPS is called AFTER the data has been fitted, a procedure invoked     -}
  30. {-    the first time one chooses selection 3 from the Main Menu.            -}
  31. {-                                                                          -}
  32. {----------------------------------------------------------------------------}
  33.  
  34. PROCEDURE SPS;
  35.  
  36.  
  37. (****************************************************************************)
  38.  
  39. IMPLEMENTATION
  40.  
  41. {----------------------------------------------------------------------------}
  42. {-                                                                          -}
  43. {-    Fitting accepts for input an array of NumPoints arbitrarily spaced    -}
  44. {-    data points, and returns as output an array of 2^n evenly spaced      -}
  45. {-    data points.  The method of interpolation used is either linear       -}
  46. {-    interpolation or natural cubic spline interpolation.                  -}
  47. {-                                                                          -}
  48. {----------------------------------------------------------------------------}
  49.  
  50.  
  51. PROCEDURE Fitting;
  52.  
  53.    VAR
  54.       NewNumPoints : INTEGER;     (* Number of points to interpolate to.    *)
  55.       i            : INTEGER;     (* Counter variable.                      *)
  56.       j            : INTEGER;     (* Counter for NewNumPoints.              *)
  57.       X1           : REAL;        (* First x value.                         *)
  58.       x            : REAL;        (* Frequencies of fitted data.            *)
  59.       increment    : REAL;        (* Distance between successive x-coor.    *)
  60.       Error        : BYTE;        (* Error code for spline routine.         *)
  61.  
  62. {----------------------------------------------------------------------------}
  63. {-                                                                          -}
  64. {-    Linear performs linear interpolation on the time domain data.         -}
  65. {-                                                                          -}
  66. {----------------------------------------------------------------------------}
  67.  
  68.    PROCEDURE Linear;
  69.  
  70.       VAR
  71.          m            : REAL;     (* Slope of interpolation line.           *)
  72.          b            : REAL;     (* Y-intercept of interpolation line.     *)
  73.          denom        : REAL;     (* Difference between succesive points.   *)
  74.          j            : INTEGER;
  75.  
  76.       BEGIN   {Linear}
  77.          i:=0;
  78.          X1:=time^[0];
  79.          TempXPtr^[0]:=time^[0];
  80.          TempYPtr^[0]:=ampl^[0];
  81.          increment:=(time^[NumPoints-1]-time^[0])/(NewNumPoints-1);
  82.          FOR j:=1 TO NewNumPoints-1 DO BEGIN
  83.             x:=X1+j*increment;
  84.             WHILE x > time^[i] DO INC (i,1);
  85.             (*** Calculate slope and y-intercept ***)
  86.             denom:=1/(time^[i]-time^[i-1]);
  87.             m:=(ampl^[i]-ampl^[i-1])*denom;
  88.             b:=(ampl^[i-1]*time^[i]-ampl^[i]*time^[i-1])*denom;
  89.             TempXPtr^[j]:=x;
  90.             TempYPtr^[j]:=m*x+b;
  91.          END;   {FOR}
  92.       END;   {Linear}
  93.  
  94. {----------------------------------------------------------------------------}
  95. {-                                                                          -}
  96. {-    CubicSpline performs a Cubic Spline interpolation on the time         -}
  97. {-    domain data. The algorithm used is a modification of the algorithm    -}
  98. {-    3.4 of Numerical Analysis, 3rd ed. by Burden & Faires, p. 122.        -}
  99. {-                                                                          -}
  100. {-    Error occurs if there is insufficient memory to perform the cubic     -}
  101. {-    spline algorithm.                                                     -}
  102. {----------------------------------------------------------------------------}
  103.  
  104.    PROCEDURE CubicSpline (Error : BYTE);
  105.  
  106.       VAR
  107.          temp         : REAL;        (* A temporary real variable.          *)
  108.          b            : TNVectorPtr; (* First spline coefficient.           *)
  109.          c            : TNVectorPtr; (* Second spline coefficient.          *)
  110.          d            : TNVectorPtr; (* Third spline coefficient.           *)
  111.          h            : TNVectorPtr; (* Dist. between successive data pts.  *)
  112.          i            : INTEGER;
  113.          j            : INTEGER;
  114.          TopOfHeap    : ^INTEGER;
  115.          dummy        : CHAR;
  116.  
  117.       BEGIN   {CubicSpline}
  118.          Error:=0;
  119.          Mark(TopOfHeap);
  120.          IF MaxAvail <= SizeOf(b)
  121.             THEN Error:=1
  122.             ELSE New (b);
  123.          IF MaxAvail <= SizeOf(c)
  124.             THEN Error:=1
  125.             ELSE New (c);
  126.          IF MaxAvail <= SizeOf(d)
  127.             THEN Error:=1
  128.             ELSE New (d);
  129.          IF MaxAvail < SizeOf(h)
  130.             THEN Error:=1
  131.             ELSE New (h);
  132.  
  133.          {--- If insufficient memory, then display error message and exit ---}
  134.          IF Error <> 0 THEN BEGIN
  135.             Release (TopOfHeap);
  136.             PrintErrorMsg ('Insufficient memory for cubic spline',
  137.                             StartColumn,23,FALSE,dummy);
  138.             Exit;
  139.          END;   {IF}
  140.  
  141.          {--- Calculate spline coefficients ---}
  142.          FOR i:=0 TO NumPoints-2 DO
  143.             h^[i]:=time^[i+1]-time^[i];
  144.          b^[0]:=1;
  145.          OutArrayX^[0]:=0;
  146.          OutArrayY^[0]:=0;
  147.          FOR i:=1 TO NumPoints-2 DO BEGIN
  148.             temp:=3*(ampl^[i+1]*h^[i-1]-ampl^[i]*(time^[i+1]-time^[i-1])
  149.                    +ampl^[i-1]*h^[i])/(h^[i-1]*h^[i]);
  150.             b^[i]:=2*(time^[i+1]-time^[i-1])-h^[i-1]*OutArrayX^[i-1];
  151.             OutArrayX^[i]:=h^[i]/b^[i];
  152.             OutArrayY^[i]:=(temp-h^[i-1]*OutArrayY^[i-1])/b^[i];
  153.          END;   {FOR}
  154.          c^[NumPoints-1]:=0;
  155.          FOR i:=NumPoints-2 DOWNTO 0 DO BEGIN
  156.             c^[i]:=OutArrayY^[i]-OutArrayX^[i]*c^[i+1];
  157.             b^[i]:=(ampl^[i+1]-ampl^[i])/h^[i]-h^[i]*(c^[i+1]+2*c^[i])/3;
  158.             d^[i]:=(c^[i+1]-c^[i])/(3*h^[i]);
  159.          END;   {FOR}
  160.  
  161.          {--- Create NewNumPoints evenly spaced points ---}
  162.          X1:=time^[0];
  163.          increment:=(time^[NumPoints-1]-time^[0])/(NewNumPoints-1);
  164.          i:=0;
  165.          TempXPtr^[0]:=time^[0];
  166.          TempYPtr^[0]:=ampl^[0];
  167.          FOR j:=1 TO NewNumPoints-1 DO BEGIN
  168.             x:=X1+j*increment;
  169.             WHILE x > time^[i] DO INC (i,1);
  170.             temp:=x-time^[i-1];
  171.             TempXPtr^[j]:=x;
  172.             TempYPtr^[j]:=ampl^[i-1]+b^[i-1]*temp+c^[i-1]*sqr(temp)
  173.                         +d^[i-1]*sqr(temp)*temp;
  174.          END;   {FOR}
  175.          Release (TopOfHeap);
  176.       END;   {CubicSpline}
  177.  
  178.    BEGIN   {Fitting}
  179.       {---  Calculate number of points to interpolate  ---}
  180.       {---  to. Must be a power of two.                ---}
  181.       i:=MinPower;
  182.       WHILE NumPoints > Pow2[i] DO INC (i,1);
  183.       NumFreqs:=Pow2[i-1];
  184.       NewNumPoints:=Pow2[i];
  185.       IF SPLINE THEN BEGIN
  186.          CubicSpline (Error);
  187.          IF Error <> 0 THEN SPLINE:=FALSE;
  188.       END;   {IF}
  189.       IF NOT SPLINE THEN Linear;
  190.       NumPoints:=NewNumPoints;
  191.       ACCEPT:=true;
  192.       time^:=TempXPtr^;
  193.       ampl^:=TempYPtr^;
  194.    END;   {Fitting}
  195.  
  196. {----------------------------------------------------------------------------}
  197. {-                                                                          -}
  198. {-    Integrate performs numerical integration on the time domain data      -}
  199. {-    via the trapezoidal rule. Note that this procedure requires           -}
  200. {-    approximately 2*m additions and 2*m multiplications.                  -}
  201. {-                                                                          -}
  202. {----------------------------------------------------------------------------}
  203.  
  204. PROCEDURE integrate;
  205.  
  206.    VAR
  207.       i   : INTEGER;         (* Counter.                                    *)
  208.       h   : REAL;            (* Distance between successive points.         *)
  209.  
  210.    BEGIN   {integrate}
  211.       h:=(time^[NumPoints-1]-time^[0])/(NumPoints);
  212.       TempXPtr^[0]:=0;
  213.       FOR i:=1 TO NumPoints-1 DO
  214.          TempXPtr^[i]:=TempXPtr^[i-1]+0.5*h*(ampl^[i-1]+ampl^[i]);
  215.       ampl^:=TempXPtr^;
  216.    END;   {integrate}
  217.  
  218. {----------------------------------------------------------------------------}
  219. {-                                                                          -}
  220. {-    Differentiate performs numerical differentiation via the so-called    -}
  221. {-    three point formulae. The routine requires approximately 2*m          -}
  222. {-    additions and 2*m multiplications.                                    -}
  223. {-                                                                          -}
  224. {----------------------------------------------------------------------------}
  225.  
  226. PROCEDURE differentiate;
  227.  
  228.    VAR
  229.       mm : INTEGER;       (* mm = 2*m-1 = Num points in input array.     *)
  230.       h  : REAL;          (* Inverse of distance between succ. points.   *)
  231.       i  : INTEGER;       (* Counter.                                    *)
  232.  
  233.    BEGIN   {differentiate}
  234.       mm:=NumPoints-1;
  235.       h:=(mm+1)/(time^[mm]-time^[0]);
  236.       TempXPtr^[0]:=0.5*(-3*ampl^[0]+4*ampl^[1]-ampl^[2])*h;
  237.       TempXPtr^[mm]:=0.5*(3*ampl^[mm]-4*ampl^[mm-1]+ampl^[mm-2])*h;
  238.       FOR i:=1 TO mm-1 DO
  239.          TempXPtr^[i]:=0.5*(ampl^[i+1]-ampl^[i-1])*h;
  240.       ampl^:=TempXPtr^;
  241.    END;   {differentiate}
  242.  
  243. {----------------------------------------------------------------------------}
  244. {-                                                                          -}
  245. {-    TransferFunction scales either the x- or y-axis by a constant. The    -}
  246. {-    first character read will decide upon which axis to act, and the      -}
  247. {-    next character read will decide the operation to be performed. Only   -}
  248. {-    the characters '*', '/', '+', and '-' are accepted for operations.    -}
  249. {-                                                                          -}
  250. {----------------------------------------------------------------------------}
  251.  
  252. PROCEDURE TransferFunction;
  253.  
  254.    VAR
  255.       i         : INTEGER;                    (* counter variable           *)
  256.       s         : string;
  257.       transfer  : REAL;                       (* scaling factor             *)
  258.       operation : char;                       (* operation to be performed  *)
  259.       X_AXIS    : BOOLEAN;                    (* apply to x or y axis?      *)
  260.  
  261.    BEGIN   {TransferFunction}
  262.       WriteXY ('Apply transfer function to x or y axis? ',StartColumn,22);
  263.       REPEAT
  264.          operation:=ReadKey;
  265.       UNTIL (UpCase(operation) in ['X','Y']);
  266.       IF UpCase(operation) = 'X'
  267.          THEN X_AXIS:=true
  268.          ELSE X_AXIS:=false;
  269.       WriteXY ('Factor: ',StartColumn,24); ClrEOL;
  270.       REPEAT
  271.          operation:=ReadKey;
  272.       UNTIL (operation in ['*','/','+','-']);
  273.       Write (operation);
  274.       ReadLn (s);
  275.       String_To_Value (s,transfer,s);
  276.       IF X_AXIS
  277.          THEN CASE operation OF
  278.             '*': FOR i:=0 TO NumPoints-1 DO time^[i]:=time^[i]*transfer;
  279.             '/': FOR i:=0 TO NumPoints-1 DO time^[i]:=time^[i]/transfer;
  280.             '+': FOR i:=0 TO NumPoints-1 DO time^[i]:=time^[i]+transfer;
  281.             '-': FOR i:=0 TO NumPoints-1 DO time^[i]:=time^[i]-transfer;
  282.          END   {CASE}
  283.          ELSE CASE operation OF
  284.             '*': FOR i:=0 TO NumPoints-1 DO ampl^[i]:=ampl^[i]*transfer;
  285.             '/': FOR i:=0 TO NumPoints-1 DO ampl^[i]:=ampl^[i]/transfer;
  286.             '+': FOR i:=0 TO NumPoints-1 DO ampl^[i]:=ampl^[i]+transfer;
  287.             '-': FOR i:=0 TO NumPoints-1 DO ampl^[i]:=ampl^[i]-transfer;
  288.          END;   {CASE}
  289.    END;   {TransferFunction}
  290.  
  291. {----------------------------------------------------------------------------}
  292. {-                                                                          -}
  293. {-    RemoveMean calculates the mean of the time domain data and            -}
  294. {-    subtracts that value from the entire array.                           -}
  295. {-                                                                          -}
  296. {----------------------------------------------------------------------------}
  297.  
  298. PROCEDURE RemoveMean;
  299.  
  300.    VAR
  301.       i   : INTEGER;              (* Counter variable.                      *)
  302.       sum : REAL;                 (* The mean of the array.                 *)
  303.  
  304.    BEGIN   {RemoveMean}
  305.       sum:=0;
  306.       FOR i:=0 TO NumPoints-1 DO
  307.          sum:=sum+ampl^[i];
  308.       sum:=sum/(NumPoints);
  309.       FOR i:=0 TO NumPoints-1 DO
  310.          ampl^[i]:=ampl^[i]-sum;
  311.    END;   {RemoveMean}
  312.  
  313. {----------------------------------------------------------------------------}
  314.  
  315. PROCEDURE SPS;
  316.  
  317.    VAR
  318.       Choice : char;
  319.  
  320.    PROCEDURE SPSMenu;
  321.  
  322.       BEGIN   {SPSMenu}
  323.          ClrScr;
  324.          WriteXY ('Signal Processing Menu',StartColumn+5,3);
  325.          WriteXY ('Select Option by typing a number: ',StartColumn,5);
  326.          WriteXY ('1. Integrate          ',StartColumn+8,8);
  327.          WriteXY ('2. Differentiate      ',StartColumn+8,10);
  328.          WriteXY ('3. FFT                ',StartColumn+8,12);
  329.          WriteXY ('4. Inverse FFT        ',StartColumn+8,14);
  330.          WriteXY ('5. Transfer Function  ',StartColumn+8,16);
  331.          WriteXY ('6. Remove DC Level    ',StartColumn+8,18);
  332.          WriteXY ('9. Exit to Main Menu  ',StartColumn+8,20);
  333.          WriteXY ('Your Choice? ',StartColumn,24);
  334.          REPEAT
  335.             Choice:=ReadKey;
  336.          UNTIL Choice IN ['1'..'6','9'];
  337.          Write (Choice);
  338.       END;   {SPSMenu}
  339.  
  340.    BEGIN   {SPS}
  341.       IF NOT ACCEPT THEN BEGIN
  342.          GotoXY (StartColumn,23);
  343.          Write ('Interpolation in progress. Please wait ... ');
  344.          Fitting;
  345.       END;   {IF}
  346.       REPEAT
  347.          SPSMenu;
  348.          CASE Choice OF
  349.             '1': Integrate;
  350.             '2': Differentiate;
  351.             '3': BEGIN   {Transform}
  352.                     Transform (Four,Fourward);
  353.                     TRANS:=TRUE;
  354.                     ORIG:=FALSE;
  355.                  END;   {Fourier Transform}
  356.             '4': IF TRANS THEN BEGIN
  357.                     Transform (Four,Inverse);
  358.                     ORIG:=TRUE;
  359.                  END;   {Inverse Transform}
  360.             '5': TransferFunction;
  361.             '6': RemoveMean;
  362.             '9': BEGIN
  363.                     ClrScr;
  364.                     Exit;
  365.                  END;
  366.          END;   {CASE}
  367.       UNTIL FALSE;
  368.    END;   {SPS}
  369.  
  370. (****************************************************************************)
  371.  
  372. BEGIN   {Initialization}
  373. END.   {UNIT SignalProcessing}