function [time, rate] = convolution_rate(spikes, sigma, dt, t_max)
% PSTH computed with convolution method. 
%
% [time, rate] = convolution_rate(spikes, sigma, dt, t_max)
%
% Arguments: 
%   spikes: 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(spikes / dt);
rate(spike_indices) = 1;
kernel = gaussKernel(sigma, dt);

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


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