clear 
clc

%% load data
load p-unit_spike_times.mat
load p-unit_stimulus.mat
sample_rate = 20000;

%% calculate STA
all_times = [];
for i = 1:length(spike_times)
    all_times = cat(1, all_times, spike_times{i});
end

[sta, sta_sd, num] = sta(stimulus_strong(:,2), all_times, 1000, sample_rate);
fig = figure();
set(fig, 'PaperUnits', 'centimeters');
set(fig, 'PaperSize', [11.7 9.0]);
set(fig, 'PaperPosition',[0.0 0.0 11.7 9.0]);
set(fig,'Color', 'white')
plot(((1:length(sta))-1000)/sample_rate, sta)
xlabel('time [s]')
ylabel('stimulus')

%% reverse reconstruction

% make binary representation of the spike times
binary_spikes = zeros(size(stimulus_strong, 1), length(spike_times));
estimated_stims = zeros(size(binary_spikes));
for i = 1:length(spike_times)
    binary_spikes(round(spike_times{i}*sample_rate), i) = 1;
    estimated_stims(:,i) = conv(binary_spikes(:,i), sta, 'same');
end

fig = figure();
set(fig, 'PaperUnits', 'centimeters');
set(fig, 'PaperSize', [11.7 9.0]);
set(fig, 'PaperPosition',[0.0 0.0 11.7 9.0]);
set(fig,'Color', 'white')
plot(stimulus_strong(:,1), stimulus_strong(:,2), 'displayname','original')
hold on
plot(stimulus_strong(:,1), mean(estimated_stims,2), 'r', 'displayname', 'reconstruction')
xlabel('time [s]')
ylabel('stimulus')
legend show

%% calculate STC

% we need to downsample the data otherwise the covariance matrixs gets too
% large 20Khz to 1kHz

downsampled_binary = zeros(size(stimulus_strong, 1)/20, length(spike_times));
downsampled_stim = zeros(size(downsampled_binary,1),1);

for i = 1:length(spike_times)
    binary_spikes(round(spike_times{i}*1000), i) = 1;
end
for i = 1:length(downsampled_stim)
    start_index = (i-1) * 1000 + 1;
    downsampled_stim(i) = mean(stimulus_strong()*20))
end