GetCov(x, y)

function result = GetCov(x, y)

%   GetCov(X,Y) normalizes by (N-1) if N>1, where N is the number of
%   observations.  This makes COV(X) the best unbiased estimate of the
%   covariance matrix if the observations are from a normal distribution.
%   For N=1, COV normalizes by N.
%
%   COV(X,Y,1) normalizes by N and produces the second
%   moment matrix of the observations about their mean.  

if nargin==0 
  error(message('GetCov:NotEnoughInputs')); 
end
if nargin>3
  error(message('GetCov:TooManyInputs')); 
end

nin = nargin;

% Check for cov(x,flag) or cov(x,y,flag)
if nin==3
  flag = varargin{end};
  if ~iscovflag(flag)
    error(message('GetCoV:notScalarFlag'));
  end   
else
  flag = 0;
end

if size(x,1) ~= size(y,1) 
   error(message('Elements not equal')); 
end

if size(x,1) <= 2
   error(message('More elements required than 2')); 
end

xbar = mean(x);
ybar = mean(y);

result = 0.0;
for i=1:size(x,1)
    result = result + (x(i) - xbar)*(y(i) - ybar);
end

if flag
   result = result/size(x,1);
else
   result = result/(size(x,1)-1);
end