Dateien hochladen nach ''

CohensD.m

@@ -0,0 +1,7 @@
+function [dret] = CohensD( distroA, distroB )
+    n1 = numel( distroA );
+    n2 = numel( distroB );
+    s = sqrt( ( (n1 - 1)*(std(distroA)^2) + (n2 - 1)*(std(distroB)^2) ) / (n1 + n2 - 2) );
+    dret = ( mean(distroA) - mean(distroB) ) / s;
+end
+

DD_OrgData.m

@@ -0,0 +1,72 @@
+%% Subscript: Organises data from datasets in multiple variables
+
+% subdivide dataset into its components
+if remI==1
+    data(remRec,:) = []; % remove certain recordings from analysis
+end
+filenames = cat(1,data{:,1});
+bsPos = cell2mat(strfind(filenames,'\')); % find positions of back slashes in filename
+startPos = bsPos(:,end); % use last back slash position as starting point
+ptPos = cell2mat(strfind(filenames,'.')); % find positions of points in filename (here only one)
+endPos = ptPos; % use point position as ending point
+recID = extractBetween(filenames,startPos,endPos,'Boundaries','Exclusive'); % extract between Positions
+switch runType
+    case 1
+        recID_oddb{c} = recID; % save recIDs for each oddball stimulus-combination
+        fileI_Oddb = 0; % file index is empty for runType 1 (oddball)
+        fileI_Ctrl = 0; % file index is empty for runType 1 (oddball)
+    case {2,3}
+        insect = intersect(recID_oddb{c},recID); % create array containing only those recordings that are available for both the oddball and the respective control paradigm
+        fileI_Oddb = ismember(recID_oddb{c},insect); % create logical array containing those oddball recordings that have a matching control recording
+        fileI_Ctrl = ismember(recID,insect);% create logical array containing those control recordings that have a matching oddball recording
+        data = data(fileI_Ctrl,:); % use only those control responses that have a matching oddball response
+        filenames = filenames(fileI_Ctrl); % use only those control filenames that have a matching oddball response
+end
+param = cat(1,data{:,2});
+seq = cat(1,data{:,4});
+rec_raw = cat(1,data(:,5));
+nFiles = size(data,1); % number of files within the variable
+%% define some additional parameters that are needed for calculations
+
+% stimulation-parameters; ATTENTION: always uses the parameters of first
+% file --> parameters must be consistent across all files!
+fs = param(1,4);
+fdown = param(1,17);
+fsDwn = fs/fdown;
+fIndx = param(1,2);
+switch runType
+    case {1,2} % Oddball, 50Per
+    Alvl = param(1,9);
+    Afrq = param(1,10);
+    Blvl = param(1,26);
+    Bfrq = param(1,25);
+    AfrqI = 999; % frequency index in MR control --> redundant but required in other runTypes
+    BfrqI = 999;
+end
+nDev = param(1,23);
+devProb = param(1,24);
+nTrials = param(1,12);
+nStd = nTrials-nDev;
+switch runType
+    case {1,2} % Oddball, 50Per
+        nBlcks = 2; % number of blocks (AB and BA)
+    case 3 % MR
+        nBlcks = 1; % only one block
+end
+stimDur = round(param(1,7),4);
+trDel = spDis/343.2; % travelling delay time of stimulus in s (343.2 is the speed of sound in air)
+stimDelay = round(param(1,13)+trDel,4); % stimulus delay in s (corrected for sound travelling delay)
+stimDelay = stimDelay+0.001; % add additional tone delay to stim delay
+pts2begin = round(stimDelay.*fsDwn); % samples before stimulus onset
+repper = round(param(1,8),4);
+recdur = repper;
+
+% additional parameters
+rawLen = zeros(1,nFiles); % preallocate
+for f = 1:nFiles
+    rawLen(f) = size(rec_raw{f},2); % number of samples in each channel of raw file
+end
+blockLenBuff = rawLen./nBlcks; % % length of one block (AB or BA) in samples, INCLUDING buffers
+blockLen = round(recdur.*nTrials.*fsDwn); % length of one block (AB or BA) in samples, WITHOUT buffers
+pSmpl = blockLen(1)/nTrials(1); % number of samples of each single response
+timetr = (0:pSmpl-1)/fsDwn(1); % time trace of one trial

DD_Plotting.m

@@ -0,0 +1,729 @@
+%% define some variables
+
+filt = 2; % 1: non, 2: 0.1 Hz, 3: 10 Hz, 4: 150 Hz, 5: 300 Hz
+cutTime = 1; % 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 = 2; % plot control data? (0: no, 1: 50 %, 2: MR, 3: both controls)
+statsType = 2; % decide which statistics to use if control is included; 1: t-test, 2: ANOVA
+plotWinI = 1; % plot RMS-windows?
+plotRespI = 1; % plot responses?
+plotDifI = 0; % plot difference curves?
+plotMtrxI = 0; % plot all responses in a matrix?
+plotHmI = 0; % plot heatmap?
+hmDataI = 1; % choose data that should be plotted as heatmap (1: Cohens D; 2: RMS of control, only useful for 50Per dataset)
+plotBpI = 0; % plot boxplots?
+stimMark = 0; % decide how to plot the stimulus (no: 0, bar: 1)
+cStart = 1; % stimulus combination to start with
+oDisEloc = 0; % plot only DisNoAM and Eloc?
+saveI = 0; % save figure?
+
+% decide which set of rms-windows to plot
+switch filt
+        case {1,2,3}
+            switch cutTime
+                case 0
+                    uWin = (1:5);
+                    boxpFigWi = 800; % width of boxplot figure
+                case 1
+                    uWin = (4:4);
+                    boxpFigWi = 300; % width of boxplot figure
+                case 2
+                    uWin = (4:5);
+                    boxpFigWi = 500; % width of boxplot figure
+            end
+    case {4,5}
+        uWin = (1:3);
+        boxpFigWi = 500; % width of boxplot figure
+end
+% uWin = uWin+5;
+
+switch filt
+    case 1
+        filtName = 'noF';
+        y_lim = [-3.5,13];
+    case 2
+        filtName = '0.1F';
+        y_lim = [-3.5,13];
+    case 3
+        filtName = '10F';
+        y_lim = [-3.5,13];
+    case 4
+        filtName = '150F';
+        y_lim = [-4.5,4.5];
+    case 5
+        filtName = '300F';
+        y_lim = [-4,4];
+end
+
+% define colors and names
+colDev = [1,0,0]; % dev: red
+colStd = [0,0,1]; % std: blue
+col50Per = [1,0,1]; % 50 %: magenta
+colMR = [0,0,0]; % MR: black
+nameDev = 'Deviant';
+nameStd = 'Standard';
+name50Per = '50 % Control';
+nameMR = 'MS Control';
+
+%% load processed data
+
+switch inclCtrl
+    case 0
+        load('DD_statsData_noCtrl.mat') % load stats data
+    case 1 % load 50Per data and rename important variables
+        load('DD_avData_50Per.mat')
+        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_cell_Ctrl = filenames_cell;
+        load('DD_statsData_50Per.mat') % load stats data
+        nameCtrl = name50Per;
+        colCtrl = col50Per; % set control color to 50 % color
+    case 2 % load MR data and rename important variables
+        load('DD_avData_MR.mat')
+        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_cell_Ctrl = filenames_cell;
+        load('DD_statsData_MR.mat') % load stats data
+        nameCtrl = nameMR;
+        colCtrl = colMR; % set control color to MR color
+    case 3
+        load('DD_avData_50Per.mat')
+        dataAv_cell_50Per = dataAv_cell; % rename
+        dataGrAv_cell_50Per = dataGrAv_cell;
+        dataDifAv_cell_50Per = dataDifAv_cell; % rename
+        dataDifGrAv_cell_50Per = dataDifGrAv_cell;
+        dataSe_cell_50Per = dataSe_cell;
+        filenames_cell_50Per = filenames_cell;
+        load('DD_avData_MR.mat')
+        dataAv_cell_MR = dataAv_cell; % rename
+        dataGrAv_cell_MR = dataGrAv_cell;
+        dataDifAv_cell_MR = dataDifAv_cell; % rename
+        dataDifGrAv_cell_MR = dataDifGrAv_cell;
+        dataSe_cell_MR = dataSe_cell;
+        filenames_cell_MR = filenames_cell;
+        load('DD_statsData_bothCtrl.mat') % load stats data
+end
+
+% load Oddball data (will overwrite Ctrl variables that were not renamed)
+load('DD_avData_Oddball.mat')
+%% do some additional calculations
+
+switch zsNormI
+    case 0
+        ylableName = 'Voltage [µV]';
+    case 1
+        ylableName = 'z-norm. voltage [a.u.]';
+end
+
+switch cutTime
+    case 1
+        switch filt
+            case {1,2,3}
+                timeStart = -5; % starting time of the plot (relative to stim onset)
+                timeEnd = 15; % ending time of the plot (relative to stim onset)
+            case {4,5}
+                timeStart = -5; % starting time of the plot (relative to stim onset)
+                timeEnd = 15; % ending time of the plot (relative to stim onset)
+        end
+    case 2
+        timeStart = -5; % starting time of the plot (relative to stim onset)
+        timeEnd = 40; % ending time of the plot (relative to stim onset)        
+end
+
+% some calculations for matrix- and heatmap-plotting
+combNameS = split(combName(cStart:end),","); % split comb names at comma
+if plotMtrxI==1||plotHmI==1
+    oDisEloc = 0;
+    cStart = 2; % skip first stimulus combination (pure tones)
+    stimOrd = ["DisNoAM","DisAM","DisMimic","Eloc","ElocMimic"]; % define order of stimuli to plot in the matrix
+    nStim = size(stimOrd,2);
+    combCoord = zeros(nComb-1,2); % preallocate
+    for c = 1:nComb-1
+        for s = 1:2
+            combCoord(c,s) = find(ismember(stimOrd,combNameS(c,s)));
+        end
+    end
+    nTiles = 1:nStim*nStim; % number of tiles in matrix
+    tileMtrx = reshape(nTiles,[nStim,nStim])'; % create matrix to use for indexing later
+end
+
+if oDisEloc==1 % plot only stim combination 7 (DisNoAM/Eloc) if oDisEloc is 1
+    cStart = 7;
+    nComb = 7;
+end
+%% plotting: responses
+
+if plotRespI==1
+    switch plotMtrxI
+        case 1
+            figure('NumberTitle','off','Name','All Stims','Position',[0,0,2000,1300],'Renderer','painters')
+            t = tiledlayout(nStim,nStim);
+    end
+    for c = cStart:nComb % run once for each stimulus combination
+
+        % extract data corresponding to current stimulation-condition from
+        % cell-arrays
+        switch plotDifI
+            case 0
+                dataAv = dataAv_cell{c};
+                dataGrAv = dataGrAv_cell{c};
+                dataSe = dataSe_cell{c};
+                switch inclCtrl
+                    case {1,2}
+                        dataAvCtrl = dataAv_cell_Ctrl{c};
+                        dataGrAvCtrl = dataGrAv_cell_Ctrl{c};
+                        dataSeCtrl = dataSe_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
+                        dataAv50Per = dataAv_cell_50Per{c};
+                        dataGrAv50Per = dataGrAv_cell_50Per{c};
+                        dataSe50Per = dataSe_cell_50Per{c};
+                        filenames_50Per = filenames_cell_50Per{c};
+                        bsPos = cell2mat(strfind(filenames_50Per,'\')); % find positions of back slashes in filename
+                        startPos = bsPos(:,end); % use last back slash position as starting point
+                        filenames_50Per = extractAfter(filenames_50Per,startPos);
+                        dataAvMR = dataAv_cell_MR{c};
+                        dataGrAvMR = dataGrAv_cell_MR{c};
+                        dataSeMR = dataSe_cell_MR{c};
+                        filenames_MR = filenames_cell_MR{c};
+                        bsPos = cell2mat(strfind(filenames_MR,'\')); % find positions of back slashes in filename
+                        startPos = bsPos(:,end); % use last back slash position as starting point
+                        filenames_MR = extractAfter(filenames_MR,startPos);
+                end
+            case 1
+                dataAv = dataDifAv_cell_Ctrl{c};
+                dataGrAv = dataDifGrAv_cell_Ctrl{c};
+                plotWinI = 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
+                        dataAv50Per = dataAv_cell_50Per{c}; % won't be plotted, only processed to keep the code simpler
+                        dataGrAv50Per = dataGrAv_cell_50Per{c};
+                        filenames_50Per = filenames_cell_50Per{c};
+                        bsPos = cell2mat(strfind(filenames_50Per,'\')); % find positions of back slashes in filename
+                        startPos = bsPos(:,end); % use last back slash position as starting point
+                        filenames_50Per = extractAfter(filenames_50Per,startPos);
+                        dataAvMR = dataAv_cell_MR{c}; % won't be plotted, only processed to keep the code simpler
+                        dataGrAvMR = dataGrAv_cell_MR{c};
+                        filenames_MR = filenames_cell_MR{c};
+                        bsPos = cell2mat(strfind(filenames_MR,'\')); % find positions of back slashes in filename
+                        startPos = bsPos(:,end); % use last back slash position as starting point
+                        filenames_MR = extractAfter(filenames_MR,startPos);
+                end
+        end
+        filenames = filenames_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(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};
+        fsDwn = fsDwn_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*fsDwn:timeWin(2)/1000*fsDwn);
+        end
+
+        % subdivide responses for simplicity
+        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};
+        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
+                AAv50Per = dataAv50Per{1};
+                BAv50Per = dataAv50Per{2};
+                AGrAv50Per = dataGrAv50Per{1};
+                BGrAv50Per = dataGrAv50Per{2};
+                ASe50Per = dataSe50Per{1};
+                BSe50Per = dataSe50Per{2};
+                AAvMR = dataAvMR{1};
+                BAvMR = dataAvMR{2};
+                AGrAvMR = dataGrAvMR{1};
+                BGrAvMR = dataGrAvMR{2};
+                ASeMR = dataSeMR{1};
+                BSeMR = dataSeMR{2};
+        end
+
+
+        % plot single channel
+        for s = 1:2 % once for A, once for B
+            switch s
+                case 1
+                    grAvStd = AStdGrAv;
+                    grAvDev = ADevGrAv;
+                    avStd = AStdAv;
+                    avDev = ADevAv;
+                    seStd = AStdSe;
+                    seDev = ADevSe;
+                    stimName = 'A';
+                    switch inclCtrl
+                        case {1,2}
+                            grAvCtrl = AGrAvCtrl;
+                            avCtrl = AAvCtrl;
+                            seCtrl = ASeCtrl;
+                        case 3
+                            grAv50Per = AGrAv50Per;
+                            av50Per = AAv50Per;
+                            se50Per = ASe50Per;
+                            grAvMR = AGrAvMR;
+                            avMR = AAvMR;
+                            seMR = ASeMR;
+                    end
+                case 2
+                    grAvStd = BStdGrAv;
+                    grAvDev = BDevGrAv;
+                    avStd = BStdAv;
+                    avDev = BDevAv;
+                    seStd = BStdSe;
+                    seDev = BDevSe;
+                    stimName = 'B';
+                    switch inclCtrl
+                        case {1,2}
+                            grAvCtrl = BGrAvCtrl;
+                            avCtrl = BAvCtrl;
+                            seCtrl = BSeCtrl;
+                        case 3
+                            grAv50Per = BGrAv50Per;
+                            av50Per = BAv50Per;
+                            se50Per = BSe50Per;
+                            grAvMR = BGrAvMR;
+                            avMR = BAvMR;
+                            seMR = BSeMR;
+                    end
+            end
+            switch plotMtrxI
+                case 0
+                    figure('NumberTitle','off','Name',[combName{c},'_',stimName],'Position',[0,0,1000,600])
+                    title(combNameS(c,s))
+                case 1
+                    switch s % mirror coordinats with s
+                        case 1
+                            nexttile(tileMtrx(combCoord(c-1,2),combCoord(c-1,1)))
+                        case 2
+                            nexttile(tileMtrx(combCoord(c-1,1),combCoord(c-1,2)))
+                    end
+            end
+            switch avType
+                case 1 % grand averages
+                    hold on
+                    switch inclCtrl
+                        case {1,2}
+                            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
+                        case 3
+                            switch plotDifI
+                                case 0 % plot control response only when difference curve is deactivated
+                                    % plot standard error
+                                    seHi50Per = grAv50Per(timeCut,filt)+se50Per(timeCut,filt);
+                                    seLo50Per = grAv50Per(timeCut,filt)-se50Per(timeCut,filt);
+                                    seX = [timetr(timeCut),fliplr(timetr(timeCut))];
+                                    seY = [seHi50Per',fliplr(seLo50Per')];
+                                    patch(seX,seY,col50Per,'FaceAlpha',0.2,'EdgeColor','none')
+                                    % plot data
+                                    plot50Per = plot(timetr(timeCut),grAv50Per(timeCut,filt),'color',col50Per,'Linewidth',2);
+                                    % plot standard error
+                                    seHiMR = grAvMR(timeCut,filt)+seMR(timeCut,filt);
+                                    seLoMR = grAvMR(timeCut,filt)-seMR(timeCut,filt);
+                                    seX = [timetr(timeCut),fliplr(timetr(timeCut))];
+                                    seY = [seHiMR',fliplr(seLoMR')];
+                                    patch(seX,seY,colMR,'FaceAlpha',0.2,'EdgeColor','none')
+                                    % plot data
+                                    plotMR = plot(timetr(timeCut),grAvMR(timeCut,filt),'color',colMR,'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):5:timeWin(2))
+                    set(gca,'XTickLabel',timeStart:5:timeEnd)
+                    switch stimMark
+                        case 1 % plot bar
+                            y_dist = y_lim(2)-y_lim(1);
+                            y_min = y_lim(1)+y_dist*0.05;
+                            x_stim = [stimWin(1);stimWin(2);stimWin(2);stimWin(1)];
+                            y_stim = [y_min+y_dist*0.05;y_min+y_dist*0.05;y_min+y_dist*0.07;y_min+y_dist*0.07];
+                            patch('XData',x_stim,'YData',y_stim,'EdgeColor','Black','EdgeAlpha',1,'FaceColor','Black','FaceAlpha',1)
+                    end
+                    switch plotMtrxI
+                        case 0
+                            xlabel('Time [ms]','FontWeight','bold','FontSize',20)
+                            ylabel(ylableName,'FontWeight','bold','FontSize',20)
+                            switch inclCtrl
+                                case 0
+                                    legend([plotDev,plotStd],nameDev,nameStd,'Location','northwest','AutoUpdate','off')
+                                case {1,2}
+                                    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
+                                case 3
+                                    switch plotDifI
+                                        case 0
+                                            legend([plotDev,plotStd,plot50Per,plotMR],nameDev,nameStd,name50Per,nameMR,'Location','northwest','AutoUpdate','off')
+                                        case 1
+                                            legend([plotDev,plotStd],nameDev,nameStd,'Location','northwest','AutoUpdate','off')
+                                    end
+                            end
+                        case 1
+                            xlabel(t,'Time [ms]','FontWeight','bold','FontSize',20)
+                            ylabel(t,ylableName,'FontWeight','bold','FontSize',20)
+                    end
+                    set(gca,'FontSize',20,'Linewidth',2,'FontName','Arial')
+                    switch plotWinI
+                        case 1
+                            x_min = timetr(winAll(:,1));
+                            x_max = timetr(winAll(:,2));
+                            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
+                                switch inclCtrl
+                                    case 0
+                                        statsV = pV(c,s,filt,w);
+                                    case {1,2,3}
+                                        switch statsType
+                                            case 1 % use t-test statistics
+                                                statsV = pV(c,s,filt,w); 
+                                            case 2 % use ANOVA statistics
+                                                statsV = multcompV{c,s,filt,w}.pValue(1); % extract pValue for comparison between dev and std condition
+                                        end
+                                end
+                                if statsV<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}
+                                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
+                            case 3
+                                switch plotDifI
+                                    case 0 % plot control response only when difference curve is deactivated
+                                        plot50Per = plot(timetr(timeCut),av50Per(timeCut,f,filt),'color',col50Per,'Linewidth',2);
+                                        plotMR = plot(timetr(timeCut),avMR(timeCut,f,filt),'color',colMR,'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 stimMark
+                            case 1 % plot bar
+                                y_dist = y_lim(2)-y_lim(1);
+                                y_min = y_lim(1)+y_dist*0.05;
+                                x_stim = [stimWin(1);stimWin(2);stimWin(2);stimWin(1)];
+                                y_stim = [y_min+y_dist*0.05;y_min+y_dist*0.05;y_min+y_dist*0.07;y_min+y_dist*0.07];
+                                patch('XData',x_stim,'YData',y_stim,'EdgeColor','Black','EdgeAlpha',1,'FaceColor','Black','FaceAlpha',1)
+                        end
+                        xlabel('Time [ms]','FontWeight','bold','FontSize',20)
+                        ylabel(ylableName,'FontWeight','bold','FontSize',20)
+                        set(gca,'FontSize',20,'Linewidth',2,'FontName','Arial')
+                        xlim([timeWin(1),timeWin(2)])
+                        if fixedY==1
+                            ylim(y_lim)
+                        end
+                    end
+            end
+            if saveI==1
+                switch plotMtrxI
+                    case 0
+                        saveas(gcf,[combName{c},'_',stimName,'_',filtName,'_cutTime',num2str(cutTime),'_ctrl',num2str(inclCtrl),'.jpg'])
+                        saveas(gcf,[combName{c},'_',stimName,'_',filtName,'_cutTime',num2str(cutTime),'_ctrl',num2str(inclCtrl),'.svg'])
+                end
+            end
+        end
+    end
+    if saveI==1
+        switch plotMtrxI
+            case 1
+                saveas(gcf,['Matrix_',filtName,'_cutTime',num2str(cutTime),'_ctrl',num2str(inclCtrl),'.jpg'])
+                saveas(gcf,['Matrix_',filtName,'_cutTime',num2str(cutTime),'_ctrl',num2str(inclCtrl),'.svg'])
+        end
+    end
+end
+%% plotting: data-heatmap
+
+if plotHmI==1
+    % do some prior calculations
+    rgb1 = [(0:0.01:1),ones(1,101)]';
+    rgb2 = [(0:0.01:1),(1:-0.01:0)]';
+    rgb3 = [ones(1,101),(1:-0.01:0)]';
+    cMap = [rgb1,rgb2,rgb3];
+    % create matrix containing the data of interest
+    hmMtrx = NaN(nStim,nStim);
+
+    for c = cStart:nComb
+        for s = 1:2
+            for w = uWin
+                switch hmDataI
+                    case 1
+                        data = round(cohensD_V(c,s,filt,w,1),2); % Cohen's D of comparison between dev and std responses
+                    case 2
+                        data = round(rmsGrAv{c}(s,3,filt,w),2); % grand average RMS values of control response
+                end
+                switch s % mirror coordinats with s
+                    case 1
+                        hmMtrx(tileMtrx(combCoord(c-1,1),combCoord(c-1,2))) = data; % ATTENTION: coordinates change with s in the opposite way as when creating a tiled layout because matrix is filled up from top to bottom while tiled layout is filled from left to right
+                    case 2
+                        hmMtrx(tileMtrx(combCoord(c-1,2),combCoord(c-1,1))) = data;
+                end
+            end
+        end
+    end
+    % create heatmap
+    xvalues = stimOrd;
+    yvalue = stimOrd;
+    figure('NumberTitle','off','Name','Heatmap','Position',[0,0,1450,1300],'Renderer','painters')
+    h = heatmap(xvalues,yvalue,hmMtrx);
+    % do calculations for heatmap
+    cmap = colormap(h);
+    sCmap = size(cmap);
+    h.Title = 'Cohens D';
+    h.XLabel = 'Active Stimulus';
+    h.YLabel = 'Modulatory Stimulus';
+    h.MissingDataColor = [0.7,0.7,0.7];
+    lims = clim;
+    hmMax = max(abs(hmMtrx),[],'all');
+    switch hmDataI
+        case 1
+            hmMin = -hmMax;
+            h.Colormap = cMap;
+        case 2
+            hmMin = 0;
+            h.Colormap = summer;
+    end
+    clim([hmMin,hmMax])
+    set(gca,'FontSize',20,'FontName','Arial')
+    if saveI==1
+        saveas(gcf,['Heatmap_',filtName,'_cutTime',num2str(cutTime),'_ctrl',num2str(inclCtrl),'.jpg'])
+        saveas(gcf,['Heatmap_',filtName,'_cutTime',num2str(cutTime),'_ctrl',num2str(inclCtrl),'.svg'])
+    end
+end
+%% plotting: boxplots
+
+if plotBpI==1
+    switch plotMtrxI
+        case 1
+            figure('NumberTitle','off','Name','Boxpl_All Stims','Position',[0,0,1500,1300],'Renderer','painters')
+            t = tiledlayout(nStim,nStim);
+            title(t,'RMS','FontWeight','bold','FontSize',25)
+    end
+    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}
+            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
+        case 3
+            xDev = [1,6,11,16,21];
+            xStd = [2,7,12,17,22];
+            x50Per = [3,8,13,18,23];
+            xMR = [4,9,14,19,24];
+            sig_groups1 = {[1,2],[6,7],[11,12],[16,17],[21,22]}; % dev vs. std
+            sig_groups2 = {[1,3],[6,8],[11,13],[16,18],[21,23]}; % dev vs. 50Per
+            sig_groups3 = {[1,4],[6,9],[11,14],[16,19],[21,24]}; % dev vs. MR
+            sig_groups4 = {[2,3],[7,8],[12,13],[17,18],[22,23]}; % 50Per vs. std
+            sig_groups5 = {[2,4],[7,9],[12,14],[17,19],[22,24]}; % MR vs. std
+            sig_groups6 = {[3,4],[8,9],[13,14],[18,19],[23,24]}; % MR vs. 50Per
+    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(rmsV{c}(:,s,:,filt,:),[],'all')); % identify minimum RMS value for y-limit
+            rmsMax = ceil(max(rmsV{c}(:,s,:,filt,:),[],'all')); % identify maximum RMS value for y-limit
+            switch plotMtrxI
+                case 0
+                    figure('NumberTitle','off','Name',['Boxpl_',combName{c},'_',stimName],'Position',[0,0,boxpFigWi,400])
+                    title('RMS')
+                case 1
+                    switch s % mirror coordinats with s
+                        case 1
+                            nexttile(tileMtrx(combCoord(c-1,2),combCoord(c-1,1)))
+                        case 2
+                            nexttile(tileMtrx(combCoord(c-1,1),combCoord(c-1,2)))
+                    end
+            end
+            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(rmsV{c}(:,s,1,filt,w))),rmsV{c}(:,s,1,filt,w),'BoxFaceColor',colDev,'MarkerColor',colDev,'Linewidth',2);
+                boxStd = boxchart(xStd(w)*ones(size(rmsV{c}(:,s,2,filt,w))),rmsV{c}(:,s,2,filt,w),'BoxFaceColor',colStd,'MarkerColor',colStd,'Linewidth',2);
+                switch inclCtrl
+                    case 0
+                        H = sigstar(sig_groups{w},pV(c,s,filt,w),cohensD_V(c,s,filt,w,1));
+                    case {1,2}
+                        boxCtrl = boxchart(xCtrl(w)*ones(size(rmsV{c}(:,s,3,filt,w))),rmsV{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
+                                statsV = 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}},statsV,eS);
+                        end
+                    case 3
+                        box50Per = boxchart(x50Per(w)*ones(size(rmsV{c}(:,s,3,filt,w))),rmsV{c}(:,s,3,filt,w),'BoxFaceColor',col50Per,'MarkerColor',col50Per,'Linewidth',2);
+                        boxMR = boxchart(xMR(w)*ones(size(rmsV{c}(:,s,4,filt,w))),rmsV{c}(:,s,4,filt,w),'BoxFaceColor',colMR,'MarkerColor',colMR,'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
+                                statsV = multcompV{c,s,filt,w}.pValue([1,2,3,5,6,9]);
+                                eS = zeros(6,1);
+                                for cmp = 1:6
+                                    eS(cmp) = cohensD_V(c,s,filt,w,cmp);
+                                end
+                                H = sigstar({sig_groups1{w},sig_groups2{w},sig_groups3{w},sig_groups4{w},sig_groups5{w},sig_groups6{w}},statsV,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}
+                    xlim([xDev(uWin(1))-1,xCtrl(uWin(end))+1])
+                    xticks(xStd(uWin))
+%                     legend([boxDev,boxStd,boxCtrl],nameDev,nameStd,nameCtrl,'Location','northwest','AutoUpdate','off')
+                case 3
+                    xlim([xDev(uWin(1))-1,xMR(uWin(end))+1])
+                    xticks((xStd(uWin)+x50Per(uWin))/2)
+%                     legend([boxDev,boxStd,box50Per,boxMR],nameDev,nameStd,name50Per,nameMR,'Location','northwest','AutoUpdate','off')
+            end
+            switch plotMtrxI
+                case 0
+                    ylabel(ylableName,'FontWeight','bold','FontSize',20)
+                    xticklabels(winTxt(uWin))
+                case 1
+                    ylabel(t,ylableName,'FontWeight','bold','FontSize',20)
+                    xticks([])
+                    ylim([0,15])
+            end
+            set(gca,'FontSize',25,'Linewidth',2,'FontName','Arial')
+            if saveI==1
+                switch plotMtrxI
+                    case 0
+                        saveas(gcf,['Boxp_',combName{c},'_',stimName,'_',filtName,'_cutTime',num2str(cutTime),'_ctrl',num2str(inclCtrl),'.jpg'])
+                        saveas(gcf,['Boxp_',combName{c},'_',stimName,'_',filtName,'_cutTime',num2str(cutTime),'_ctrl',num2str(inclCtrl),'.svg'])
+                end
+            end
+        end
+    end
+    if saveI==1
+        switch plotMtrxI
+            case 1
+                saveas(gcf,['Boxpl_Matrix_',filtName,'_cutTime',num2str(cutTime),'_ctrl',num2str(inclCtrl),'.jpg'])
+                saveas(gcf,['Boxpl_Matrix_',filtName,'_cutTime',num2str(cutTime),'_ctrl',num2str(inclCtrl),'.svg'])
+        end
+    end
+end

DD_Processing.m

@@ -0,0 +1,490 @@
+%% define some variables
+
+% remove recordings
+remI = 0; % remove certain recordings from analysis? (0: no, 1: yes)
+remRec = [1,2,3,4,5,6,13]; % define recordings to be removed from analysis (if remI==1)
+
+% select processing options
+detrendI = 0; % detrend data? (1: yes, 0: no)
+zsNormI = 0; % z-normalise data data? (1: yes, 0: no)
+demeanI = 0; % demean data? (this is independent from baseline correction!)
+subsampleI = 1; % calculate averages from all responses (= 0) or only from dev after std and std before dev (= 1)?
+
+% filter settings
+type = 1; % 1: Butterworth, 2: IFFT
+order = 4;
+lowc = [0.1,40,150,300]; % low cut frequency in Hz
+hic = 2500; % high cut frequency in Hz
+nFilt = size(lowc,2)+1; % number of filters used (+1 because of unfiltered data)
+
+% rejection criterion
+artRejI = 0; % turn artefact rejection on (1) or off (0)
+rejThr = 150; % uV threshold (every trial with min/max value exceeding this will be rejected before averaging)
+
+% speaker distance
+spDis = 0.15; % distance from ear to speaker in m (to correct for sound-travel delay)
+
+% baseline correction settings
+blcWin = [-0.001,0]; % window used for baseline correction in s; postive values: time after stim onset, negative values: time before stim onset
+%% data processing
+
+for runType = 1:3 % run once for Oddball (1), once for 50 % (2) and once for MR dataset (3)
+    
+    combName = ["20000,60000";"DisAM,DisNoAM";"DisAM,DisMimic";"DisAM,Eloc";"DisAM,ElocMimic";...
+        "DisNoAM,DisMimic";"DisNoAM,Eloc";"DisNoAM,ElocMimic";...
+        "DisMimic,Eloc";"DisMimic,ElocMimic";...
+        "Eloc,ElocMimic"];
+    % load data
+    switch runType
+        case 1
+            load('DD_Data_Oddball.mat')
+            fileName = "DD_avData_Oddball";
+        case 2
+            load('DD_Data_50Per.mat')
+            fileName = "DD_avData_50Per";
+        case 3
+            load('DD_Data_MR.mat')
+            fileName = "DD_avData_MR";
+    end
+    nComb = size(combName,1); % number of different stimulus combination to be processed
+    
+    dataAv_cell = cell(1,nComb); % preallocate
+    dataGrAv_cell = cell(1,nComb); % preallocate
+    dataDifAv_cell = cell(1,nComb); % preallocate
+    dataDifGrAv_cell = cell(1,nComb); % preallocate
+    dataSe_cell = cell(1,nComb); % preallocate
+    dataTt_cell = cell(1,nComb); % preallocate
+    filenames_cell = cell(1,nComb); % preallocate
+    recID_cell = cell(1,nComb); % preallocate
+    fileI_Oddb_cell = cell(1,nComb); % preallocate
+    fileI_Ctrl_cell = cell(1,nComb); % preallocate
+    timetr_cell = cell(1,nComb); % preallocate
+    pts2begin_cell = cell(1,nComb); % preallocate
+    nBlcks_cell = cell(1,nComb); % preallocate
+    nFiles_cell = cell(1,nComb); % preallocate
+    nFilt_cell = cell(1,nComb); % preallocate
+    nDevUsed_cell = cell(1,nComb); % preallocate
+    nStdUsed_cell = cell(1,nComb); % preallocate
+    stimDur_cell = cell(1,nComb); % preallocate
+    stimDelay_cell = cell(1,nComb); % preallocate
+    fsDwn_cell = cell(1,nComb); % preallocate
+
+    for c = 1:nComb % run analysis once for each stimulation-condition
+    % for c = 7:7 % run analysis once for each stimulation-condition    
+        switch runType
+            case 1
+                switch c % switch dataset depending on loop iteration
+                    case 1
+                        data = Data_Oddball_20_60;
+                    case 2
+                        data = Data_Oddball_DisAM_DisNoAM;
+                    case 3
+                        data = Data_Oddball_DisAM_DisMimic;
+                    case 4
+                        data = Data_Oddball_DisAM_Eloc;
+                    case 5
+                        data = Data_Oddball_DisAM_ElocMimic;
+                    case 6
+                        data = Data_Oddball_DisNoAM_DisMimic;
+                    case 7
+                        data = Data_Oddball_DisNoAM_Eloc;
+                    case 8
+                        data = Data_Oddball_DisNoAM_ElocMimic;
+                    case 9
+                        data = Data_Oddball_DisMimic_Eloc;
+                    case 10
+                        data = Data_Oddball_DisMimic_ElocMimic;
+                    case 11
+                        data = Data_Oddball_Eloc_ElocMimic;
+                end
+            case 2
+                switch c % switch dataset depending on loop iteration
+                    case 1
+                        data = Data_50Per_20_60;
+                    case 2
+                        data = Data_50Per_DisAM_DisNoAM;
+                    case 3
+                        data = Data_50Per_DisAM_DisMimic;
+                    case 4
+                        data = Data_50Per_DisAM_Eloc;
+                    case 5
+                        data = Data_50Per_DisAM_ElocMimic;
+                    case 6
+                        data = Data_50Per_DisNoAM_DisMimic;
+                    case 7
+                        data = Data_50Per_DisNoAM_Eloc;
+                    case 8
+                        data = Data_50Per_DisNoAM_ElocMimic;
+                    case 9
+                        data = Data_50Per_DisMimic_Eloc;
+                    case 10
+                        data = Data_50Per_DisMimic_ElocMimic;
+                    case 11
+                        data = Data_50Per_Eloc_ElocMimic;
+                end
+            case 3
+                switch c
+                    case 1 % 20000,60000
+                        data = Data_MR_PureTones;
+                        AfrqI = 1; % frequency index of A in MR control 
+                        BfrqI = 7; % frequency index of B in MR control
+                    case 2 % DisAM,DisNoAM
+                        data = Data_MR_Vocs;
+                        AfrqI = 2; % frequency index of DisAM-Response in MR control 
+                        BfrqI = 3; % frequency index of DisNoAM-Response in MR control
+                    case 3 % DisAM,DisMimic
+                        data = Data_MR_Vocs;
+                        AfrqI = 2; % frequency index of DisAM-Response in MR control 
+                        BfrqI = 10; % frequency index of DisMimic-Response in MR control
+                    case 4 % DisAM,Eloc
+                        data = Data_MR_Vocs;
+                        AfrqI = 2; % frequency index of DisAM-Response in MR control 
+                        BfrqI = 1; % frequency index of Eloc-Response in MR control
+                    case 5 % DisAM,ElocMimic
+                        data = Data_MR_Vocs;
+                        AfrqI = 2; % frequency index of DisAM-Response in MR control 
+                        BfrqI = 9; % frequency index of ElocMimic-Response in MR control
+                    case 6 % DisNoAM,DisMimic
+                        data = Data_MR_Vocs;
+                        AfrqI = 3; % frequency index of DisNoAM-Response in MR control 
+                        BfrqI = 10; % frequency index of DisMimic-Response in MR control
+                    case 7 % DisNoAM,Eloc
+                        data = Data_MR_Vocs;
+                        AfrqI = 3; % frequency index of DisNoAM-Response in MR control 
+                        BfrqI = 1; % frequency index of Eloc-Response in MR control
+                    case 8 % DisNoAM,ElocMimic
+                        data = Data_MR_Vocs;
+                        AfrqI = 3; % frequency index of DisNoAM-Response in MR control 
+                        BfrqI = 9; % frequency index of ElocMimic-Response in MR control
+                    case 9 % DisMimic,Eloc
+                        data = Data_MR_Vocs;
+                        AfrqI = 10; % frequency index of DisMimic-Response in MR control 
+                        BfrqI = 1; % frequency index of Eloc-Response in MR control
+                    case 10 % DisMimic,ElocMimic
+                        data = Data_MR_Vocs;
+                        AfrqI = 10; % frequency index of DisMimic-Response in MR control 
+                        BfrqI = 9; % frequency index of ElocMimic-Response in MR control
+                    case 11 % Eloc,ElocMimic
+                        data = Data_MR_Vocs;
+                        AfrqI = 1; % frequency index of Eloc-Response in MR control 
+                        BfrqI = 9; % frequency index of ElocMimic-Response in MR control
+                end
+        end
+        
+        
+        % Organise dataset in multiple variables (in separate script)
+        DD_OrgData
+        
+        % cut out buffers from raw data and split data into step-blocks
+        recRawCut = zeros(blockLen,nBlcks,nFiles); % preallocate
+        for f = 1:nFiles
+            for b = 1:nBlcks
+                recRawCut(:,b,f) = rec_raw{f}(1+(b-1)*blockLenBuff(f):blockLen+(b-1)*blockLenBuff(f));
+                recRawCut(:,b,f) =recRawCut(:,b,f)*96; % convert values from points to µV 
+            end
+        end
+        
+        % artefact rejection
+        recRawCut = reshape(recRawCut,pSmpl,nTrials,nBlcks,nFiles); % reshape to have one trial in each column
+        nTrialsAcc = zeros(nBlcks,nFiles); % preallocate
+        seqAcc = zeros(nTrials,nBlcks,nFiles); % preallocate
+        switch artRejI
+            case 0 % if artefact rejection is turned off
+                nTrialsAcc(:) = nTrials;
+                for f = 1:nFiles
+                    for b = 1:nBlcks
+                        seqAcc(:,b,f) = seq(f,:);
+                    end
+                end
+            case 1 % run if artefact rejection is turned on
+                % preallocate
+                maxV = zeros(1,nTrials,nBlcks,nFiles);
+                rejIdx = zeros(1,nTrials,nBlcks,nFiles);
+                for f = 1:nFiles
+                    for b = 1:nBlcks
+                        for t = 1:nTrials
+                            maxV(1,t,b,f) = max(abs(recRawCut(:,t,b,f))); % find (absolute) max value in each trial for unfiltered data (hence "1" in 5th dimension)
+                            rejIdx(1,t,b,f) = maxV(:,t,b,f)<rejThr; % compare max value of trial with rejection threshold, "1" if max is smaller than threshold -> trial will be accepted later
+                        end % repeat for every trial in current file and (level-)step
+                        nTrialsAcc(b,f) = sum(rejIdx(:,:,b,f)); % count number of accepted trials for current file and step
+                        recRawCut(:,1:nTrialsAcc(b,f),b,f) = recRawCut(:,logical(rejIdx(:,:,b,f)),b,f); % apply previously created rejection index-array as logical array on 2nd dimension (trials) of filtered data. Indexing is applied on all filters (5th dim) simultaneously; array is filled only until number of accepted responses is reached (3rd dim), the rest is still filled with zeros
+                        seqAcc(1:nTrialsAcc(b,f),b,f) = seq(f,logical(rejIdx(:,:,b,f)));
+                    end
+                end
+        end
+        
+        % process data
+        recPro = recRawCut; % rename
+        recFilt = zeros(pSmpl*nTrials,nBlcks,nFiles,nFilt-1); % preallocate
+        for f = 1:nFiles
+            for b = 1:nBlcks
+                % detrending
+                if detrendI==1
+                    recPro(:,1:nTrialsAcc(b,f),b,f) = detrend(recPro(:,1:nTrialsAcc(b,f),b,f));
+                end
+                % z score normalisation
+                if zsNormI==1
+                    recPro(:,1:nTrialsAcc(b,f),b,f) = zscoreJo(recPro(:,1:nTrialsAcc(b,f),b,f),pSmpl);
+                end
+                % demeaning
+                if demeanI==1
+                    recPro(:,1:nTrialsAcc(b,f),b,f) = demean(recPro(:,1:nTrialsAcc(b,f),b,f));
+                end
+                % filtering
+                recProRS = reshape(recPro,nTrials*pSmpl,nBlcks,nFiles);
+                for ft = 1:nFilt-1 % "-1" because "nFilt" includes raw-array which must be excluded here
+                    recFilt(1:nTrialsAcc(b,f)*pSmpl,b,f,ft) = Filter_Butter(recProRS(1:nTrialsAcc(b,f)*pSmpl,b,f),order,lowc(ft),hic,"bandpass",fsDwn);
+                end
+            end
+        end
+        recProFilt = cat(4,recProRS,recFilt); % concatenate unfiltered and filtered data in 4th dimension
+        recProFilt = reshape(recProFilt,pSmpl,nTrials,nBlcks,nFiles,nFilt); % reshape to split raw-data into separate trials (2nd dimension)
+        
+        % baseline correction
+        blcStartPts = round(pts2begin(1)+blcWin(1)*fsDwn); % transform time window for baseline correction into points
+        blcEndPts = round(pts2begin(1)+blcWin(2)*fsDwn);
+        for f = 1:nFiles
+            for b = 1:nBlcks
+                for ft = 1:nFilt
+                    recProFilt(:,1:nTrialsAcc(b,f),b,f,ft) = blCorr(recProFilt(:,1:nTrialsAcc(b,f),b,f,ft),blcStartPts,blcEndPts);
+                end
+            end
+        end
+        
+        % split blocks into deviant and standard responses
+        [devR,stdR,nDevAcc,nStdAcc] = sortBlck(recProFilt,seqAcc,nTrialsAcc,AfrqI,BfrqI,runType);
+        nDevUsed = nDevAcc; % define number of accepted deviants as number of used deviants --> changes later again if subsampling is on!
+        nStdUsed = nStdAcc; % define number of accepted standards as number of used standards --> changes later again if subsampling is on!
+        
+        % prepare data for time trace analysis (different script)
+        devAll = devR; % rename all accepted dev and std responses
+        stdAll = stdR;
+        devTt = zeros(pSmpl,100,nBlcks,nFiles,nFilt); % preallocate
+        stdTt = zeros(pSmpl,100,nBlcks,nFiles,nFilt); % preallocate
+        stdTtRs = zeros(pSmpl,9,100,nBlcks,nFiles,nFilt); % preallocate
+%         stdTtMean = zeros(pSmpl,1,100,nBlcks,nFiles,nFilt); % preallocate
+        stdTt2 = zeros(pSmpl,100,nBlcks,nFiles,nFilt); % preallocate
+        switch runType % change depending on Oddball/Ctrl
+            case 1 % oddball
+                for f = 1:nFiles
+                    for b = 1:nBlcks
+                        devTt(:,:,b,f,:) = devAll(:,1:100,b,f,:); % use all dev responses
+                        stdTt(:,:,b,f,:) = stdAll(:,1:9:900,b,f,:); % use only every 9th std response (to have the same number as dev responses)
+                        stdTtRs(:,:,:,b,f,:) = reshape(stdAll(:,1:900,b,f,:),pSmpl,9,100,1,1,nFilt); % reshape to have 9 consecutive std responses in one dimension
+                        stdTt2(:,:,b,f,:) = mean(stdTtRs(:,:,:,b,f,:),2); % average 9 consecutive std responses and save as stdTt2
+                    end
+                end
+            case 2 % 50Per
+                for f = 1:nFiles
+                    for b = 1:nBlcks
+                        devTt(:,:,b,f,:) = devAll(:,1:5:500,b,f,:); % use only every 5th 50Per response (to have the same number als dev responses)
+                        stdTt(:,:,b,f,:) = stdAll(:,1:5:500,b,f,:);
+                    end
+                end
+            case 3 % MR
+                for f = 1:nFiles
+                    for b = 1:nBlcks
+                        devTt(:,:,b,f,:) = devAll(:,1:100,b,f,:); % use 100 MR responses; ATTENTION: Number of MR responses varies between 95 and 105 -> use only 90 responses in final analysis
+                        stdTt(:,:,b,f,:) = stdAll(:,1:100,b,f,:);
+                    end
+                end
+        end
+        devTtAv = reshape(mean(devTt,4),pSmpl,100,nBlcks,nFilt); % average responses across individuals to be saved and further processed in another script
+        stdTtAv = reshape(mean(stdTt,4),pSmpl,100,nBlcks,nFilt);
+        stdTt2Av = reshape(mean(stdTt2,4),pSmpl,100,nBlcks,nFilt);
+        
+        % subsample oddball responses (extract dev after std & std before
+        % dev) and reassign variables "nDevUsed" & "nStdUsed"
+        % ATTENTION: this section only works correctly if the same files
+        % are available for Oddball and Control conditions!
+        switch runType
+            case 1 % if oddball run is processed
+                switch subsampleI
+                    case 1
+                        [devR,stdR,nTrialsSub] = subsample(devR,stdR,seqAcc,nTrialsAcc); % subsample oddball responses --> number of trials will be much less afterwards (= nTrialsSub)
+                        nDevUsed = nTrialsSub; % define number of subsampled deviants as number of used deviants
+                        nStdUsed = nTrialsSub; % define number of subsampled standards as number of used standards
+                        nDevUsed_oddb{c} = nDevUsed; % assign number of used deviants to separate variable. (std not necessary bc ctrls will only need number of deviants. Additionally, nStdAcc and nDevAcc should be equal anyway if subsampling is on)
+                end
+            case {2,3} % if control run is processed
+                switch subsampleI
+                    case 1
+                        switch runType
+                            case 2 % 50Per
+                                nDevUsed = nDevUsed_oddb{c}; % use respective number of used oddball deviants as number of used control deviants
+                                nStdUsed = nDevUsed_oddb{c}; % ATTENTION: use respective number of used oddball deviants as number of used control standards
+                        end
+                        % in case of MR, use all accepted dev and std
+                        % responses
+                end
+                % run random perutation to shuffle responses --> important for 50Per
+                % control: to make sure not only the earliest responses are used
+                % for averaging
+                for f = 1:nFiles
+                    for b = 1:nBlcks
+                        rDev = randperm(nDevAcc(b,f)); % create random vector containing all ACCEPTED dev responses (they cover the whole sequence)
+                        devR(:,1:nDevAcc(b,f),b,f,:) = devR(:,rDev,b,f,:); % shuffle those vectors that contain an actual response and no zeros --> timing within sequence is now irrelevant
+                        rStd = randperm(nStdAcc(b,f)); % repeat the same for std responses
+                        stdR(:,1:nStdAcc(b,f),b,f,:) = stdR(:,rStd,b,f,:);
+                    end
+                end
+        end
+
+        % calculate average, grand-average & standard error of sorted responses
+        devAv = zeros(pSmpl,nBlcks,nFiles,nFilt);
+        stdAv = zeros(pSmpl,nBlcks,nFiles,nFilt);
+        for f = 1:nFiles
+            for b = 1:nBlcks
+                devAv(:,b,f,:) = mean(devR(:,1:nDevUsed(b,f),b,f,:),2); % averaging is performed only on those vectors that contain an actual response and not only zeros (2nd dim) --> only those (meaningful) vectors have been mixed before
+                stdAv(:,b,f,:) = mean(stdR(:,1:nStdUsed(b,f),b,f,:),2);
+            end
+        end
+        % calculate grand average
+        devGrAv = reshape(mean(devAv,3),pSmpl,nBlcks,nFilt);
+        stdGrAv = reshape(mean(stdAv,3),pSmpl,nBlcks,nFilt);
+        % calculate standard deviation
+        devStdD = reshape(std(devAv,0,3),pSmpl,nBlcks,nFilt);
+        stdStdD = reshape(std(stdAv,0,3),pSmpl,nBlcks,nFilt);
+        % calculate standard error
+        devSe = devStdD/sqrt(nFiles);
+        stdSe = stdStdD/sqrt(nFiles);
+        
+        % split dataset into several sub-variables for simplisity
+        switch runType
+            case 1
+                % sorted responses
+                ADev{c} = devAv(:,1,:,:); % A is deviant in block 1 (AB)
+                ADev{c} = reshape(ADev{c},pSmpl,nFiles,nFilt); % reshape av data into 3D array
+                AStd{c} = stdAv(:,2,:,:); % A is standard in block 2 (BA)
+                AStd{c} = reshape(AStd{c},pSmpl,nFiles,nFilt); % reshape av data into 3D array
+                BDev{c} = devAv(:,2,:,:); % B is deviant in block 2 (BA)
+                BDev{c} = reshape(BDev{c},pSmpl,nFiles,nFilt); % reshape av data into 3D array
+                BStd{c} = stdAv(:,1,:,:); % B is standard in block 1 (AB)
+                BStd{c} = reshape(BStd{c},pSmpl,nFiles,nFilt); % reshape av data into 3D array
+                ADevGrAv{c} = devGrAv(:,1,:);
+                ADevGrAv{c} = reshape(ADevGrAv{c},pSmpl,nFilt); % reshape grav into 2D array
+                AStdGrAv{c} = stdGrAv(:,2,:);
+                AStdGrAv{c} = reshape(AStdGrAv{c},pSmpl,nFilt); % reshape grav data into 2D array
+                BDevGrAv{c} = devGrAv(:,2,:);
+                BDevGrAv{c} = reshape(BDevGrAv{c},pSmpl,nFilt); % reshape grav data into 2D array
+                BStdGrAv{c} = stdGrAv(:,1,:);
+                BStdGrAv{c} = reshape(BStdGrAv{c},pSmpl,nFilt); % reshape grav data into 2D array
+                ADevSe{c} = devSe(:,1,:);
+                ADevSe{c} = reshape(ADevSe{c},pSmpl,nFilt); % reshape standard error into 2D array
+                AStdSe{c} = stdSe(:,2,:);
+                AStdSe{c} = reshape(AStdSe{c},pSmpl,nFilt); % reshape standard error into 2D array
+                BDevSe{c} = devSe(:,2,:);
+                BDevSe{c} = reshape(BDevSe{c},pSmpl,nFilt); % reshape standard error into 2D array
+                BStdSe{c} = stdSe(:,1,:);
+                BStdSe{c} = reshape(BStdSe{c},pSmpl,nFilt); % reshape standard error into 2D array
+                ADevTt{c} = devTtAv(:,:,1,:); % A is deviant in block 1 (AB)
+                ADevTt{c} = reshape(ADevTt{c},pSmpl,100,nFilt); % reshape av data into 3D array
+                AStdTt{c} = stdTtAv(:,:,2,:); % A is standard in block 2 (BA)
+                AStdTt{c} = reshape(AStdTt{c},pSmpl,100,nFilt); % reshape av data into 3D array
+                BDevTt{c} = devTtAv(:,:,2,:); % B is deviant in block 2 (BA)
+                BDevTt{c} = reshape(BDevTt{c},pSmpl,100,nFilt); % reshape av data into 3D array
+                BStdTt{c} = stdTtAv(:,:,1,:); % B is standard in block 1 (AB)
+                BStdTt{c} = reshape(BStdTt{c},pSmpl,100,nFilt); % reshape av data into 3D array
+                AStdTt2{c} = stdTt2Av(:,:,2,:); % A is standard in block 2 (BA)
+                AStdTt2{c} = reshape(AStdTt2{c},pSmpl,100,nFilt); % reshape av data into 3D array
+                BStdTt2{c} = stdTt2Av(:,:,1,:); % B is standard in block 1 (AB)
+                BStdTt2{c} = reshape(BStdTt2{c},pSmpl,100,nFilt); % reshape av data into 3D array
+            case {2,3}
+                switch runType
+                    case 2
+                        % sorted responses
+                        ACtrl{c} = devAv(:,1,:,:); % A is deviant in block 1 (AB) --> in this case (50Per), it does not matter which condition is used, as long as it is the right stimulus
+                        ACtrl{c} = reshape(ACtrl{c},pSmpl,nFiles,nFilt); % reshape av data into 3D array
+                        BCtrl{c} = devAv(:,2,:,:); % B is deviant in block 2 (BA) --> in this case (50Per), it does not matter which condition is used, as long as it is the right stimulus
+                        BCtrl{c} = reshape(BCtrl{c},pSmpl,nFiles,nFilt); % reshape av data into 3D array
+                        ACtrlGrAv{c} = devGrAv(:,1,:);
+                        ACtrlGrAv{c} = reshape(ACtrlGrAv{c},pSmpl,nFilt); % reshape grav into 2D array
+                        BCtrlGrAv{c} = devGrAv(:,2,:);
+                        BCtrlGrAv{c} = reshape(BCtrlGrAv{c},pSmpl,nFilt); % reshape grav data into 2D array
+                        ACtrlSe{c} = devSe(:,1,:);
+                        ACtrlSe{c} = reshape(ACtrlSe{c},pSmpl,nFilt); % reshape standard error into 2D array
+                        BCtrlSe{c} = devSe(:,2,:);
+                        BCtrlSe{c} = reshape(BCtrlSe{c},pSmpl,nFilt); % reshape standard error into 2D array
+                        ACtrlTt{c} = devTtAv(:,:,1,:); % A is deviant in block 1 (AB) --> in this case (50Per), it does not matter which condition is used, as long as it is the right stimulus
+                        ACtrlTt{c} = reshape(ACtrlTt{c},pSmpl,100,nFilt); % reshape av data into 3D array
+                        BCtrlTt{c} = devTtAv(:,:,2,:); % B is deviant in block 2 (BA) --> in this case (50Per), it does not matter which condition is used, as long as it is the right stimulus
+                        BCtrlTt{c} = reshape(BCtrlTt{c},pSmpl,100,nFilt); % reshape av data into 3D array
+                    case 3
+                        % sorted responses
+                        ACtrl{c} = devAv; % A was treated as deviant for calculations
+                        ACtrl{c} = reshape(ACtrl{c},pSmpl,nFiles,nFilt); % reshape av data into 3D array
+                        BCtrl{c} = stdAv; % B was treated as standard for calculations
+                        BCtrl{c} = reshape(BCtrl{c},pSmpl,nFiles,nFilt); % reshape av data into 3D array
+                        ACtrlGrAv{c} = devGrAv;
+                        ACtrlGrAv{c} = reshape(ACtrlGrAv{c},pSmpl,nFilt); % reshape grav into 2D array
+                        BCtrlGrAv{c} = stdGrAv;
+                        BCtrlGrAv{c} = reshape(BCtrlGrAv{c},pSmpl,nFilt); % reshape grav into 2D array
+                        ACtrlSe{c} = devSe;
+                        ACtrlSe{c} = reshape(ACtrlSe{c},pSmpl,nFilt); % reshape standard error into 2D array
+                        BCtrlSe{c} = stdSe;
+                        BCtrlSe{c} = reshape(BCtrlSe{c},pSmpl,nFilt); % reshape standard error into 2D array
+                        ACtrlTt{c} = devTtAv; % A was treated as deviant for calculations
+                        ACtrlTt{c} = reshape(ACtrlTt{c},pSmpl,100,nFilt); % reshape av data into 3D array
+                        BCtrlTt{c} = stdTtAv; % B was treated as standard for calculations
+                        BCtrlTt{c} = reshape(BCtrlTt{c},pSmpl,100,nFilt); % reshape av data into 3D array
+                end
+                % calculate difference curves
+                ADifDev{c} = ADev{c}-ACtrl{c};
+                ADifStd{c} = AStd{c}-ACtrl{c};
+                BDifDev{c} = BDev{c}-BCtrl{c};
+                BDifStd{c} = BStd{c}-BCtrl{c};
+                ADifDevGrAv{c} = ADevGrAv{c}-ACtrlGrAv{c};
+                ADifStdGrAv{c} = AStdGrAv{c}-ACtrlGrAv{c};
+                BDifDevGrAv{c} = BDevGrAv{c}-BCtrlGrAv{c};
+                BDifStdGrAv{c} = BStdGrAv{c}-BCtrlGrAv{c};
+        end
+        
+        
+        % organise processed dev and std responses in cell array
+        switch runType
+            case 1
+                dataAv= {ADev{c},AStd{c},BDev{c},BStd{c}};
+                dataGrAv = {ADevGrAv{c},AStdGrAv{c},BDevGrAv{c},BStdGrAv{c}};
+                dataDifAv = [];
+                dataDifGrAv = [];
+                dataSe = {ADevSe{c},AStdSe{c},BDevSe{c},BStdSe{c}};
+                dataTt = {ADevTt{c},AStdTt{c},BDevTt{c},BStdTt{c},AStdTt2{c},BStdTt2{c}};
+            case {2,3}
+                dataAv = {ACtrl{c},BCtrl{c}};
+                dataGrAv = {ACtrlGrAv{c},BCtrlGrAv{c}};
+                dataDifAv = {ADifDev{c},ADifStd{c},BDifDev{c},BDifStd{c}};
+                dataDifGrAv = {ADifDevGrAv{c},ADifStdGrAv{c},BDifDevGrAv{c},BDifStdGrAv{c}};
+                dataSe = {ACtrlSe{c},BCtrlSe{c}};
+                dataTt = {ACtrlTt{c},BCtrlTt{c}};
+        end
+        
+        % concatenate important data in another cell array with one cell
+        % for each stimulation-condition
+        dataAv_cell{c} = dataAv;
+        dataGrAv_cell{c} = dataGrAv;
+        dataDifAv_cell{c} = dataDifAv;
+        dataDifGrAv_cell{c} = dataDifGrAv;
+        dataSe_cell{c} = dataSe;
+        dataTt_cell{c} = dataTt;
+        filenames_cell{c} = filenames;
+        recID_cell{c} = recID;
+        fileI_Oddb_cell{c} = fileI_Oddb;
+        fileI_Ctrl_cell{c} = fileI_Ctrl;
+        timetr_cell{c} = timetr;
+        pts2begin_cell{c} = pts2begin;
+        nBlcks_cell{c} = nBlcks;
+        nFiles_cell{c} = nFiles;
+        nFilt_cell{c} = nFilt;
+        stimDur_cell{c} = stimDur;
+        stimDelay_cell{c} = stimDelay;
+        fsDwn_cell{c} = fsDwn;
+        nDevUsed_cell{c} = nDevUsed;
+        nStdUsed_cell{c} = nStdUsed;
+    end
+    save(fileName,'dataAv_cell','dataGrAv_cell','dataDifAv_cell',...
+        'dataDifGrAv_cell','dataSe_cell','dataTt_cell','filenames_cell','recID_cell',...
+        'fileI_Oddb_cell','fileI_Ctrl_cell','timetr_cell','pts2begin_cell','combName',...
+        'nComb','nBlcks_cell','nFiles_cell','nFilt_cell','nDevUsed_cell','nStdUsed_cell',...
+        'stimDur_cell','stimDelay_cell','fsDwn_cell','subsampleI',...
+        'artRejI','zsNormI','detrendI','demeanI')
+end

DD_Stats.m

@@ -0,0 +1,178 @@
+%% define some variables
+
+% define time windows to be analysed (in s); ATTENTION: compared to
+% previous paper, all values are -0.4 ms, since now a correction for sound
+% travelling delay has been performed (as opposed to before)
+
+% Eloc
+win1 = [0,0.001]; % ABR (i)
+win2 = [0.001,0.0022]; % ABR (ii/iii)
+win3 = [0.0022,0.0034]; % ABR (iv)
+win4 = [0,0.01]; % ABR (whole)
+win5 = [0.01,0.02]; % ABR (late)
+% DisNoAM
+win6 = [0,0.0013]; % ABR (i)
+win7 = [0.0013,0.0025]; % ABR (ii/iii)
+win8 = [0.0025,0.0041]; % ABR (iv)
+win9 = [0,0.01]; % ABR (whole)
+win10 = [0.01,0.024]; % ABR (late)
+
+winTxt = ["i","ii/iii","iv","ABR (whole)","ABR (late)"];
+
+inclCtrl = 3; % include control data? (0: no, 1: 50 %, 2: MR, 3: both controls)
+%% load datasets
+
+switch inclCtrl
+    case {1,2} % load 50Per data and rename important variables
+        switch inclCtrl
+            case 1
+                load('DD_avData_50Per.mat')
+            case 2
+                load('DD_avData_MR.mat')
+        end
+        dataAv_cell_ctrl = dataAv_cell; % rename
+        fileI_Oddb_cell_ctrl = fileI_Oddb_cell; % rename
+        fileI_Ctrl_cell_ctrl = fileI_Ctrl_cell; % rename
+        dataTt_cell_ctrl = dataTt_cell; % rename
+    case 3
+        load('DD_avData_50Per.mat')
+        dataAv_cell_50Per = dataAv_cell; % rename
+        fileI_Oddb_cell_50Per = fileI_Oddb_cell; % rename
+        fileI_Ctrl_cell_50Per = fileI_Ctrl_cell; % rename
+        load('DD_avData_MR.mat')
+        dataAv_cell_MR = dataAv_cell; % rename
+        fileI_Oddb_cell_MR = fileI_Oddb_cell; % rename
+        fileI_Ctrl_cell_MR = fileI_Ctrl_cell; % rename
+end
+
+% load oddball data
+load('DD_avData_Oddball.mat')
+%% do some additional calulations that are required
+
+fs = fsDwn_cell{1}(1);
+pts2begin = pts2begin_cell{1}(1);
+winAll = round([win1;win2;win3;win4;win5;win6;win7;win8;win9;win10]*fs); % combine time windows and convert them into points
+winAll = winAll+pts2begin; % correct for stim delay
+nWin = size(winAll,1);
+
+nStims = 2; % number of different stimuli per combination
+nFilt = nFilt_cell{1};
+
+switch inclCtrl % number of different stimulus condition (depends on whether MR is included or not)
+    case 0
+        filename = 'DD_statsData_noCtrl';
+        nCond = 2; % 2 conditions: dev and std
+        dataAv = dataAv_cell; % only oddball data
+    case {1,2}
+        switch inclCtrl
+            case 1
+                filename = 'DD_statsData_50Per';
+            case 2
+                filename = 'DD_statsData_MR';
+        end
+        nCond = 3; % 3 conditions: dev, std and control
+        dataAv = cell(size(dataAv_cell)); % preallocate
+        for c = 1:nComb % once for each stimulus combination
+            for s = 1:nStims % once for each stimulus frequency
+                for cn = 1:nCond % once for each stimulus condition (dev, std, ctrl)
+                    switch cn % change dataset depending on condition (Oddball vs. control)
+                        case {1,2} % if oddball
+                            dataAv{c}{(s-1)*nCond+cn} = dataAv_cell{c}{(s-1)*2+cn}; % concatenate oddball data of respective stimulus & condition (dev/std)
+                        case 3 % if control
+                            dataAv{c}{(s-1)*nCond+cn} = dataAv_cell_ctrl{c}{s}; % concatenate control data of respective stimulus
+                    end
+                end
+            end
+        end
+    case 3
+        filename = 'DD_statsData_bothCtrl';
+        nCond = 4; % 4 conditions: dev, std, 50Per and MR
+        dataAv = cell(size(dataAv_cell)); % preallocate
+        for c = 1:nComb % once for each stimulus combination
+            for s = 1:nStims % once for each stimulus frequency
+                for cn = 1:nCond % once for each stimulus condition (dev, std, ctrl)
+                    switch cn % change dataset depending on condition (Oddball vs. control)
+                        case {1,2} % if oddball
+                            dataAv{c}{(s-1)*nCond+cn} = dataAv_cell{c}{(s-1)*2+cn}; % concatenate oddball data of respective stimulus & condition (dev/std)
+                        case 3 % if control
+                            dataAv{c}{(s-1)*nCond+cn} = dataAv_cell_50Per{c}{s}; % concatenate control data of respective stimulus
+                        case 4
+                            dataAv{c}{(s-1)*nCond+cn} = dataAv_cell_MR{c}{s}; % concatenate control data of respective stimulus
+                    end
+                end
+            end
+        end
+end
+%% calulate RMS values within defined time windows
+
+rmsV = cell(1,nComb); % preallocate
+rmsGrAv = cell(1,nComb); % preallocate
+for c = 1:nComb % once for each stimulus combination
+    nFiles = nFiles_cell{c};
+    rmsV{c} = zeros(nFiles,nStims,nCond,nFilt,nWin); % preallocate
+    for s = 1:nStims % once for each stimulus frequency
+        for cn = 1:nCond % once for each stimulus condition (dev, std, ctrl)
+            for ft = 1:nFilt % once for each filter
+                for w = 1:nWin % once for each time window
+                    rmsV{c}(:,s,cn,ft,w) = rms(dataAv{c}{(s-1)*nCond+cn}(winAll(w,1):winAll(w,2),:,ft)); % calculate rms
+                    rmsGrAv{c}(s,cn,ft,w) = mean(rmsV{c}(:,s,cn,ft,w),1); % grand average of RMS data
+                end
+            end
+        end
+    end
+end
+%% run t-test or ANOVA on each time window
+
+cohensD_V = zeros(nComb,nStims,nFilt,nWin,3); % preallocate
+hV = zeros(nComb,nStims,nFilt,nWin); % preallocate
+pV = zeros(nComb,nStims,nFilt,nWin); % preallocate
+ciV = zeros(nComb,nStims,nFilt,nWin,2); % preallocate
+statsV = cell(nComb,nStims,nFilt,nWin); % preallocate
+switch inclCtrl
+    case {1,2}
+        anovaV = cell(nComb,nStims,nFilt,nWin); % preallocate
+        multcompV = cell(nComb,nStims,nFilt,nWin); % preallocate
+        within = table([1,2,3]','VariableNames',{'stim_prob'}); % define within model
+    case 3
+        cohensD_V = zeros(nComb,nStims,nFilt,nWin,6); % preallocate
+        anovaV = cell(nComb,nStims,nFilt,nWin); % preallocate
+        multcompV = cell(nComb,nStims,nFilt,nWin); % preallocate
+        within = table([1,2,3,4]','VariableNames',{'stim_prob'}); % define within model
+end
+for c = 1:nComb % once for each stimulus combination
+    for s = 1:nStims % once for each stimulus frequency
+        for ft = 1:nFilt % once for each filter
+            for w = 1:nWin % once for each time window
+                cohensD_V(c,s,ft,w,1) = CohensD(rmsV{c}(:,s,1,ft,w),rmsV{c}(:,s,2,ft,w)); % calculate Cohens D; dev vs. std
+                [hV(c,s,ft,w),pV(c,s,ft,w),ciV(c,s,ft,w,:),statsV{c,s,ft,w}] = ttest(rmsV{c}(:,s,1,ft,w),rmsV{c}(:,s,2,ft,w)); % run t-test
+                switch inclCtrl
+                    case {1,2} % with control data
+                        t = table(rmsV{c}(:,s,1,ft,w),rmsV{c}(:,s,2,ft,w),rmsV{c}(:,s,3,ft,w),'VariableNames',{'dev','std','ctrl'}); % required for ANOVA
+                        rm = fitrm(t,'dev-ctrl~1','WithinDesign',within); % required for ANOVA
+                        anovaV{c,s,ft,w} = ranova(rm); % run ANOVA
+                        multcompV{c,s,ft,w} = multcompare(rm,'stim_prob','ComparisonType','bonferroni'); % run multiple comparison to find differences between responses
+                        cohensD_V(c,s,ft,w,2) = CohensD(rmsV{c}(:,s,1,ft,w),rmsV{c}(:,s,3,ft,w)); % calculate Cohens D; dev vs. ctrl
+                        cohensD_V(c,s,ft,w,3) = CohensD(rmsV{c}(:,s,3,ft,w),rmsV{c}(:,s,2,ft,w)); % ctrl vs. std
+                    case 3
+                        t = table(rmsV{c}(:,s,1,ft,w),rmsV{c}(:,s,2,ft,w),rmsV{c}(:,s,3,ft,w),rmsV{c}(:,s,4,ft,w),'VariableNames',{'dev','std','50Per','MR'}); % required for ANOVA
+                        rm = fitrm(t,'dev-MR~1','WithinDesign',within); % required for ANOVA
+                        anovaV{c,s,ft,w} = ranova(rm); % run ANOVA
+                        multcompV{c,s,ft,w} = multcompare(rm,'stim_prob','ComparisonType','bonferroni'); % run multiple comparison to find differences between responses
+                        cohensD_V(c,s,ft,w,2) = CohensD(rmsV{c}(:,s,1,ft,w),rmsV{c}(:,s,3,ft,w)); % calculate Cohens D; dev vs. 50Per
+                        cohensD_V(c,s,ft,w,3) = CohensD(rmsV{c}(:,s,1,ft,w),rmsV{c}(:,s,4,ft,w)); % calculate Cohens D; dev vs. MR
+                        cohensD_V(c,s,ft,w,4) = CohensD(rmsV{c}(:,s,3,ft,w),rmsV{c}(:,s,2,ft,w)); % 50Per vs. std
+                        cohensD_V(c,s,ft,w,5) = CohensD(rmsV{c}(:,s,4,ft,w),rmsV{c}(:,s,2,ft,w)); % MR vs. std
+                        cohensD_V(c,s,ft,w,6) = CohensD(rmsV{c}(:,s,4,ft,w),rmsV{c}(:,s,3,ft,w)); % MR vs. 50Per
+                end
+            end
+        end
+    end
+end
+%% save data
+
+switch inclCtrl
+    case 0
+        save(filename,'winAll','winTxt','rmsV','rmsGrAv','hV','pV','ciV','statsV','cohensD_V')
+    case {1,2,3}
+        save(filename,'winAll','winTxt','rmsV','rmsGrAv','hV','pV','ciV','statsV','anovaV','multcompV','cohensD_V')
+end

DD_UniteDatapoints.m

@@ -0,0 +1,66 @@
+projectdir = 'C:\Users\johan\OneDrive\Uni\PhD\Matlab Rawdata\ABR_Dataset_bats\Deviance detection\MR\Vocs\';
+mFileName = 'DD_Data_MR.mat';
+mVarName = 'Data_MR_Vocs';
+dinfo = dir(fullfile(projectdir));
+dinfo([dinfo.isdir]) = [];     %get rid of all directories including . and ..
+nfiles = length(dinfo);
+Data_MR_Vocs = [];
+for j = 1:nfiles
+    filename = fullfile(projectdir, dinfo(j).name);
+    fid = fopen(filename, 'r');
+    Parameters = fread(fid,60,'float32')';
+    sf = Parameters(4);
+    fdown = Parameters(17);
+    FileIndex = Parameters(2);
+    datale = Parameters(44);  % length of datablock
+    blocklength = Parameters(44);
+    nstims = Parameters(12);
+    xs = (0:datale-1)/sf*fdown;
+    % AB locked runs
+    if FileIndex == 8
+        avdata = fread(fid,8*blocklength,'float32')';
+        seq = fread(fid,nstims,'float32')';
+        rawdata = fread(fid,'float32')';
+        fclose(fid);
+    end
+    % Std/Dev only runs & MR runs
+    if FileIndex==6||FileIndex==7||FileIndex==9||FileIndex==10||FileIndex==11||FileIndex==12
+        avdata = fread(fid,4*blocklength,'float32')';
+        seq = fread(fid,nstims,'float32')';
+        rawdata = fread(fid,'float32')';
+        fclose(fid);
+    end
+    % MR runs
+    if FileIndex==13
+        avdata = fread(fid,4*blocklength,'float32')';
+        seq = fread(fid,nstims,'float32')';
+        rest = fread(fid,'float32')'; % "rest" contains rawdata and MR sequence
+        rawdata = rest(1:end-nstims);
+        seq_MR = rest(end-nstims+1:end);
+        fclose(fid);
+    end
+    % Rep Suppr runs
+    if FileIndex==14
+        dur_unit = Parameters(52);
+        recpts = floor(dur_unit*sf/fdown);
+        av_unit_NF = fread(fid,recpts,'float32')';
+        av_unit = fread(fid,recpts,'float32')';
+        rawdata = fread(fid,'float32')';
+    end
+    % Save files
+    File = convertCharsToStrings(filename);
+    switch FileIndex
+        case {1,2,3,4,5,6,7,8,9,10,11,12}
+            Data_temp = {File,Parameters,avdata,seq,rawdata};
+        case 13
+            Data_temp = {File,Parameters,avdata,seq_MR,rawdata};
+        case 14
+            Data_temp = {File,Parameters,av_unit_NF,av_unit,rawdata};
+    end
+    Data_MR_Vocs = [Data_MR_Vocs;Data_temp];
+end
+if isfile(mFileName)==0
+    save(mFileName,mVarName,'-v7.3')
+else
+    save(mFileName,mVarName,'-append')
+end

DD_avData_50Per.mat

@@ -0,0 +1 @@
+/annex/objects/MD5-s113753369--c2a9e1002bf0208c1a9e460d9466c4d2

DD_avData_MR.mat

@@ -0,0 +1 @@
+/annex/objects/MD5-s113063327--87aee43e801f13b8fd9e38f87433b5f1

DD_avData_Oddball.mat

@@ -0,0 +1 @@
+/annex/objects/MD5-s262722686--d992c479ef29a3d4dfd3a088778f32ce

Filter_Butter.m

@@ -0,0 +1,11 @@
+%% Butterworth_Fco
+
+function [outtrace] = Filter_Butter(signal,order,lf,hf,type,fs)
+if nargin~= 6
+    error('You should have six input parameters. Check this');
+end
+Wn = [lf,hf]./(fs/2); % non-dimensional scale frequency
+[z,p,k] = butter(order,Wn,type);
+[sos,g] = zp2sos(z,p,k);
+outtrace = filtfilt(sos,g,signal);
+end

Filter_IFFT.m

@@ -0,0 +1,28 @@
+function [outtrace] = Filter_IFFT(intrace, frq1, frq2, sfrq, no)
+% Calc FFT-like square------
+fmax=sfrq/2;
+ll = length(intrace);
+if mod(length(intrace),2) > 0
+    intrace = intrace(1:ll-1);
+end
+Fpts = length(intrace)/2;
+lpfrqs=round(Fpts*frq1/fmax);
+if lpfrqs < 1
+    lpfrqs = 1;
+end
+hpfrqs=round(Fpts*frq2/fmax);
+Fwindow = ones(1,Fpts)*0.0000001;
+Fwindow(lpfrqs:hpfrqs)=1;
+Fwindow = [Fwindow fliplr(Fwindow)]';
+if no == 1
+    Fwindow = Fwindow +1;
+    Fwindow = mod(Fwindow,2);
+end
+%----------------------------
+FFT = fft(intrace);
+mag = abs(FFT);
+mag = (mag.*Fwindow);  % .* elementwise multiplication and exch row/columns
+phase = angle(FFT);    % exch row/columns
+% % ifft---------------------
+filteredSpec =mag.*exp(1i.*phase); % Lu .* array multiplication
+outtrace=real(ifft(filteredSpec))';  % % exch row/columns    Lu

blCorr.m

@@ -0,0 +1,14 @@
+function [outData] = blCorr(inData,timeWinStart,timeWinEnd)
+%DEMEAN performs baseline correction using the given time window
+%   input stream is divided into trials and average voltage within the set
+%   time window is calculated for each trial. Afterwards this average 
+%   activity is subtracted from the whole trial to normalise activity in
+%   set time window to 0.
+
+% calculate average voltage within time window
+avV = mean(inData(timeWinStart:timeWinEnd,:));
+
+% subtract average voltage from trials
+outData = inData-avV;
+end
+

cutblockTrial.m

@@ -0,0 +1,11 @@
+%% Function: cut data into blocks, based on trial-triggers
+
+function [blocks] = cutblockTrial(data,trigs_all,blocklen,nBlck)
+nTrials = size(trigs_all,3);
+pSmpl = blocklen/nTrials;
+blocks = zeros(blocklen,nBlck);
+for b = 1:nBlck
+    for s = 0:nTrials-1
+        blocks(s*pSmpl+1:s*pSmpl+pSmpl,b) = data(trigs_all(1,b,s+1):trigs_all(1,b,s+1)+pSmpl-1); % cut block at each trigger point, creating "s" sections of pSmpl per block
+    end
+end

demean.m

@@ -0,0 +1,14 @@
+function [outData] = demean(inData)
+%DEMEAN performs baseline correction using the given time window
+%   input stream is divided into trials and average voltage within the set
+%   time window is calculated for each trial. Afterwards this average 
+%   activity is subtracted from the whole trial to normalise activity in
+%   set time window to 0.
+
+% calculate average voltage within time window
+avV = mean(inData);
+
+% subtract average voltage from trials
+outData = inData-avV;
+end
+

sigstar.m

@@ -0,0 +1,312 @@
+function varargout=sigstar(groups,stats,esize,nosort)
+    % sigstar - Add significance stars to bar charts, boxplots, line charts, etc,
+    %
+    % H = sigstar(groups,stats,nsort)
+    %
+    % Purpose
+    % Add stars and lines highlighting significant differences between pairs of groups. 
+    % The user specifies the groups and associated p-values. The function handles much of 
+    % the placement and drawing of the highlighting. Stars are drawn according to:
+    %   * represents p<=0.05
+    %  ** represents p<=1E-2
+    % *** represents p<=1E-3
+    %
+    %
+    % Inputs
+    % groups - a cell array defining the pairs of groups to compare. Groups defined 
+    %          either as pairs of scalars indicating locations along the X axis or as 
+    %          strings corresponding to X-tick labels. Groups can be a mixture of both 
+    %          definition types.
+    % stats -  a vector of p-values the same length as groups. If empty or missing it's 
+    %          assumed to be a vector of 0.05s the same length as groups. Nans are treated
+    %          as indicating non-significance.
+    % nsort -  optional, 0 by default. If 1, then significance markers are plotted in 
+    %          the order found in groups. If 0, then they're sorted by the length of the 
+    %          bar.
+    %
+    % Outputs
+    % H - optionally return handles for significance highlights. Each row is a different
+    %     highlight bar. The first column is the line. The second column is the text (stars).
+    %     
+    %
+    % Examples
+    % 1. 
+    % bar([5,2,1.5])
+    % sigstar({[1,2], [1,3]})
+    %
+    % 2. 
+    % bar([5,2,1.5])
+    % sigstar({[2,3],[1,2], [1,3]},[nan,0.05,0.05])
+    %
+    % 3.  **DOESN'T WORK IN 2014b**
+    % R=randn(30,2);
+    % R(:,1)=R(:,1)+3;
+    % boxplot(R)
+    % set(gca,'XTick',1:2,'XTickLabel',{'A','B'})
+    % H=sigstar({{'A','B'}},0.01);
+    % ylim([-3,6.5])
+    % set(H,'color','r')
+    %
+    % 4. Note the difference in the order with which we define the groups in the 
+    %    following two cases. 
+    % x=[1,2,3,2,1];
+    % subplot(1,2,1)
+    % bar(x)
+    % sigstar({[1,2], [2,3], [4,5]})
+    % subplot(1,2,2)
+    % bar(x)
+    % sigstar({[2,3],[1,2], [4,5]})
+    %
+    % ALSO SEE: demo_sigstar
+    %
+    % KNOWN ISSUES:
+    % 1. Algorithm for identifying whether significance bar will overlap with 
+    %    existing plot elements may not work in some cases (see line 277)
+    % 2. Bars may not look good on exported graphics with small page sizes.
+    %    Simply increasing the width and height of the graph with the 
+    %    PaperPosition property of the current figure should fix things.
+    %
+    % Rob Campbell - CSHL 2013
+
+
+
+    %Input argument error checking
+
+    %If the user entered just one group pair and forgot to wrap it in a cell array 
+    %then we'll go easy on them and wrap it here rather then generate an error
+    if ~iscell(groups) & length(groups)==2
+        groups={groups};
+    end
+
+    if nargin<3
+        stats=repmat(0.05,1,length(groups));
+    end
+    if isempty(stats)
+        stats=repmat(0.05,1,length(groups));
+    end
+    if nargin<4
+        nosort=0;
+    end
+
+
+
+
+    %Check the inputs are of the right sort
+    if ~iscell(groups)
+        error('groups must be a cell array')
+    end
+
+    if ~isvector(stats)
+        error('stats must be a vector')
+    end
+
+    if length(stats)~=length(groups)
+        error('groups and stats must be the same length')
+    end
+
+
+
+
+
+
+    %Each member of the cell array groups may be one of three things:
+    %1. A pair of indices.
+    %2. A pair of strings (in cell array) referring to X-Tick labels
+    %3. A cell array containing one index and one string
+    %
+    % For our function to run, we will need to convert all of these into pairs of
+    % indices. Here we loop through groups and do this. 
+
+    xlocs=nan(length(groups),2); %matrix that will store the indices 
+    xtl=get(gca,'XTickLabel');  
+
+    for ii=1:length(groups)
+        grp=groups{ii};
+
+        if isnumeric(grp)
+            xlocs(ii,:)=grp; %Just store the indices if they're the right format already
+
+        elseif iscell(grp) %Handle string pairs or string/index pairs
+
+            if isstr(grp{1})
+                a=strmatch(grp{1},xtl);
+            elseif isnumeric(grp{1})
+                a=grp{1};
+            end
+            if isstr(grp{2})
+                b=strmatch(grp{2},xtl);
+            elseif isnumeric(grp{2})
+                b=grp{2};
+            end
+
+            xlocs(ii,:)=[a,b];
+        end
+
+        %Ensure that the first column is always smaller number than the second
+        xlocs(ii,:)=sort(xlocs(ii,:));
+
+    end
+
+    %If there are any NaNs we have messed up. 
+    if any(isnan(xlocs(:)))
+        error('Some groups were not found')
+    end
+
+
+
+
+
+
+    %Optionally sort sig bars from shortest to longest so we plot the shorter ones first
+    %in the loop below. Usually this will result in the neatest plot. If we waned to 
+    %optimise the order the sig bars are plotted to produce the neatest plot, then this 
+    %is where we'd do it. Not really worth the effort, though, as few plots are complicated
+    %enough to need this and the user can define the order very easily at the command line. 
+    if ~nosort
+        [~,ind]=sort(xlocs(:,2)-xlocs(:,1),'ascend');
+        xlocs=xlocs(ind,:);groups=groups(ind);
+        stats=stats(ind);
+        esize=esize(ind);
+    end
+
+
+
+    %-----------------------------------------------------
+    %Add the sig bar lines and asterisks 
+    holdstate=ishold;
+    hold on
+
+    H=ones(length(groups),2); %The handles will be stored here
+
+    y=ylim;
+    yd=myRange(y)*0.1; %separate sig bars vertically by 5% 
+
+    for ii=1:length(groups)
+        thisY=findMinY(xlocs(ii,:))+yd;
+        H(ii,:)=makeSignificanceBar(xlocs(ii,:),thisY,stats(ii),round(esize(ii),2));
+    end
+    %-----------------------------------------------------
+
+
+
+
+    %Now we can add the little downward ticks on the ends of each line. We are
+    %being extra cautious and leaving this it to the end just in case the y limits
+    %of the graph have changed as we add the highlights. The ticks are set as a
+    %proportion of the y axis range and we want them all to be the same the same
+    %for all bars.
+    yd=myRange(ylim)*0.01; %Ticks are 1% of the y axis range
+    for ii=1:length(groups)
+        y=get(H(ii,1),'YData');
+        y(1)=y(1)-yd;
+        y(4)=y(4)-yd;   
+        set(H(ii,1),'YData',y)
+    end
+
+
+
+
+    %Be neat and return hold state to whatever it was before we started
+    if ~holdstate
+        hold off
+    elseif holdstate
+        hold on
+    end
+
+
+    %Optionally return the handles to the plotted significance bars (first column of H)
+    %and asterisks (second column of H).
+    if nargout>0
+        varargout{1}=H;
+    end
+
+
+end %close sigstar
+
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%Internal functions
+
+function H=makeSignificanceBar(x,y,p,es)
+    %makeSignificanceBar produces the bar and defines how many asterisks we get for a 
+    %given p-value
+
+
+    if p<=1E-3
+        stars = {num2str(es),'***'};
+        offset = 0.03;
+    elseif p<=1E-2
+        stars = {num2str(es),'**'};
+        offset = 0.03;
+    elseif p<=0.05
+        stars = {num2str(es),'*'};
+        offset = 0.03;
+    else
+        stars='n.s.';
+        offset = 0.03;
+    end
+            
+    x=repmat(x,2,1);
+    y=repmat(y,4,1);
+
+    H(1)=plot(x(:),y,'-k','Linewidth',2,'Tag','sigstar_bar');
+
+    %Increase offset between line and text if we will print "n.s."
+    %instead of a star. 
+%     offset = 0.025;
+%     switch stars
+%         case {'*','**','***'}
+%             offset = 0.015;
+%         case {'n.s.'}
+%             offset = 0.055;
+%     end
+
+    starY=mean(y)+myRange(ylim)*offset;
+    H(2)=text(mean(x(:)),starY,stars,...
+        'HorizontalAlignment','center',...
+        'BackGroundColor','none',...
+        'Tag','sigstar_stars');
+
+    Y=ylim;
+    if Y(2)<starY*1.1
+        ylim([Y(1),starY+myRange(Y)*0.1])
+    end
+
+
+end %close makeSignificanceBar
+
+
+
+function Y=findMinY(x)
+    % The significance bar needs to be plotted a reasonable distance above all the data points
+    % found over a particular range of X values. So we need to find these data and calculat the 
+    % the minimum y value needed to clear all the plotted data present over this given range of 
+    % x values. 
+    %
+    % This version of the function is a fix from Evan Remington
+    oldXLim = get(gca,'XLim');
+    oldYLim = get(gca,'YLim');
+
+    axis(gca,'tight')
+    
+    %increase range of x values by 0.1 to ensure correct y max is used
+    x(1)=x(1)-0.1;
+    x(2)=x(2)+0.1;
+    
+    set(gca,'xlim',x) %Matlab automatically re-tightens y-axis
+
+    yLim = get(gca,'YLim'); %Now have max y value of all elements within range.
+    Y = max(yLim);
+
+    axis(gca,'normal')
+    set(gca,'XLim',oldXLim,'YLim',oldYLim)
+
+end %close findMinY
+
+
+function rng=myRange(x)
+    %replacement for stats toolbox range function
+    rng = max(x) - min(x);
+end %close myRange
+

sortBlck.m

@@ -0,0 +1,35 @@
+function [outDev,outStd,nDevAcc,nStdAcc] = sortBlck(inData,seqAcc,nTrialsAcc,AfrqI,BfrqI,runType)
+%SORTBLCK Divides blocks into deviant and standard responses
+%   Detailed explanation goes here
+
+% define some variables
+pSmpl = size(inData,1);
+nTrials = size(inData,2);
+nBlcks = size(inData,3);
+nFiles = size(inData,4);
+nFilt = size(inData,5);
+
+% apply sequence on each block
+nDevAcc = zeros(nBlcks,nFiles); % preallocate
+nStdAcc = zeros(nBlcks,nFiles); % preallocate
+outDev = zeros(pSmpl,nTrials,nBlcks,nFiles,nFilt); % preallocate
+outStd = zeros(pSmpl,nTrials,nBlcks,nFiles,nFilt); % preallocate
+for f = 1:nFiles
+    for b = 1:nBlcks
+        % create sequences
+        switch runType
+            case {1,2}
+                seqDev = logical(seqAcc(1:nTrialsAcc(b,f),b,f));
+                seqStd = ~seqAcc(1:nTrialsAcc(b,f),b,f);
+            case 3
+                seqDev = seqAcc(1:nTrialsAcc(b,f),b,f)==AfrqI;
+                seqStd = seqAcc(1:nTrialsAcc(b,f),b,f)==BfrqI;
+        end
+        nDevAcc(b,f) = sum(seqDev);
+        nStdAcc(b,f) = sum(seqStd);
+        % apply sequences on data
+        outDev(:,1:nDevAcc(b,f),b,f,:) = inData(:,seqDev,b,f,:);
+        outStd(:,1:nStdAcc(b,f),b,f,:) = inData(:,seqStd,b,f,:);
+    end
+end
+end

subsample.m

@@ -0,0 +1,34 @@
+function [outDev,outStd,nTrialsSub] = subsample(inDev,inStd,seqAcc,nTrialsAcc)
+%SUBSAMPLE Extracts deviants after standards and standards before deviants
+%   Detailed explanation goes here
+
+% define some variables
+pSmpl = size(inDev,1);
+nTrials = size(inDev,2)*0.2;
+nBlcks = size(inDev,3);
+nFiles = size(inDev,4);
+nFilt = size(inDev,5);
+
+% perform the subsampling
+nTrialsSub = zeros(nBlcks,nFiles); % preallocate
+outDev = zeros(pSmpl,nTrials,nBlcks,nFiles,nFilt); % preallocate
+outStd = zeros(pSmpl,nTrials,nBlcks,nFiles,nFilt); % preallocate
+for f = 1:nFiles
+    for b = 1:nBlcks
+        % create sequences
+        seqDev = logical(seqAcc(1:nTrialsAcc(b,f),b,f));
+        seqStd = ~seqAcc(1:nTrialsAcc(b,f),b,f);
+        
+        % shift dev and std sequence by one point and index into it with original std and dev sequence to obtain stds before devs and devs after stds
+        seqDev_as = logical([zeros(1,1);seqStd(1:nTrialsAcc(b,f)-1)]); % shift original std sequences by 1 point (backwards)
+        seqDev_as = seqDev_as(seqDev); % index into shifted standard sequences with original dev sequences to obtain dev after stds
+        seqStd_bd = logical([seqDev(2:nTrialsAcc(b,f));zeros(1,1)]); % shift original dev sequences by 1 point (forwards)
+        seqStd_bd = seqStd_bd(seqStd); % index into shifted deviant sequences with original std sequences to obtain stds before devs
+        nTrialsSub(b,f) = sum(seqStd_bd);
+        % apply newly generated sequences on dev and std datasets
+        outDev(:,1:nTrialsSub(b,f),b,f,:) = inDev(:,seqDev_as,b,f,:);
+        outStd(:,1:nTrialsSub(b,f),b,f,:) = inStd(:,seqStd_bd,b,f,:);
+    end
+end
+end
+

zscoreJo.m

@@ -0,0 +1,15 @@
+function [outData] = zscoreJo(inData,pSmpl)
+%ZSCOREJO zscore function adapted to ABR data-structure
+%   basically the default zscore function of Matlab but adapted to the
+%   structure of ABR data
+
+% divide input data into trials
+inData = reshape(inData,[],1);
+
+% perform z score normalisation
+outData = zscore(inData);
+
+% restore original shape of input data
+outData = reshape(outData,pSmpl,[]);
+end
+