DevianceDetectionHumanABR/Deviance Detection_Humans/DD_Processing.m

%% define some variables
tic
clear
% remove recordings
remI = 0; % remove certain recordings from analysis? (0: no, 1: yes)
% define recordings to be removed from analysis (if remI==1)
remRec{1} = []; % in dataset 1 (25 ms)
remRec{2} = []; % in dataset 2 (30 ms)
remRec{3} = []; % in dataset 3 (50 ms)

% select processing options
dwnSmpleI = 0; % down sample data? (1: yes, 0: no)
dwnSmplRate = 5000; % which should be the new sampling rate if signal is down sampled?
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!)
shCtrlI = 0; % shuffle control trials before averaging?
subsampleI = 1; % calculate averages from all responses (= 0) or only from dev after std and std before dev (= 1)?
% subStd = 2; % decide which std responses are used when subsampling; 1: Std after Dev, 2: Std before Dev (does not affect omission responses)
recDurMult = 1; % multiplier of recdur (1 = original rec dur)

% define channels to be processed; type "all" to process all channels
channels = ["Plus"];

% filter settings
type = 1; % 1: Butterworth, 2: IFFT
order = 4;
lowc = [10,45,50,55]; % low cut frequency in Hz
hic = 1500; % high cut frequency in Hz
nFilt = size(lowc,2)+1; % number of filters used (+1 because of unfiltered data)

% rejection criterion
artRejI = 1; % turn artifact rejection on (1) or off (0)
rejThr = 200; % uV threshold (every trial with min/max value exceeding this will be rejected before averaging)
remTrArtI = 1; % should trigger artifact be removed from the data?
trArtDate = '20231220'; % day before which trigger artifact existed

% baseline correction settings
blcWin = 0.001; % length of time window for baseline correction in s
chirpPeak = [0.00648,0.00842]; % amplitude peaks of chirps (first: low, second: high)
blcEnd = chirpPeak; % end of baseline correction window in s (first: low freq chirp, second: high freq chirp)


% speaker distance
spDis = 0.26; % distance from ear to speaker in m (to correct for sound-travel delay)

% extract consecutive standards
nConS = 8; % number of consecutive standards to be extracted from oddball sequences
exactNConSI = 1; % should only the exact number (and not less) of consecutive standards be used? (= 1)

% permutation
nPerm = 1; % number of permutations (dev and std)
%% data processing
for subStd = 2:2 % once for Std before Dev, once for Std after Dev
    for repDaSe = 1:3 % once for each repetition rate
        switch repDaSe
            case 1 % 25 ms
                repRateName = '_25ms';
                nRunType = [1,3,4,5,6]; % Oddball, Low only, High only, Low OM, High OM (50Per missing)
            case 2 % 30 ms
                repRateName = '_30ms';
                nRunType = [1,3,4,5,6]; % Oddball, Low only, High only, Low OM, High OM (50Per missing)
            case 3 % 50 ms
                repRateName = '_50ms';
                nRunType = (1:6); % Oddball, 50Per, Low only, High only, Low OM, High OM
        end
        % name combinations and load data
        nRunType = 1; % ONLY FOR TESTING REASONS - REMOVE!!!
        for runType = nRunType % run once for each condition
            combName = ["Low,High"];
            switch runType
                case 1
                    load(['DD_Data_Oddball',repRateName,'.mat'])
                    fileName = append("DD_avData_Oddball",repRateName);
                case 2
                    load(['DD_Data_50Per',repRateName,'.mat'])
                    fileName = append("DD_avData_50Per",repRateName);
                case {3,4,5,6}
                    load(['DD_Data_OR',repRateName,'.mat'])
                    switch runType
                        case 3
                            fileName = append("DD_avData_LowOnly",repRateName);
                        case 4
                            fileName = append("DD_avData_HighOnly",repRateName);
                        case 5
                            fileName = append("DD_avData_HighOm",repRateName);
                        case 6
                            fileName = append("DD_avData_LowOm",repRateName);
                    end
            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
            dataAvP_cell = cell(1,nComb); % preallocate
            dataGrAvP_cell = cell(1,nComb); % preallocate
            dataDifAvP_cell = cell(1,nComb); % preallocate
            dataDifGrAvP_cell = cell(1,nComb); % preallocate
            dataSeP_cell = cell(1,nComb); % preallocate
            filenames_stim_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
            nCh_s_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
            fs_cell = cell(1,nComb); % preallocate
            ch_idx_cell = cell(1,nComb); % preallocate
            nameCh_cell = cell(1,nComb); % preallocate
            chLbl_cell = cell(1,nComb); % preallocate
            devR_cell = cell(1,nComb); % preallocate
            stdR_cell = cell(1,nComb); % preallocate
            conSAv = cell(1,nComb); % preallocate
            conSGrAv = cell(1,nComb); % preallocate
            conSSe = cell(1,nComb); % preallocate
            conST = cell(1,nComb); % preallocate
    
            for c = 1:nComb % run analysis once for each stimulation-combination
                switch runType
                    case 1 % Oddball
                        data = Data_Oddball_Low_High;
                    case 2 % 50Per
                        data = Data_50Per_Low_High;
                    case 3 % Low only
                        data = Data_OR_Low_HighOm;
                    case 4 % High only
                        data = Data_OR_LowOm_High;
                    case 5 % High OM (Low present)
                        data = Data_OR_Low_HighOm;
                    case 6 % Low OM (High present)
                        data = Data_OR_LowOm_High;
                end
    
                % Organise dataset in multiple variables (in separate script)
                DD_OrgData
    
                % extract data of channel to be processed
                ch_idx = zeros(nCh_s(1),nFiles,'single'); % preallocate
                rec_raw_ch = cell(nFiles,1); % preallocate
                recRawCut = zeros(blockLen,nBlcks,nFiles); % preallocate
                for f = 1:nFiles
                    ch_idx(:,f) = find(ismember(chLbl(f,:),nameCh)); % calculate indices of channels to be processed
                    rec_raw_ch{f} = rec_raw{f}(ch_idx(:,f),:); % pick data from channels to be processed
                    recRawCut(:,:,f) = cutblockTrial(rec_raw_ch{f},trigs_all(f,:,:),blockLen,nBlcks); % cut out buffers and overlengths from raw data and split data into step-blocks
                    %             figure; plot(recRawCut(:,:,f)) % plot raw trace to identify artifacts
                end
    
                % artifact 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 artifact 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 artifact 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
                                    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
                
                % remove trigger artifact
                if remTrArtI==1
                    for f = 1:nFiles
                        switch artDateI(f)
                            case 0 % if no Stimtrak artifact exists
                                nArt = 1; % remove only artifact at the beginning
                            case 1 % if Stimtrak artifact exists
                                nArt = 2; % remove both artifacts, at the beginning and 11 ms later
                        end
                        for art = 1:nArt % once for each artifact per stimulus
                            switch art
                                case 1
                                    artDur = 0.0015; % duration of trigger artifact in ms
                                    artStartPts = 1; % for first artifact, startimg point is always 1
                                case 2
                                    artStart = 0.0107; % start of trigger artifact in ms
                                    artDur = 0.0015; % duration of trigger artifact in ms
                                    artStartPts = round(artStart*fs); % convert start of trigger artifact  to points
                            end
                            artDurPts = round(artDur*fs); % convert duration of trigger artifact to points
                            artEndPts = artStartPts+artDurPts; % calculate end of trigger artifact in points
                            for b = 1:nBlcks
                                for t = 1:nTrials
                                    recRawCut(artStartPts:artEndPts,t,b,f) = linspace(recRawCut(artStartPts,t,b,f),recRawCut(artEndPts,t,b,f),artDurPts+1); % replace artifact by linear vector from start to end
                                end
                            end
                        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",fs);
                        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)
    
                % split blocks into deviant and standard responses
                [devR,stdR,nDevAcc,nStdAcc] = sortBlck(recProFilt,seqAcc,nTrialsAcc);
                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!

                % extract sequences of consecutive standards
                avNStd = countConS(recProFilt,seqAcc,nTrialsAcc); % extract clusters of consecutive standards
                [conS,conSI] = findConS(recProFilt,seqAcc,nTrialsAcc,nConS,exactNConSI); % extract clusters of consecutive standards
                conST{c} = zeros(nBlcks,nFiles);
                for f = 1:nFiles
                    for b = 1:nBlcks
                        conST{c}(b,f) = size(conSI{b,f},2);
                    end
                end
    
                % perform baseline correction: after Dev/Std split because it
                % depends on processed stimulus (different to bat analysis)
                blcStartPts = round(pts2begin+(blcEnd-blcWin)*fs); % transform time window for baseline correction into points
                blcEndPts = round(pts2begin+blcEnd*fs);
                conSRs = cell(nBlcks,nFiles); % preallocate
                for f = 1:nFiles
                    for b = 1:nBlcks
                        switch runType % blc changes for omission responses
                            case {1,2,3,4}
                                switch b % change time window for baseline correction depending on processed block (ADev/BStd & BDev/AStd) and stim combination (c)
                                    case 1
                                        blcStartPtsDev = blcStartPts(1);
                                        blcStartPtsStd = blcStartPts(2);
                                        blcEndPtsDev = blcEndPts(1);
                                        blcEndPtsStd = blcEndPts(2);
                                    case 2
                                        blcStartPtsDev = blcStartPts(2);
                                        blcStartPtsStd = blcStartPts(1);
                                        blcEndPtsDev = blcEndPts(2);
                                        blcEndPtsStd = blcEndPts(1);
                                end
                            case 5 % High OM: Although high was omitted, we care about a potential "Low" omission response, triggered by the presented low freq chirps. Hence bcl only for low-timing
                                blcStartPtsDev = blcStartPts(1);
                                blcStartPtsStd = blcStartPts(1);
                                blcEndPtsDev = blcEndPts(1);
                                blcEndPtsStd = blcEndPts(1);
                            case 6 % Low OM: Although low was omitted, we care about a potential "High" omission response, triggered by the presented high freq chirps. Hence bcl only for high-timing
                                blcStartPtsDev = blcStartPts(2);
                                blcStartPtsStd = blcStartPts(2);
                                blcEndPtsDev = blcEndPts(2);
                                blcEndPtsStd = blcEndPts(2);
                        end
                        for ft = 1:nFilt
                            devR(:,1:nDevAcc(b,f),b,f,ft) = blCorr(devR(:,1:nDevAcc(b,f),b,f,ft),blcStartPtsDev,blcEndPtsDev);
                            stdR(:,1:nStdAcc(b,f),b,f,ft) = blCorr(stdR(:,1:nStdAcc(b,f),b,f,ft),blcStartPtsStd,blcEndPtsStd);
                            % TEST normalise each cluster instead of each
                            % response in cluster
                            conSRs{b,f} = zeros(pSmpl*nConS,size(conS{b,f},3),nFilt); % preallocate
                            for cs = 1:nConS
                                conSRs{b,f}(1+pSmpl*(cs-1):pSmpl*cs,:,ft) = conS{b,f}(:,cs,:,ft); % concatenate consecutive responses of one cluster (required for blc)
                            end
                            conSRs{b,f}(:,:,ft) = blCorr(reshape(conSRs{b,f}(:,:,ft),pSmpl*nConS,[]),blcStartPtsStd,blcEndPtsStd); % perform blc
                            for cs = 1:nConS
                                conS{b,f}(:,cs,:,ft) = conSRs{b,f}(1+pSmpl*(cs-1):pSmpl*cs,:,ft); % split responses again
                            end
                            % ALTERNATIVE - normalise each response in
                            % cluster
                            % for cs = 1:nConS
                            %     conS{b,f}(:,cs,:,ft) = blCorr(reshape(conS{b,f}(:,cs,:,ft),pSmpl,[]),blcStartPtsStd,blcEndPtsStd);
                            % end
                        end
                    end
                end
    
                % 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 subsampleI
                    case 1
                        switch runType
                            case {1,2,3,4} % Dev after Std, Std before Dev
                                switch subStd
                                    case 1 % std after dev
                                        [devR,stdR,nTrialsDev,nTrialsStd] = subsampleOm(devR,stdR,seqAcc,nTrialsAcc); % subsample oddball responses --> number of trials will be much less afterwards (= nTrialsSub)
                                        nDevUsed = nTrialsDev; % define number of subsampled deviants as number of used deviants
                                        nStdUsed = nTrialsStd; % define number of subsampled standards as number of used standards
                                    case 2 % std before dev
                                        [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
                                end
                                switch runType
                                    case 1
                                        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 {5,6} % Dev after Std, Std after Dev
                                [devR,stdR,nTrialsDev,nTrialsStd] = subsampleOm(devR,stdR,seqAcc,nTrialsAcc); % subsample oddball responses --> number of trials will be much less afterwards (= nTrialsSub)
                                nDevUsed = nTrialsDev; % define number of subsampled deviants as number of used deviants
                                nStdUsed = nTrialsStd; % define number of subsampled standards as number of used standards
                        end
                end
    

                % if activated, run random perutation to shuffle 50Per and MR responses --> important for 50Per
                % control: to make sure not only the earliest responses are used
                % for averaging
                if shCtrlI==1
                    switch runType
                        case 2
                            for f = 1:nFiles
                                for b = 1:nBlcks
                                    rDev = randperm(nDevUsed(b,f)); % create random vector containing all used dev responses (they cover the whole sequence)
                                    devR(:,1:nDevUsed(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(nStdUsed(b,f)); % repeat the same for std responses
                                    stdR(:,1:nStdUsed(b,f),b,f,:) = stdR(:,rStd,b,f,:);
                                end
                            end
                    end
                end
    
                % calculate average & grand-average of sorted responses
                switch runType
                    case 2 % 50 Per
                        nDevUsed = nDevUsed_oddb{c}; % change nDevUsed to trial number from oddball recording to end up with the same number of trials for averaging
                        nStdUsed = nDevUsed_oddb{c}; % do the same for nStdUsed
                end
                devAv = zeros(pSmpl,nBlcks,nFiles,nFilt);
                stdAv = zeros(pSmpl,nBlcks,nFiles,nFilt);
                conSAv{c} = zeros(pSmpl,nBlcks,nFiles,nConS,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)
                        stdAv(:,b,f,:) = mean(stdR(:,1:nStdUsed(b,f),b,f,:),2);
                        for cs = 1:nConS
                            conSAv{c}(:,b,f,cs,:) = mean(conS{b,f}(:,cs,conSI{b,f}(cs,:),:),3);
                        end
                    end
                end
                
                % perform permutations
                devAvP = zeros(pSmpl,nBlcks,nFiles,nFilt,nPerm);
                stdAvP = zeros(pSmpl,nBlcks,nFiles,nFilt,nPerm);
                for f = 1:nFiles
                    for b = 1:nBlcks
                        devRU = devR(:,1:nDevUsed(b,f),b,f,:); % exctract meaningful data
                        stdRU = stdR(:,1:nStdUsed(b,f),b,f,:);
                        trialsRU = cat(2,devRU,stdRU); % concatenate dev and std data
                        for p = 1:nPerm
                            trialsRP = trialsRU(:,(randperm(size(trialsRU,2))),:,:,:); % shuffle dev and std trials
                            devRP = trialsRP(:,(1:nDevUsed(b,f)),:,:,:); % use first half of shuffled data as dev permutation
                            stdRP = trialsRP(:,(nDevUsed(b,f)+1:end),:,:,:); % use second half of shuffled data as std permutation
                            devAvP(:,b,f,:,p) = mean(devRP,2); % calculate average
                            stdAvP(:,b,f,:,p) = mean(stdRP,2);
                        end
                    end
                end
                
                % calculate grand average
                devGrAv = reshape(mean(devAv,3),pSmpl,nBlcks,nFilt);
                stdGrAv = reshape(mean(stdAv,3),pSmpl,nBlcks,nFilt);
                conSGrAv{c} = reshape(mean(conSAv{c},3),pSmpl,nBlcks,nConS,nFilt);
                devGrAvP = reshape(mean(devAvP,3),pSmpl,nBlcks,nFilt,nPerm);
                stdGrAvP = reshape(mean(stdAvP,3),pSmpl,nBlcks,nFilt,nPerm);
                % calculate standard deviation
                devStdD = reshape(std(devAv,0,3),pSmpl,nBlcks,nFilt);
                stdStdD = reshape(std(stdAv,0,3),pSmpl,nBlcks,nFilt);
                conSStdD = reshape(std(conSAv{c},0,3),pSmpl,nBlcks,nConS,nFilt);
                devStdDP = reshape(std(devAvP,0,3),pSmpl,nBlcks,nFilt,nPerm);
                stdStdDP = reshape(std(stdAvP,0,3),pSmpl,nBlcks,nFilt,nPerm);
                % calculate standard error
                devSe = devStdD/sqrt(nFiles);
                stdSe = stdStdD/sqrt(nFiles);
                conSSe{c} = conSStdD/sqrt(nFiles);
                devSeP = devStdDP/sqrt(nFiles);
                stdSeP = stdStdDP/sqrt(nFiles);
    
                % split dataset into several sub-variables for simplisity
                switch runType
                    case {1,3,4,5,6}
                        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 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
                        % do the same for permutation data
                        ADevP{c} = devAvP(:,1,:,:,:); % A is deviant in block 1 (AB)
                        ADevP{c} = reshape(ADevP{c},pSmpl,nFiles,nFilt,nPerm); % reshape av data into 4D array
                        AStdP{c} = stdAvP(:,2,:,:,:); % A is standard in block 2 (BA)
                        AStdP{c} = reshape(AStdP{c},pSmpl,nFiles,nFilt,nPerm); % reshape av data into 4D array
                        BDevP{c} = devAvP(:,2,:,:,:); % B is deviant in block 2 (BA)
                        BDevP{c} = reshape(BDevP{c},pSmpl,nFiles,nFilt,nPerm); % reshape av data into 4D array
                        BStdP{c} = stdAvP(:,1,:,:,:); % B is standard in block 1 (AB)
                        BStdP{c} = reshape(BStdP{c},pSmpl,nFiles,nFilt,nPerm); % reshape av data into 4D array
                        ADevGrAvP{c} = devGrAvP(:,1,:,:);
                        ADevGrAvP{c} = reshape(ADevGrAvP{c},pSmpl,nFilt,nPerm); % reshape grav into 3D array
                        AStdGrAvP{c} = stdGrAvP(:,2,:,:);
                        AStdGrAvP{c} = reshape(AStdGrAvP{c},pSmpl,nFilt,nPerm); % reshape grav data into 3D array
                        BDevGrAvP{c} = devGrAvP(:,2,:,:);
                        BDevGrAvP{c} = reshape(BDevGrAvP{c},pSmpl,nFilt,nPerm); % reshape grav data into 3D array
                        BStdGrAvP{c} = stdGrAvP(:,1,:,:);
                        BStdGrAvP{c} = reshape(BStdGrAvP{c},pSmpl,nFilt,nPerm); % reshape grav data into 3D array
                        ADevSeP{c} = devSeP(:,1,:,:);
                        ADevSeP{c} = reshape(ADevSeP{c},pSmpl,nFilt,nPerm); % reshape standard error 3D array
                        AStdSeP{c} = stdSeP(:,2,:,:);
                        AStdSeP{c} = reshape(AStdSeP{c},pSmpl,nFilt,nPerm); % reshape standard error into 3D array
                        BDevSeP{c} = devSeP(:,2,:,:);
                        BDevSeP{c} = reshape(BDevSeP{c},pSmpl,nFilt,nPerm); % reshape standard error into 3D array
                        BStdSeP{c} = stdSeP(:,1,:,:);
                        BStdSeP{c} = reshape(BStdSeP{c},pSmpl,nFilt,nPerm); % reshape standard error into 3D array
                    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
                        % do the same for permutation data
                        ACtrlP{c} = devAvP(:,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
                        ACtrlP{c} = reshape(ACtrlP{c},pSmpl,nFiles,nFilt,nPerm); % reshape av data into 4D array
                        BCtrlP{c} = devAvP(:,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
                        BCtrlP{c} = reshape(BCtrlP{c},pSmpl,nFiles,nFilt,nPerm); % reshape av data into 4D array
                        ACtrlGrAvP{c} = devGrAvP(:,1,:,:);
                        ACtrlGrAvP{c} = reshape(ACtrlGrAvP{c},pSmpl,nFilt,nPerm); % reshape grav into 3D array
                        BCtrlGrAvP{c} = devGrAvP(:,2,:,:);
                        BCtrlGrAvP{c} = reshape(BCtrlGrAvP{c},pSmpl,nFilt,nPerm); % reshape grav data into 3D array
                        ACtrlSeP{c} = devSeP(:,1,:,:);
                        ACtrlSeP{c} = reshape(ACtrlSeP{c},pSmpl,nFilt,nPerm); % reshape standard error into 3D array
                        BCtrlSeP{c} = devSeP(:,2,:,:);
                        BCtrlSeP{c} = reshape(BCtrlSeP{c},pSmpl,nFilt,nPerm); % reshape standard error into 3D array
                end

                % organise processed dev and std responses in cell array
                switch runType
                    case {1,3,4,5,6}
                        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}};
                        % permutation data
                        dataAvP= {ADevP{c},AStdP{c},BDevP{c},BStdP{c}};
                        dataGrAvP = {ADevGrAvP{c},AStdGrAvP{c},BDevGrAvP{c},BStdGrAvP{c}};
                        dataDifAvP = [];
                        dataDifGrAvP = [];
                        dataSeP = {ADevSeP{c},AStdSeP{c},BDevSeP{c},BStdSeP{c}};
                    case 2
                        dataAv = {ACtrl{c},BCtrl{c}};
                        dataGrAv = {ACtrlGrAv{c},BCtrlGrAv{c}};
                        dataDifAv = [];
                        dataDifGrAv = []; % empty bc unused (use line below if needed)
                        % dataDifAv = {ADifDev{c},ADifStd{c},BDifDev{c},BDifStd{c}};
                        % dataDifGrAv = {ADifDevGrAv{c},ADifStdGrAv{c},BDifDevGrAv{c},BDifStdGrAv{c}};
                        dataSe = {ACtrlSe{c},BCtrlSe{c}};
                        % permutation data
                        dataAvP = {ACtrlP{c},BCtrlP{c}};
                        dataGrAvP = {ACtrlGrAvP{c},BCtrlGrAvP{c}};
                        dataDifAvP = [];
                        dataDifGrAvP = []; % empty bc unused (use line below if needed)
                        dataSeP = {ACtrlSeP{c},BCtrlSeP{c}};
                end
                % for consecutive std clusters: swap blocks to have A in
                % block 1 and B in block 2 - to be consistent with the
                % rest
                conST{c}([1,2],:) = conST{c}([2,1],:); % swap blocks
                conSAv{c}(:,[1,2],:,:,:) = conSAv{c}(:,[2,1],:,:,:); % swap blocks
                conSGrAv{c}(:,[1,2],:,:) = conSGrAv{c}(:,[2,1],:,:); % swap blocks
                conSSe{c}(:,[1,2],:,:) = conSSe{c}(:,[2,1],:,:); % swap blocks

                % 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;
                dataAvP_cell{c} = dataAvP;
                dataGrAvP_cell{c} = dataGrAvP;
                dataDifAvP_cell{c} = dataDifAvP;
                dataDifGrAvP_cell{c} = dataDifGrAvP;
                dataSeP_cell{c} = dataSeP;
                filenames_stim_cell{c} = filenames_stim;
                recID_cell{c} = recID;
                fileI_Oddb_cell{c} = fileI_Oddb;
                fileI_Ctrl_cell{c} = fileI_Ctrl;
                timetr_cell{c} = timetr;
                pts2begin_cell{c} = pts2begin;
                nCh_s_cell{c} = nCh_s;
                nBlcks_cell{c} = nBlcks;
                nFiles_cell{c} = nFiles;
                nFilt_cell{c} = nFilt;
                stimDur_cell{c} = stimDur;
                stimDelay_cell{c} = stimDelay;
                fs_cell{c} = fs;
                ch_idx_cell{c} = ch_idx;
                nameCh_cell{c} = nameCh;
                chLbl_cell{c} = chLbl;
                nDevUsed_cell{c} = nDevUsed;
                nStdUsed_cell{c} = nStdUsed;
                devR_cell{c} = devR;
                stdR_cell{c} = stdR;
            end
            switch runType
                case {1,2,3,4}
                    switch subStd
                        case 1 % std after dev
                            fileName = append(fileName,"_SaD");
                        case 2
                            fileName = append(fileName,"_SbD");
                    end
            end
            save(fileName,'dataAv_cell','dataGrAv_cell','dataDifAv_cell','dataDifGrAv_cell','dataSe_cell',...
                'filenames_stim_cell','recID_cell',...
                'fileI_Oddb_cell','fileI_Ctrl_cell','timetr_cell','pts2begin_cell','combName',...
                'nComb','nCh_s_cell','nBlcks_cell','nFiles_cell','nFilt_cell','nDevUsed_cell','nStdUsed_cell',...
                'stimDur_cell','stimDelay_cell','fs_cell','ch_idx_cell','nameCh_cell',...
                'chLbl_cell','subsampleI','artRejI','zsNormI','detrendI','demeanI','blcEnd','nPerm','pSmpl',...
                'conSAv','conSGrAv','conSSe','nConS','conST')
            save(append(fileName,'_perm'),'dataAvP_cell','dataGrAvP_cell',...
                'dataDifAvP_cell','dataDifGrAvP_cell','dataSeP_cell','nPerm') % save trial data for permutation in separate file
        end
    end
end
% copy relevant files to Rep Suppression folder
% curDir = pwd;
% savDir = 'C:\Users\johan\OneDrive\Uni\PhD\MATLAB\ABR2022_Joh\Analysis\Repetition Suppression_Humans';
% copyfile([curDir,'\DD_avData_Oddball_25ms_SbD.mat'],savdir)
% copyfile([curDir,'\DD_avData_Oddball_50ms_SbD.mat'],savdir)
toc
beep