function polezero(B, A, N, FS)
% polezero(B, A, N, FS)
%
% Draws a pole-zero plot and frequency spectrum for the following filter:
%
%           B(1) + B(2)z^(-1) + B(3)z^(-2) + ...
% H(z) = -------------------------------------------
%           A(1) + A(2)z^(-1) + A(3)z^(-2) + ...
%
% Poles and zeros are drawn in the complex Z-plane.
% Poles are drawn as 'x', zeros are drawn as 'o'.
% The unit circle is drawn with a solid line.
% Only the positive-imaginary part of the Z-plane is plotted; 
%  we assume that all roots are either real or in complex-conjugate pairs.
%
% The frequency spectrum is drawn using the freqz function.
%

% Calculate the spectrum
if nargin>3,
 H = freqz(B, A, N, FS);
else
 H = freqz(B, A, N, FS);
end;
clf;

% Plot the unit circle
angles = [0:0.01:(2*pi)];
subplot(2,2,1);
plot(cos(angles),sin(angles));
hold on;

% Plot the Re(z) and Im(z) axis
ticks = [-1.5:0.01:1.5];
plot(ticks, zeros(size(ticks)));
plot(zeros(size(ticks)), ticks);
ylabel('Imaginary Part of Z');
xlabel('Real Part of Z');

% Plot the poles and zeros
zeros = roots(B);
poles = roots(A);
plot(real(zeros),imag(zeros),'ro');
plot(real(poles),imag(poles),'rx');
axis([-1.1 1.1 -1.1 1.1]);
hold off;

% Plot the magnitude spectrum
subplot(2,1,2);
f = [0:(N-1)]*FS/(2*N);
plot(f, 20*log10(abs(H)));
xlabel('Frequency (Hertz)');
ylabel('Magnitude Response (dB)');
