function [acf, C] = estacf(sig, maxlag)

% Laboriously find acf from data
acf = zeros(maxlag+1,1);
for i=1:length(acf)
  
  % Current lag distance
  lag = i-1;
  
  % Current sum and number of contributing points
  tsum = 0.0;
  nump = 0;
  
  % Starting indices for pair of samples
  si1 = 1;
  si2 = si1 + lag;
  
  while si2<length(sig)
    
    % Accumulate pairwise products for lag over signal
    tsum = tsum + sig(si1)*sig(si2);
    nump = nump + 1;
    
    % Indices for next sample pair
    si1 = si1 + 1;
    si2 = si2 + 1;
    
  end
  
  % Fill in autocorrelation function - note that Matlab vector indices
  % start at 1, so acf(1) correponds to r(0)
  acf(i) = tsum/nump;
  
end

% Build covariance matrix 
C = zeros(maxlag+1, maxlag+1);
for r=1:size(C,1)
  for c=1:size(C,2)
    
    C(r,c) = acf(abs(r-c)+1);
    
  end
end
