home *** CD-ROM | disk | FTP | other *** search
/ HAM Radio 1 / HamRadio.cdr / tech / design2 / chebfil.pas < prev    next >
Encoding:
Pascal/Delphi Source File  |  1986-07-27  |  10.3 KB  |  390 lines
  1. program fildisp;
  2.  
  3. {               Filter Response & Stores Data to Disk        }
  4.  
  5. type     XYData = record
  6.                         XX,YY: real;
  7.                   end;
  8.          vector = array[1..151] of real;
  9.          vector1= array[0..20] of real;
  10.          vector20 = array[0..20] of real;
  11.          Plotarray = array[1..100,1..2] of real;
  12.  
  13. var      A              : Plotarray;
  14.          graph_title    : string[15];
  15.          N              : integer;
  16.          Fr,mag         : vector;
  17.          wp_lo,wp_hi    : real;
  18.          wo_geom        : real;
  19.          w_beg,w_end    : real;
  20.          order          : integer;
  21.          eta            : real;
  22.          w_inc          : real;
  23.          km             : integer;
  24.          w              : real;
  25.          i,j,k,m        : integer;
  26.          wp             : real;
  27.          p,p_2          : real;
  28.          B,r            : real;
  29.          ak,bk,g        : vector1;
  30.          gk_sum         : real;
  31.          L_dB           : real;
  32.          Qu             : real;
  33.          T_do           : real;
  34.          Td_sum         : real;
  35.          Peak           : real;
  36.          Loss,Loss_wp   : real;
  37.          phi            : real;
  38.          Anz,Ans1,Ans2  : char;
  39.          T_Rhs,T_Lhs    : real;
  40.          Value          : string[20];
  41.          Const_2,Const_1: real;
  42.          Pole_Re        : vector20;
  43.          Pole_Im        : vector20;
  44.          Time           : vector;
  45.          t1,t2,t3       : real;
  46.          Odd_Order      : Boolean;
  47.          Ang,Ang_1      : real;
  48.          xmax,xmin,ymax,ymin: real;
  49.          Title              : string[15];
  50.          file_name          : file;
  51.  
  52. label 500;
  53. {                              Math Library Functions                       }
  54. {                                   Version  2.0                            }
  55. {                                                                           }
  56. Function Tan( X: Real): real;
  57. begin
  58.     Tan:= Sin(X)/Cos(X);
  59. end;
  60. { ///////////////////////////}
  61.  
  62. Function Sign ( X:Real ): real;
  63.     var y: real;
  64.     begin
  65.          if X < 0 then Sign:= -1.0 else Sign:= 1.0;
  66.     end;
  67. { ///////////////////////////}
  68.  
  69. Function Atan4( Y,X:Real ): Real ;  {  4 Quadrant ArcTan  }
  70.     const Pi = 3.141592654;
  71.     var T1 : real;
  72.     begin
  73.          if abs(X) > 1.0E-8 * ABS(Y) then
  74.               begin
  75.                    T1:= ArcTan( Y/X );
  76.                    if X < 0 then T1:= T1 + Pi;
  77.               end
  78.               else T1:= 0.5*Pi*Sign ( Y );
  79.          Atan4:= T1;
  80.     end;
  81. { ///////////////////////////}
  82.  
  83. Function Log10( X: real): real;
  84.     begin
  85.          Log10:= 0.434294482*Ln( X );
  86.     end;
  87. { /////////////////////////// }
  88.  
  89. Function Cosh ( X: real ): real;
  90.     var temp: real;
  91.     begin
  92.          temp:= exp( X );
  93.          Cosh:= 0.5* (temp + 1/temp );
  94.     end;
  95. { /////////////////////////// }
  96.  
  97. Function Sinh ( X:real ): real;
  98.     var temp : real;
  99.     begin
  100.          temp:= exp( X );
  101.          Sinh:= 0.5* (temp - 1/temp);
  102.     end;
  103. { /////////////////////////// }
  104.  
  105. Function ArcSinh( X: real) : real;
  106.     begin
  107.          ArcSinh:= Ln ( X + sqrt( sqr(X) + 1 ));
  108.     end;
  109. { /////////////////////////// }
  110.  
  111. Function XtoY ( X,Y: real): real;
  112.     begin
  113.          XtoY:= exp( Y*Ln(X));
  114.     end;
  115. { /////////////////////////// }
  116.  
  117. Function Tanh( X: real): real;
  118.     var temp,temp2: real;
  119.     begin
  120.          temp:= exp( X );
  121.          temp2:= 1/temp;
  122.          Tanh:= ( temp -temp2 )/( temp + temp2 );
  123.     end;
  124. { /////////////////////////// }
  125.  
  126. Function ArcCosh( X: real): real;
  127.     begin
  128.          ArcCosh:= Ln ( X + sqrt( sqr(X) - 1));
  129.     end;
  130. { /////////////////////////// }
  131.  
  132. Function ArcCos( X: real): real;
  133.     begin
  134.          ArcCos:= arctan( sqrt( 1 - sqr(X))/X );
  135.     end;
  136. { /////////////////////////// }
  137.  
  138. Procedure Write_to_File( X,Y: vector);
  139.  
  140. var j              : integer;
  141.     write_file     : file of XYData;
  142.     read_file      : file of XYData;
  143.     file_rec       : XYData;
  144.     file_name      : string[10];
  145.     file_size      : integer;
  146.     OK             : Boolean;
  147.     xyz            : boolean;
  148.  
  149. label 1;
  150.  
  151. begin
  152. clrscr;
  153. gotoxy(10,10);
  154. write('Store Inputed Data to Disk  y / n  ?  ');
  155. readln(Ans1);
  156. if Ans1 = 'Y' then
  157.     begin
  158.     clrscr;
  159.     OK:= False;
  160.     gotoxy(10,10);
  161.     write('Write Data to What File  ?  ');
  162.     readln(file_name);
  163.     Assign(write_file,file_name);
  164.     {$I-} Reset( write_file)  {$I+};
  165.     OK:= (IOResult = 0);   { Determine if Write-to File Already Exists }
  166.     if OK then
  167.          begin
  168.          clrscr;
  169.          gotoxy(10,10);
  170.          write('File Already Exists . . .  Overwrite  y / n  ? ');
  171.          readln(Ans2);
  172.          if Ans2 <> 'y' then goto 1;
  173.          end;
  174.     Assign(write_file,file_name);
  175.     Rewrite(write_file);
  176.     with file_rec do
  177.          begin
  178.          for j:= 1 to km do   { Writes Data to Disk File }
  179.               begin
  180.               XX:= X[j];
  181.               YY:= Y[j];
  182.               write(write_file,file_rec);
  183.               end;
  184.          end;
  185.     close(write_file);
  186.  end;
  187.  1: delay(10);
  188.  end;
  189.  
  190.  
  191. begin
  192. 500: clrscr;
  193. gotoxy(1,5);
  194. write('     Enter    Low,    Hi    Ripple Freqs,   MHz  ');
  195. readln(wp_lo,wp_hi);
  196. write('     Enter    Start   Stop   Freqs    MHz  ');
  197. readln(w_beg,w_end);
  198. write('     Enter   Freq   Increment    MHz  ');
  199. readln(w_inc);
  200. write('     Enter   dB   Passband   Ripple  ');
  201. readln(r);
  202. write('     Enter   Order  of  Filter;    Unloaded Q   ');
  203. readln(order,Qu);
  204. eta:= sqrt(XtoY(10.0,r/10.0) - 1.0);
  205. B:= Ln(1/tanh(r/17.37) );
  206. p:= sinh(B/(2*order));
  207. p_2:= sqr(p);
  208. for k:= 1 to order do
  209.     begin
  210.     ak[k]:= sin((2*k-1)*Pi/(2*order));
  211.     bk[k]:= p_2 + sqr(sin(k*Pi/order));
  212.     end;
  213. g[1]:= 2*ak[1]/p;
  214. for k:= 2 to order do
  215.     begin
  216.     g[k]:= 4*ak[k-1]*ak[k]/(bk[k-1]*g[k-1]);
  217.     end;
  218. if ( order mod 2 = 1 ) then
  219.     g[order+1]:= 1.0 else g[order+1]:= sqr(1.0/tanh(0.25*B));
  220.     gk_sum:= 0.0;
  221. for k:= 1 to order do
  222.     begin
  223.     gk_sum:= gk_sum + g[k];
  224.     end;
  225. L_dB:= 4.343*gk_sum*sqrt(wp_lo*wp_hi)/( (wp_hi-wp_lo)*Qu );
  226. T_do:= gk_sum/( 2*Pi*0.001*(wp_hi-wp_lo) );
  227. phi:= 1/order*Arcsinh(1/eta);
  228. if order mod 2 = 0 then Odd_order:= False else Odd_order:= True;
  229. Const_2:= phi;
  230. Const_1:= sinh(Const_2);
  231. Const_2:= cosh(Const_2);
  232.  
  233. Peak:= 1.0;
  234. if order mod 2 = 0 then
  235.     begin
  236.          for i:= 1 to order div 2 do
  237.               begin
  238.                    Pole_Re[i]:= -Const_1*sin(Pi*(2*i-1)/(2*order) );
  239.                    Pole_Im[i]:=  Const_2*Cos(Pi*(2*i-1)/(2*order) );
  240.                    Peak:= Peak * sqr( sqr(Pole_Re[i] ) + sqr(Pole_Im[i] ) );
  241.               end;
  242.          end
  243.     else
  244.          begin
  245.               for i:= 1 to (order-1) div 2 do
  246.                    begin
  247.                    Pole_Re[i]:= -Const_1*sin(Pi*(2*i-1)/(2*order) );
  248.                    Pole_Im[i]:=  Const_2*Cos(Pi*(2*i-1)/(2*order) );
  249.                    Peak:= Peak * sqr( sqr(Pole_Re[i] ) + sqr(Pole_Im[i] ) );
  250.               end;
  251.          Pole_Re[i+1]:= -Const_1;
  252.          Pole_Im[i+1]:= 0.0;
  253.          Peak:= Peak * sqr( Const_1 );
  254.     end;
  255.  
  256. w:= w_beg;
  257. m:= 1;
  258. wo_geom:= sqrt( wp_lo*wp_hi);
  259. T_Lhs:= -1.0;
  260. T_Rhs:= -1.0;
  261. while w <= w_end do
  262.     begin
  263.     Time[m]:= 0.0;
  264.     wp:= wo_geom/(wp_hi-wp_lo);
  265.     wp:= wp*(w/wo_geom - wo_geom/w );
  266.     Fr[m]:= w;
  267.     if not Odd_order then
  268.          begin
  269.          Ang:= 0.0;
  270.          for i:= 1 to order div 2 do
  271.               begin
  272.               t1:= wp - Pole_Im[i];
  273.               t2:= wp + Pole_im[i];
  274.               t3:= Pole_Re[i];
  275.               Ang_1:= -Atan4(t1,t3);
  276.               Time[m]:= Time[m] - sqr(Cos(Ang_1) ) /t3 ;
  277.               Ang:= Ang + Ang_1;
  278.               Ang_1:= -Atan4( t2,t3);
  279.               Time[m]:= Time[m] - sqr(Cos(Ang_1) )/t3;
  280.               Ang:= Ang + Ang_1;
  281.               end;
  282.          end;
  283.     if Odd_order or (order = 1) then
  284.     begin
  285.          Ang:= 0.0;
  286.          for i:= 1 to (order-1) div 2 do
  287.               begin
  288.               t1:= wp - Pole_Im[i];
  289.               t2:= wp + Pole_im[i];
  290.               t3:= Pole_Re[i];
  291.               Ang_1:= -Atan4(t1,t3);
  292.               Time[m]:= Time[m] - sqr(Cos(Ang_1) ) /t3 ;
  293.               Ang:= Ang + Ang_1;
  294.               Ang_1:= -Atan4( t2,t3);
  295.               Time[m]:= Time[m] - sqr(Cos(Ang_1) )/t3;
  296.               Ang:= Ang + Ang_1;
  297.               end;
  298.          t3:= Pole_Re[(order+1) div 2];
  299.          Ang_1:= -Atan4(wp,t3);
  300.          Time[m]:= Time[m] - sqr( Cos( Ang_1) ) /t3;
  301.          Ang:= Ang + Ang_1;
  302.          end;
  303. Time[m]:= Time[m] * 1000 / ( Pi * (wp_hi-wp_lo) );
  304. if w < wo_geom then
  305.     begin
  306.     if Time[m] > T_Lhs then T_Lhs:= Time[m];
  307.     end;
  308. if w > wo_geom then
  309.     begin
  310.     if Time[m] > T_Rhs then T_Rhs:= Time[m];
  311.     end;
  312. w:= w + w_inc;
  313. m:= m + 1;
  314. end;
  315. w:= w_beg;
  316. m:= 1;
  317. clrscr;
  318. while w < w_end do
  319.     begin
  320.     wp:= wo_geom/(wp_hi - wp_lo);
  321.     wp:= abs( wp*(w/wo_geom - wo_geom/w ) );
  322.     wp:= wp*1.0E10;
  323.     wp:= Int(wp)*1.0E-10;
  324.     if ( (wp > 1.0) and ( w < wp_lo) ) then
  325.          begin
  326.          mag[m]:= -10*Log10(1 + sqr(eta*Cosh(order*Arccosh(wp)))) - L_db*T_Lhs/T_do;
  327.          end;
  328.     if ( (wp > 1.0) and ( w > wp_hi) ) then
  329.          begin
  330.          mag[m]:= -10*Log10(1 + sqr(eta*Cosh(order*Arccosh(wp)))) - L_db*T_Rhs/T_do;
  331.          end;
  332.     if  wp <= 1.0 then
  333.          begin
  334.          mag[m]:= -10*Log10(1 + sqr(eta*Cos(order*Arccos(wp)))) - L_db*Time[m]/T_do;
  335.          end;
  336. if m <= 23 then
  337.     begin
  338.     gotoxy(1,m+1);
  339.     write(Fr[m]:9:4,'  ',mag[m]:7:2);
  340.     end;
  341. if ( ( m > 23 ) and ( m <= 46) ) then
  342.     begin
  343.     gotoxy(20,m-22);
  344.     write('  ',Fr[m]:9:4,'  ',mag[m]:7:2);
  345.     end;
  346. if ( ( m > 46 ) and ( m <= 69) ) then
  347.     begin
  348.     gotoxy(40,m-45);
  349.     write('  ',Fr[m]:9:4,'  ',mag[m]:7:2);
  350.     end;
  351. if ( ( m > 69 ) and ( m <= 92) ) then
  352.     begin
  353.     gotoxy(60,m-68);
  354.     write('  ',Fr[m]:9:4,' ',mag[m]:7:2);
  355.     end;
  356. w:= w + w_inc;
  357. m:= m + 1;
  358. end;
  359. repeat until keypressed;
  360. clrscr;
  361. gotoxy(10,5);
  362. write('Enter Title for Graph  ');
  363. readln(Graph_Title);
  364. for i:= 1 to m-1 do
  365.     begin
  366.     A[i,1]:= Fr[i];
  367.     A[i,2]:= mag[i];
  368.     end;
  369. N:= m - 1;
  370. assign(file_name,'PlotGr.Chn');
  371. Chain(file_name);
  372. readln(N);
  373. clrscr;
  374. gotoxy(10,5);
  375. write('Enter Title for Graph for Group Delay ');
  376. readln(Graph_Title);
  377. for i:= 1 to m-1 do
  378.     begin
  379.     A[i,1]:= Fr[i];
  380.     A[i,2]:= Time[i];
  381.     end;
  382. N:= m - 1;
  383. assign(file_name,'PlotGr.Chn');
  384. Chain(file_name);
  385. Write_to_File( Fr,mag);
  386. goto 500;
  387. end.
  388.  
  389.  
  390.