function [time, rate] = convolution_rate(spike_times, sigma, dt, t_max)
% Calculates the firing rate with the convolution method. 
% Arguments: 
%           spike_times, a vector containing the spike times.
%           sigma, the standard deviation of the Gaussian kernel in seconds.
%           dt, the temporal resolution in seconds.
%           t_max, the trial duration in seconds.
%
% 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