function [time, rate] = instantaneousRate(spike_times, dt, t_max)
% Function calculates the firing rate as the inverse of the interspike
% interval.
% Arguments:
%            spike_times, vector containing the times of the spikes.
%            dt, the temporal resolutions of the recording.
%            t_max, the duration of the trial.
%
% Returns two vectors containing the  time and the rate.

time = 0:dt:t_max-dt;
rate = zeros(size(time));

isis = diff([0 spike_times]);
inst_rate = 1 ./ isis;
spike_indices = [1 round(spike_times ./ dt)];

for i = 2:length(spike_indices)
    rate(spike_indices(i - 1):spike_indices(i)) = inst_rate(i - 1);
end