%% Using the spike-triggered average to characterize a neural receptive field

close all;
clearvars;
load('c1p8.mat');

%% parameters

Fs = 500;                       % sampling rate in Hertz
dt = 1/Fs * 1000;               % sampling period in msec, this is the answer to Q2

Hours_recording = round(numel(rho)/Fs/60); % this is the answer to Q3

STA_length = 300;               % spike-triggered average (STA) length in msec

STA_samples = STA_length/dt;    % samples for STA, this is the answer to Q4

%% retrieve spikes that occur after 300 msec

Spike_indices = find(rho);      % find the indices of all of the spikes
% select spikes that occur after 300 msec into the recording
Spike_indices = Spike_indices(Spike_indices>STA_samples);
% number of valid spikes (i.e after 300 msec into the recording)
Num_valid_spikes = numel(Spike_indices); % this is the answer to Q5 
% we exclude spikes occuring before 300 msec into the recording as for those 
% spikes we do not have stimulus spanning 300 msec before the spike to average  

%% compute the STA

STA = zeros(STA_samples,1);     % create vector to store STA

    for idx = 1 : Num_valid_spikes  % loop over all spikes 
        % add corresponding stimulus vector
        STA = STA + stim(Spike_indices(idx)-STA_samples+1:Spike_indices(idx)); 
    end
    
% divide by the number of valid spikes to get the average stimulus
STA = STA / Num_valid_spikes;

%% plot the STA

t_vector = dt*(-STA_samples+1:0);   % vector of time samples in msec

figure(1);
plot(t_vector, STA, 'LineWidth',2);
xlabel('Time before a spike (ms)');
ylabel('Bar Velocity (arbitrary units)');
title('Spike-Triggered Average of H1 neuron');
prettyfigure;

% Q6: The neuron prefers a sudden increase in pattern velocity and then a slow
% (exponential like) decay. The STA looks like the response of a leaky
% integrate and fire neuron. Therefore we can think that it integrates the
% stimulus and then forgets about it with a time constant of roughly 20 msec.
