DevianceDetectionHumanABR/Repetition Suppression_Humans/RS_Plotting.m

%% define some variables

filt = 6; % 1: non, 2: 0.1 Hz, 3: 10 Hz, 4: 40 Hz, 5: 45 Hz, 6: 50 Hz
cutTime = 0; % plot only certain time window? (0: no, 1: short, 2: long)
avType = 1; % plot grand average (1) or individual averages (2)
varType = 1; % plot stats of amplitude or timing variable? (1: amplitude, 2: timimg)
repRateI = 1; % which rep rate to plot (1: 25 ms, 2: 50 ms)
fixedY = 0; % use dynamic (0) or fixed (1) y-axis?
plotRespI = 0;
plotRespOsciI = 1; % plot oscillogram of the responses?
plotRespSpecI = 0; % plot spectrogram of the responses?
plotWinI = 1; % plot stat windows?
plotSeI = 0; % plot data with or without Standard error?
win = 1; % stats window to be plotted (it's always 2 if varType is 2)
plotCurveI = 1; % plot fitted curve?
curveInt = 1; % which interval should be plotted on x? (1: Rep#; 2: time)
saveI = 1; % save figure?

chirpPeak = [0.00648,0.00842];
%% load processed data

switch repRateI
    case 1
        load('RS_avData_25ms.mat') % load preprocessed data
        load('RS_statsData_25ms.mat') % load stats data
    case 2
        load('RS_avData_50ms.mat') % load data
        load('RS_statsData_50ms.mat') % load stats data
end
col = 'k';
%% do some additional calculations

% decide which set of rms-windows to plot
if varType==2
    win = 2; % always use second window if plotting peak timing
end
uWin = (win:win);

% define y-label name
switch zsNormI
    case 0
        ylabelName = 'Voltage [µV]';
    case 1
        ylabelName = 'z-norm. Voltage [a.u.]';
end

% define x-Axis
switch cutTime
    case 0
        xtickSteps = 100; % ms-steps on x-axis that should be plotted
    case 1
        xtickSteps = 5; % ms-steps on x-axis that should be plotted
        switch filt
            case {1,2,3}
                timeStart = 0; % starting time of the plot (relative to stim onset)
                timeEnd = 20; % ending time of the plot (relative to stim onset)
            case {4,5}
                timeStart = 0; % starting time of the plot (relative to stim onset)
                timeEnd = 20; % ending time of the plot (relative to stim onset)
        end
    case 2
        xtickSteps = 50; % ms-steps on x-axis that should be plotted
        timeStart = 0; % starting time of the plot (relative to stim onset)
        switch repRateI
            case 1
                timeEnd = 350; % ending time of the plot (relative to stim onset)
            case 2
                timeEnd = 450; % ending time of the plot (relative to stim onset)
        end
end

% define some names
combNameS = split(combName,","); % split comb names at comma
if filt==1
    filtName = 'noF';
else
    filtName = num2str(lowc(filt-1));
end

% define y-axis depending on filter setting and plotted time frame
switch filt
    case 1
        switch cutTime
            case 0
                y_lim = [-0.8,1];
            case 1
                y_lim = [-0.8,1];
            case 2
                y_lim = [-0.8,1];
        end
    case 2
        switch cutTime
            case 0
                y_lim = [-3,2];
            case 1
                y_lim = [-3,2];
            case 2
                y_lim = [-3,2];
        end
    case 3
        switch cutTime
            case 0
                y_lim = [-0.9,1];
            case 1
                y_lim = [-0.3,0.55];
            case 2
                y_lim = [-0.9,1];
        end
    case 4
        switch cutTime
            case 0
                y_lim = [-0.3,0.7];
            case 1
                y_lim = [-0.3,0.7];
            case 2
                y_lim = [-0.3,0.7];
        end
    case 5
        switch cutTime
            case 0
                y_lim = [-0.3,0.7];
            case 1
                y_lim = [-0.3,0.7];
            case 2
                y_lim = [-0.3,0.7];
        end
    case 6
        switch cutTime
            case 0
                y_lim = [-0.3,0.7];
            case 1
                y_lim = [-0.3,0.7];
            case 2
                y_lim = [-0.3,0.7];
        end
end

% do some prior calculations for colormap
rgb1 = [ones(1,100),ones(1,100),(1:-0.01:0.01)]';
rgb2 = [ones(1,100),(1:-0.01:0.01),zeros(1,100)]';
rgb3 = [(1:-0.01:0.01),zeros(1,100),zeros(1,100)]';
cMap = [rgb1,rgb2,rgb3];

% do some prior calculations for tiled plots
yMax = 1.5; % y-axis maximum in kHz for spectrogram

% change stats variables if timing is selected
if varType==2
    varS = timS;
    varAvS = timAvS;
    varSeS = timSeS;
end
%% plotting: responses

if plotRespI==1
    switch avType
        case 1
            for f = 1:2 % once for each stimulation type (block-stim and o1st-stim)
                for s = 1:2 % run once for each stimulus combination
                    figure('NumberTitle','off','Name',[num2str(f),'_',SOAName,],'Position',[0,0,1500,1300])
                    set(gca,'FontSize',20,'Linewidth',2,'FontName','Arial')
                    c = (f-1)*2+s; % index of dataset to be used
                    nResp = nResp_cell{c};
                    % extract data corresponding to current stimulation-condition from
                    % cell-arrays
                    dataAv = dataAv_cell{c};
                    dataGrAv = dataGrAv_cell{c};
                    dataSe = dataSe_cell{c};
                    filenames = filenames_stim_cell{c};
                    recID = recID_cell{c};
                    bsPos = cell2mat(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(timetrUnit_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_stim = fs_stim_cell{c};
                    fs = fs_cell{c};
                    stimCat = stimCat_cell{c};
                    pts2begin = pts2begin_cell{c};

                    % 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
                    dataPlot = dataGrAv(timeCut,filt);
                    % oscillogram of responses
                    if plotRespOsciI==1
                        hold on
                        box on
                        % plot standard error
                        seHiCtrl = dataPlot+dataSe(timeCut,filt);
                        seLoCtrl = dataPlot-dataSe(timeCut,filt);
                        seX = [timetr(timeCut),fliplr(timetr(timeCut))];
                        seY = [seHiCtrl',fliplr(seLoCtrl')];
                        patch(seX,seY,col,'FaceAlpha',0.2,'EdgeColor','none')
                        % plot data
                        plotGrAv = plot(timetr(timeCut),dataPlot,col,'Linewidth',2);
                        xticks(timeWin(1):xtickSteps:timeWin(2))
                        set(gca,'XTickLabel',timeStart:xtickSteps:timeEnd)
                        xlabel('Time [ms]','FontWeight','bold','FontSize',20)
                        ylabel(ylabelName,'FontWeight','bold','FontSize',20)
                        set(gca,'FontSize',20,'Linewidth',2,'FontName','Arial')
                        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);
                                ylim(y_lim)
                        end
                        xlim([timeWin(1),timeWin(2)])
                        switch plotWinI
                            case {1,2}
                                x_min = timetr(winAllSti{c}(:,:,1))';
                                x_max = timetr(winAllSti{c}(:,:,2))';
                                x = cat(3,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 r = 1:nResp
                                    for w = uWin
                                        statsV = pV(c,filt,r,w);
                                        % markWin = patch('XData',x(1,w,:),'YData',y(:,w),'EdgeColor','black','EdgeAlpha',1,'FaceColor','none','Linewidth',2,'LineStyle','--'); % window for first response, to indicate area all other responses were compared to
                                        % uistack(markWin,'bottom')
                                        switch plotWinI
                                            case 1
                                                if statsV<0.05
                                                    sigWin = patch('XData',x(r,w,:),'YData',y(:,w),'EdgeColor','black','EdgeAlpha',1,'FaceColor','black','FaceAlpha',.2,'Linewidth',2);
                                                else
                                                    sigWin = patch('XData',x(r,w,:),'YData',y(:,w),'EdgeColor','black','EdgeAlpha',1,'FaceColor','none','Linewidth',2);
                                                end
                                            case 2
                                                sigWin = patch('XData',x(r,w,:),'YData',y(:,w),'EdgeColor','red','EdgeAlpha',1,'FaceColor','none','Linewidth',2,'LineStyle','--');
                                        end
                                        uistack(sigWin,'bottom')
                                    end
                                end
                        end
                    end
                    % spectrogram of responses
                    if plotRespSpecI==1
                        spectrogram(dataPlot,56,52,[],fs,'yaxis',"MinThreshold",-75)
                        colormap(cMap)
                        xLimAuto = xlim;
                        xlim([xLimAuto(1),timeEnd])
                        xticks(0:xtickSteps:timeEnd)
                        ylim([0,yMax])
                        colorbar('off')
                        ylabel('Frequency [kHz]','FontWeight','bold')
                        xlabel('Time [ms]','FontWeight','bold')
                        box on
                        set(gca,'FontSize',20,'Linewidth',2,'FontName','Arial')
                    end
                    hold on
                    title([convertStringsToChars(combNameS(s)),'-Freq. Chirp'])
                    if saveI==1
                        saveas(gcf,[num2str(c),'_',filtName,'F_',SOAName,'_cutTime',num2str(cutTime),'_to',num2str(timeEnd),'_win',num2str(win),'.jpg'])
                        saveas(gcf,[num2str(c),'_',filtName,'F_',SOAName,'_cutTime',num2str(cutTime),'_to',num2str(timeEnd),'_win',num2str(win),'.svg'])
                        % saveas(gcf,[num2str(f),'_',filtName,'_',SOAName,'_cutTime',num2str(cutTime),'_to',num2str(timeEnd),'.fig'])
                    end
                end
            end
        case 2
            for f = 1:1 % once for each stimulation type (block-stim and o1st-stim)
                for s = 1:2 % run once for each stimulus combination
                    c = (f-1)*2+s; % index of dataset to be used
                    nResp = nResp_cell{c};
                    % extract data corresponding to current stimulation-condition from
                    % cell-arrays
                    dataAv = dataAv_cell{c};
                    dataGrAv = dataGrAv_cell{c};
                    dataSe = dataSe_cell{c};
                    filenames = filenames_stim_cell{c};
                    recID = recID_cell{c};
                    bsPos = cell2mat(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(timetrUnit_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_stim = fs_stim_cell{c};
                    fs = fs_cell{c};
                    stimCat = stimCat_cell{c};
                    pts2begin = pts2begin_cell{c};

                    % 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,4),...
                                round((stimDelay)+timeEnd,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
                    for fi = 1:nFiles
                        figure('NumberTitle','off','Name',[num2str(f),'_',SOAName,'_',convertStringsToChars(filenames(fi))],'Position',[0,0,1500,1300])
                        set(gca,'FontSize',20,'Linewidth',2,'FontName','Arial')
                        dataPlot = dataAv(timeCut,fi,filt);
                        % oscillogram of responses
                        hold on
                        box on
                        % plot data
                        plotAv = plot(timetr(timeCut),dataPlot,col,'Linewidth',2);
                        xticks(timeWin(1):xtickSteps:timeWin(2))
                        set(gca,'XTickLabel',timeStart:xtickSteps:timeEnd)
                        xlabel('Time [ms]','FontWeight','bold','FontSize',20)
                        ylabel(ylabelName,'FontWeight','bold','FontSize',20)
                        set(gca,'FontSize',20,'Linewidth',2,'FontName','Arial')
                        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);
                                ylim(y_lim)
                        end
                        xlim([timeWin(1),timeWin(2)])
                        switch plotWinI
                            case {1,2}
                                x_min = timetr(winAllSti{c}(:,:,1))';
                                x_max = timetr(winAllSti{c}(:,:,2))';
                                x = cat(3,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 r = 1:nResp
                                    for w = uWin
                                        sigWin = patch('XData',x(r,w,:),'YData',y(:,w),'EdgeColor','red','EdgeAlpha',1,'FaceColor','none','Linewidth',2,'LineStyle','--');
                                        uistack(sigWin,'bottom')
                                    end
                                end
                        end
                        hold on
                        title([convertStringsToChars(combNameS(s)),'-Freq. Chirp'])
                    end
                end
            end
    end
end
%% plotting: curve fit

switch varType
    case 1
        titleC = 'ABR Amplitude';
        yLabelNameLC = 'ABR p-p Amp. [µV]';
        yLabelNameRC = 'Rel. ABR p-p Amp.';
    case 2
        titleC = 'Peak V Latency';
        yLabelNameLC = 'Peak V Latency [ms]';
        yLabelNameRC = 'Rel. Peak V Latency';
end

xMin = 0;
xMax = 10;
uWin = win;
if plotCurveI==1
    for f = 1:1 % run only for block-stimulation
        for s = 1:2 % run once for each stimulus combination
            c = (f-1)*2+s; % index of dataset to be used
            nResp = nResp_cell{c};
            varMax = max(abs(varAvS{c}(filt,1:end-1,uWin)));
            varPlot = reshape(varAvS{c}(filt,1:end-1,uWin),1,nResp-1);
            sePlot = reshape(varSeS{c}(filt,1:end-1,uWin),1,nResp-1);
            figure('NumberTitle','off','Name',[convertStringsToChars(combNameS(c)),'_',SOAName,'_Win',num2str(uWin)],'Position',[0,0,1500,500])
            hold on
            box on
            % plotFit = plot(fV{c,filt,uWin},'k');
            plotReg = plot(mdl{varType,c,filt,uWin});
            col = plotReg.Color;
            plotReg(1).Color = [0,0,0];
            plotReg(2).Color = [0,0,0];
            plotReg(3).Color = [0,0,0];
            plotReg(4).Color = [0,0,0];
            mdl{varType,c,filt,uWin}
            % plot absolute and relative y axis
            yyaxis left % absolute
            switch plotSeI
                case 0
                    plotVar = plot(varPlot,'k','LineStyle','none','Marker','.','MarkerSize',30,'Linewidth',2);
                case 1
                    plotVar = errorbar(varPlot,sePlot,'k','LineStyle','none','Marker','.','MarkerSize',30,'Linewidth',2);
            end
            % adjust ylim
            yLimL = ylim;
            switch fixedY
                case 0
                    yMinL = yLimL(1)-diff(yLimL)*0.1;
                    yMaxL = yLimL(2)+diff(yLimL)*0.1;
                case 1
                    yMinL = yMinAbs;
                    yMaxL = yMaxAbs;
            end
            ylim([yMinL,yMaxL]);
            yLimL = ylim;
            ylabel(yLabelNameLC,'FontWeight','bold','FontSize',30)
            yyaxis right % relative
            ylim(yLimL/varMax)
            yLimR = ylim;
            ylabel(yLabelNameRC,'FontWeight','bold','FontSize',30)
            % fV{varType,c,filt,uWin}
            % gofV{varType,c,filt,uWin}
            % cosmetics
            set(plotReg,'Linewidth',2);
            legend([plotVar,plotReg(2),plotReg(3)],'Data','Lin. Reg. Model','95 % Conf. Int.','Location','northeast')
            xticks(1:1:nResp-1)
            xticklabels(1:int(curveInt)*1:int(curveInt)*nResp-1)
            xlim([0,nResp])
            xLim = xlim;
            xlabel('Repetition in Stim. Train','FontWeight','bold','FontSize',30)
            % create text
            p = round(mdl{varType,c,filt,uWin}.Coefficients.pValue(2),3);
            if p<0.001
                ast = '***';
            elseif p<0.01
                ast = '**';
            elseif p<0.05
                ast = '*';
            else
                ast = '';
            end
            Rs = round(mdl{varType,c,filt,uWin}.Rsquared.Ordinary,3);
            txt = ['\itp',ast,'\rm = ',num2str(p),', \itR^2\rm = ',num2str(Rs)];
            text((xLim(2)*0.4),yLimR(2)-diff(yLimR)*0.1,txt,'FontSize',30)
            title([convertStringsToChars(combNameS(c)),'-Freq. Chirp'])
            set(gca,'FontSize',30,'Linewidth',2,'FontName','Arial','YColor','k')
            if saveI==1
                saveas(gcf,['Reg_',num2str(c),'_',filtName,'F_',SOAName,'_win',num2str(win),'_varType',num2str(varType),'.jpg'])
                saveas(gcf,['Reg_',num2str(c),'_',filtName,'F_',SOAName,'_win',num2str(win),'_varType',num2str(varType),'.svg'])
            end
        end
    end
end