function Vcovar=yulewalk(a,sigmae,N)
% Computes the first N covariance samples for the AR process
%  A(z) y(k) =  e(k)
% The coefficients of A(z) are stored in a (row vector)
% sigmae=variance of white noise e(k)

n=length(a)-1;
M1=fliplr(toeplitz([a(n+1); zeros(n,1)],fliplr(a)));
M2=[zeros(1,n+1); zeros(n,1) toeplitz(a(1:n)',[a(1) zeros(1,n-1)])];
Vcovar=inv(M1+M2)*[sigmae; zeros(n,1)];
for i=n+2:N+1,
Vcovar(i)=-a(2:n+1)*flipud(Vcovar(i-n:i-1));
end
Vcovar=Vcovar(1:N+1);
