function [pacf,phi,phiE,v] = pacf(acf,n,fig,pacf,phi,v)
% COMPUTING THE PARTIAL AUTOCORRELATION FUNCTION AND BEST LINEAR PREDICTION
%                              COEFFICIENTS
%                       (Durbin-Levinson algorithm)
% function [pacf,phi,phiE,v] = pacf(acf,n,fig,pacf,phi,v)
% acf ... estimate of the autocovariance function obtained via function 'acf'.
%         of a stationary zero-mean random process {X(j)}, p=length(acf)-1 >= 1
% n   ... number of data samples used for acf estimate
% fig ... figure handle for the pacf plot (optional: default=[]=no plot)
%         at most 50 pacf values plotted.
% pacf... estimate of the partial autocorrelation function
%         input : previous estimate of length(pacf)<length(acf)
%                 default: pacf=[] <=> pacf=1
%         output: updated up-to-date estimate of length(pacf)=length(acf)
%                 (column vector)
% phi ... coefficients of the best linear prediction of X(m+1) via X(j), j<=m
%         X(m+1)=phi(1)*X(m)+phi(2)*X(m-1)+...+phi(m)*X(1) where m=length(phi)
%         It holds: pacf(m+1)=phi(m), length(pacf)=m+1.
%         If x is AR(p) process then phi is estimate of its AR coefficients.
%         input : previous estimate   for m<p
%                 default: phi=[], m=0
%         output: up-to-date estimate for m=p (column vector)
% phiE... phiE(j)=asymptotic 95% confidence halfwidth for phi(j)
% v   ... estimate of the mean square error (MSE): E(Xest(m+1)-X(m+1))^2 (scalar)
%         If x is AR(p) process then v is the estimate of the variance of the
%         generating noise.
%         input : previous estimate   for m<p
%                 default: v=[]
%         output: up-to-date estimate for m=p

% (C) 1996, V. Vesely, Masaryk University of Brno (Czech Republic)
%
error(nargchk(2,6,nargin));
if nargin<6, v    = []; end
if nargin<5, phi  = []; end
if nargin<4, pacf = []; end
if nargin<3, fig  = []; end
if isempty(pacf), pacf = 1 ; end
p = min(length(acf),ceil(n/2))-1;
m = length(phi);
if m>=p, error('Submitted acf is too short'); end
if length(pacf)~=m+1,
   error('Submitted pacf has to be longer then phi exactly by one');
end
if xor(isempty(v),isempty(phi)), error('phi or v estimate missing'); end
phi = [reshape(phi ,m  ,1); zeros(p-m,1)];
pacf= [reshape(pacf,m+1,1); zeros(p-m,1)];
if isempty(v), v=acf(1); end
for k=m:p-1
    if k>0
       phi(k+1)  = ( acf(k+2)-dot(phi(k:-1:1),acf(2:k+1)) )./v;
    else
       phi(k+1)  = acf(k+2)./v;
    end
    phi(1:k)  =   phi(1:k)-phi(k+1)*conj(phi(k:-1:1));
    v         =   v*(1-abs(phi(k+1))^2);
    pacf(k+2) =   phi(k+1);
end
if nargout>2
   phiE  = 1.96*sqrt(v*diag(inv(toeplitz(acf(1:p))))./n);
   if ~isempty(fig)
      figure(fig);
      errorbar(1:length(phi),phi,phiE,phiE,'o');
      title('Autoregression coefficients');
      fig = figure;
      sz = get(0,'ScreenSize');
      figp = get(gcf,'Position');
      figp(2) = max(1,figp(2)-figp(4));
      set(gcf,'Position',figp);
   end
end
pacfE = 1.96*sqrt(1./n);
% Plot pacf estimate:
  if ~isempty(fig)
     figure(fig)
     npl = min(p+1,50);
     t=(1:npl)';
     stem(t-1,pacf(t)); hold on;
     plot([1;npl]-1,[1;1]*[pacfE,-pacfE],'w--');
     title('Partial autocorrelation function');
     hold off;
end
