function y = ach1(x, ns, sr)
% function y = ach1(x, ns, sr) - calculate autocorrelation histogram 
% from spike time train.  First bin is spike count (as in true ACF).

x = diff(x);		% intervals
mask = x >= 0;		% remove negative intervals (between presentations)
x = x .* mask;

y = zeros(1, ns);
N = size(x, 2) - ns;

% scale to express intervals as number of histogram bins
x  = x .* sr;
n = N;
for i = 1:N

	% List all all-order intervals starting from x(i)
	% To save time, use sub-array of length n (adjusted later in loop)
	bump = i+n;
	if bump > N
		bump = N;
	end
	z = cumsum(x(i:bump));	
	z = ceil(z);
	
	% restrict to histogram range;
	mask = z < ns-1;
	z = z .* mask;
	z = nonzeros(z)' + 1;
	
	% adjust n 
	n = round(size(z, 2) * 1.1);

	% for each interval, increment corresponding bin
	y(z) = y(z) + 1;
end
y(1) = N;