function acv = nanacv(X,N)
%{
    Computes the autocovariance for series X over N lags.
    (Alas, very slowly...)
    Ignores NaNs in X.

    Jan 2016 BWB
%}

acv = zeros(2*N+1,1);
out = zeros(N+1,1);
X = X - nanmean1(X);
out(1) = nanvar1(X); % zero lag value is the variance

% iterate for N lags, rotate time series one point at each step
% Xs = X;
% for i=2:N+1
%     Xn = Xs(N);	% save last pnt
%     Xs(2:end) = Xs(1:end-1); Xs(1) = Xn;    % shift series one point
%     aa = nancov1(X,Xs,'pairwise'); out(i) = aa(1,2);    % compute cov 
% end;

% This method appears closer to the xcov(X, 'unbiased') result.
% But it is one-sided with N+1 points.  xcov() output is 2N+1 points.
for i=2:N-1     % just shift and compute over overlapping segment
    aa = nancov1(X(i:end),X(1:end-i+1),'pairwise'); out(i) = aa(1,2);
end;

% figure; plot(0:N,out); grid; title('Y = nanacv(X,N), overlap mode');

% create 2-sided output array
acv(N+1:end) = out(1:end);
acv(1:N) = flipud(out(2:end));

end