function [lagcovs, lags, nobs] = lagcov(a,b,maxlag)
  % [lagcovs, lags, nobs] = lagcov(a,b,maxlag)
  % finds the lagged covariance of the vectors a and b,
  % lagged from a minimum of -maxlag to a maximum of maxlag
  % indices.  lags contains the vector of the number of lags for
  % the covariances in lagcovs.  A negative lag means b lags behind a.
    
  N=min(length(a),length(b));
  maxlag=min(maxlag,N);

  lagcovs = zeros(2*maxlag+1,1);
  nobs    = zeros(2*maxlag+1,1);
  lags    = zeros(2*maxlag+1,1);
  
  for lag = -maxlag:0;
    % lag<0 ==> b lags a
    lagcovs(maxlag+1+lag) = nanmean( (a(1:N+lag)-nanmean(a(1:N+lag))).*(b(1-lag:N)-nanmean(b(1-lag:N))) );
    nobs(maxlag+1+lag)    = sum(isfinite(a(1:N+lag)) & isfinite(b(1-lag:N)));
    lags(maxlag+1+lag)    = lag;
  end
  for lag = 1:maxlag;
    lagcovs(maxlag+1+lag) = nanmean( (a(1+lag:N-maxlag+lag)-nanmean(a(1+lag:N-maxlag+lag))).*(b(1:N-maxlag)-nanmean(b(1:N-maxlag))) );
    nobs(maxlag+1+lag)    = sum(isfinite(a(1+lag:N-maxlag+lag)) & isfinite(b(1:N-maxlag)));
    lags(maxlag+1+lag)    = lag;
  end