% E03_5.m            20oct14
% tf2zp,  zp2tf,  zplane, impz, zp2sos
clear all; clf;

%%  slide 3.5-50
num = [2  16  44  56  32];  %  numerator coefficients:  z^(-n), n=0...N
den = [3  3  -15  18  -12];   % denominator coefficients
% num = [2  6  4  6  3];  %  numerator coefficients:  z^(-n), n=0...N
% den = [3  -2  1  1  -2];   % denominator coefficients
[z,p,k] = tf2zp(num,den);
disp('Zeros are at');  disp(z);   disp('Poles are at');  disp(p);
disp('Gain constant');  disp(k);   disp('Radius of poles');  disp(abs(p));
sos = zp2sos(z,p,k);   disp('Second-order sections');  disp(real(sos));
figure(1);  zplane(num,den);  

[b, a] = zp2tf(z,p,k);  b, a     % num and den are scaled by d0=3

figure(2);
subplot(211); zplane(sos(1,1:3),sos(1,4:6));
subplot(212); zplane(sos(2,1:3),sos(2,4:6));

% fvtool(num,den);   % (unstable filter for causal signals; check impulse response)
% figure(3);  impz(num,den,20);   % inverse transform by long division

return;

%%  slide 03.5-60
num1 = [1  2];    den1 = [1  .4  -.12];   % slide 03.5-60
num2 = impz(num1,den1,6)',   den2 = 1;   % FIR approximation
figure(2);
subplot(221); zplane(num1,den1);
subplot(222); zplane(num2,den2);
subplot(223); impz(num1,den1);
subplot(224); impz(num2,den2);

figure(3);
[re,im] = meshgrid(-3:.05:2,-1.5:.05:1.5);
z = re + j*im;
X1 = (num1(1).*z.^0 + num1(2).*z.^-1) ./   ...
    (den1(1).*z.^0 +den1(2).*z.^-1 +den1(3).*z.^-2);
X1m = log(abs(X1)+.01);
mesh(re,im,X1m); axis([-3 2 -1.5 1.5 -5 5]);

figure(4);
X2 = zeros(size(z));
for k = 1:length(num2)
    X2 = X2 + num2(k) .* z.^(-k+1);
end
X2m = log(abs(X2)+.01);
mesh(re,im,X2m); axis([-3 2 -1.5 1.5 -5 5]);

%%
 fvtool(num1,den1);
 fvtool(num2,den2);


