Add a useful package (from Source forge) for octave

[CreaPhase.git] / octave_packages / nan-2.5.5 / xcovf.m

diff --git a/octave_packages/nan-2.5.5/xcovf.m b/octave_packages/nan-2.5.5/xcovf.m
new file mode 100644 (file)
index 0000000..0cc1512
--- /dev/null
+++ b/octave_packages/nan-2.5.5/xcovf.m
@@ -0,0 +1,118 @@
+function [C,N,LAGS] = xcovf(X,Y,MAXLAG,SCALEOPT)
+% XCOVF generates cross-covariance function. 
+% XCOVF is the same as XCORR except 
+%   X and Y can contain missing values encoded with NaN.
+%   NaN's are skipped, NaN do not result in a NaN output. 
+%   The output gives NaN only if there are insufficient input data
+%
+% [C,N,LAGS] = xcovf(X,MAXLAG,SCALEOPT);
+%      calculates the (auto-)correlation function of X
+% [C,N,LAGS] = xcovf(X,Y,MAXLAG,SCALEOPT);
+%      calculates the crosscorrelation function between X and Y
+%
+%  SCALEOPT   [character string] specifies the type of scaling applied
+%          to the correlation vector (or matrix). is one of:
+%    'none'      return the unscaled correlation, R,
+%    'biased'    return the biased average, R/N, 
+%    'unbiased'  return the unbiassed average, R(k)/(N-|k|), 
+%    'coeff'     return the correlation coefficient, R/(rms(x).rms(y)),
+%          where "k" is the lag, and "N" is the length of X.
+%          If omitted, the default value is "none".
+%          If Y is supplied but does not have the ame length as X,
+%          scale must be "none".
+%
+%
+% see also: COVM, XCORR
+
+%      $Id: xcovf.m 9608 2012-02-10 09:56:25Z schloegl $
+%      Copyright (C) 2005,2010,2011 by Alois Schloegl <alois.schloegl@gmail.com>       
+%       This function is part of the NaN-toolbox
+%       http://pub.ist.ac.at/~schloegl/matlab/NaN/
+
+%    This program is free software; you can redistribute it and/or modify
+%    it under the terms of the GNU General Public License as published by
+%    the Free Software Foundation; either version 3 of the License, or
+%    (at your option) any later version.
+%
+%    This program is distributed in the hope that it will be useful,
+%    but WITHOUT ANY WARRANTY; without even the implied warranty of
+%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+%    GNU General Public License for more details.
+%
+%    You should have received a copy of the GNU General Public License
+%    along with this program; If not, see <http://www.gnu.org/licenses/>.
+
+if nargin<2,
+        Y = [];
+        MAXLAG = [];
+        SCALEOPT = 'none';
+elseif ischar(Y),
+        MAXLAG = Y;
+        SCALEOPT=MAXLAG;
+        Y=[];
+elseif all(size(Y)==1),
+        if nargin<3
+                SCALEOPT = 'none';
+        else
+                SCALEOPT = MAXLAG;
+        end;
+        MAXLAG = Y; 
+        Y = [];
+end;
+
+if 0,
+
+elseif isempty(Y) && isempty(MAXLAG)
+        NX = isnan(X);
+        X(NX) = 0;
+        [C,LAGS] = xcorr(X,'none');
+        [N,LAGS] = xcorr(1-NX,'none');
+elseif ~isempty(Y) && isempty(MAXLAG)
+        NX = isnan(X);
+        NY = isnan(Y);
+        X(NX) = 0;
+        Y(NY) = 0;
+        [C,LAGS] = xcorr(X,Y,'none');
+        [N,LAGS] = xcorr(1-NX,1-NY,'none');
+elseif isempty(Y) && ~isempty(MAXLAG)
+        NX = isnan(X);
+        X(NX) = 0;
+        [C,LAGS] = xcorr(X,MAXLAG,'none');
+        [N,LAGS] = xcorr(1-NX,MAXLAG,'none');
+elseif ~isempty(Y) && ~isempty(MAXLAG)
+        NX = isnan(X);
+        NY = isnan(Y);
+        X(NX) = 0;
+        Y(NY) = 0;
+        [C,LAGS] = xcorr(X,Y,MAXLAG,'none');
+        [N,LAGS] = xcorr(1-NX,1-NY,MAXLAG,'none');
+end;
+
+if 0,
+
+elseif strcmp(SCALEOPT,'none')
+       % done
+
+elseif strcmp(SCALEOPT,'coeff')
+       ix = find(LAGS==0);
+       if ~any(size(X)==1), %% ~isvector(X)    
+               c  = C(ix,1:size(X,2)+1:end);   %% diagonal elements
+               v  = c.^-0.5; % sqrt(1./c(:));
+               v  = v'*v;
+               C  = C.*repmat(v(:).',size(C,1),1);
+       elseif isempty(Y)
+               C = C/C(ix);
+       else 
+               C = C/sqrt(sumsq(X)*sumsq(Y));
+       end;            
+
+elseif strcmp(SCALEOPT,'biased')
+       C = C./repmat(max(N),size(C,1),1);
+
+elseif strcmp(SCALEOPT,'unbiased')
+       C = C./(repmat(max(N),size(C,1),1)-repmat(LAGS,1,size(C,2)));
+
+else
+        warning('invalid SCALEOPT - not supported');
+end;
+