bachelor_thesis/DataAnalysisToolbox/plot_FI.m

function plot_FI(files, varargin)
%% Create a classic FI curve (with spike count instead of spike frequency)
% Create a plot which shows the spike count of the T cell and the
% interneuron during a stimulation plotted against the stimulation current.
% -------------------------------------------------------------------------
% Input     [additional arguments are positional, 
% -----      optional arguments are name-value pairs]
%
%   files:      ID-number(s) of the '.mat' data-file(s), 
%                   Type: integer or vector
%   trials:     number(s) of the trials to be analysed, type: integer or
%   (additional)    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)
%   ErrorBars:      set whether error bars (standard deviation) 
%   (optional)          are plotted (1 or 0). Default: 1 (yes)
% ------
% Output
% ------
%   visually output of the plot(s)
% -------------------------------------------------------------------------
% Program by: Bjarne Schultze [last modified 02.03.2022]
% -------------------------------------------------------------------------

%% Preperations
% parse input using a function (checking for vaild inputs included)
[file_num, trials, filter_choice, tol, MinPeakPromT, stimCell, ...
    error_bars, ~,~] = parse_input_peakfun(files, varargin{:});

% access the data and applying filters if required
recData = extract_data(file_num, filter_choice);

% iterate through all requested files
for file = 1:length(files)
    % get the number of trials
    if isequal(trials,'all')
        trials = 1:size(recData{file,1}.recordingT,2);   
    end
    % differentiate between a stimulated T-cell or interneuron
    if isequal(stimCell, "INT")
        stimulated = setdiff(trials, recData{file,1}.stimulatedT(:),...
            'stable');
        plot_title = sprintf("File: %02.f - Interneuron stimulated", ...
            file_num(file));
    else
        stimulated = recData{file,1}.stimulatedT(:);
        plot_title = sprintf("File: %02.f - T cell stimulated", ...
            file_num(file));
        if ~isequal(trials,'all')
            stimulated = intersect(stimulated,trials);
        end
    end
    
    %% Calculation of spiking frequencies
    % find the points, where the positive stimuli start
    [stim_peaks, stim_locs] = findpeaks(recData{file,1}.stimulus, 10000, ...
        'SortStr', 'ascend');
    
    % length of a stimulus (the same for each stimulus)
    stim_len = 0.5;
    
    % pre-define vectors for the calculation results
    peak_countT = zeros(length(stim_peaks), length(stimulated));
    peak_countINT = peak_countT;
    
    % iterate through the trails calculating the spiking frequency
    for exp = 1:length(stimulated)
        % determine the locations of spikes in T cell and interneuron
        [~, spike_locsT] = ...
            findpeaks(recData{file,1}.recordingT(:,stimulated(exp)), ...
            10000, 'MinPeakProminence', MinPeakPromT);
        [~, spike_locsINT] = ...
            findIntSpikes(recData{file,1}.recordingINT(:,stimulated(exp)));
        spike_locsINT = spike_locsINT/10000; 
        
        % iterate through the stimuli counting the spikes for both cells
        for stimulus = 1:length(stim_peaks)
            % select the time frame (plus tolerance) where one cell was stimulated 
            search_indexT = spike_locsT >= stim_locs(stimulus)- tol & ...
                            spike_locsT <= stim_locs(stimulus) + ...
                            stim_len + tol;
            % calculate the spiking frequency in the stimulus time frame
            peak_countT(stimulus,exp) = ...
                length(spike_locsT(search_indexT));
            
            % the same calculations for the interneuron
            search_indexINT = spike_locsINT >= stim_locs(stimulus) - ...
                tol & spike_locsINT <= stim_locs(stimulus) + ...
                stim_len + tol;
            peak_countINT(stimulus,exp) = ...
                length(spike_locsINT(search_indexINT));
        end
    end
    
    %% Averaging of the spiking frequencies
    % find the average and standard deviation of the spiking frequency for
    % each stimulus
    m_peak_countT = mean(peak_countT, 2);
    sd_peak_countT = std(peak_countT, 0, 2);
    m_peak_countINT = mean(peak_countINT, 2);
    sd_peak_countINT = std(peak_countINT, 0, 2);
    
    %% Plotting the results
    % creating a figure handle
    figure();
    % Plotting the spiking frequency against the current
    if isequal(error_bars, 1)                       
        % plot with error bars
        errorbar(stim_peaks, m_peak_countT, sd_peak_countT);
        axis padded;
        hold on;
        errorbar(stim_peaks, m_peak_countINT, sd_peak_countINT);
        hold off;
        title(plot_title);
        xlabel("current [nA]"); ylabel("spike count");
        legend("T cell", "Interneuron", 'Location', 'northwest', ...
            'Autoupdate', 'off');
        subtitle("- error bars representing standard deviation " + ...
            "between trials -");
    else                                            
        % plot without error bars
        plot(stim_peaks, m_peak_countT, 'ro-', ...
            stim_peaks, m_peak_countINT, 'bo-','MarkerSize', 4);
        axis padded;
        % axis padded
        legend("T cell", "Interneuron", 'Location','northwest',...
            'AutoUpdate','off');
        title(plot_title);
        xlabel("stimulation current [nA]"); ylabel("spike count");
    end
    % add a referece line at y=0
    yline(0, 'Color', [0.6 0.6 0.6])

% end of file iteration
end
% end of function definition
end