home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
HAM Radio 1
/
HamRadio.cdr
/
tech
/
eepub05
/
chebfil.pas
< prev
next >
Wrap
Pascal/Delphi Source File
|
1986-07-27
|
11KB
|
390 lines
program fildisp;
{ Filter Response & Stores Data to Disk }
type XYData = record
XX,YY: real;
end;
vector = array[1..151] of real;
vector1= array[0..20] of real;
vector20 = array[0..20] of real;
Plotarray = array[1..100,1..2] of real;
var A : Plotarray;
graph_title : string[15];
N : integer;
Fr,mag : vector;
wp_lo,wp_hi : real;
wo_geom : real;
w_beg,w_end : real;
order : integer;
eta : real;
w_inc : real;
km : integer;
w : real;
i,j,k,m : integer;
wp : real;
p,p_2 : real;
B,r : real;
ak,bk,g : vector1;
gk_sum : real;
L_dB : real;
Qu : real;
T_do : real;
Td_sum : real;
Peak : real;
Loss,Loss_wp : real;
phi : real;
Anz,Ans1,Ans2 : char;
T_Rhs,T_Lhs : real;
Value : string[20];
Const_2,Const_1: real;
Pole_Re : vector20;
Pole_Im : vector20;
Time : vector;
t1,t2,t3 : real;
Odd_Order : Boolean;
Ang,Ang_1 : real;
xmax,xmin,ymax,ymin: real;
Title : string[15];
file_name : file;
label 500;
{ Math Library Functions }
{ Version 2.0 }
{ }
Function Tan( X: Real): real;
begin
Tan:= Sin(X)/Cos(X);
end;
{ ///////////////////////////}
Function Sign ( X:Real ): real;
var y: real;
begin
if X < 0 then Sign:= -1.0 else Sign:= 1.0;
end;
{ ///////////////////////////}
Function Atan4( Y,X:Real ): Real ; { 4 Quadrant ArcTan }
const Pi = 3.141592654;
var T1 : real;
begin
if abs(X) > 1.0E-8 * ABS(Y) then
begin
T1:= ArcTan( Y/X );
if X < 0 then T1:= T1 + Pi;
end
else T1:= 0.5*Pi*Sign ( Y );
Atan4:= T1;
end;
{ ///////////////////////////}
Function Log10( X: real): real;
begin
Log10:= 0.434294482*Ln( X );
end;
{ /////////////////////////// }
Function Cosh ( X: real ): real;
var temp: real;
begin
temp:= exp( X );
Cosh:= 0.5* (temp + 1/temp );
end;
{ /////////////////////////// }
Function Sinh ( X:real ): real;
var temp : real;
begin
temp:= exp( X );
Sinh:= 0.5* (temp - 1/temp);
end;
{ /////////////////////////// }
Function ArcSinh( X: real) : real;
begin
ArcSinh:= Ln ( X + sqrt( sqr(X) + 1 ));
end;
{ /////////////////////////// }
Function XtoY ( X,Y: real): real;
begin
XtoY:= exp( Y*Ln(X));
end;
{ /////////////////////////// }
Function Tanh( X: real): real;
var temp,temp2: real;
begin
temp:= exp( X );
temp2:= 1/temp;
Tanh:= ( temp -temp2 )/( temp + temp2 );
end;
{ /////////////////////////// }
Function ArcCosh( X: real): real;
begin
ArcCosh:= Ln ( X + sqrt( sqr(X) - 1));
end;
{ /////////////////////////// }
Function ArcCos( X: real): real;
begin
ArcCos:= arctan( sqrt( 1 - sqr(X))/X );
end;
{ /////////////////////////// }
Procedure Write_to_File( X,Y: vector);
var j : integer;
write_file : file of XYData;
read_file : file of XYData;
file_rec : XYData;
file_name : string[10];
file_size : integer;
OK : Boolean;
xyz : boolean;
label 1;
begin
clrscr;
gotoxy(10,10);
write('Store Inputed Data to Disk y / n ? ');
readln(Ans1);
if Ans1 = 'Y' then
begin
clrscr;
OK:= False;
gotoxy(10,10);
write('Write Data to What File ? ');
readln(file_name);
Assign(write_file,file_name);
{$I-} Reset( write_file) {$I+};
OK:= (IOResult = 0); { Determine if Write-to File Already Exists }
if OK then
begin
clrscr;
gotoxy(10,10);
write('File Already Exists . . . Overwrite y / n ? ');
readln(Ans2);
if Ans2 <> 'y' then goto 1;
end;
Assign(write_file,file_name);
Rewrite(write_file);
with file_rec do
begin
for j:= 1 to km do { Writes Data to Disk File }
begin
XX:= X[j];
YY:= Y[j];
write(write_file,file_rec);
end;
end;
close(write_file);
end;
1: delay(10);
end;
begin
500: clrscr;
gotoxy(1,5);
write(' Enter Low, Hi Ripple Freqs, MHz ');
readln(wp_lo,wp_hi);
write(' Enter Start Stop Freqs MHz ');
readln(w_beg,w_end);
write(' Enter Freq Increment MHz ');
readln(w_inc);
write(' Enter dB Passband Ripple ');
readln(r);
write(' Enter Order of Filter; Unloaded Q ');
readln(order,Qu);
eta:= sqrt(XtoY(10.0,r/10.0) - 1.0);
B:= Ln(1/tanh(r/17.37) );
p:= sinh(B/(2*order));
p_2:= sqr(p);
for k:= 1 to order do
begin
ak[k]:= sin((2*k-1)*Pi/(2*order));
bk[k]:= p_2 + sqr(sin(k*Pi/order));
end;
g[1]:= 2*ak[1]/p;
for k:= 2 to order do
begin
g[k]:= 4*ak[k-1]*ak[k]/(bk[k-1]*g[k-1]);
end;
if ( order mod 2 = 1 ) then
g[order+1]:= 1.0 else g[order+1]:= sqr(1.0/tanh(0.25*B));
gk_sum:= 0.0;
for k:= 1 to order do
begin
gk_sum:= gk_sum + g[k];
end;
L_dB:= 4.343*gk_sum*sqrt(wp_lo*wp_hi)/( (wp_hi-wp_lo)*Qu );
T_do:= gk_sum/( 2*Pi*0.001*(wp_hi-wp_lo) );
phi:= 1/order*Arcsinh(1/eta);
if order mod 2 = 0 then Odd_order:= False else Odd_order:= True;
Const_2:= phi;
Const_1:= sinh(Const_2);
Const_2:= cosh(Const_2);
Peak:= 1.0;
if order mod 2 = 0 then
begin
for i:= 1 to order div 2 do
begin
Pole_Re[i]:= -Const_1*sin(Pi*(2*i-1)/(2*order) );
Pole_Im[i]:= Const_2*Cos(Pi*(2*i-1)/(2*order) );
Peak:= Peak * sqr( sqr(Pole_Re[i] ) + sqr(Pole_Im[i] ) );
end;
end
else
begin
for i:= 1 to (order-1) div 2 do
begin
Pole_Re[i]:= -Const_1*sin(Pi*(2*i-1)/(2*order) );
Pole_Im[i]:= Const_2*Cos(Pi*(2*i-1)/(2*order) );
Peak:= Peak * sqr( sqr(Pole_Re[i] ) + sqr(Pole_Im[i] ) );
end;
Pole_Re[i+1]:= -Const_1;
Pole_Im[i+1]:= 0.0;
Peak:= Peak * sqr( Const_1 );
end;
w:= w_beg;
m:= 1;
wo_geom:= sqrt( wp_lo*wp_hi);
T_Lhs:= -1.0;
T_Rhs:= -1.0;
while w <= w_end do
begin
Time[m]:= 0.0;
wp:= wo_geom/(wp_hi-wp_lo);
wp:= wp*(w/wo_geom - wo_geom/w );
Fr[m]:= w;
if not Odd_order then
begin
Ang:= 0.0;
for i:= 1 to order div 2 do
begin
t1:= wp - Pole_Im[i];
t2:= wp + Pole_im[i];
t3:= Pole_Re[i];
Ang_1:= -Atan4(t1,t3);
Time[m]:= Time[m] - sqr(Cos(Ang_1) ) /t3 ;
Ang:= Ang + Ang_1;
Ang_1:= -Atan4( t2,t3);
Time[m]:= Time[m] - sqr(Cos(Ang_1) )/t3;
Ang:= Ang + Ang_1;
end;
end;
if Odd_order or (order = 1) then
begin
Ang:= 0.0;
for i:= 1 to (order-1) div 2 do
begin
t1:= wp - Pole_Im[i];
t2:= wp + Pole_im[i];
t3:= Pole_Re[i];
Ang_1:= -Atan4(t1,t3);
Time[m]:= Time[m] - sqr(Cos(Ang_1) ) /t3 ;
Ang:= Ang + Ang_1;
Ang_1:= -Atan4( t2,t3);
Time[m]:= Time[m] - sqr(Cos(Ang_1) )/t3;
Ang:= Ang + Ang_1;
end;
t3:= Pole_Re[(order+1) div 2];
Ang_1:= -Atan4(wp,t3);
Time[m]:= Time[m] - sqr( Cos( Ang_1) ) /t3;
Ang:= Ang + Ang_1;
end;
Time[m]:= Time[m] * 1000 / ( Pi * (wp_hi-wp_lo) );
if w < wo_geom then
begin
if Time[m] > T_Lhs then T_Lhs:= Time[m];
end;
if w > wo_geom then
begin
if Time[m] > T_Rhs then T_Rhs:= Time[m];
end;
w:= w + w_inc;
m:= m + 1;
end;
w:= w_beg;
m:= 1;
clrscr;
while w < w_end do
begin
wp:= wo_geom/(wp_hi - wp_lo);
wp:= abs( wp*(w/wo_geom - wo_geom/w ) );
wp:= wp*1.0E10;
wp:= Int(wp)*1.0E-10;
if ( (wp > 1.0) and ( w < wp_lo) ) then
begin
mag[m]:= -10*Log10(1 + sqr(eta*Cosh(order*Arccosh(wp)))) - L_db*T_Lhs/T_do;
end;
if ( (wp > 1.0) and ( w > wp_hi) ) then
begin
mag[m]:= -10*Log10(1 + sqr(eta*Cosh(order*Arccosh(wp)))) - L_db*T_Rhs/T_do;
end;
if wp <= 1.0 then
begin
mag[m]:= -10*Log10(1 + sqr(eta*Cos(order*Arccos(wp)))) - L_db*Time[m]/T_do;
end;
if m <= 23 then
begin
gotoxy(1,m+1);
write(Fr[m]:9:4,' ',mag[m]:7:2);
end;
if ( ( m > 23 ) and ( m <= 46) ) then
begin
gotoxy(20,m-22);
write(' ',Fr[m]:9:4,' ',mag[m]:7:2);
end;
if ( ( m > 46 ) and ( m <= 69) ) then
begin
gotoxy(40,m-45);
write(' ',Fr[m]:9:4,' ',mag[m]:7:2);
end;
if ( ( m > 69 ) and ( m <= 92) ) then
begin
gotoxy(60,m-68);
write(' ',Fr[m]:9:4,' ',mag[m]:7:2);
end;
w:= w + w_inc;
m:= m + 1;
end;
repeat until keypressed;
clrscr;
gotoxy(10,5);
write('Enter Title for Graph ');
readln(Graph_Title);
for i:= 1 to m-1 do
begin
A[i,1]:= Fr[i];
A[i,2]:= mag[i];
end;
N:= m - 1;
assign(file_name,'PlotGr.Chn');
Chain(file_name);
readln(N);
clrscr;
gotoxy(10,5);
write('Enter Title for Graph for Group Delay ');
readln(Graph_Title);
for i:= 1 to m-1 do
begin
A[i,1]:= Fr[i];
A[i,2]:= Time[i];
end;
N:= m - 1;
assign(file_name,'PlotGr.Chn');
Chain(file_name);
Write_to_File( Fr,mag);
goto 500;
end.