DevianceDetectionHumanABR/Deviance Detection_Humans/DD_Plotting.m

%% define some variables

filt = 4; % 1: non, 2: 10 Hz, 3: 45 Hz, 4: 50 Hz, 5: 55 Hz
cutTime = 2; % plot only certain time window? (0: no, 1: short, 2: long)
avType = 1; % plot grand average (1) or individual averages (2)
fixedY = 1; % use dynamic (0) or fixed (1) y-axis?
inclCtrl = 0; % plot control data? (0: no, 1: 50 %, 2: MR, 3: DevOnly)
statsType = 1; % decide which statistics to use if control is included; 1: t-test, 2: ANOVA
plotStatsI = 1; % plot stats-windows?
varType = 1; % plot stats of amplitude or timing variable? (1: amplitude, 2: timimg)
uWin = (1); % used windows for stats (1:4) (it's always 4 if varType is 2)
repRateI = 1; % decide which rep rate to plot (1: 25 ms, 2: 30 ms, 3: 50 ms)
subStdI = 2; % decide which dataset to plot; 1: Std after Dev, 2: Std before Dev (does not affect omission responses)
anaType = 1; % 1: window-based analysis; 2: point-based analysis
plotRespI = 1; % plot responses?
plotDifI = 0; % plot difference curves?
plotOmI = 0; % plot omission responses?
plotBpI = 1; % plot boxplots?
plotCsResI = 0; % plot consecutive standard responses?
plotCsValI = 0; % plot response strengths of consecutive standard responses
clMark = 0; % plot significant clusters as bars? (no: 0, bar: 1)
cStart = 1; % stimulus combination to start with
oLowHigh = 0; % plot only Low and High chirp?
saveI = 0; % save figure?

chirpPeak = [0.00648,0.00842];

% decide which set of rms-windows to plot
% boxpFigWi = 800; % width of boxplot figure
switch size(uWin,2)
    case 1
        boxpFigWi = 300; % width of boxplot figure
    case 2
        boxpFigWi = 500; % width of boxplot figure
end

% switch plotOmI
%     case 1
%         uWin = 1;
% end



switch filt
    case 1
        switch plotOmI
            case 0
                filtName = 'noF';
                y_lim = [-0.8,1];
            case 1
                filtName = 'noF';
                y_lim = [-0.8,1];
        end
    case 2
        switch plotOmI
            case 0
                filtName = '10F';
                y_lim = [-0.8,1];
            case 1
                filtName = '10F';
                y_lim = [-0.8,1];
        end
    case 3
        switch plotOmI
            case 0
                filtName = '45F';
                y_lim = [-0.8,1];
            case 1
                filtName = '45F';
                y_lim = [-0.8,1];
        end
    case 4
        switch plotOmI
            case 0
                filtName = '50F';
                y_lim = [-0.35,0.55];
            case 1
                filtName = '50F';
                y_lim = [-0.35,0.55];
        end
    case 5
        switch plotOmI
            case 0
                filtName = '55F';
                y_lim = [-0.3,0.55];
            case 1
                filtName = '55F';
                y_lim = [-0.3,0.55];
        end

end

% define colors and names
colDev = [1,0,0]; % dev: red
colStd = [0,0,1]; % std: blue
col50Per = [1,0,1]; % 50Per: magenta
colMR = [0,0,0]; % MR: black
colDevO = [0,0,0]; % Dev only: black
switch plotOmI
    case 0
        nameDev = 'Deviant';
        nameStd = 'Standard';
    case 1
        nameDev = 'Omitted deviant';
        nameStd = 'Omitted standard';
end
name50Per = '50 % Control';
nameMR = 'MS Control';
nameDevO = 'Dev Only';

%% load processed data

switch repRateI
    case 1
        repRate = '25ms';
    case 2
        repRate = '30ms';
    case 3
        repRate = '50ms';
end
switch subStdI
    case 1
        subStd = 'SaD';
    case 2
        subStd = 'SbD';
end

switch inclCtrl
    case 0
        load(['DD_CbpData_noCtrl_',repRate,'_',subStd,'.mat']) % load stats data
        load(['DD_statsData_noCtrl_',repRate,'_',subStd,'.mat']) % load stats data
    case 1 % load 50Per data and rename important variables
        load(['DD_avData_50Per_',repRate,'_',subStd,'.mat'])
        load('DD_CbpData_50Per_',repRate,'_',subStd,'.mat') % load stats data
        load('DD_statsData_50Per_',repRate,'_',subStd,'.mat') % load stats data
        dataAv_cell_Ctrl = dataAv_cell; % rename
        dataGrAv_cell_Ctrl = dataGrAv_cell;
        dataDifAv_cell_Ctrl = dataDifAv_cell; % rename
        dataDifGrAv_cell_Ctrl = dataDifGrAv_cell;
        dataSe_cell_Ctrl = dataSe_cell;
        filenames_stim_cell_Ctrl = filenames_stim_cell;
        nameCtrl = name50Per;
        colCtrl = col50Per; % set control color to 50 % color
    case 2 % load MR data and rename important variables
        load(['DD_avData_MR_',repRate,'_',subStd,'.mat'])
        load('DD_CbpData_MR_',repRate,'_',subStd,'.mat') % load stats data
        load('DD_statsData_MR_',repRate,'_',subStd,'.mat') % load stats data
        dataAv_cell_Ctrl = dataAv_cell; % rename
        dataGrAv_cell_Ctrl = dataGrAv_cell;
        dataDifAv_cell_Ctrl = dataDifAv_cell; % rename
        dataDifGrAv_cell_Ctrl = dataDifGrAv_cell;
        dataSe_cell_Ctrl = dataSe_cell;
        filenames_stim_cell_Ctrl = filenames_stim_cell;
        nameCtrl = nameMR;
        colCtrl = colMR; % set control color to MR color
    case 3
        load(['DD_avData_LowOnly_',repRate,'_',subStd,'.mat'])
        load('DD_CbpData_DevO_',repRate,'_',subStd,'.mat') % load stats data
        load('DD_statsData_DevO_',repRate,'_',subStd,'.mat') % load stats data
        dataAv_cell_CoLow = dataAv_cell; % rename
        dataGrAv_cell_CoLow = dataGrAv_cell;
        dataSe_cell_CoLow = dataSe_cell;
        filenames_cell_CoLow = filenames_stim_cell;
        load(['DD_avData_HighOnly_',repRate,'_',subStd,'.mat'])
        dataAv_cell_CoHigh = dataAv_cell; % rename
        dataGrAv_cell_CoHigh = dataGrAv_cell;
        dataSe_cell_CoHigh = dataSe_cell;
        filenames_cell_CoHigh = filenames_stim_cell;
        nameCtrl = nameDevO;
        colCtrl = colDevO; % set control color to MR color
end
switch plotOmI
    case 0
        % load Oddball data (will overwrite Ctrl variables that were not renamed)
        load(['DD_avData_Oddball_',repRate,'_',subStd,'.mat'])
    case 1
        load(['DD_avData_HighOm_',repRate,'.mat'])
        % load Om data
        dataAv_cell_HighOm = dataAv_cell; % rename
        dataGrAv_cell_HighOm = dataGrAv_cell;
        dataSe_cell_HighOm = dataSe_cell;
        filenames_cell_HighOm = filenames_stim_cell;
        load(['DD_avData_LowOm_',repRate,'.mat'])
        dataAv_cell_LowOm = dataAv_cell; % rename
        dataGrAv_cell_LowOm = dataGrAv_cell;
        dataSe_cell_LowOm = dataSe_cell;
        filenames_cell_LowOm = filenames_stim_cell;
        load(['DD_CbpData_Om_',repRate,'.mat']) % load stats data
        load(['DD_statsData_Om_',repRate,'.mat']) % load stats data
end
%% do some additional calculations

switch zsNormI
    case 0
        yLabelName = 'Voltage [µV]';
    case 1
        yLabelName = 'z-norm. voltage [a.u.]';
end

% change stats variables if timing is selected
if varType==2
    varS = timS;
    pV = pTV;
    cohensD_V = cohensDT_V;
    statsV = statsTV;
    uWin = 2;
end

switch cutTime
    case 1
        switch filt
            case {1,2,3}
                timeStart = -2; % starting time of the plot (relative to stim onset)
                timeEnd = 20; % ending time of the plot (relative to stim onset)
            case {4,5}
                timeStart = -2; % starting time of the plot (relative to stim onset)
                timeEnd = 20; % ending time of the plot (relative to stim onset)
        end
    case 2
        switch repRateI
            case 1
                timeStart = -2; % starting time of the plot (relative to stim onset)
                timeEnd = 10.5; % ending time of the plot (relative to stim onset)
            case 2
                timeStart = -2; % starting time of the plot (relative to stim onset)
                timeEnd = 10.5; % ending time of the plot (relative to stim onset)
            case 3
                timeStart = -2; % starting time of the plot (relative to stim onset)
                timeEnd = 10.5; % ending time of the plot (relative to stim onset)
        end
end

if oLowHigh==1 % plot only stim combination 7 (DisNoAM/Eloc) if oDisEloc is 1
    cStart = 1;
    nComb = 1;
end

combNameS = split(combName(cStart:end),","); % split comb names at comma
%% plotting: responses

if plotRespI==1
    for c = cStart:nComb % run once for each stimulus combination

        % extract data corresponding to current stimulation-condition from
        % cell-arrays
        switch plotDifI
            case 0
                switch plotOmI
                    case 0
                        dataAv = dataAv_cell{c};
                        dataGrAv = dataGrAv_cell{c};
                        dataSe = dataSe_cell{c};
                    case 1
                        dataAvHighOm = dataAv_cell_HighOm{c};
                        dataGrAvHighOm = dataGrAv_cell_HighOm{c};
                        dataSeHighOm = dataSe_cell_HighOm{c};
                        dataAvLowOm = dataAv_cell_LowOm{c};
                        dataGrAvLowOm = dataGrAv_cell_LowOm{c};
                        dataSeLowOm = dataSe_cell_LowOm{c};
                        combName = append("OR_",combName);
                end
                switch inclCtrl
                    case {1,2}
                        dataAvCtrl = dataAv_cell_Ctrl{c};
                        dataGrAvCtrl = dataGrAv_cell_Ctrl{c};
                        dataSeCtrl = dataSe_cell_Ctrl{c};
                        filenames_Ctrl = filenames_stim_cell_Ctrl{c};
                        bsPos = cell2mat(strfind(filenames_Ctrl,'\')); % find positions of back slashes in filename
                        % bsPos = strfind(filenames_Ctrl,'\'); % find positions of back slashes in filename
                        startPos = bsPos(:,end); % use last back slash position as starting point
                        filenames_Ctrl = extractAfter(filenames_Ctrl,startPos);
                    case 3
                        dataAvCoLow = dataAv_cell_CoLow{c};
                        dataGrAvCoLow = dataGrAv_cell_CoLow{c};
                        dataSeCoLow = dataSe_cell_CoLow{c};
                        filenames_CoLow = filenames_cell_CoLow{c};
                        bsPos = cell2mat(strfind(filenames_CoLow,'\')); % find positions of back slashes in filename
                        startPos = bsPos(:,end); % use last back slash position as starting point
                        filenames_CoLow = extractAfter(filenames_CoLow,startPos);
                        dataAvCoHigh = dataAv_cell_CoHigh{c};
                        dataGrAvCoHigh = dataGrAv_cell_CoHigh{c};
                        dataSeCoHigh = dataSe_cell_CoHigh{c};
                        filenames_CoHigh = filenames_cell_CoHigh{c};
                        bsPos = cell2mat(strfind(filenames_CoHigh,'\')); % find positions of back slashes in filename
                        startPos = bsPos(:,end); % use last back slash position as starting point
                        filenames_CoHigh = extractAfter(filenames_CoHigh,startPos);
                end
            case 1
                dataAv = dataDifAv_cell_Ctrl{c};
                dataGrAv = dataDifGrAv_cell_Ctrl{c};
                plotStatsI = 0; % plot no significance windows when plotting diference curves
                switch inclCtrl
                    case 0
                        error('To plot difference curves, activate inclCtrl')
                    case {1,2}
                        dataAvCtrl = dataAv_cell_Ctrl{c}; % won't be plotted, only processed to keep the code simpler
                        dataGrAvCtrl = dataGrAv_cell_Ctrl{c};
                        filenames_Ctrl = filenames_cell_Ctrl{c};
                        bsPos = cell2mat(strfind(filenames_Ctrl,'\')); % find positions of back slashes in filename
                        startPos = bsPos(:,end); % use last back slash position as starting point
                        filenames_Ctrl = extractAfter(filenames_Ctrl,startPos);
                    case 3
                        dataAvCoLow = dataAv_cell_CoLow{c}; % won't be plotted, only processed to keep the code simpler
                        dataGrAvCoLow = dataGrAv_cell_CoLow{c};
                        filenames_CoLow = filenames_cell_CoLow{c};
                        bsPos = cell2mat(strfind(filenames_CoLow,'\')); % find positions of back slashes in filename
                        startPos = bsPos(:,end); % use last back slash position as starting point
                        filenames_CoLow = extractAfter(filenames_CoLow,startPos);
                        dataAvCoHigh = dataAv_cell_CoHigh{c}; % won't be plotted, only processed to keep the code simpler
                        dataGrAvCoHigh = dataGrAv_cell_CoHigh{c};
                        filenames_CoHigh = filenames_cell_CoHigh{c};
                        bsPos = cell2mat(strfind(filenames_CoHigh,'\')); % find positions of back slashes in filename
                        startPos = bsPos(:,end); % use last back slash position as starting point
                        filenames_CoHigh = extractAfter(filenames_CoHigh,startPos);
                end
        end
        filenames = filenames_stim_cell{c};
        recID = recID_cell{c};
        bsPos = cell2mat(strfind(filenames,'\')); % find positions of back slashes in filename
        % bsPos = strfind(filenames,'\'); % find positions of back slashes in filename
        startPos = bsPos(:,end); % use last back slash position as starting point
        filenames = extractAfter(filenames,startPos);
        timetr = round(timetr_cell{c}*1000,4);
        stimDur = stimDur_cell{c}(1)*1000;
        stimDelay = stimDelay_cell{c}(1)*1000;
        stimWin = [stimDelay,stimDelay+stimDur];
        nFiles = nFiles_cell{c};
        fs = fs_cell{c};

        % subdivide responses for simplicity
        switch plotOmI
            case 0
                ADevAv = dataAv{1};
                AStdAv = dataAv{2};
                BDevAv = dataAv{3};
                BStdAv = dataAv{4};
                ADevGrAv = dataGrAv{1};
                AStdGrAv = dataGrAv{2};
                BDevGrAv = dataGrAv{3};
                BStdGrAv = dataGrAv{4};
                ADevSe = dataSe{1};
                AStdSe = dataSe{2};
                BDevSe = dataSe{3};
                BStdSe = dataSe{4};
            case 1
                ADevAv = dataAvHighOm{3};
                AStdAv = dataAvHighOm{4};
                BDevAv = dataAvLowOm{1};
                BStdAv = dataAvLowOm{2};
                ADevGrAv = dataGrAvHighOm{3};
                AStdGrAv = dataGrAvHighOm{4};
                BDevGrAv = dataGrAvLowOm{1};
                BStdGrAv = dataGrAvLowOm{2};
                ADevSe = dataSeHighOm{3};
                AStdSe = dataSeHighOm{4};
                BDevSe = dataSeLowOm{1};
                BStdSe = dataSeLowOm{2};
        end
        switch inclCtrl
            case {1,2}
                AAvCtrl = dataAvCtrl{1};
                BAvCtrl = dataAvCtrl{2};
                AGrAvCtrl = dataGrAvCtrl{1};
                BGrAvCtrl = dataGrAvCtrl{2};
                ASeCtrl = dataSeCtrl{1};
                BSeCtrl = dataSeCtrl{2};
            case 3
                AAvCtrl = dataAvCoLow{1};
                BAvCtrl = dataAvCoHigh{3};
                AGrAvCtrl = dataGrAvCoLow{1};
                BGrAvCtrl = dataGrAvCoHigh{3};
                % TEST (plots std only when activated)
%                 AGrAvCtrl = dataGrAvCoLow{2};
%                 BGrAvCtrl = dataGrAvCoHigh{4};
                ASeCtrl = dataSeCoLow{1};
                BSeCtrl = dataSeCoHigh{3};
        end
            

        % plot single channel
        for s = 1:2 % once for A, once for B
            % process timetrace
            switch cutTime
                case 0
                    timeWin = floor([timetr(1),timetr(end)]); % data in this window will be plotted (relative to recording onset)
                    timeStart = timeWin(1)-(stimDelay);
                    timeEnd = timeWin(2)-(stimDelay);
                    timeCut = (1:size(timetr,2));
                case {1,2}
                    timeWin = [round((stimDelay)+timeStart+chirpPeak(s)*1000,4),...
                        round((stimDelay)+timeEnd+chirpPeak(s)*1000,4)]; % data in this window will be plotted (relative to recording onset)
                    timeCut = round(timeWin(1)/1000*fs:timeWin(2)/1000*fs+1); % "+1" to reach the proper length, otherwise the index will be 1 point too short since timetr starts at 0 ms
            end
            switch s
                case 1
                    grAvStd = AStdGrAv;
                    grAvDev = ADevGrAv;
                    avStd = AStdAv;
                    avDev = ADevAv;
                    seStd = AStdSe;
                    seDev = ADevSe;
                    stimName = 'A';
                    switch inclCtrl
                        case {1,2,3}
                            grAvCtrl = AGrAvCtrl;
                            avCtrl = AAvCtrl;
                            seCtrl = ASeCtrl;
                    end
                case 2
                    grAvStd = BStdGrAv;
                    grAvDev = BDevGrAv;
                    avStd = BStdAv;
                    avDev = BDevAv;
                    seStd = BStdSe;
                    seDev = BDevSe;
                    stimName = 'B';
                    switch inclCtrl
                        case {1,2,3}
                            grAvCtrl = BGrAvCtrl;
                            avCtrl = BAvCtrl;
                            seCtrl = BSeCtrl;
                    end
            end
            figure('NumberTitle','off','Name',[combName{c},'_',stimName],'Position',[0,0,1000,600])
            title([convertStringsToChars(combNameS(s,c)),'-Freq. Chirp'])
            switch avType
                case 1 % grand averages
                    hold on
                    switch inclCtrl
                        case {1,2,3}
                            switch plotDifI
                                case 0 % plot control response only when difference curve is deactivated
                                    % plot standard error
                                    seHiCtrl = grAvCtrl(timeCut,filt)+seCtrl(timeCut,filt);
                                    seLoCtrl = grAvCtrl(timeCut,filt)-seCtrl(timeCut,filt);
                                    seX = [timetr(timeCut),fliplr(timetr(timeCut))];
                                    seY = [seHiCtrl',fliplr(seLoCtrl')];
                                    patch(seX,seY,colCtrl,'FaceAlpha',0.2,'EdgeColor','none')
                                    % plot data
                                    plotCtrl = plot(timetr(timeCut),grAvCtrl(timeCut,filt),'color',colCtrl,'Linewidth',2);
                            end
                    end
                    % plot standard error
                    seHiDev = grAvDev(timeCut,filt)+seDev(timeCut,filt);
                    seLoDev = grAvDev(timeCut,filt)-seDev(timeCut,filt);
                    seX = [timetr(timeCut),fliplr(timetr(timeCut))];
                    seY = [seHiDev',fliplr(seLoDev')];
                    patch(seX,seY,colDev,'FaceAlpha',0.2,'EdgeColor','none')
                    % plot data
                    plotDev = plot(timetr(timeCut),grAvDev(timeCut,filt),'color',colDev,'Linewidth',2);
                    % plot standard error
                    seHiStd = grAvStd(timeCut,filt)+seStd(timeCut,filt);
                    seLoStd = grAvStd(timeCut,filt)-seStd(timeCut,filt);
                    seX = [timetr(timeCut),fliplr(timetr(timeCut))];
                    seY = [seHiStd',fliplr(seLoStd')];
                    patch(seX,seY,colStd,'FaceAlpha',0.2,'EdgeColor','none')
                    % plot data
                    plotStd = plot(timetr(timeCut),grAvStd(timeCut,filt),'color',colStd,'Linewidth',2);
                    xticks(timeWin(1):2:timeWin(2))
                    set(gca,'XTickLabel',timeStart:2:timeEnd)
                    switch clMark
                        case 1 % plot bars respresenting clusters
                            yDist = y_lim(2)-y_lim(1);
                            yMin = y_lim(1)+yDist*0.05;
                            yCl = [yMin+yDist*0.05;yMin+yDist*0.05;yMin+yDist*0.07;yMin+yDist*0.07];
                            % count number of significant clusters in vector
                            pSmpl = length(grAvStd);
                            vecSigClIs = [0,reshape(vecSigClI(c,s,filt,:),1,pSmpl)]; % add 0 at start to make sure possible cluster starting at point 1 is considered
                            vecSigClIe = [reshape(vecSigClI(c,s,filt,:),1,pSmpl),0]; % add 0 at the end to make sure possible cluster lasting until the end is considered
                            startSigCl = find(diff(vecSigClIs)==1); % find all points where difference vector is 1 (starting points of clusters in original vector)
                            endSigCl = find(diff(vecSigClIe)==-1); % find all points where difference vector is -1 (ending points of clusters in original vector)
                            nCl = numel(startSigCl); % return number of cluster starting points = number of clusters
                            xCl  = zeros(4,nCl);
                            for cl = 1:nCl
                                xCl(:,cl) = ([startSigCl(cl);endSigCl(cl);endSigCl(cl);startSigCl(cl)])./fs*1000;
                                patch('XData',xCl(:,cl),'YData',yCl,'EdgeColor','Black','EdgeAlpha',1,'FaceColor','Black','FaceAlpha',1)
                            end
                    end
                    xlabel('Time [ms]','FontWeight','bold','FontSize',20)
                    ylabel(yLabelName,'FontWeight','bold','FontSize',20)
                    switch inclCtrl
                        case 0
                            legend([plotDev,plotStd],nameDev,nameStd,'Location','northwest','AutoUpdate','off')
                        case {1,2,3}
                            switch plotDifI
                                case 0
                                    legend([plotDev,plotStd,plotCtrl],nameDev,nameStd,nameCtrl,'Location','northwest','AutoUpdate','off')
                                case 1
                                    legend([plotDev,plotStd],nameDev,nameStd,'Location','northwest','AutoUpdate','off')
                            end
                    end
                    set(gca,'FontSize',20,'Linewidth',2,'FontName','Arial')
                    switch plotStatsI
                        case 1
                            x_min = timetr(winAll{s}(:,1));
                            x_max = timetr(winAll{s}(:,2));
                            x = [x_min;x_max;x_max;x_min];
                            switch fixedY
                                case 0
                                    y_limAuto = ylim;
                                    yMin = y_limAuto(1);
                                    y_max = y_limAuto(2);
                                case 1
                                    yMin = y_lim(1);
                                    y_max = y_lim(2);
                            end
                            y_temp = [yMin;yMin;y_max;y_max];
                            y = repmat(y_temp,1,uWin(end));
                            for w = uWin
                                switch inclCtrl
                                    case 0
                                        pVTemp = pV(c,s,filt,w);
                                        % pVTemp
                                        % statsV{c,s,filt,w}
                                        % cohensD_V(c,s,filt,w,1)
                                    case {1,2,3}
                                        switch statsType
                                            case 1 % use t-test statistics
                                                pVTemp = pV(c,s,filt,w);
                                            case 2 % use ANOVA statistics
                                                pVTemp = multcompV{c,s,filt,w}.pValue(1); % extract pValue for comparison between dev and std condition
                                                % pVTemp
                                                % statsV{c,s,filt,w}
                                                % cohensD_V(c,s,filt,w,1)
                                        end
                                end
                                if pVTemp<0.05
                                    sigWin = patch('XData',x(:,w),'YData',y(:,w),'EdgeColor','black','EdgeAlpha',1,'FaceColor','black','FaceAlpha',.2,'Linewidth',2);
                                else
                                    sigWin = patch('XData',x(:,w),'YData',y(:,w),'EdgeColor','black','EdgeAlpha',1,'FaceColor','none','Linewidth',2);
                                end
                                uistack(sigWin,'bottom')
                            end
                    end
                    xlim([timeWin(1),timeWin(2)])
                    if fixedY==1
                        ylim(y_lim)
                    end
                case 2 % individual averages
                    for f = 1:nFiles
                        figure ('NumberTitle','off','Name',[combName{c},'_',stimName,'_',convertStringsToChars(recID(f))],'Position',[0,0,1000,500])
                        hold on
                        switch inclCtrl
                            case {1,2,3}
                                switch plotDifI
                                    case 0 % plot control response only when difference curve is deactivated
                                        plotCtrl = plot(timetr(timeCut),avCtrl(timeCut,f,filt),'color',colCtrl,'Linewidth',2);
                                end
                        end
                        plotStd = plot(timetr(timeCut),avStd(timeCut,f,filt),'color',colStd,'Linewidth',2);
                        plotDev = plot(timetr(timeCut),avDev(timeCut,f,filt),'color',colDev,'Linewidth',2);
                        xticks(timeWin(1):5:timeWin(2))
                        set(gca,'XTickLabel',timeStart:5:timeEnd)
                        switch clMark
                            case 1 % plot bar
                                yDist = y_lim(2)-y_lim(1);
                                yMin = y_lim(1)+yDist*0.05;
                                x_stim = [stimWin(1);stimWin(2);stimWin(2);stimWin(1)];
                                yCl = [yMin+yDist*0.05;yMin+yDist*0.05;yMin+yDist*0.07;yMin+yDist*0.07];
                                patch('XData',x_stim,'YData',yCl,'EdgeColor','Black','EdgeAlpha',1,'FaceColor','Black','FaceAlpha',1)
                        end
                        xlabel('Time [ms]','FontWeight','bold','FontSize',20)
                        ylabel(yLabelName,'FontWeight','bold','FontSize',20)
                        set(gca,'FontSize',20,'Linewidth',2,'FontName','Arial')
                        xlim([timeWin(1),timeWin(2)])
                        if fixedY==1
                            ylim(y_lim)
                        end
                        switch plotStatsI
                            case 1
                                x_min = timetr(winAll{s}(:,1)+blcEnd(s)*fs(1));
                                x_max = timetr(winAll{s}(:,2)+blcEnd(s)*fs(1));
                                x = [x_min;x_max;x_max;x_min];
                                switch fixedY
                                    case 0
                                        y_limAuto = ylim;
                                        yMin = y_limAuto(1);
                                        y_max = y_limAuto(2);
                                    case 1
                                        yMin = y_lim(1);
                                        y_max = y_lim(2);
                                end
                                y_temp = [yMin;yMin;y_max;y_max];
                                y = repmat(y_temp,1,uWin(end));
                                for w = uWin
                                    sigWin = patch('XData',x(:,w),'YData',y(:,w),'EdgeColor','black','EdgeAlpha',1,'FaceColor','none','Linewidth',2);
                                    uistack(sigWin,'bottom')
                                    txt = {['Dev: ',num2str(varS{c}(f,s,1,filt,w))],['Std: ',num2str(varS{c}(f,s,2,filt,w))]};
                                    text(x(1,w)-x(1,w)*0.2,y(3,w)-y(3,w)*0.2,txt,'FontSize',20,'Linewidth',2,'FontName','Arial')
                                end
                        end
                    end
            end
            %             xline(((pts2begin_cell{1}(1)/fs(1)))*1000,'Linewidth',2)
            % xline((chirpPeak(s)+(pts2begin_cell{1}(1)/fs(1)))*1000,'Linewidth',2) % mark blc end
            if saveI==1
                saveas(gcf,[combName{c},'_',stimName,'_',filtName,'_cutTime',num2str(cutTime),'_ctrl',num2str(inclCtrl),'_uWin',num2str(uWin),'_varType',num2str(varType),'_repRate',['_',repRate],'_subStd',['_',subStd],'.jpg'])
                saveas(gcf,[combName{c},'_',stimName,'_',filtName,'_cutTime',num2str(cutTime),'_ctrl',num2str(inclCtrl),'_uWin',num2str(uWin),'_varType',num2str(varType),'_repRate',['_',repRate],'_subStd',['_',subStd],'.svg'])
                saveas(gcf,[combName{c},'_',stimName,'_',filtName,'_cutTime',num2str(cutTime),'_ctrl',num2str(inclCtrl),'_uWin',num2str(uWin),'_varType',num2str(varType),'_repRate',['_',repRate],'_subStd',['_',subStd],'.fig'])
            end
        end
    end
end
%% plotting: boxplots

switch varType
    case 1
        titleBp = 'ABR Amplitude';
        yLabelNameBp = 'ABR p-p Amp. [µV]';
        % yLabelNameBp = yLabelName;
    case 2
        titleBp = 'Peak V Latency';
        yLabelNameBp = 'Peak V Latency [ms]';
end

if plotBpI==1
    switch inclCtrl
        case 0
            xDev = [1,4,7,10,13];
            xStd = [2,5,8,11,14];
            sig_groups = {[1,2],[4,5],[7,8],[10,11],[13,14]};
        case {1,2,3}
            xDev = [1,5,9,13,17];
            xStd = [2,6,10,14,18];
            xCtrl = [3,7,11,15,19];
            sig_groups1 = {[1,2],[5,6],[9,10],[13,14],[17,18]}; % dev vs. std
            sig_groups2 = {[1,3],[5,7],[9,11],[13,15],[17,19]}; % dev vs. ctrl
            sig_groups3 = {[2,3],[6,7],[10,11],[14,15],[18,19]}; % std vs. ctrl
    end

    for c = cStart:nComb % run once for each stimulus combination
        for s = 1:2 % once for A, once for B
            switch s
                case 1
                    stimName = 'A';
                case 2
                    stimName = 'B';
            end
            rmsMin = floor(min(varS{c}(:,s,:,filt,:),[],'all')); % identify minimum RMS value for y-limit
            rmsMax = ceil(max(varS{c}(:,s,:,filt,:),[],'all')); % identify maximum RMS value for y-limit
            figure('NumberTitle','off','Name',['Boxpl_',combName{c},'_',stimName],'Position',[0,0,boxpFigWi,400])
            % title(titleBp)
            ylim([rmsMin-0.5,rmsMax+1]) % doesn't set the actual y-lim, is just required to get an equal distance between lines and asterisks of sigstar. Sigstar sets the actual y-lim itself
            for w = uWin
                hold on
                boxDev = boxchart(xDev(w)*ones(size(varS{c}(:,s,1,filt,w))),varS{c}(:,s,1,filt,w),'BoxFaceColor',colDev,'MarkerColor',colDev,'Linewidth',2);
                boxStd = boxchart(xStd(w)*ones(size(varS{c}(:,s,2,filt,w))),varS{c}(:,s,2,filt,w),'BoxFaceColor',colStd,'MarkerColor',colStd,'Linewidth',2);
                varSdif = mean(varS{c}(:,s,2,filt,w)-varS{c}(:,s,1,filt,w))
                varSdifStdd = std(varS{c}(:,s,2,filt,w)-varS{c}(:,s,1,filt,w))
                switch inclCtrl
                    case 0
                        H = sigstar(sig_groups{w},pV(c,s,filt,w),cohensD_V(c,s,filt,w,1));
                    case {1,2,3}
                        boxCtrl = boxchart(xCtrl(w)*ones(size(varS{c}(:,s,3,filt,w))),varS{c}(:,s,3,filt,w),'BoxFaceColor',colCtrl,'MarkerColor',colCtrl,'Linewidth',2);
                        switch statsType
                            case 1
                                H = sigstar(sig_groups1{w},pV(c,s,filt,w),cohensD_V(c,s,filt,w,1));
                            case 2
                                pVTemp = multcompV{c,s,filt,w}.pValue([1,2,4]);
                                eS = zeros(3,1);
                                for cmp = 1:3
                                    eS(cmp) = cohensD_V(c,s,filt,w,cmp);
                                end
                                H = sigstar({sig_groups1{w},sig_groups2{w},sig_groups3{w}},pVTemp,eS);
                        end
                end
                set(H(:,2),'FontSize',25,'FontName','Arial')
            end
            switch inclCtrl
                case 0
                    xlim([xDev(uWin(1))-1,xStd(uWin(end))+1])
                    xticks((xDev(uWin)+xStd(uWin))/2)
%                     legend([boxDev,boxStd],nameDev,nameStd,'Location','northwest','AutoUpdate','off')
                case {1,2,3}
                    xlim([xDev(uWin(1))-1,xCtrl(uWin(end))+1])
                    xticks(xStd(uWin))
%                     legend([boxDev,boxStd,boxCtrl],nameDev,nameStd,nameCtrl,'Location','northwest','AutoUpdate','off')
            end
            ylabel(yLabelNameBp,'FontWeight','bold','FontSize',20)
            % xticklabels(winTxt(uWin))
            xticks([])
            set(gca,'FontSize',25,'Linewidth',2,'FontName','Arial')
            if saveI==1
                saveas(gcf,['Boxp_',combName{c},'_',stimName,'_',filtName,'_cutTime',num2str(cutTime),'_ctrl',num2str(inclCtrl),'_uWin','_varType',num2str(varType),'_uWin',num2str(uWin),'_repRate',['_',repRate],'_subStd',['_',subStd],'.jpg'])
                saveas(gcf,['Boxp_',combName{c},'_',stimName,'_',filtName,'_cutTime',num2str(cutTime),'_ctrl',num2str(inclCtrl),'_uWin','_varType',num2str(varType),'_uWin',num2str(uWin),'_repRate',['_',repRate],'_subStd',['_',subStd],'.svg'])
                saveas(gcf,['Boxp_',combName{c},'_',stimName,'_',filtName,'_cutTime',num2str(cutTime),'_ctrl',num2str(inclCtrl),'_uWin','_varType',num2str(varType),'_uWin',num2str(uWin),'_repRate',['_',repRate],'_subStd',['_',subStd],'.fig'])
            end
        end
    end
end
%% consecutive stds
if plotCsValI==1   
    for c = 1:nComb % once for each stimulus combination
        for s = 1:2 % once for A, once for B
            switch s
                case 1
                    stimName = 'A'; % ConS starts with B (High), since B is std in oddball block 1 
                case 2
                    stimName = 'B';
            end
            switch varType
                case 1
                    varMax = max(abs(varAvCsS{c}(s,filt,uWin,:)));
                    varPlot = reshape(varAvCsS{c}(s,filt,uWin,:),1,nConS)/varMax; % plot amplitude variable
                    sePlot = reshape(varSeCsS{c}(s,filt,uWin,:),1,nConS);
                case 2
                    varMax = max(abs(timAvCsS{c}(s,filt,uWin,:)));
                    varPlot = reshape(timAvCsS{c}(s,filt,uWin,:),1,nConS)/varMax; % plot timing variable
                    sePlot = reshape(timSeCsS{c}(s,filt,uWin,:),1,nConS);
            end
            figure('NumberTitle','off','Name',[combName{c},'_',stimName],'Position',[0,0,1000,600])
            hold on
            title([convertStringsToChars(combNameS(s,c)),'-Freq. Chirp'])
            errConSV = errorbar(varPlot,sePlot,'k','LineStyle','none','Marker','.','MarkerSize',30,'Linewidth',2);
            % plotConSV = plot(varPlot,'k','LineStyle','none','Marker','.','MarkerSize',30,'Linewidth',2);
        end
    end
end
if plotCsResI==1
    for c = 1:nComb % once for each stimulus combination
        pSmpl = size(conSGrAv{c},1);
        fs = fs_cell{c};
        timetrConS = (0:nConS*pSmpl-1)/fs(1); % time trace of one trial
        filenames = filenames_stim_cell{c};
        recID = recID_cell{c};
        bsPos = cell2mat(strfind(filenames,'\')); % find positions of back slashes in filename
        % bsPos = strfind(filenames,'\'); % find positions of back slashes in filename
        startPos = bsPos(:,end); % use last back slash position as starting point
        filenames = extractAfter(filenames,startPos);
        nFiles = nFiles_cell{c};
        for s = 1:2 % once for A, once for B
            switch s
                case 1
                    stimName = 'A'; % ConS starts with B (High), since B is std in oddball block 1 
                case 2
                    stimName = 'B';
            end
            conSPlot = zeros(pSmpl*nConS,1);
            switch avType
                case 1 % grand average
                    for cs = 1:nConS
                        conSPlot(1+pSmpl*(cs-1):pSmpl*cs) = conSGrAv{c}(:,s,cs,filt);
                        % figure
                        % plot(conSGrAv{c}(:,s,cs,filt),'k','Linewidth',2)
                    end
                    figure('NumberTitle','off','Name',[combName{c},'_',stimName],'Position',[0,0,1000,600])
                    title([convertStringsToChars(combNameS(s,c)),'-Freq. Chirp'])
                    hold on
                    plotConS = plot(timetrConS,conSPlot,'k','Linewidth',2);
                    switch plotStatsI
                        case {1,2}
                            for cn = 1:nConS
                                x_min = timetrConS(winAll{s}(:,1)+round(blcEnd(s)*fs(1)+pSmpl*(cn-1)+1));
                                x_max = timetrConS(winAll{s}(:,2)+round(blcEnd(s)*fs(1)+pSmpl*(cn-1)+1));
                                x = [x_min;x_max;x_max;x_min];
                                switch fixedY
                                    case 0
                                        y_limAuto = ylim;
                                        y_min = y_limAuto(1);
                                        y_max = y_limAuto(2);
                                    case 1
                                        y_min = y_lim(1);
                                        y_max = y_lim(2);
                                end
                                y_temp = [y_min;y_min;y_max;y_max];
                                y = repmat(y_temp,1,uWin(end));
                                for w = uWin
                                    sigWin = patch('XData',x(:,w),'YData',y(:,w),'EdgeColor','red','EdgeAlpha',1,'FaceColor','none','Linewidth',2,'LineStyle','--');
                                    uistack(sigWin,'bottom')
                                end
                            end
                    end
                case 2 % individual average
                    for f = 1:nFiles
                        for cs = 1:nConS
                            conSPlot(1+pSmpl*(cs-1):pSmpl*cs) = conSAv{c}(:,s,f,cs,filt);
                        end
                        figure('NumberTitle','off','Name',[combName{c},'_',stimName,'_',convertStringsToChars(filenames(f))],'Position',[0,0,1000,600])
                        title([convertStringsToChars(combNameS(s,c)),'-Freq. Chirp'])
                        hold on
                        plotConS = plot(timetrConS,conSPlot,'k','Linewidth',2);
                        switch plotStatsI
                            case {1,2}
                                for cn = 1:nConS
                                    x_min = timetrConS(winAll{s}(:,1)+round(blcEnd(s)*fs(1)+pSmpl*(cn-1)+1));
                                    x_max = timetrConS(winAll{s}(:,2)+round(blcEnd(s)*fs(1)+pSmpl*(cn-1)+1));
                                    x = [x_min;x_max;x_max;x_min];
                                    switch fixedY
                                        case 0
                                            y_limAuto = ylim;
                                            y_min = y_limAuto(1);
                                            y_max = y_limAuto(2);
                                        case 1
                                            y_min = y_lim(1);
                                            y_max = y_lim(2);
                                    end
                                    y_temp = [y_min;y_min;y_max;y_max];
                                    y = repmat(y_temp,1,uWin(end));
                                    for w = uWin
                                        sigWin = patch('XData',x(:,w),'YData',y(:,w),'EdgeColor','red','EdgeAlpha',1,'FaceColor','none','Linewidth',2,'LineStyle','--');
                                        uistack(sigWin,'bottom')
                                    end
                                end
                        end
                    end
            end
        end
    end
end