bachelor_thesis/DataAnalysisToolbox/spikeTrigAvg_presyn.m

function spikeTrigAvg_presyn(files, varargin)
%% Calculate and plot the spike triggered average 
% Spike locations for the not stimulated cell are determined within the 
% time frames of the stimuli. Starting at this spike locations the 
% membrane potential of the stimulated cell is gathered backwards for a 
% given length. This data is summed up for each spike and then divided by 
% the number of spikes in total. The result is plotted to see, whether
% there is a distinct behavior in the stimulated cell that is correlated
% with the occurrence of a spike in the not stimulated cell.
% -------------------------------------------------------------------------
% Input     [additional arguments are positional, 
% -----      optional arguments are name-value pairs]
%
%   files:          ID(s) of the datasets in the 'CleanDatasets' directory,
%                       type: integer or vector
%   trials:         number(s) of the trials to be analysed, type: integer
%   (additional)        or vector, default: 'all'
%   Filter:         select between three filter options: 
%   (optional)          'hn'    applys a hum-noise-remove filter
%                       'hnhfn' additonally applys a low pass filter to
%                              remove high frequency noise
%                       If not set, no filter will be applyed (default).
%   Tolerance:      set tolerance in seconds that is added to the stimulus
%   (optional)          time frame on both ends. Default: 0.05 s)
%   MinPeakPromT:   set the value for 'MinPeakProminence' which is used in
%   (optional)          the 'findpeaks' function to find T-cell peaks. 
%                       Default: 17
%   StimCell:       set whether the data for the stimulated T-cell or the 
%   (optional)          stimulated interneuron should be used.
%                       'INT' for the interneuron
%                       'T'   for the T-cell (default)
% ------
% Output
% ------
%   visually output of the plot(s), one per dataset
% -------------------------------------------------------------------------
% Program by: Bjarne Schultze [last modified 01.03.2022]
% -------------------------------------------------------------------------

%% Preperations
% parse the user input
[file_num, trials, filter_choice, tol, MinPeakPromT, stimCell, ...
    ~,~,~] = parse_input_peakfun(files, varargin{:});
% extract the requested data
recData = extract_data(file_num, filter_choice);

% length of a stimulus (the same for each stimulus)
stim_len = 0.5;
% set the sampling rate
samp_rate = 10000;
% set how long the time frame should be to gather the triggering data
trig_len = 0.03;                     % in seconds
trig_len_dp = trig_len * samp_rate;      % in datapoints

% find the points, where the positive stimuli start
[~, stim_locs] = findpeaks(recData{1,1}.stimulus, samp_rate, ...
    'SortStr', 'ascend');

% get the amplitudes for the cells in the stimulus time frame
amplitudes = ampAtStim(file_num);

for file = 1:length(file_num)
    % get the number or trials 
    trial_num = size(recData{file,1}.recordingT,2);
    % distinguish between a stimulated T-cell and a stimulated interneuron
    if isequal(stimCell, "INT")
        % get the columns with data for a stimulated interneuron
        stimulated = setdiff(1:trial_num, recData{file,1}.stimulatedT(:), ...
            'stable');
        % set an appropriate plot title
        plot_title = ...
            sprintf("Spike Triggered Average - File: %02.f - " + ...
            "Interneuron stimulated", file_num(file));
        % assign the data to working variables
        recordingDataNSTIM = recData{file,1}.recordingT(:,stimulated);
        recordingDataSTIM = recData{file,1}.recordingINT(:,stimulated);
        MinPeakPromNSTIM = MinPeakPromT;
        ampSTIM = amplitudes{file,1}.AmplitudeINT(:,stimulated);
    else
        % get the columns with data for a stimulated T-cell
        stimulated = recData{file,1}.stimulatedT(:);
        % set an appropriate plot title
        plot_title = ...
            sprintf("Spike Triggered Average - File: %02.f - T cell " + ...
            "stimulated", file_num(file));
        % assign the data to working variables
        recordingDataSTIM = recData{file,1}.recordingT(:,stimulated);
        recordingDataNSTIM = recData{file,1}.recordingINT(:,stimulated);
        ampSTIM = amplitudes{file,1}.AmplitudeT(:,stimulated);
    end
    if ~isequal(trials,'all')
        % select single trials if requested 
        [~,experiments,~] = intersect(stimulated,trials);
    else
        experiments = 1:length(stimulated);
    end
    
    %% Main calculations 
    % pre-define a variable to sum up the triggering membrane potentials
    trigger_sum = 0;                                         
    spike_counter = 0;
    % iterate through the trials 
    for exp = experiments
        % search for spikes in the not-stimulated cell 
        if isequal(stimCell,'T')
            % use custom function to find interneuron spikes
            [~, spike_locsNSTIM] = findIntSpikes(recordingDataNSTIM(:,exp)); 
            spike_locsNSTIM = spike_locsNSTIM/samp_rate; 
        else
            % use standard findpeaks function for T cell spikes
            [~, spike_locsNSTIM] = findpeaks(recordingDataNSTIM(:,exp), ...
                samp_rate, 'MinPeakProminence', MinPeakPromNSTIM);
        end
        
        for stimulus = 1:length(stim_locs)
            % search for spikes in the stimulus time frame (plus tolerance)
            % Notice: Not using the whole stimulus time frame to not take 
            %         the increase of membrane potentail at the stimulus
            %         onset into account.
            search_indexNSTIM = spike_locsNSTIM >= stim_locs(stimulus) + ...
                tol & spike_locsNSTIM <= stim_locs(stimulus) + stim_len;
            % select the spikes in the stimulus time frame
            spike_locsNSTIM_select = spike_locsNSTIM(search_indexNSTIM);
            
            if isempty(spike_locsNSTIM_select)
                continue
            else
                % iterate through the spikes creating a search box for each
                for spike = 1:length(spike_locsNSTIM_select)
                    % get the index of the spike in the time vector
                    spike_ind = uint32(spike_locsNSTIM_select(spike) * ...
                        samp_rate);
                    % create the index frame of the length of 'trig_len_dp' 
                    % datapoints to select the data (+1 to get the requested 
                    % search frame)
                    search_box = (spike_ind - trig_len_dp + 1):spike_ind;
                    % select the data of the stimulated cell and subtract the
                    % passive response (stimulus+3, because here only the
                    % positive stimuli are examined)
                    trigger_data = recordingDataSTIM(search_box,exp) ...
                        - ampSTIM(stimulus+3,exp);
                    
                    % select the data of the stimulated cell and sum it up
                    trigger_sum = trigger_sum + trigger_data;
                end
            end
        % add the number of spikes to the spike counter
        spike_counter = spike_counter + spike;
        end
    end
    
    % calculate the spike triggered average
    spTrAvg = trigger_sum ./ spike_counter;
    
    % give back the number of spikes used for the calculations
    fprintf("Number of spikes used: %0.f\n\n", spike_counter);
    
    %% Plotting the result
    if spike_counter > 0
        % create a time vector for the trigger data
        timeVector = ...
            sort(-recData{file,1}.timeVector(1:length(trigger_sum)));
        timeVector = timeVector * 1000;     % in ms
        
        % plot the spike triggered average
        figure();
        plot(timeVector, spTrAvg);
        title(plot_title);
        xlabel("time before spike [ms]"); 
        ylabel("average trigger stimulus [mV]")
    else
        fprintf("No spikes were found for the not stimulated cell in " + ...
            "dataset %g!\n", file_num(file));
    end
% end of iteration through the files
end
% end of function definition
end