%--------------------------------------------------
% Tony Hyun Kim
% 6.003: PS#9, Matlab exercise
%--------------------------------------------------
function lab9

%--------------------------------------------------
% Part (a)
%--------------------------------------------------
b = [1];
a1 = [1 0 1];
a2 = [1 1/2 1];
a3 = [1 2 1];
a4 = [1 4 1];

subplot(221);
plotpz(b,a1)
title('[1 0 1]');

subplot(222);
plotpz(b,a2);
title('[1 1/2 1]');

subplot(223);
plotpz(b,a3);
title('[1 2 1]');

subplot(224);
plotpz(b,a4);
title('[1 4 1]');
axis([-5 5 -5 5]);

%--------------------------------------------------
% Part (b)
%--------------------------------------------------
omega = [-5:0.1:5];
figure;

H1 = abs(freqs(b,a1,omega));
H2 = abs(freqs(b,a2,omega));
H3 = abs(freqs(b,a3,omega));
H4 = abs(freqs(b,a4,omega));

subplot(221);
plot(omega,H1);
title('[1 0 1]');

subplot(222);
plot(omega,H2);
title('[1 1/2 1]');

subplot(223);
plot(omega,H3);
title('[1 2 1]');

subplot(224);
plot(omega,H4);
title('[1 4 1]');

%--------------------------------------------------
% Part (c)
%--------------------------------------------------
figure;
zetarange = [0 logspace(-1,1,99)];
azeta = zeros(3,100);
zetapoles = zeros(2,100);
for k = 1:100
    azeta(:,k) = [1 2*zetarange(k) 1]';
    r = roots(azeta(:,k));
    if (length(r) > 1)
        zetapoles(:,k) = r;
    else
        zetapoles(:,k) = [r r];
    end
end
plot(real(zetapoles(1,:)),imag(zetapoles(1,:)),'x','MarkerSize',20);
hold on;
plot(real(zetapoles(2,:)),imag(zetapoles(2,:)),'x','MarkerSize',20);
plotpz(b,a1);
plotpz(b,a2);
plotpz(b,a3);
plotpz(b,a4);
axis('equal');
axis([-4 0 -2 2]);

%--------------------------------------------------
% Plot the pole zero plot. Coefficient vectors
% 'b' and 'a' are numerator and denominator resp.
%--------------------------------------------------
function [ps zs] = plotpz(b,a)
zs = roots(b);
ps = roots(a);
plot(real(zs),imag(zs),'ro','MarkerSize',20);
hold on;
plot(real(ps),imag(ps),'rx','MarkerSize',20);
grid;
axis([-3 3 -3 3]);