Ntheta = 300; htheta = 2*pi/Ntheta;
theta = 0:htheta:2*pi;
z = exp(1i*theta);
rz = real(z);
iz = imag(z);

% poles

z1p = 0.9*exp(1i*pi/2);
z2p = 0.2*exp(1i*29*pi/32);
z3p = conj(z1p);
z4p = conj(z2p);
P = [z1p z2p z3p z4p];
RP = real(P);
IP = imag(P);

% zeros

z1z = 0.96*exp(1i*3*pi/8);
z2z = 0.94*exp(1i*11*pi/16);
z3z = 0.90*exp(1i*pi/16);
z4z = conj(z1z);
z5z = conj(z2z);
z6z = conj(z3z);
Z = [z1z z2z z3z z4z z5z z6z];
RZ = real(Z);
IZ = imag(Z);

% the transfer function evaluated on the unit circle

p = (z-z1z).*(z-z2z).*(z-z3z).*(z-z4z).*(z-z5z).*(z-z6z);
q = (z-z1p).*(z-z2p).*(z-z3p).*(z-z4p);
H = p./q;
rH = real(H);
iH = imag(H);
aH = abs(H);

% plot to MATLAB figure file

L = 30;
LL = 16;
figure(1)
subplot(1,3,1)
plot(rz,iz,'b.')
hold on
plot(RP,IP,'ro')
plot(RZ,IZ,'bd')
hold off
title('Zeros and Poles','fontsize',12)
grid on, box on
axis([-1.1 1.1 -1.1 1.1]), set(gca,'xtick',-1:1:1,'ytick',-1:1:1)
xlabel Re(z), ylabel Im(z)
subplot(1,3,2)
plot(theta,aH)
title('|H(e^{i\theta})|','fontsize',12)
grid on, box on
axis([0 2*pi 0 LL]), set(gca,'xtick',0:pi/2:2*pi,'ytick',0:LL/4:LL)
xlabel \theta, ylabel |H(e^{i\theta})|
subplot(1,3,3)
plot(rH,iH,'bo') 
title('H(e^{i\theta})','fontsize',12)
grid on, box on
axis([-L L -L L]), set(gca,'xtick',-L:L:L,'ytick',-L:L:L) 
xlabel Re(z), ylabel Im(z)