A = [1 -1 0.5]; B = 1;
sig2 = 1; N = 256;
Nspec = 128;

f = (0:Nspec-1)/Nspec;
g = abs(freqz(B,A,Nspec)).^2;
count = 100;

avsx = zeros(Nspec,1);
for ii = 1:count
  x = genarma(A,B,sig2,N);
  sx = abs(freqz(x,1,Nspec)).^2/N;
  avsx = avsx + sx;
  plot(f,[avsx/ii g]); drawnow;
end;
avsx = avsx / count;
plot(f,[avsx g]);
