function [time, rate] = convolution_rate(spike_times, sigma, dt, t_max)
% Calculates the firing rate with the convolution method. Function takes
% the spike times, the width of the Gaussian kernel (sigma, i.e. the
% standard deviation), the temporal resolution, and the the duration of 
% the trial. 
%
% Returns two vectors containing the time and the rate.

time = 0:dt:t_max-dt;
rate = zeros(size(time));
spike_indices = round(spike_times/dt);
rate(spike_indices) = 1;
kernel = gauss_kernel(sigma, dt);

rate = conv(rate, kernel, 'same');
end


function y = gauss_kernel(s, step)
    x = -4*s:step:4*s;
    y = exp(-0.5*(x./s).^2)/sqrt(2*pi)/s;
end