function [xy, nobs] = nancov(x,varargin)
%   [xy, nobs] = nancov(x, [y, flag])
%   
%   NANCOV Covariance matrix.
%   Works exactly like COV, but treats NaNs as missing data. The
%   number of observations NOBS is optionally returned. Because it
%   treats NaNs specially, it runs slower than COV.
%
%   See also COV
  
%   Simon de Szoeke 1-8-02
%   adapted from COV:
%   J. Little 5-5-86
%   Revised 6-9-88 LS 3-10-94 BJ
%   Copyright (c) 1984-98 by The MathWorks, Inc.
%   $Revision: 5.12 $  $Date: 1997/11/21 23:23:35 $

if nargin==0, error('Not enough input arguments.'); end
if nargin>3, error('Too many input arguments.'); end
if ndims(x)>2, error('Inputs must be 2-D.'); end

nin = nargin;

% Check for cov(x,flag) or cov(x,y,flag)
if (nin==3) | ((nin==2) & (length(varargin{end})==1));
  flag = varargin{end};
  nin = nin - 1;
else
  flag = 0;
end

if nin == 2,
  x = x(:);
  y = varargin{1}(:);
  if length(x) ~= length(y), 
    error('The lengths of x and y must match.');
  end
  x = [x y];
end

if length(x)==prod(size(x))
  x = x(:);
end

[m,n] = size(x);

if m==1,  % Handle special case
  xy = 0;

else
  xc = x - repmat(nanmean(x),m,1);   % Remove mean
  nobs = double(isfinite(xc)')*double(isfinite(xc));
                                     % number of observations in each covariance
  xc(find(isnan(xc))) = 0;                 % zero out NaNs so they don't
                                     % contribute to the covariance
  if flag
    xy = xc' * xc ./ nobs;
  else
    xy = xc' * xc ./ (nobs-1);
  end
end
