function [pdf, centers] = isihist(isis, binwidth)
% Compute normalized histogram of interspike intervals.
%
% [pdf, centers] = isihist(isis, binwidth)
%
% Arguments:
%   isis: vector of interspike intervals in seconds
%   binwidth: optional width in seconds to be used for the isi bins
%
% Returns:
%   pdf: vector with probability density of interspike intervals in Hz
%   centers: vector with centers of interspikeintervalls in seconds

    if nargin < 2
        % compute good binwidth:
        nperbin = 200;     % average number of data points per bin
        bins = length(isis)/nperbin;   % number of bins
        binwidth = max(isis)/bins;
        if binwidth < 5e-4             % half a millisecond
            binwidth = 5e-4;
        end
    end
    bins = 0.5*binwidth:binwidth:max(isis);
    % histogram data:
    [nelements, centers] = hist(isis, bins);
    % normalization (integral = 1):
    pdf = nelements / sum(nelements) / binwidth;
end