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

% First consider a pole at z = .9
figure
ps = [.9];
zs = [0];
zplane(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
% First consider a pole at z = .9
ps = [.9];
zs = [0];
zplane(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,

hold off
a = levinson(R,p);
ps = roots(a);
a = poly(ps);
zs = [];
b = poly(zs);

sp = num2str(p);
hold off
zplane(zs,ps);
title(['Poles of model of Order ',sp],'fontsize',18);
axis([-1.25,1.25,-1.25,1.25])
pause
hold off
%close all
[H2, W] = freqz(1, [1 -alpha]);
plot(W/pi,20*log10(abs(H2)),'r','linewidth',2);	% Plot ideal curve and hold
hold on
% freqz(1,a);
[H,W] = freqz(1,a);
plot(W/pi,20*log10(abs(H)),'b','linewidth',2);
title(['Response of model of Order ',sp],'fontsize',18);
xlabel('\omega/\pi','fontsize',18)
ylabel('20*log10(|H|)','fontsize',18)
axis([0,1,-40,20])
pause
hold off

end




% Now we'll repeat for a zero at z = .9
ps = [0];
zs = [.9];
hold off
zplane(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,

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


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

sp = num2str(p);
hold off
zplane(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','linewidth',2);	% Plot ideal curve and hold
hold on

[H,W] = freqz(1,a);
plot(W/pi,20*log10(abs(H)),'b','linewidth',2);
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','linewidth',2);	% Plot ideal curve and hold
hold on

a = levinson(R,p);
ps = roots(a);
zs = [];
hold off
sp = num2str(p);
zplane(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','linewidth',2);	% Plot ideal curve and hold
hold on

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

end
