% AR_LEVINSON.M
%
% P. Flandrin, 30 sept. 2003
%
% identifie un mod�le AR par l'algorithme
% de Levinson
%
% inputs  = - x    : signal � identifier
%           - pmax : ordre max du mod�le
%
% outputs = - AA   : matrice [pmax x pmax] des coeffs AR
%           - P    : puissance de l'erreur de pr�diction
%           - FPE  : Final Prediction Error (Akaike)
%           - popt : ordre optimal au sens FPE

function [AA,P,FPE,popt] = AR_levinson(x,pmax);

N = length(x);
c = xcorr(x);

ceff = c(N:N+pmax+1);
ceff = ceff/ceff(1);

AA = [];
P = [];
Q = [];
k = [];

AA(1,1) = -ceff(2)/ceff(1);
P(1) = ceff(1) + AA(1,1)*ceff(2);

for m = 1:pmax-1
	
	Q(m) = ceff(m+2) + AA(m,1:m)*ceff(m+1:-1:2)';
	k(m) = -Q(m)/P(m);
	P(m+1) = P(m)*(1-k(m)^2);
	
	for n = 1:m
		
		AA(m+1,n) = AA(m,n) + k(m)*AA(m,m+1-n);
		
	end

	AA(m+1,m+1) = k(m);

end

FPE = P.*(N+1+(1:pmax))./(N+1-(1:pmax));
popt = find(FPE == min(FPE));