


PXCORR Efficient cross-correlation for point process data.
C = PXCORR(A,B,Fs) returns the cross-correlation between A and B at
sampling rate Fs for the point process whose times are given in the
length N(>1) vector A and the length M (>1) vector B. C will be a row
vector with length given by
2*Fs*max(max(A)-min(B),max(B)-min(A))+1
i.e., twice the maximum time difference between events in A and B
measured in units of 1/Fs (plus 1 for 0 lag). The vector C estimates
C(t) = E[A(x-t)*B(x)]
The event times listed in A and B are assumed to be sorted and will be
rounded to bins of duration 1/Fs.
The computation here is O(N*J), where 0<=J<=M is the mean number of
events in the vector B that are within MAXLAG (see below) of events in
the vector A (when MAXLAG is not specified, J=M). Matlab's XCORR is
O(K log K), where K = max(max(A),max(B))*Fs -- note that if MAXLAG is
large and the mean interevent interval is small, XCORR can be faster.
C = PXCORR(A,Fs) returns the auto-correlation of the events listed in A
and is equivalent to PXCORR(A,A,Fs)
C = PXCORR(...,MAXLAG) only considers lags upto MAXLAG; the length of C
will then be 2*Fs*MAXLAG + 1.
PXCORR(...,'sort') indicates the inputs are not yet sorted.
[C,LAGS] = PXCORR(...) also returns a vector of lag indices.

0001 function [C, lags] = pxcorr(x, varargin) 0002 %PXCORR Efficient cross-correlation for point process data. 0003 % C = PXCORR(A,B,Fs) returns the cross-correlation between A and B at 0004 % sampling rate Fs for the point process whose times are given in the 0005 % length N(>1) vector A and the length M (>1) vector B. C will be a row 0006 % vector with length given by 0007 % 2*Fs*max(max(A)-min(B),max(B)-min(A))+1 0008 % i.e., twice the maximum time difference between events in A and B 0009 % measured in units of 1/Fs (plus 1 for 0 lag). The vector C estimates 0010 % C(t) = E[A(x-t)*B(x)] 0011 % The event times listed in A and B are assumed to be sorted and will be 0012 % rounded to bins of duration 1/Fs. 0013 % 0014 % The computation here is O(N*J), where 0<=J<=M is the mean number of 0015 % events in the vector B that are within MAXLAG (see below) of events in 0016 % the vector A (when MAXLAG is not specified, J=M). Matlab's XCORR is 0017 % O(K log K), where K = max(max(A),max(B))*Fs -- note that if MAXLAG is 0018 % large and the mean interevent interval is small, XCORR can be faster. 0019 % 0020 % C = PXCORR(A,Fs) returns the auto-correlation of the events listed in A 0021 % and is equivalent to PXCORR(A,A,Fs) 0022 % 0023 % C = PXCORR(...,MAXLAG) only considers lags upto MAXLAG; the length of C 0024 % will then be 2*Fs*MAXLAG + 1. 0025 % 0026 % PXCORR(...,'sort') indicates the inputs are not yet sorted. 0027 % 0028 % [C,LAGS] = PXCORR(...) also returns a vector of lag indices. 0029 0030 0031 %%%%%%%%%%%%%%%%%%%%%%%%%% Argument Parsing %%%%%%%%%%%%%%%%%%%%%%%%% 0032 if (nargin<2), error('Not enough arguments.'); end; 0033 if (length(x)<=1), error('Point processes must consist of >1 event.'); end; 0034 if (~isvectord(x)), error('Matrix data are not supported.'); end; 0035 0036 presorted = 1; % already sorted events? 0037 if (strcmpi(varargin{end},'sort')), 0038 presorted = 0; 0039 varargin = varargin(1:end-1); 0040 if(isempty(varargin)), error('Not enough arguments.'); end; 0041 end; 0042 0043 if (length(varargin{1}) == 1), y = x; % make auto-corr case look like X 0044 else 0045 y = varargin{1}; 0046 varargin = varargin(2:end); 0047 if(isempty(varargin)), error('Not enough arguments.'); end; 0048 end 0049 if (length(y)<=1), error('Point processes must consist of >1 event.'); end; 0050 if (~isvectord(x)), error('Matrix data are not supported.'); end; 0051 0052 x = x(:)'; y = y(:)'; % enforce row vectors 0053 if (~presorted), x = sort(x); y = sort(y); end; % only do this if needed 0054 0055 nargs = length(varargin); 0056 Fs = varargin{1}; 0057 if(length(Fs)~=1), error('Fs must be a scalar.'); end; 0058 0059 if (nargs == 1) % limiting the lag saves time 0060 maxlag = max(x(end)-y(1),y(end)-x(1)); 0061 elseif (nargs == 2) 0062 maxlag = varargin{2}; 0063 if(length(maxlag)~=1 || isinf(maxlag)), error('MAXLAG must be a finite scalar.'); end; 0064 elseif ((nargs == 3) && (length(varargin{3})<=1)) 0065 error('Point processes must consist of > 1 event.'); 0066 else 0067 error('Invalid syntax.'); 0068 end 0069 0070 x = round(x*Fs); y = round(y*Fs); 0071 maxlag = ceil(maxlag * Fs); 0072 0073 if (nargout == 2), lags = (-maxlag:maxlag)./Fs; end 0074 0075 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Correlate %%%%%%%%%%%%%%%%%%%%%%%%%%%% 0076 C = CORE_pxcorr(x,y,maxlag);