Use of TFESTIMATE to plot the frequency response of a system from its I/O data

Input and output data with frequency response estimate using tfestimate the implementation of pwelch for cross power spectum. Zooming in on both is interesting:

The m-file:

N=2^16;n=0:N-1;x=randn(size(n));
b=[0 1];a=[1 -1.8 .9];
y=filter(b,a,x);
subplot(1,2,1)
plot(n,y,n,x)
title('Input and output data; N=2^1^6;')
legend('output, y','input, x')
[H,w]=freqz(b,a);
[Pxyh,F]=tfestimate(x,y,hamming(N),0,N,1);
[Pxyr,F]=tfestimate(x,y,rectwin(N),0,N,1);
subplot(1,2,2)
plot(F,abs(Pxyr),F,abs(Pxyh),w/2/pi,abs(H))
title('Frequency response and estimates using tfestimate')
legend('tfe:Rectwin','tfe:Hamming','Freq Resp of filter')