function [PSTH, bin, ST] = psth(T, stim_times, bounds,dt)
%PSTH Peri-stimulus spiking histogram 

% Initialization
bin = bounds(1):dt:bounds(2);
PSTH = zeros(size(bin,2),1);

% Get bins for each stimulus
bins = bin(ones(size(stim_times)),:);
bins = bins + stim_times(:,ones(size(bin)));

ST = cell(numel(stim_times),1);

% For all spikes occuring in at times T
for i = 1:numel(stim_times)

    % Current window
    PSTH = PSTH + histc(T, bins(i,:)); 
    
    % Gather relative spike times
    ST{i} = T(T >=  bins(i,1) & T <=  bins(i,end)) - bins(i,1) +  bounds(1);

end

% Normalize by number of stimuli and time bin with to get rate
PSTH = PSTH./numel(stim_times)./dt;

end

