Home > chronux_1_1 > 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.

   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.

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 %   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 (~isvector(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 (~isvector(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);

Generated on Sun 13-Aug-2006 11:49:44 by m2html © 2003