bachelor_thesis/DataAnalysisToolbox/intSpikeHeight.m

function varargout = intSpikeHeight(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, ~,~,~,~,~,~, plotflag] = ...
    parse_input_peakfun(files, varargin{:});

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

% preallocate vectors to gather the results
m_prom = zeros(length(fileNums),1);
sd_prom = zeros(length(fileNums),1);

%% 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;
        peaks_subset1 = peak_values(select_peaks);
        % select the corresponding values/properties
        locs_subset1 = locations(select_peaks);
        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;
        peaks_subset2 = peaks_subset1(select_peaks);
        % select the coorresponding values/properties
        locs_subset2 = locs_subset1(select_peaks);
        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;
        peaks_subset3 = peaks_subset2(select_peaks);
        % select the corresponding values/properties
        locs_subset3 = locs_subset2(select_peaks);
        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
    
        % plot the peaks that are found plus their prominence and width
        if plotflag
            % create appropriate title
            plot_title = sprintf("File %g - trial %g", fileNums(file), ...
                trial);
            figure
            findpeaks(dataINT, 10000,...
                'Annotate','extents', 'WidthReference','halfprom',...
                'MinPeakProminence', 1,'MaxPeakWidth',200);
            locs_subset3_plot = locs_subset3 / 10000;
            hold on; plot(locs_subset3_plot, peaks_subset3, 'og'); hold off;
            title(plot_title)
            xlabel("time [s]"); ylabel("membrane potential [mV]")
            legend('signal','peaks from findpeaks function','prominence',...
                'width (half-prominence)',['peaks identified as ' ...
                'action potentials'], 'Location', 'southeast')
        end
    % end of iteration through the trials
    end
    
    % calculate mean and sd over all spikes from all trials
    m_prom(file,1) = mean(all_proms, 'omitnan');
    sd_prom(file,1) = std(all_proms, 'omitnan');

% end of iteration through the files
end

% return mean and sd values
varargout{1} = m_prom;
if nargout > 1
    varargout{2} = sd_prom;
end

% end of function definition
end