Implementation of ARYule

The graphs show an expt testing  ARYule 100 times on a 2-pole system with 2^8 data points:

>> for i=1:100
AAaryule
drawnow
end
>>

% Aryule method as implemented in AA Book(ref)

N=2^7;n=0:N-1;g=randn(size(n));%3.5*(rand(size(n))-0.5);
%b=1;a=[1 -.9 .3];
b=1;a=conv([1 -.9 .8],[1 .3 .8]);
x=filter(b,a,g);

maxlag=15;
p=4;
[rx,lags]=xcorr(x,x,maxlag,'biased');
rx=rx(maxlag+1:maxlag+1+p);
R=toeplitz(rx);
R1=inv(R);
b2=1/R1(1,1);b0=sqrt(b2);            %find b(0) squared
a1=b2*R1(2:p+1,1);                      %fin AR coefficients a(1) - a(5)

% Obtain a coeffs using 'inner' matrix
RR=R(1:p,1:p);
a2=inv(RR)*-rx(2:p+1)';
[a1(1:p-1) a2(1:p-1)];
[H,w]=freqz(b,a);[H1,w]=freqz(b,[1 a1']);
figure(1)
plot(w,20*log10(abs(H)),w,20*log10(abs(H1)));hold on
axis([0 pi -20 20]);

figure(2)
th=0:.1:2*pi;plot(cos(th),sin(th));hold on
plot(roots([1 a1']),'xr')
%zplane(1,[1 a1']);hold on