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

clear;

p1=-0.993551; 
p2= 0.977820;  
p3=-0.952320;   
p4= 0.906600;    
p5=-0.824390;   
p6= 0.683398;    
p7=-0.463090;  
p8= 0.165550;   

printf("pole position\n");
printf("%10.6f\n%10.6f\n%10.6f\n%10.6f\n", p1, p2, p3, p4);
printf("%10.6f\n%10.6f\n%10.6f\n%10.6f\n", p5, p6, p7, p8);
   
z=poly(0, 'z');
zi=(z-1/p1)*(z-1/p2)*(z-1/p3)*(z-1/p4)*(z-1/p5)*(z-1/p6)*(z-1/p7)*(z-1/p8);
zi=zi/(z-p1)/(z-p2)/(z-p3)/(z-p4)/(z-p5)/(z-p6)/(z-p7)/(z-p8);
zq=(z+1/p1)*(z+1/p2)*(z+1/p3)*(z+1/p4)*(z+1/p5)*(z+1/p6)*(z+1/p7)*(z+1/p8);
zq=zq/(z+p1)/(z+p2)/(z+p3)/(z+p4)/(z+p5)/(z+p6)/(z+p7)/(z+p8);

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, -1440, 0.5, 0], axesflag=1, style=2, nax=[1, 3, 0, 9]);  // blue 
plot2d(fr2, phi2, rect=[0, -1440, 0.5, 0], axesflag=1, style=5, nax=[1, 3, 0, 9]);  // 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);
zc=(z+1/p5)*(z+1/p6)*(z-p5)*(z-p6)/(z-1/p5)/(z-1/p6)/(z+p5)/(z+p6);
zd=(z+1/p7)*(z+1/p8)*(z-p7)*(z-p8)/(z-1/p7)/(z-1/p8)/(z+p7)/(z+p8);

ha=syslin('d', za);
hb=syslin('d', zb);
hc=syslin('d', zc);
hd=syslin('d', zd);

[fra, repfa]=repfreq(ha, 0, 0.5, 0.0002);
[frb, repfb]=repfreq(hb, 0, 0.5, 0.0002);
[frc, repfc]=repfreq(hc, 0, 0.5, 0.0002);
[frd, repfd]=repfreq(hd, 0, 0.5, 0.0002);

[dba, phia]=dbphi(repfa);
[dbb, phib]=dbphi(repfb);
[dbc, phic]=dbphi(repfc);
[dbd, phid]=dbphi(repfd);

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(frc, phic, rect=[0, -100, 0.5, 0], axesflag=1, style=6);   // pink  
plot2d(frd, phid, rect=[0, -100, 0.5, 0], axesflag=1, style=13);  // dark green
plot2d(fr3, phi3, rect=[0, -100, 0.5, 0], axesflag=1, style=19);  // brown 
xtitle("phase of sub-system", "normalized frequency", "phase [degree]");