load('ampullary.mat')
load('electroreceptor_stimulus.mat')
sample_rate = 20000; % Hz
max_time = 10;

%% create instantaneous firing rate on the basis of the interspike intervals
t = times{1};
firing_rate = [0 1./diff(t)];
start = 1;
response = zeros(1, round(max_time * sample_rate));
for i = 1:length(t)
    response(1,start:round(t(i) * sample_rate)) = firing_rate(i);
    start = round(t(i) * sample_rate);
end
fig = figure();
set(gcf, 'PaperUnits', 'centimeters');
set(gcf, 'PaperSize', [11.7 9.0]);
set(gcf, 'PaperPosition',[0.0 0.0 11.7 9.0]);
set(gcf,'Color', 'white')
plot((1/sample_rate:1/sample_rate:max_time), response)
xlabel('time [s]')
ylabel('firing rate [Hz]')
ylim([0 300])
xlim([0 1])
title('instanataneous firing rate')
saveas(fig, 'isi.pdf','pdf')


%% create PSTH using the binning method
bin_width = 0.0125; % s
edges = 0:bin_width:max_time;
firing_rate = [];
for i = 1:size(times,2)
    t = times{i};
    [n, time] = hist(t, edges);
    if isempty(firing_rate)
        firing_rate = n / bin_width / size(times,2);
    else
        firing_rate = firing_rate + (n / bin_width / size(times,2));
    end
end

response = zeros(1, round(max_time * sample_rate));
start_index = 1;
for i = 1:length(edges)
    end_index = round(edges(i) * sample_rate);
    response(start_index:end_index) = firing_rate(i);
    start_index = end_index +1; 
end

fig = figure();
set(gcf, 'PaperUnits', 'centimeters');
set(gcf, 'PaperSize', [11.7 9.0]);
set(gcf, 'PaperPosition',[0.0 0.0 11.7 9.0]);
set(gcf,'Color', 'white')
plot(0:1/sample_rate:max_time-1/sample_rate, response)
ylim([0 300])
xlim([0 1])
xlabel('time [s]')
ylabel('firing rate [Hz]')
title('binning method')
saveas(fig, 'binning.pdf', 'pdf')


%% create PSTH using the kernel-convolution method
kernel_width = 0.0125; %s
binary_spikes = zeros(size(times,2), round(max_time*sample_rate));
responses = zeros(size(binary_spikes));
window = hann(kernel_width*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;
    responses(i,:) = conv(binary_spikes(i,:), window, 'same')*sample_rate;
end
fig = figure();
set(gcf, 'PaperUnits', 'centimeters');
set(gcf, 'PaperSize', [11.7 9.0]);
set(gcf, 'PaperPosition',[0.0 0.0 11.7 9.0]);
set(gcf,'Color', 'white')

subplot(2,1,1)
plot((1/sample_rate:1/sample_rate:max_time), mean(responses,1))
ylim([0 300])
xlim([0 1])
xlabel('time [s]')
ylabel('firing rate [Hz]')
title('convolution method')

subplot(2,1,2)
plot(stimulus_strong(:,1), stimulus_strong(:,2))
ylim([-1 1])
xlim([0 1])
xlabel('time [s]')
ylabel('stimulus intensity [arb. units]')