nanxcov

 NANXCOV   Cross-covariance function estimates. Skips NaNs.
   XCOV(A,B), where A and B are length M vectors, returns the
   length 2*M-1 cross-covariance sequence in a column vector.

   XCOV(A), when A is a vector, is the auto-covariance sequence.
   XCOV(A), when A is an M-by-N matrix, is a large matrix with
   2*M-1 rows whose N^2 columns contain the cross-covariance
   sequences for all combinations of the columns of A.
   The zeroth lag of the output covariance is in the middle of the
   sequence, at element or row M.

   The cross-covariance is the cross-correlation function of
   two sequences with their means removed:
        C(m) = E[(A(n+m)-MA)*conj(B(n)-MB)]
   where MA and MB are the means of A and B respectively.

   XCOV(...,MAXLAG) computes the (auto/cross) covariance over the
   range of lags:  -MAXLAG to MAXLAG, i.e., 2*MAXLAG+1 lags.
   If missing, default is MAXLAG = M-1.

   [C,LAGS] = XCOV(...) returns a vector of lag indices (LAGS).

   XCOV(...,SCALEOPT), normalizes the covariance according to SCALEOPT:
       biased   - scales the raw cross-covariance by 1/M.
       unbiased - scales the raw covariance by 1/(M-abs(k)), where k
                  is the index into the result.
       coeff    - normalizes the sequence so that the covariances at
                  zero lag are identically 1.0.
       none     - no scaling (this is the default).

   See also XCORR, CORRCOEF, CONV, COV and XCORR2.

0001 function [xycov,lags,nanp] = nanxcov(x,y,option1,option2)
0002 
0003 % NANXCOV   Cross-covariance function estimates. Skips NaNs.
0004 %   XCOV(A,B), where A and B are length M vectors, returns the
0005 %   length 2*M-1 cross-covariance sequence in a column vector.
0006 %
0007 %   XCOV(A), when A is a vector, is the auto-covariance sequence.
0008 %   XCOV(A), when A is an M-by-N matrix, is a large matrix with
0009 %   2*M-1 rows whose N^2 columns contain the cross-covariance
0010 %   sequences for all combinations of the columns of A.
0011 %   The zeroth lag of the output covariance is in the middle of the
0012 %   sequence, at element or row M.
0013 %
0014 %   The cross-covariance is the cross-correlation function of
0015 %   two sequences with their means removed:
0016 %        C(m) = E[(A(n+m)-MA)*conj(B(n)-MB)]
0017 %   where MA and MB are the means of A and B respectively.
0018 %
0019 %   XCOV(...,MAXLAG) computes the (auto/cross) covariance over the
0020 %   range of lags:  -MAXLAG to MAXLAG, i.e., 2*MAXLAG+1 lags.
0021 %   If missing, default is MAXLAG = M-1.
0022 %
0023 %   [C,LAGS] = XCOV(...) returns a vector of lag indices (LAGS).
0024 %
0025 %   XCOV(...,SCALEOPT), normalizes the covariance according to SCALEOPT:
0026 %       biased   - scales the raw cross-covariance by 1/M.
0027 %       unbiased - scales the raw covariance by 1/(M-abs(k)), where k
0028 %                  is the index into the result.
0029 %       coeff    - normalizes the sequence so that the covariances at
0030 %                  zero lag are identically 1.0.
0031 %       none     - no scaling (this is the default).
0032 %
0033 %   See also XCORR, CORRCOEF, CONV, COV and XCORR2.
0034 
0035 %   Author(s): L. Shure, 1-9-88
0036 
0037 
0038 
0039 % skips NaNs, P. Sturm, 10-Nov-2008
0040 % outputs isnan vector of skipped NaNs, P. Sturm, 9-Nov-2009
0041 
0042 %   References:
0043 %     [1] J.S. Bendat and A.G. Piersol, "Random Data:
0044 %         Analysis and Measurement Procedures", John Wiley
0045 %         and Sons, 1971, p.332.
0046 %     [2] A.V. Oppenheim and R.W. Schafer, Digital Signal
0047 %         Processing, Prentice-Hall, 1975, pg 539.
0048 
0049 nanx = isnan(x);
0050 nany = isnan(y);
0051 if nargin == 1
0052     x(nanx) = [];
0053     mx = size(x, 1);
0054     [xycov,l] = xcorr(x-ones(mx,1)*mean(x));
0055     nanp = nanx;
0056 elseif nargin == 2
0057     x(nanx|nany) = [];
0058     y(nanx|nany) = [];
0059     mx = size(x, 1);
0060     if ischar(y)||(~ischar(y)&&length(y)==1)
0061         [xycov,l] = xcorr(x-ones(mx,1)*mean(x),y);
0062     else
0063         my = size(y, 1);
0064         [xycov,l] = xcorr(x-ones(mx,1)*mean(x),y-ones(my,1)*mean(y));
0065     end
0066     nanp = (nanx|nany);
0067 elseif nargin == 3
0068     x(nanx|nany) = [];
0069     y(nanx|nany) = [];
0070     mx = size(x, 1);
0071     my = size(y, 1);
0072     if length(y)==1
0073         [xycov,l] = xcorr(x-ones(mx,1)*mean(x),y,option1);
0074     else
0075         [xycov,l] = xcorr(x-ones(mx,1)*mean(x),y-ones(my,1)*mean(y),option1);
0076     end
0077     nanp = (nanx|nany);
0078 elseif nargin == 4
0079     x(nanx|nany) = [];
0080     y(nanx|nany) = [];
0081     mx = size(x, 1);
0082     my = size(y, 1);
0083     [xycov,l] = xcorr(x-ones(mx,1)*mean(x),y-ones(my,1)*mean(y),...
0084                          option1,option2);
0085     nanp = (nanx|nany);
0086 end
0087 if nargout > 1
0088     lags = l;
0089 end