bachelor_thesis/DataAnalysisToolbox/spikeTrigAvg_postsyn_funtest.m

function spikeTrigAvg_postsyn_funtest(file, varargin)
%% TEST-FUNCTION: Visual verification of the calculations of spikeTrigAvg_postsyn
% In addition to the original function this function plots the data of each
% trial and marks the times for spikes in the other cell as well as the
% data gathered for the calculations. Green dot: start of time frame for
% data selection. Orange line: selected data. Grey line: time of a spike in
% the other cell. 
% Only one file is ananlysed at a time and no trials can be excluded.
% -------------------------------------------------------------------------
% Input     [optional arguments are name-value pairs]
% -----
%
%   file:           number of the '.mat' file in the directory,
%                       type: integer
%   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).
%   MinPeakPromT:   set the value for 'MinPeakProminence' which is used in
%   (optional)          the 'findpeaks' function to find T cell peaks. 
%                       Type: float. Default: 17
%   StimCell:       set whether the data for a stimulated T cell or a 
%   (optional)          stimulated interneuron should be used.
%                       'INT' for the interneuron
%                       'T'   for the T cell (default)
% ------
% Output
% ------
%   visually output of the plots (one plot per trial + one for the result)
% -------------------------------------------------------------------------
% Program by: Bjarne Schultze [last modified 02.03.2022]
% -------------------------------------------------------------------------

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

% select whether or not to plot the rebounce spikes seperately 
% ('Rebounce' for yes)
select_spikes = "Rebounce";

% get the number or trials 
trial_count = size(recData.recordingT,2);
% set how long the time frame should be to gather the triggered data
trig_len = 0.03;                     % in seconds
trig_len_dp = trig_len * 10000;      % in datapoints

% 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_count, recData.stimulatedT(:), 'stable');
    % set an appropriate plot title
    plot_title = ...
        sprintf("Postsynaptic Spike Response - File: %02.f - " + ...
        "Interneuron stimulated", file_num);
    % assign the data to working variables
    recDataNSTIM = recData.recordingT(:,stimulated);
    recDataSTIM = recData.recordingINT(:,stimulated);
else
    % get the columns with data for a stimulated T-cell
    stimulated = recData.stimulatedT(:);
    % set an appropriate plot title
    plot_title = ...
        sprintf("Postsynaptic Spike Response - File: %02.f - " + ...
        "T-cell stimulated", file_num);
    % assign the data to working variables
    recDataSTIM = recData.recordingT(:,stimulated);
    recDataNSTIM = recData.recordingINT(:,stimulated);
    MinPeakPromSTIM = MinPeakPromT;
end

%% Main calculations
% pre-define a variable to sum up the triggered membrane potentials
triggered_sum = 0; rbs_triggered_sum = 0;                                         
spike_counter = 0; rbSpike_counter = 0;
% iterate through the trials
for exp = 1:length(stimulated)
    % search for spikes in the stimulated cell
    if isequal(stimCell, 'INT')
        [~,spike_indSTIM_all] = findIntSpikes(recDataSTIM(:,exp));
    else
        [~, spike_indSTIM_all] = findpeaks(recDataSTIM(:,exp), ...
            'MinPeakProminence', MinPeakPromSTIM);
    end
    % exclude spikes at the very beginning and end to ensure a complete
    % time frame for the triggered data
    spike_indSTIM_sel = spike_indSTIM_all > 50 & ...
        spike_indSTIM_all < (size(recDataSTIM,1) - trig_len_dp);
    spike_indSTIM = spike_indSTIM_all(spike_indSTIM_sel);
    
    % if no spikes are found the next trial is examined
    if isempty(spike_indSTIM)
        continue
    else
        % -----------------------------------------------------------------
        % TEST PART:
        % plot the spike times over the membrane potentail of the 
        % not-stimulated cell
        figure();
        plot(recData.timeVector,recDataNSTIM(:,exp), '-');
        xline((spike_indSTIM./10000), 'Color', [0.3,0.3,0.3]);
        title(plot_title); subtitle(sprintf("Trial: %02.f",stimulated(exp))); 
        xlabel("time [s]"); ylabel("membrane potential [mV]");
        % -----------------------------------------------------------------

        % select the membrane potentail of the not-stimulated cell in 
        % an index frame of the length 'trig_len_dp' datapoints for 
        % each spike of the stimulated cell
        for spike = 1:length(spike_indSTIM)
            % select the data (triggered by a spike)
            triggered_data = ...
                recDataNSTIM((spike_indSTIM(spike)-50):...
                (spike_indSTIM(spike) + trig_len_dp - 1), exp);
            % -------------------------------------------------------------
            % TEST PART:
            % select an apropriate time vector for this data
            triggered_time = recData.timeVector((spike_indSTIM(spike)-50):...
                (spike_indSTIM(spike) + trig_len_dp - 1), 1);
            % add the selected data to the plot
            hold on;
            plot(triggered_time, triggered_data, '-', 'Color', ...
                [0.85,0.33,0.1]);
            plot(triggered_time(1), triggered_data(1), 'o', 'Color',...
                [0.4,0.83,0.07])
            hold off;
            % -------------------------------------------------------------
            % sum up the data
            triggered_sum = triggered_sum + triggered_data;
        end
        % count the spikes
        spike_counter = spike_counter + spike;
        
        % -----------------------------------------------------------------
        % TEST PART:
        % only analyse the data for rebounce spikes
        if isequal(select_spikes, "Rebounce")
            % select the spikes that are not on top of a passive response
            % (rebounce spikes)
            rebSpikes_search = ...
                (spike_indSTIM < 42000 & spike_indSTIM > 35000) | ...
                (spike_indSTIM < 122000 & spike_indSTIM > 115000) | ...
                (spike_indSTIM < 242000 & spike_indSTIM > 235000);
            rebounceSpikes = spike_indSTIM(rebSpikes_search);
             if isempty(rebounceSpikes)
                continue
            else
                % select the membrane potentail of the not-stimulated cell 
                % in an index frame of the length 'trig_len_dp' datapoints 
                % (-1 to get the requested search frame) for each spike of 
                % the stimulated cell
                for rbSpike = 1:length(rebounceSpikes)
                    % select the data (triggered by a spike)
                    rbs_triggered_data = ...
                        recDataNSTIM((rebounceSpikes(rbSpike)-50):...
                        (rebounceSpikes(rbSpike) + trig_len_dp - 1), exp);
                    % select an apropriate time vector for this data
                    rb_triggered_time = ...
                        recData.timeVector((rebounceSpikes(rbSpike)-50):...
                        (rebounceSpikes(rbSpike) + ...
                        trig_len_dp - 1), 1);
                    % add the selected data to the plot
                    hold on;
                    plot(rb_triggered_time, rbs_triggered_data, '-m');
                    hold off; 
                    % sum up the data
                    rbs_triggered_sum = rbs_triggered_sum + ...
                        rbs_triggered_data;
                end
                % count the spikes
                rbSpike_counter = rbSpike_counter + rbSpike;
             end
        end
        % -----------------------------------------------------------------
    end
end

% feedback about the number of spikes used for the calculations
fprintf("\nSpikes in total: %0.f\nthereof rebounce spikes: %0.f\n",...
    spike_counter, rbSpike_counter);

% calculate the averaged postsynaptic membrane potential
afterSpkAvg = triggered_sum ./ spike_counter;
afterRbSpkAvg = rbs_triggered_sum ./ rbSpike_counter;

%% Plotting the result
% create a time vector for the selected data
timeVector = [sort(-recData.timeVector(1:50));...
    recData.timeVector(1:trig_len_dp)];
timeVector = timeVector * 1000;     % in ms

% plot the averaged postsynaptic membrane potential
figure();
plot(timeVector, afterSpkAvg, timeVector, afterRbSpkAvg);
legend("All spikes", "Only rebounce spikes", 'AutoUpdate', 'off');
title(plot_title);
xline(0, 'Color', [0.6,0.6,0.6]);
xlabel("time after spike [ms]"); 
ylabel("averaged membrane potentail [mV]");

% end of function definition
end