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

clear;

p1=-0.975418; 
p2= 0.917498;    
p3=-0.829167;
p4= 0.685127;    
p5=-0.463441;   
p6= 0.165534;   

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

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

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-180, rect=[0, -1080, 0.5, 0], axesflag=1, style=2, nax=[1, 3, 0, 7]);  // blue
plot2d(fr2, phi2-180, rect=[0, -1080, 0.5, 0], axesflag=1, style=5, nax=[1, 3, 0, 7]);  // 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);

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

[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);

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

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(fr3, phi3, rect=[0, -100, 0.5, 0], axesflag=1, style=19);  // brown 
xtitle("phase of sub-system", "normalized frequency", "phase [degree]");