home *** CD-ROM | disk | FTP | other *** search
/ Programmer 7500 / MAX_PROGRAMMERS.iso / PASCAL / NUMPAS.ZIP / FREQRES.PAS < prev    next >
Encoding:
Pascal/Delphi Source File  |  1985-06-25  |  3.8 KB  |  141 lines
  1. program freqres(input,output);
  2. (* Program evaluates magnitude in decibels and angle in degrees of a *)
  3. (* rational function  K*N(s)/D(s) for s = jw over a desired range of w. *)
  4.  
  5. type vector = array[0..10] of real;
  6. var coefn,coefd: vector;
  7.     gain,freqlo,freqhi,freq,rmag,jmag,
  8.     rmagn,jmagn,rmagd,jmagd,cmag,angle,spacing: real;
  9.     ndegree,ddegree,npts: integer;
  10.  
  11. procedure readin(var a1: vector; n1: integer);
  12. (* Reads in polynomial coefficients from the keyboard. *)
  13.  
  14. var i: integer;
  15.  
  16. begin
  17.    for i := 0 to n1 do
  18.       begin
  19.          write('Input the ',i,' coefficent: ');
  20.          readln(a1[i]);
  21.       end;
  22.    writeln('The polynomial coefficients are:');
  23.    writeln;
  24.    for i := 0 to n1 do
  25.       begin
  26.          write(a1[i]:6:3,'  ');
  27.       end;
  28.    writeln;
  29. end;
  30.  
  31. procedure decade(npts1: integer; var delta1: real);
  32. (* Sets frequency increment so as to get uniform spacing of points on *)
  33. (* a log(frequency) scale for  npts1 = 2^i points per decade. *)
  34.  
  35. var i,j,n1: integer;
  36.  
  37. begin
  38.    i := 0;
  39.    n1 := npts1;
  40.    while n1 >= 2 do
  41.       begin
  42.          n1 := n1 div 2;
  43.          i := i + 1;
  44.       end;
  45.    delta1 := 10.0;
  46.    for j := 1 to i do
  47.       begin
  48.          delta1 := sqrt(delta1);
  49.       end;
  50. end;
  51.  
  52. procedure evalpoly(coef: vector; w1,w2: real; var b1,b2: real; n1: integer);
  53. (* Evaluates a polynomial N(s) in the form *)
  54. (* A(0) + A(1)s + A(2)s^2 +...+A(n1)s^n1  for s = w1 + jw2.  In this *)
  55. (* program w1 = 0. *)
  56.  
  57. var a1,a2: vector;
  58.     i: integer;
  59.  
  60. begin
  61.    a1[0] := 1.0;
  62.    a2[0] := 0.0;
  63.    b1 := a1[0]*coef[0];
  64.    b2 := a2[0]*coef[0];
  65.    for i := 1 to n1 do
  66.       begin
  67.          a1[i] := w1*a1[i-1] - w2*a2[i-1];  (* real part of w^i *)
  68.          a2[i] := w1*a2[i-1] + w2*a1[i-1];  (* imag part of w^i *)
  69.          b1 := b1 + a1[i]*coef[i];   (* real part of trunc polynomial *)
  70.          b2 := b2 + a2[i]*coef[i];   (* imag part of trunc polynomial *)
  71.       end;
  72. end;
  73.  
  74. procedure evaltf(a1,a2,b1,b2: real; var d1,d2: real);
  75. (* Evaluates a complex quotient  d1 + jd2 = (a1 + ja2)/(b1 + jb2). *)
  76.  
  77. var d: real;
  78.  
  79. begin
  80.    d := b1*b1 + b2*b2;
  81.    d1 := (a1*b1 + a2*b2)/d;
  82.    d2 := (a2*b1 - a1*b2)/d;
  83. end;
  84.  
  85. procedure rectpol(a1,a2: real; var cmag1,angle1: real);
  86. (* Converts rectangular form a1 + ja2  to polar form. *)
  87.  
  88. var eta: real;
  89.  
  90. begin
  91.    cmag1 := sqrt(a1*a1 + a2*a2);
  92.    if (a1 = 0.0) and (a2 = 0.0) then
  93.       angle1 := 0.0;
  94.    if (a1 = 0.0) and (a2 > 0.0) then
  95.       angle1 := 90.0;
  96.    if (a1 = 0.0) and (a2 < 0.0) then
  97.       angle1 := -90.0;
  98.    if abs(a1) > 0.0 then
  99.       eta := arctan(a2/a1)*57.296;
  100.    if a1 > 0 then
  101.       angle1 := eta;
  102.    if (a1 < 0.0) and (a2 >= 0.0) then
  103.       angle1 := eta + 180.0;
  104.    if (a1 < 0.0) and (a2 < 0.0) then
  105.       angle1 := eta - 180.0;
  106. end;
  107.  
  108. begin  (*  main  *)
  109.    write('Input gain constant: ');
  110.    readln(gain);
  111.    writeln;
  112.    write('Input numerator and denominator degree: ');
  113.    readln(ndegree,ddegree);
  114.    writeln;
  115.    write('Input lowest and highest frequencies: ');
  116.    readln(freqlo,freqhi);
  117.    writeln;
  118.    writeln('Input number of points per decade (powers of 2): ');
  119.    readln(npts);
  120.    writeln;
  121.    decade(npts,spacing);
  122.    writeln('Input numerator coefficients in response to prompt:');
  123.    writeln;
  124.    readin(coefn,ndegree);
  125.    writeln('Input denominator coefficients in response to prompt:');
  126.    writeln;
  127.    readin(coefd,ddegree);
  128.    writeln(' rad/s ',' decibels ',' degrees ');
  129.    writeln;
  130.    freq := freqlo;
  131.    while freq < freqhi do
  132.       begin
  133.          evalpoly(coefn,0,freq,rmagn,jmagn,ndegree);
  134.          evalpoly(coefd,0,freq,rmagd,jmagd,ddegree);
  135.          evaltf(rmagn,jmagn,rmagd,jmagd,rmag,jmag);
  136.          rectpol(rmag,jmag,cmag,angle);
  137.          writeln(freq:6:1,'  ',20*ln(gain*cmag)/ln(10):8:3,'  ',angle:5:2);
  138.          freq := freq*spacing;
  139.       end;
  140. end.
  141.