% Quick demo of LD recursion ....
% Last revised 11/18/97 - rms

% First consider a pole at z = .9
ps = [.9];
zs = [0];
zplane0(zs,ps);
title('Poles and zeros of original system:','fontsize',18)
pause
alpha = .9;	
n = 0:100;
h = alpha.^n;	% impulse response
stem(n,h)
xlabel('n','fontsize',18)
ylabel('h[n]','fontsize',18)
axis([0 30 0 1])
pause
R = conv(h,fliplr(h));
n = -100:100;
stem(n,R)
xlabel('n','fontsize',18)
ylabel('R[n]','fontsize',18)	% correlation function
axis([-30 30 0 inf])
pause
R = R(101:201);

for p = 1:4,


a = levinson(R,p);
ps = roots(a);
zs = [];

sp = num2str(p);
zplane0(zs,ps);
title(['Poles of model of Order ',sp],'fontsize',18);
pause

[H2, W] = freqz(1, [1 -alpha]);
plot(W/pi,20*log10(abs(H2)),'r');	% Plot ideal curve and hold
hold on

[H,W] = FREQZ(1,a);
plot(W/pi,20*log10(abs(H)));
title(['Response of model of Order ',sp],'fontsize',18);
xlabel('\omega/\pi','fontsize',18)
ylabel('20*log10(|H|)','fontsize',18)
pause
hold off

end




% Now we'll repeat for a zero at z = .9
ps = [0];
zs = [.9];
zplane0(zs,ps);
title('Poles and zeros of original system:','fontsize',18)
pause
alpha = .9;	
n = 0:100;
h = [1 -1*alpha zeros(1,99)];
stem(n,h)
xlabel('n','fontsize',18)
ylabel('h[n]','fontsize',18)
axis([0 30 -1 1])
pause
R = conv(h,fliplr(h));
n = -100:100;
stem(n,R)
xlabel('n','fontsize',18)
ylabel('R[n]','fontsize',18)
axis([-30 30 -inf inf])
pause
R = R(101:201);

for p = 1:10,

[H2] = freqz([1 -alpha], 1);
plot(W/pi,20*log10(abs(H2)),'r');	% Plot ideal curve and hold
hold on


a = levinson(R,p);
ps = roots(a);
zs = [];

sp = num2str(p);
zplane0(zs,ps);
title(['Poles of model of Order ',sp],'fontsize',18);
pause

[H2] = freqz([1 -alpha], 1);
plot(W/pi,20*log10(abs(H2)),'r');	% Plot ideal curve and hold
hold on

[H,W] = freqz(1,a);
plot(W/pi,20*log10(abs(H)));
title(['Response of model of Order ',sp],'fontsize',18);
xlabel('\omega/\pi','fontsize',18)
ylabel('20*log10(|H|)','fontsize',18)
pause
hold off

end


for p = 20:10:50,

[H2, W] = freqz(1,[1 -alpha]);
plot(W/pi,20*log10(abs(H2)),'r');	% Plot ideal curve and hold
hold on

a = levinson(R,p);
ps = roots(a);
zs = [];

sp = num2str(p);
zplane0(zs,ps);
title(['Response of model of Order ',sp],'fontsize',18);
xlabel('\omega/\pi','fontsize',18)
ylabel('20*log10(|H|)','fontsize',18)
pause

[H2, W] = freqz([1 -alpha], 1);
plot(W/pi,20*log10(abs(H2)),'r');	% Plot ideal curve and hold
hold on

[H,W] = FREQZ(1,a);
plot(W/pi,20*log10(abs(H)));
title(['Poles of model of Order ',sp],'fontsize',18);
pause
hold off

end