% This is a calling script for the z-plane routines
% Last revised 9/20/99 - rms

% First, a simple pole and zero pair

newplot;
zeroz = [0 1]';
poles = [.8*exp(j*pi/4) .8*exp(-j*pi/4)]';

zplane0(zeroz,poles);
pause
zplane3(zeroz,poles);
pause
freqz2(zeroz,poles);
pause

% Moving the poles closer to the unit circle:

newplot;
zeroz = [0 1]';
poles = -[.9*exp(j*pi/4) .9*exp(-j*pi/4)]';

zplane0(zeroz,poles);
pause
zplane3(zeroz,poles);
pause
freqz2(zeroz,poles);
pause


% A notch filter

newplot;
zeroz = [exp(j*pi/4) exp(-j*pi/4)]';
poles = [.9*exp(j*pi/4) .9*exp(-j*pi/4)]';

zplane0(zeroz,poles);
pause
zplane3(zeroz,poles);
pause
freqz2(zeroz,poles);
pause


% Now, some real filters .... specs [Wp Ws Rp Rs] = [0.30    0.40    1.00   60.00]

% First, using ellip, with [N Wn] = [6 .3]

Wp = .3; Ws = .4; Rp = 1; Rs = 60;
[N Wn] = ellipord(Wp, Ws, Rp, Rs);
[B,A] = ELLIP(N,Rp,Rs,Wn);

zeroz = roots(B);
poles = roots(A);

newplot;

zplane0(zeroz,poles);
pause
zplane3(zeroz,poles);
pause
freqz2(zeroz,poles);
pause
freqz(B,A);
pause


% Now repeating, but using Parks-McClellan

f = [0 .3 .4 1];
a = [1 1 0 0];

N = 65;
B=REMEZ(N,f,a);	% N = 65 obtained empirically ...Èmax(20*log10(abs(FREQZ(B,A,W)))) = ~-.60

zeroz = roots(B);
A = 1;

newplot;

zplane0(zeroz,[]);
pause
zplane3(zeroz,[]);
pause
freqz2(zeroz,[]);
pause
freqz(B,A);