



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.


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