function [xro,xroE,mi,miE,xacf] = acf(x,fig)
% DIRECT COMPUTATION OF THE AUTOCORRELATION (AUTOCOVARIANCE) FUNCTION ESTIMATE
%                 (suitable for short sequences)
% function [xro,xroE,mi,miE,xacf] = acf(x,fig)
% x ... time series data vector of length n
% fig ... figure handle for the acf plot (optional: default=[]=no plot)
%         at most 50 acf values plotted.
% xro ... estimate of the autocorrelation function of length n/4 (column vector)
%         (leading significant part only)
% xroE ... xroE(q+1)=asymptotic 95% confidence halfwidth for xro(j+1), j>q
%          (q=0,1,...n/2-1) under the hypothesis that x is an MA(q) process.
%          (column vector)
% mi   ... estimate of the mean
% miE  ... = [miE1,miE2] = 95% asymptotic confidence halfwidths of the mean
%          miE2 is under assumption of normality of the generating noise
%          miE2<miE1.
% xacf ... full estimate of the autocovariance function of length n
%          (column vector)

% (C) 1996, V. Vesely, Masaryk University of Brno (Czech Republic)
%
error(nargchk(1,2,nargin));
if size(x,1)==1, x=x.'; end
if size(x,2)>1 , error('x is not a vector'); end
n  = length(x);
if n<50
   disp('WARNING: few samples (less than 50) to get a reliable estimate');
end
mi  = mean(x);
x0  = x-mi;
xacf = conv(x0,flipud(x0));
xacf = xacf(n:2*n-1)./n;
n2   = ceil(n/2);
xro  = xacf(1:n2)./xacf(1);
%
% Computing for q=0,1,2,...,n2-1
% xroE(q+1) = 1.96*sqrt(1+2*(|xro(1+1)|^2+...+|xro(q+1)|^2)./n)
xroE = [1;2*abs(xro(2:n2)).^2];
xroE = 1.96*sqrt(sum(triu(xroE*ones(1,n2)))./n).';
%
% Computing
% miE1 = 1.96*sqrt(xacf(0+1)*(1+2*(sum{k=1,..,n2-1}{rxro(k+1)}))./n)
% miE2 = 1.96*sqrt(xacf(0+1)*(1+2*(sum{k=1,..,n2-1}{(1-(k./n))*rxro(k+1)}))./n),
% kde rxro(:) = real(xro(:))
  v1 = xacf(1)*(2*sum(real(xro))-1);
  v2 = xacf(1)*(2*sum( (1-(0:n2-1).'./n).*real(xro) )-1);
  if (v1<0 | v2<0)
     disp(['WARNING: The process behaves like non-stationary => ',...
           'the estimate miE may not be correct']);
  end
  miE1 = 1.96*sqrt(abs(v1./n));
  miE2 = 1.96*sqrt(abs(v2./n));
  miE = [miE1,miE2];
if nargin==2
   figure(fig);
% Plot data with the mean estimate:
   t=(1:n).';
   plot(t,x,'y',[1;n],[mi;mi],'w',[1;n],[1;1]*[mi+miE,mi-miE],'w:');
   title('Data with the mean estimate');
% Plot acf estimate:
   figure
   sz = get(0,'ScreenSize');
   figp = get(gcf,'Position');
   figp(2) = max(1,figp(2)-figp(4));
   set(gcf,'Position',figp);
   npl = min(n2,50);
   t=(1:npl)';
   stem(t-1,xro(t)); hold on;
   plot(t-1,[xroE(t),-xroE(t)],'w--',[1;npl]-1,[1;1]*[xroE(1),-xroE(1)],'w:');
   title('Autocorrelation function');
   hold off;
end
