octave_packages/nan-2.5.5/xcovf.m

   1 function [C,N,LAGS] = xcovf(X,Y,MAXLAG,SCALEOPT)
   2 % XCOVF generates cross-covariance function.
   3 % XCOVF is the same as XCORR except
   4 %   X and Y can contain missing values encoded with NaN.
   5 %   NaN's are skipped, NaN do not result in a NaN output.
   6 %   The output gives NaN only if there are insufficient input data
   7 %
   8 % [C,N,LAGS] = xcovf(X,MAXLAG,SCALEOPT);
   9 %      calculates the (auto-)correlation function of X
  10 % [C,N,LAGS] = xcovf(X,Y,MAXLAG,SCALEOPT);
  11 %      calculates the crosscorrelation function between X and Y
  12 %
  13 %  SCALEOPT   [character string] specifies the type of scaling applied
  14 %          to the correlation vector (or matrix). is one of:
  15 %    'none'      return the unscaled correlation, R,
  16 %    'biased'    return the biased average, R/N,
  17 %    'unbiased'  return the unbiassed average, R(k)/(N-|k|),
  18 %    'coeff'     return the correlation coefficient, R/(rms(x).rms(y)),
  19 %          where "k" is the lag, and "N" is the length of X.
  20 %          If omitted, the default value is "none".
  21 %          If Y is supplied but does not have the ame length as X,
  22 %          scale must be "none".
  23 %
  24 %
  25 % see also: COVM, XCORR
  26
  27 %       $Id: xcovf.m 9608 2012-02-10 09:56:25Z schloegl $
  28 %       Copyright (C) 2005,2010,2011 by Alois Schloegl <alois.schloegl@gmail.com>
  29 %       This function is part of the NaN-toolbox
  30 %       http://pub.ist.ac.at/~schloegl/matlab/NaN/
  31
  32 %    This program is free software; you can redistribute it and/or modify
  33 %    it under the terms of the GNU General Public License as published by
  34 %    the Free Software Foundation; either version 3 of the License, or
  35 %    (at your option) any later version.
  36 %
  37 %    This program is distributed in the hope that it will be useful,
  38 %    but WITHOUT ANY WARRANTY; without even the implied warranty of
  39 %    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  40 %    GNU General Public License for more details.
  41 %
  42 %    You should have received a copy of the GNU General Public License
  43 %    along with this program; If not, see <http://www.gnu.org/licenses/>.
  44
  45 if nargin<2,
  46         Y = [];
  47         MAXLAG = [];
  48         SCALEOPT = 'none';
  49 elseif ischar(Y),
  50         MAXLAG = Y;
  51         SCALEOPT=MAXLAG;
  52         Y=[];
  53 elseif all(size(Y)==1),
  54         if nargin<3
  55                 SCALEOPT = 'none';
  56         else
  57                 SCALEOPT = MAXLAG;
  58         end;
  59         MAXLAG = Y;
  60         Y = [];
  61 end;
  62
  63 if 0,
  64
  65 elseif isempty(Y) && isempty(MAXLAG)
  66         NX = isnan(X);
  67         X(NX) = 0;
  68         [C,LAGS] = xcorr(X,'none');
  69         [N,LAGS] = xcorr(1-NX,'none');
  70 elseif ~isempty(Y) && isempty(MAXLAG)
  71         NX = isnan(X);
  72         NY = isnan(Y);
  73         X(NX) = 0;
  74         Y(NY) = 0;
  75         [C,LAGS] = xcorr(X,Y,'none');
  76         [N,LAGS] = xcorr(1-NX,1-NY,'none');
  77 elseif isempty(Y) && ~isempty(MAXLAG)
  78         NX = isnan(X);
  79         X(NX) = 0;
  80         [C,LAGS] = xcorr(X,MAXLAG,'none');
  81         [N,LAGS] = xcorr(1-NX,MAXLAG,'none');
  82 elseif ~isempty(Y) && ~isempty(MAXLAG)
  83         NX = isnan(X);
  84         NY = isnan(Y);
  85         X(NX) = 0;
  86         Y(NY) = 0;
  87         [C,LAGS] = xcorr(X,Y,MAXLAG,'none');
  88         [N,LAGS] = xcorr(1-NX,1-NY,MAXLAG,'none');
  89 end;
  90
  91 if 0,
  92
  93 elseif strcmp(SCALEOPT,'none')
  94         % done
  95
  96 elseif strcmp(SCALEOPT,'coeff')
  97         ix = find(LAGS==0);
  98         if ~any(size(X)==1), %% ~isvector(X)
  99                 c  = C(ix,1:size(X,2)+1:end);   %% diagonal elements
 100                 v  = c.^-0.5; % sqrt(1./c(:));
 101                 v  = v'*v;
 102                 C  = C.*repmat(v(:).',size(C,1),1);
 103         elseif isempty(Y)
 104                 C = C/C(ix);
 105         else
 106                 C = C/sqrt(sumsq(X)*sumsq(Y));
 107         end;
 108
 109 elseif strcmp(SCALEOPT,'biased')
 110         C = C./repmat(max(N),size(C,1),1);
 111
 112 elseif strcmp(SCALEOPT,'unbiased')
 113         C = C./(repmat(max(N),size(C,1),1)-repmat(LAGS,1,size(C,2)));
 114
 115 else
 116         warning('invalid SCALEOPT - not supported');
 117 end;
 118