gistfile1.matlab

function [tsOffsets, instantRate] = calcCrossCorrelogram(ts1, ts2, window)
    % INPUT:
    % ts1 - timestamp array, center of cross-correlogram
    % ts2 - timestamp array
    % window - the window of the timestamps to compute, e.g. [-0.1 0.1] for
    %           0.1 seconds around each spike in ts1.
    %
    % OUTPUT:
    % tsOffsets - the offsets from each spike in ts1 that has spikes nearby
    % instantRate - 1/ISI of each spike in ts1 that has spikes nearby
    % 
    % If you want to make a cross-correlogram histogram, run the following:
    % hist(tsOffsets, 50); 
    
    startWindow = 1;
    spikeIdx = 2;
    ccInfo = zeros(length(ts1), 4)+NaN; % startIdx, endIdx, spikeTime, 1/ISI
    while spikeIdx <= length(ts1)

        % Seek to the beginning of the current spike's window
        i = startWindow;
        while ts2(i) <= (ts1(spikeIdx) + window(1)) && i < length(ts2)
            i = i+1;
        end;
        if ts2(i) > ( ts1(spikeIdx) + window(2))
            spikeIdx = spikeIdx+1;
            continue;
        end
        startWindow = i; % save the location for later
        
        % Find all the spike indices that fall within the window
        while ts2(i) <= (ts1(spikeIdx) + window(2)) && i < length(ts2)
            i = i+1;
        end
        endWindow = i-1;
        
        ccInfo(spikeIdx,1) = startWindow;
        ccInfo(spikeIdx,2) = endWindow;
        ccInfo(spikeIdx,3) = ts1(spikeIdx);
        ccInfo(spikeIdx,4) = 1./(ts1(spikeIdx) - ts1(spikeIdx-1));
        spikeIdx = spikeIdx+1;
    end
    
    ccInfo = ccInfo(~isnan(ccInfo(:,1)),:);
    
    % Now I've got all of the indices into spikes, and their offset
    % from the center-spike is easily calculable.
    diffArray = diff(ccInfo (:,1:2),1,2); 
    done = 0; tsOffsets = []; instantRate = [];
    incr = 0;
    
    while ~done
        tmp = find(diffArray >= incr);
        if isempty(tmp); done = 1; end
        idx = ccInfo(tmp,1)+incr; % brilliant!
        centerTimes = ccInfo(tmp,3);
        
        tsOffsets = [tsOffsets ts2(idx)-centerTimes'];
        instantRate = [instantRate ccInfo(tmp,4)'];
        incr = incr+1;
    end
 
    
end