function [counts, bins] = counthist(spikes, w)
% computes count histogram and compare with Poisson distribution
%
% [counts, bins] = counthist(spikes, w)
%
% Arguments:
%   spikes: a cell array of vectors of spike times in seconds
%   w: observation window duration in seconds for computing the counts
%
% Returns:
%   counts: the histogram of counts normalized to probabilities
%   bins: the bin centers for the histogram

    % collect spike counts:
    tmax = spikes{1}(end);
    n = [];
    r = [];
    for k = 1:length(spikes)
        times = spikes{k};
%       alternative 1: count the number of spikes in each window:
%         for tk = 0:w:tmax-w
%             nn = sum( ( times >= tk ) & ( times < tk+w ) );
%             %nn = length( find( ( times >= tk ) & ( times < tk+w ) ) );
%             n = [ n nn ];
%         end
%     alternative 2: use the hist function to do that!
        tbins = 0.5*w:w:tmax-0.5*w;
        nn = hist(times, tbins);
        n = [ n nn ];
        % the rate of the spikes:
        rate = (length(times)-1)/(times(end) - times(1));
        r = [ r rate ];
    end

    % histogram of spike counts:
    maxn = max( n );
    [counts, bins ] = hist( n, 0:1:maxn+10 );
    % normalize to probabilities:
    counts = counts / sum( counts );

    % plot:
    if nargout == 0
        bar( bins, counts );
        hold on;
        % Poisson distribution:
        rate = mean( r );
        x = 0:1:maxn+10;
        a = rate*w;
        y = a.^x.*exp(-a)./factorial(x);
        plot( x, y, 'r', 'LineWidth', 3 );
        hold off;
        xlabel( 'counts k' );
        ylabel( 'P(k)' );
    end
end