// conventional IIR Hilbert transform filter pair example 

clear;

a0=0.6923878;
a1=0.9360654;
a2=0.9882295;
a3=0.9987488;

b0=0.4021921;
b1=0.8561711;
b2=0.9722910;
b3=0.9952885;
 
z=poly(0, 'z');

z1=   ((a0^2)*(z^2)-1)/(z^2-a0^2);
z1=z1*((a1^2)*(z^2)-1)/(z^2-a1^2);
z1=z1*((a2^2)*(z^2)-1)/(z^2-a2^2);
z1=z1*((a3^2)*(z^2)-1)/(z^2-a3^2);

z2=   ((b0^2)*(z^2)-1)/(z^2-b0^2);
z2=z2*((b1^2)*(z^2)-1)/(z^2-b1^2);
z2=z2*((b2^2)*(z^2)-1)/(z^2-b2^2);
z2=z2*((b3^2)*(z^2)-1)/(z^2-b3^2);


h1=syslin('d', z1);
h2=syslin('d', z1/z);
h3=syslin('d', z2);
h4=syslin('d', z1/z/z2);

scf(1);  clf;  plzr(h1);  
scf(2);  clf;  plzr(h2);
scf(3);  clf;  plzr(h3);
//scf(4);  clf;  plzr(h4);

[fr1, repf1]=repfreq(h1, 0, 0.5, 0.0001);
[fr2, repf2]=repfreq(h2, 0, 0.5, 0.0001);
[fr3, repf3]=repfreq(h3, 0, 0.5, 0.0001);
[fr4, repf4]=repfreq(h4, 0, 0.5, 0.0001);
[db1, phi1]=dbphi(repf1);
[db2, phi2]=dbphi(repf2);
[db3, phi3]=dbphi(repf3);
[db3, phi4]=dbphi(repf4);

scf(11);
clf;
plot2d(fr1, phi1, rect=[0, -1440-180, 0.5, 0], axesflag=1, style=3, nax=[1, 3, 0, 10]);  // green
plot2d(fr2, phi2, rect=[0, -1440-180, 0.5, 0], axesflag=1, style=2, nax=[1, 3, 0, 10]);  // blue
plot2d(fr3, phi3, rect=[0, -1440-180, 0.5, 0], axesflag=1, style=5, nax=[1, 3, 0, 10]);  // red

scf(12);
clf;
plot2d(fr4, abs(phi4+90), rect=[0, 0, 0.5, 1], axesflag=1);
xtitle("phase error", "normalized frequency", "phase [degree]");