


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.
PXCORR(A,Fs) returns the auto-correlation of the events listed in A
and is equivalent to PXCORR(A,A,Fs)
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 % PXCORR(A,Fs) returns the auto-correlation of the events listed in A 0021 % and is equivalent to PXCORR(A,A,Fs) 0022 % 0023 % 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 %%%%%%%%%%%%%%%%%%%%%%%%%% Argument Parsing %%%%%%%%%%%%%%%%%%%%%%%%% 0031 if (nargin<2), error('Not enough arguments.'); end; 0032 if (length(x)<=1), error('Point processes must consist of >1 event.'); end; 0033 if (~isavector(x)), error('Matrix data are not supported.'); end; 0034 0035 presorted = 1; % already sorted events? 0036 if (strcmpi(varargin{end},'sort')), 0037 presorted = 0; 0038 varargin = varargin(1:end-1); 0039 if(length(varargin)==0), error('Not enough arguments.'); end; 0040 end; 0041 0042 if (length(varargin{1}) == 1), y = x; % make auto-corr case look like X 0043 else, 0044 y = varargin{1}; 0045 varargin = varargin(2:end); 0046 if(length(varargin)==0), error('Not enough arguments.'); end; 0047 end 0048 if (length(y)<=1), error('Point processes must consist of >1 event.'); end; 0049 if (~isavector(x)), error('Matrix data are not supported.'); end; 0050 0051 x = x(:)'; y = y(:)'; % enforce row vectors 0052 if (~presorted), x = sort(x); y = sort(y); end; % only do this if needed 0053 0054 nargs = length(varargin); 0055 Fs = varargin{1}; 0056 if(length(Fs)~=1), error('Fs must be a scalar.'); end; 0057 0058 if (nargs == 1) % limiting the lag saves time 0059 maxlag = max(x(end)-y(1),y(end)-x(1)); 0060 elseif (nargs == 2) 0061 maxlag = varargin{2}; 0062 if(length(maxlag)~=1 || isinf(maxlag)), error('MAXLAG must be a finite scalar.'); end; 0063 elseif ((nargs == 3) && (length(varargin{3})<=1)) 0064 error('Point processes must consist of > 1 event.'); 0065 else 0066 error('Invalid syntax.'); 0067 end 0068 0069 x = round(x*Fs); y = round(y*Fs); 0070 maxlag = ceil(maxlag * Fs); 0071 0072 if (nargout == 2), lags = [-maxlag:maxlag]./Fs; end 0073 0074 0075 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Correlate %%%%%%%%%%%%%%%%%%%%%%%%%%%% 0076 C = CORE_pxcorr(x,y,maxlag);