load('ampullary.mat')
sample_rate = 20000; % Hz
max_time = 0;
for i = 1:size(times,2)
    max_time = max([max_time, max(times{i})]);
end
fig = figure();
set(gcf,'Color', 'white')
%% create PSTH on the basis of the interspike intervals
subplot(3,1,1)
hold on
% 1. get the interspike intervals for each trial
for i = 1:size(times,2)
    t = times{i};
    isi = diff(t);
    plot(t(2:end), 1./isi)
end

xlabel('time [s]')
ylabel('firing rate [Hz]')
box('off')
title('instanataneous firing rate')

%% create PSTH using the binning method
subplot(3,1,2)
box('off')
bin_width = 0.02; % s
edges = 0:bin_width:max_time;
firing_rate = [];
for i = 1:size(times,2)
    t = times{i};
    [n, t] = hist(t, edges);
    if isempty(firing_rate)
        firing_rate = n / bin_width;
    else
        firing_rate = firing_rate + (n / bin_width / size(times,2));
    end
end
plot(t,firing_rate)
xlabel('time [s]')
ylabel('firing rate [Hz]')
title('binning method')

%% create PSTH using the kernel-convolution method
subplot(3,1,3)
binary_spikes = zeros(size(times,2), round(max_time*sample_rate));
resps = zeros(size(binary_spikes));
window = hann(bin_width/4*sample_rate,'symmetric');
window = window/sum(window);

for i = 1:size(times,2)
    t = times{i};
    temp = round(t*sample_rate);
    if temp(1) <= 0 
      temp(1) = 1;
    end
    binary_spikes(i, temp) = 1;
    resps(i,:) = conv(binary_spikes(i,:), window, 'same');
end
plot((0:1/sample_rate:max_time), mean(resps,2))