XCorr1

XCorr1 - Compute cross-correlograms and their modes for pairs of curves.

  USAGE

    [xc,mode] = XCorr1(x,y,<options>)

    x              first set of curves (MxN matrix: M curves, N bins)
    y              second set of curves
    <options>      optional list of property-value pairs (see table below)

    =========================================================================
     Properties    Values
    -------------------------------------------------------------------------
     'type'        'linear' for linear curves (default), 'circular' otherwise
    =========================================================================

0001 function [xc,mode] = XCorr1(x,y,varargin)
0002 
0003 %XCorr1 - Compute cross-correlograms and their modes for pairs of curves.
0004 %
0005 %  USAGE
0006 %
0007 %    [xc,mode] = XCorr1(x,y,<options>)
0008 %
0009 %    x              first set of curves (MxN matrix: M curves, N bins)
0010 %    y              second set of curves
0011 %    <options>      optional list of property-value pairs (see table below)
0012 %
0013 %    =========================================================================
0014 %     Properties    Values
0015 %    -------------------------------------------------------------------------
0016 %     'type'        'linear' for linear curves (default), 'circular' otherwise
0017 %    =========================================================================
0018 %
0019 
0020 % Copyright (C) 2012 by Michaël Zugaro
0021 %
0022 % This program is free software; you can redistribute it and/or modify
0023 % it under the terms of the GNU General Public License as published by
0024 % the Free Software Foundation; either version 3 of the License, or
0025 % (at your option) any later version.
0026 
0027 % Defaults
0028 type = 'linear';
0029 
0030 % Check number of parameters
0031 if nargin < 2 | mod(length(varargin),2) ~= 0,
0032   error('Incorrect number of parameters (type ''help <a href="matlab:help XCorr1">XCorr1</a>'' for details).');
0033 end
0034 
0035 % Check parameter sizes
0036 if ~isdmatrix(x) || ~isdmatrix(y),
0037     error('Both sets of curves should be MxN matrices (type ''help <a href="matlab:help XCorr1">XCorr1</a>'' for details).');
0038 end
0039 if ~all(size(x)==size(y)),
0040     error('Both sets of curves should have the same sizes (type ''help <a href="matlab:help XCorr1">XCorr1</a>'' for details).');
0041 end
0042 
0043 % Parse parameter list
0044 for i = 1:2:length(varargin),
0045     if ~ischar(varargin{i}),
0046         error(['Parameter ' num2str(i+2) ' is not a property (type ''help <a href="matlab:help XCorr1">XCorr1</a>'' for details).']);
0047     end
0048     switch(lower(varargin{i})),
0049         case 'type',
0050             type = varargin{i+1};
0051             if ~isastring(type,'linear','circular'),
0052                 error('Incorrect value for property ''type'' (type ''help <a href="matlab:help XCorr1">XCorr1</a>'' for details).');
0053             end
0054         otherwise,
0055             error(['Unknown property ''' num2str(varargin{i}) ''' (type ''help <a href="matlab:help XCorr1">XCorr1</a>'' for details).']);
0056     end
0057 end
0058 
0059 [m,n] = size(x);
0060 
0061 if strcmp(type,'linear'),
0062     xc = nan(m,2*n-1);
0063     for i = 1:m,
0064         % Cross-correlagram i
0065         xci = xcorr(x(i,:),y(i,:),'coeff');
0066         % mode of cross-correlogram i (in [-1..1])
0067         mi = round(mean(find(xci==max(xci))));
0068         mi = mi/n-1;
0069         % Store in matrices
0070         xc(i,:) = xci;
0071         mode(i,1) = mi;
0072     end
0073 else
0074     xc = nan(m,2*n-1);
0075     for i = 1:m,
0076         % Cross-correlagram i
0077         xci = xcorr(x(i,:),y(i,:),'coeff');
0078         % mode of cross-correlogram i (in [-1..1])
0079         mi = round(mean(find(xci==max(xci))));
0080         mi = mi/n-1;
0081         % For circular curves, very large modes must be corrected
0082         if abs(mi) > 0.5,
0083             if mi > 0,
0084                 xci = [xci(n+1:end) xci(1:n)];
0085                 mi = mi - 1;
0086             else
0087                 xci = [xci(n:end) xci(1:n-1)];
0088                 mi = mi +1;
0089             end
0090         end
0091         xc(i,:) = xci;
0092         mode(i,1) = mi;
0093     end
0094 end
0095 
0096 
0097 %      xc = nan(m,2*n-1);
0098 %      for i = 1:m,
0099 %          % Cross-correlagram i
0100 %          X = repmat(x(i,:),n,1);
0101 %          Y = repmat(y(i,:),1,n+1);
0102 %          Y(n+1:n+1:end) = [];
0103 %          Y = reshape(Y,n,[]);
0104 %          xci = sum(X.*Y,2)';
0105 %          xci = [xci(2:end) xci];
0106 %          % mode of cross-correlogram i (in [0..1])
0107 %          % (we need to compute mean values of x, but x is circular, so we transform it to angles)
0108 %          a = find(xci==max(xci));
0109 %          a = exp(j*a*2*pi/n);
0110 %          mi = round(wrap(angle(mean(a)),2)/(2*pi)*n)/n-1;
0111 %          if abs(mi) > 0.5,
0112 %              if mi > 0,
0113 %                  mi = mi - 1;
0114 %              else
0115 %                  mi = mi +1;
0116 %              end
0117 %          end
0118 %          % Zero outside primary field
0119 %          m = round((mi+1)*n);
0120 %          h = ceil(n/2);
0121 %          xci(m+h:end) = 0;
0122 %          xci(1:m-h) = 0;
0123 %          % Store in matrices
0124 %          xc(i,:) = xci;
0125 %          mode(i,1) = mi;
0126 %      end
0127 %      % Scale so that autocorrelation at shift 0 is 1
0128 %      x0 = sum(x.^2,2);
0129 %      y0 = sum(y.^2,2);
0130 %      scale = sqrt(x0.*y0);
0131 %      xc = xc./repmat(scale,1,2*n-1);
0132 %
0133 %      end