bachelor_thesis/DataAnalysisToolbox/intSpikeHeight_trend.m

function intSpikeHeight_trend(files, varargin)
%% Function to measure the amplitude of interneuron spikes
% This function is written to detect spikes in interneuron recordings and
% obtain their mean amplitude. A three step mechanism is used to 
% distinguish the spikes from noise or other rapid changes in the membrane 
% potential. If requested a plot is generated showing the selected peaks
% (set the plotflag to 1). The function can handle multiple files at a
% time.
% -------------------------------------------------------------------------
% Input     [additional arguments are positional, 
% -----      optional arguments are name-value pairs]
%
%   file:           ID(s) of the file(s) in the 'CleanDatasets' folder 
%                       type: integer or vector
%   trials:         select which trials should be used, default: 'all'
%   (additional)        type: integer or vector
%   PlotFlag:       select if a plot should be created for each trial
%   (optional)          showing the detected peaks. 0 (no plots, default), 
%                       1 (plots)
% ------
% Output
% ------
%   varargout{1}:   measured average spike amplitude
%   varargout{2}:   standard deviation across all spikes of one cell
% -------------------------------------------------------------------------
% Program by:   Bjarne Schultze     [last modified: 20.01.2022]
% -------------------------------------------------------------------------

%% Preparations
% parse the user input
[fileNums, trials, ~,~,~,~,~,~, ~] = ...
    parse_input_peakfun(files, varargin{:});

% access the data using the function 'extract_data', filtering hum noise
recData = extract_data(fileNums,'hn');

%% Main calculations
% iterate through the requested files
for file = 1:length(fileNums)
    % get the number of trials where the interneuron was stimulated
    trialNum = size(recData{file,1}.recordingINT,2);
    stimulated = setdiff(1:trialNum, recData{file,1}.stimulatedT(:), ...
        'stable');
    % select single trials if requested
    if ~isequal(trials, 'all')
        stimulated = intersect(stimulated,trials);
    end

    % iterate through the trials computing the mean spike height/prominence
    for trial_index = 1:length(stimulated)
        % access the data of the current trial
        trial = stimulated(trial_index);
        dataINT = recData{file,1}.recordingINT(:,trial);
        % find the peaks in the data
        [peak_values, locations, ~, prom] = findpeaks(dataINT,...
            'Annotate','extents', 'WidthReference','halfprom',...
            'MinPeakProminence', 1,'MaxPeakWidth',200);
    
        % preassign a vector to gather the quater prominence width
        qprom_width = zeros(length(peak_values),1);
        % measure the peak width at one quarter of the prominence for each
        % spike (default would be half prom)
        for spike = 1:length(peak_values)
            % set the quarter prom as reference value for width measurement 
            measure_level = peak_values(spike) - 0.25*prom(spike);
            % select a window of 100 datapoints around the current spike if
            % possible, otherwise a one-sided smaller window
            if locations(spike) < 51
                spikeData = dataINT(1:locations(spike)+50);
                spikeLoc = locations(spike);
            elseif locations(spike)+50 > locations(end)
                spikeData = dataINT(locations(spike)-50:end);
                spikeLoc = 51;
            else
                spikeData = ...
                    dataINT(locations(spike)-50:locations(spike)+50);
                spikeLoc = 51;
            end
            % devide the points in left and right of the peak
            dataL = spikeData(1:spikeLoc);
            dataR = spikeData(spikeLoc+1:end);
            % take the points (left and right) that are closest to the
            % measure level of a quarter of the prominence
            [~,minL] = min((dataL - measure_level), [], ...
                'ComparisonMethod','abs');
            [~,minR] = min((dataR - measure_level), [], ...
                'ComparisonMethod','abs');
            % calculate the final width at quarter peak prominence
            qprom_width(spike) = minR + spikeLoc - minL;
        end
        % sort and select the all found peaks in a 3-step algorithm 
        % 1st step: exclude every peak beneath the overall mean +1 mV 
        select_peaks = peak_values > mean(dataINT)+1;
        % select the corresponding values/properties
        qprom_width_subset1 = qprom_width(select_peaks);
        prom_subset1 = prom(select_peaks);
        
        % 2nd step: select only peaks with a quarter peak prom over 4 and
        % under 45
        select_peaks = 4 < qprom_width_subset1 & qprom_width_subset1 < 45;
        % select the coorresponding values/properties
        prom_subset2 = prom_subset1(select_peaks);
        
        % 3rd step: exclude the peaks smaller than 20% of the largest peak
        select_peaks = prom_subset2 > max(prom_subset2)*0.2;
        % select the corresponding values/properties
        prom_subset3 = prom_subset2(select_peaks);
        
        % add the prominence values to the all_proms vector
        if eq(trial_index,1)
            all_proms = prom_subset3;
        else
            all_proms = [all_proms;prom_subset3];
        end
    
    % end of iteration through the trials
    end
    
    % calculate mean and sd over all spikes from all trials
    m_prom = mean(all_proms, 'omitnan');
    sd_prom = std(all_proms, 'omitnan');

    % define borders (mean +- sd) to fill
    fillX = [1, length(all_proms), length(all_proms), 1, 1];
    fillY = [(m_prom+sd_prom), (m_prom+sd_prom),(m_prom-sd_prom), ...
        (m_prom-sd_prom),(m_prom+sd_prom)];

    % plot the spike amplitudes over the index
    figure();
    plot(1:length(all_proms), all_proms)
    yline(m_prom, 'Color', [0 0.55 0], 'LineWidth', 2)
    yline((m_prom+sd_prom), 'Color', [1 0.67 0])
    yline((m_prom-sd_prom), 'Color', [1 0.67 0])
    % add an filled panel around the mean indicating one sd
    hold on;
    fill(fillX, fillY, [1 0.67 0], 'FaceAlpha', 0.2, 'LineStyle', 'none')
    hold off;
    % adapt the x limit
    xlim([0,length(all_proms)])
    % labelling
    xlabel("spike index over all trials");
    ylabel("spike amplitude [mV]")
    legend("spike amplitudes", "mean", sprintf("mean \x00B1 sd"))
    title(sprintf("Variation of the measured spike amplitude - " + ...
        "File: %g", fileNums(file)))

% end of iteration through the files
end

% end of function definition
end