% AUTOCORR calculates the autocorrelation function of a spike train
% (taken from Dayan and Abbott).
function [Q,Q_norm,s] = autocorr(sX,ds,smin,smax)

% INPUT 
% sX: spiketrain 
% ds: bin width of the autocorrelogram
% smin, smax: minimum and maximum of the range in which the autocorr is
%             calculated
%
% OUTPUT  
% Q:  autocorrelogram (raw)
% Q_norm:  autocorrelogram (normalized)


n = length(sX);
T = sX(n) - sX(1);

k = 0;
% create a list of all interspike intervals within [smin,smax]
for i = 1:n
    for j = 1:n
      sij = sX(i)-sX(j);
%     if ((sij>smin) && (sij <smax))  % in octave
      if ((sij>smin) & (sij <smax))  % in matlab
        k = k+1;
        Q1(k) = sij;
      end
    end
end

s = smin-ds/2:ds:smax+ds/2;

% Q = histc(Q1,s);     % matlab            (bin edges   specified by 's')
[Q,s] = hist(Q1,s); % octave and matlab (bin centers specified by 's')

Q_norm = Q/T - n*n*ds/(T*T);