Home > chronux_2_00 > spikesort > utility > datatools > pxcorr.m

pxcorr

PURPOSE ^

PXCORR Efficient cross-correlation for point process data.

SYNOPSIS ^

function [C, lags] = pxcorr(x, varargin)

DESCRIPTION ^

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.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

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);

Generated on Fri 15-Aug-2008 11:35:42 by m2html © 2003