b=1;
a=[1 -1.2 .9];
% zplane(b,a);
n=1024;
x=randn(n,1);
y=filter(b,a,x);
pery=abs(fft(y));
pery=(pery.*pery)/(2*pi*n);
[h,w]=freqz(b,a,n/2);
% the following scale factor makes \int_0^{2\pi}h^2 df = 1
h=h/(sqrt(2*pi));
win=hanning(32)/(sum(hanning(32)));
pery=conv(pery,win);
pery=pery(16:length(pery)-16);

subplot(211);
plot(w,log10(pery(1:n/2)),'-',w,2*log10(abs(h)),'--');
v=axis;
axis([w(1) w(length(w)) v(3) v(4)]);
title('Smoothed Periodogram and theoretical density');
ylabel('LOG10');
xlabel('Frequency in radians');
gtext('Dashed curve is theoretical density');
gtext(['Hanning window width = ',int2str(length(win))]);

subplot(212);
j1=floor(n*.5/(2*pi));
j2=ceil(n*1.5/(2*pi));
plot(w(j1:j2),log10(pery(j1:j2)),'-',w(j1:j2),2*log10(abs(h(j1:j2))),'--'); 
v=axis;
axis([w(j1) w(j2) v(3) v(4)]);
title('Smoothed Periodogram and theoretical density');
xlabel('Frequency in radians');
ylabel('LOG10');
print -deps pgram9
 
% figure out the scale related to variance of y or maybe pi
% look in B&D ;  compute the variance theor from the coefs
% check what fft does and also freqz