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