Home > chronux_2_00 > spikesort > utility > datatools > private > CORE_pxcorr.m

CORE_pxcorr   Windows

PURPOSE ^

CORE_PXCORR Core computational routine for PXCORR.

SYNOPSIS ^

function Z = CORE_pxcorr(x,y,maxlag)

DESCRIPTION ^

CORE_PXCORR       Core computational routine for PXCORR.
   Z = CORE_PXCORR(X,Y,MAXLAG), where X and Y are length M and N
   vectors respectively, returns the length (2*MAXLAG+1) vector
   Z such that Z(i) = #(X(j)-Y(k) == i), |i| <= MAXLAG.

   CONDITIONS (including requirements for corresponding MEX file)
   ----------
   X and Y must contain sorted integer values.
   X and Y must be row vectors of type DOUBLE.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function Z = CORE_pxcorr(x,y,maxlag)
0002 %CORE_PXCORR       Core computational routine for PXCORR.
0003 %   Z = CORE_PXCORR(X,Y,MAXLAG), where X and Y are length M and N
0004 %   vectors respectively, returns the length (2*MAXLAG+1) vector
0005 %   Z such that Z(i) = #(X(j)-Y(k) == i), |i| <= MAXLAG.
0006 %
0007 %   CONDITIONS (including requirements for corresponding MEX file)
0008 %   ----------
0009 %   X and Y must contain sorted integer values.
0010 %   X and Y must be row vectors of type DOUBLE.
0011 
0012 Z = zeros(1,2*maxlag+1);
0013 
0014 % We do this efficiently by stepping along X and keeping track of those
0015 % indices in Y that are near by (i.e., within a distance of MAXLAG).
0016 limit = length(x);
0017 a = 1;  c = 1;   
0018 
0019 for b = 1:length(y)
0020     while((y(b)-x(a)) > maxlag),        % move left bracket until less than MAXLAG before x(b)
0021             a=a+1;   if (a > limit), return; end;
0022     end
0023     if (c < a), c = a; end;             % catch up the right bracket if it was passed by
0024     if (c <= limit)
0025         while((x(c)-y(b)) <= maxlag),   % move right bracket until more than MAXLAG after x(b)
0026             c=c+1;   if (c > limit), break; end;
0027         end
0028     end
0029 
0030     offset = -y(b)+maxlag+1;            % add 'em up
0031     for bb = a:(c-1)
0032         ind = x(bb)+offset;
0033         Z(ind) = Z(ind) + 1;
0034     end
0035 end
0036 
0037 
0038 return;
0039 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0040 %%%%%%%%%%%%%%%%%%%%%%%%%%%%% TEST CODE %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0041 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0042 % Fs = 1000; N = 5000;  mu = 0.05;  maxlag = 1;
0043 % intervals = exprnd(mu,2,N) + 2/Fs;     lags = [-maxlag:1/Fs:maxlag];
0044 % events1 = cumsum(intervals(1,:),2);  events1(events1>(mu*N)) = [];  events1 = round(events1*Fs);
0045 % events2 = cumsum(intervals(2,:),2);  events2(events2>(mu*N)) = [];  events2 = round(events2*Fs);
0046 % tic;  cross  = CORE_pxcorr(events1,events2,ceil(maxlag*Fs));  t(1) = toc;
0047 % tic;  cross2 = xcorr(double(rasterize(events1)), double(rasterize(events2)),ceil(maxlag*Fs)); t(2) = toc;
0048 % printf('\nCORE_pxcorr took %5.3f sec and equivalent native code took %5.3f sec.', t(1), t(2));
0049 % printf('The MSE between the two results was %6.4f.', mean((cross-cross2).^2));

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