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

function [AA,KK,CC,P,FPE,popt] = AR_levinson2(x,pmax,cmax);

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

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

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

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

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);
	KK(m+1) = k(m);
	CC(m+1,1) = - AA(m+1,1);
	
	for nn = 2:m+1
		
		CC(m+1,nn) = - AA(m+1,nn);
	
		for l = 1:nn-1
		
			CC(m+1,nn) = CC(m+1,nn) - (l/nn)*CC(m+1,l)*AA(m+1,nn-l);
			
		end
	
	end
	
	for nn = m+2:cmax
		
		CC(m+1,nn) = 0;
		
		for l = 1:m+1
			
			CC(m+1,nn) = CC(m+1,nn) - (1-l/nn)*CC(m+1,nn-l)*AA(m+1,l);
			
		end
	
	end

end

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