% Estimate the ISI density of a spike train P
% dt is the binsize (in whatiever units you're using)
% of the input process P.
% n is the number of bins to be used in the density.
% The density is returned in H and the domain
% of the density in X so that the jth element
% of H represents the density at the jth element of 
% X.  So that you might make the following call:
%        [h,x]=isi_density(p,100,.1);
%        plot(x,h);
% To get the density and plot it with the axes 
% correctly marked.
% NOTE: If you would rather pass a vector of bin centers for
% instead of the number of bins, simply pass it in n.  See
% the help documentation for hist.

function [H,X]=isi_density(P,n,dt)

% Get a representation of P
% as a vector of spike times
s=find(P);


% Initialize the vector that will
% hold the interspike intervals.
% Note that there is one less
% ISI than there are spikes
isis=zeros(1,numel(s)-1);



% Fill the isi vector
% Note that it is in units of 1/dt.
for i=1:numel(s)-1
    isis(i)=s(i+1)-s(i);
end

% Put isis into the correct units
isis=isis*dt;

% Create a histogram that is normalized so as
% to sum to 1 (by dividing by numel(s)).
% Note that, by multiplying by dt, we
% get the correct units.

% Use Matlab's built in function to get 
% an ISI histogram
[H,X]=hist(isis,n);

% Normalize the histogram so that it sums to 1
H=H/(numel(s)*dt);