function Vcovar=yulewalk(a,b,sigmae,N)
% Computes the first N covariance samples for the ARMA process
%  A(z) y(k) = B(z) e(k)
% The coefficients of A(z) and B(z) are stored in row vectors A and B in ascending powers of z^-1
% sigmae=variance of white noise e(k)

% NOTA: c'e' un metodo molto piu' furbo! (anche se approssimato)
%		G=tf(b,[a zeros(m-n)],1);     
%		g=impulse(G,1000);
%		c=xcorr(g,g,N);
%

n=length(a)-1;
m=length(b)-1;

if m==0,
	% AR process (b is a scalar)
	a=a/b;
	rcol=[sigmae; zeros(n,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)*rcol;
	for i=n+2:N+1,
		Vcovar(i)=-a(2:n+1)*flipud(Vcovar(i-n:i-1));
	end
	Vcovar=Vcovar(1:N+1);
else
	aext=[a zeros(1,m-n)];
	Ma=toeplitz(aext(1:m+1)',[aext(1) zeros(1,m)]);
	Vcovey=inv(Ma)*sigmae*b';
	Mey=toeplitz([Vcovey(1); zeros(n,1)],Vcovey');
	rcol=Mey*b';
	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)*rcol;
	Vcoveyext=[zeros(N,1); Vcovey];
	for i=n+2:N+1,
		Vcovar(i)=-a(2:n+1)*flipud(Vcovar(i-n:i-1))+b*Vcoveyext(N-i+2:N-i+2+m);
	end
	Vcovar=Vcovar(1:N+1);
end

