DevianceDetectionNaturalVocalisationsBats/DD_Processing.m

%% 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