// phase and pole/zero plot of 4th order Hilbert transform
// IIR filter pair

clear;

p1=-0.947997;  
p2= 0.814444;    
p3=-0.583754;   
p4= 0.218006;   

printf("pole position\n%10.6f\n%10.6f\n%10.6f\n%10.6f\n", p1, p2, p3, p4);

z=poly(0, 'z');
zi=((z-1/p1)*(z-1/p2)*(z-1/p3)*(z-1/p4))/((z-p1)*(z-p2)*(z-p3)*(z-p4));
zq=((z+1/p1)*(z+1/p2)*(z+1/p3)*(z+1/p4))/((z+p1)*(z+p2)*(z+p3)*(z+p4));

h1=syslin('d', zi);
h2=syslin('d', zq);
h3=syslin('d', zq/zi);

scf(1);
clf;
plzr(h1);
xtitle("Hi(z)");

scf(2);
clf;
plzr(h2);
xtitle("Hq(z)");

scf(3);
clf;
plzr(h3);
xtitle("Hp(z)=Hq(z)/Hi(z)");

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

scf(11);
clf;
plot2d(fr1, phi1, rect=[0, -720, 0.5, 0], axesflag=1, style=2, nax=[1, 3, 0, 5]);  // blue
plot2d(fr2, phi2, rect=[0, -720, 0.5, 0], axesflag=1, style=5, nax=[1, 3, 0, 5]);  // red
xtitle("phase of Hi(z), Hq(z)", "normalized frequency", "phase [degree]");

scf(12);
clf;
plot2d(fr3, phi3, rect=[0, -100, 0.5, 0], axesflag=1);  
xtitle("phase of Hp(z)=Hq(z)/Hi(z)", "normalized frequency", "phase [degree]");

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

//------------
// sub-system phase plot

za=(z+1/p1)*(z+1/p2)*(z-p1)*(z-p2)/(z-1/p1)/(z-1/p2)/(z+p1)/(z+p2);
zb=(z+1/p3)*(z+1/p4)*(z-p3)*(z-p4)/(z-1/p3)/(z-1/p4)/(z+p3)/(z+p4);

ha=syslin('d', za);
hb=syslin('d', zb);

[fra, repfa]=repfreq(ha, 0, 0.5, 0.0002);
[frb, repfb]=repfreq(hb, 0, 0.5, 0.0002);

[dba, phia]=dbphi(repfa);
[dbb, phib]=dbphi(repfb);

scf(21);
clf;
plot2d(fra, phia, rect=[0, -100, 0.5, 0], axesflag=1, style=2);   // blue  
plot2d(frb, phib, rect=[0, -100, 0.5, 0], axesflag=1, style=5);   // red
plot2d(fr3, phi3, rect=[0, -100, 0.5, 0], axesflag=1, style=19);  // brown 
xtitle("phase of sub-system", "normalized frequency", "phase [degree]");