%% define some variables

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 (50 ms)

% 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!)
recDurMult = 1; % multiplier of recdur (1 = original rec dur)

% subsample blocks
subSConS = 0; % decide whether to subsample number of blocks in order to have the same number of trials as in oddball condition

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

% filter settings
type = 1; % 1: Butterworth, 2: IFFT
order = 4;
lowc = [0.1,10,40,45,50]; % 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 artefact 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
% TEST
blcWin = 0.005; % 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)
%% data processing


% name combinations and load data
for repDaSe = 1:2 % run once for each condition

    % load DD data (for number of conS trials)
    switch repDaSe
        case 1
            load("DD_avData_Oddball_25ms_SbD.mat")
        case 2
            load("DD_avData_Oddball_50ms_SbD.mat")
    end
    % calculate mean number of trials
    conST = round(mean(mean(conST{1},1),2));
    
    switch repDaSe
        case 1 % 25 ms
            repRateName = '_25ms';
        case 2 % 50 ms
            repRateName = '_50ms';
    end

    load(['RS_Data',repRateName,'.mat'])
    fileName = append("RS_avData",repRateName);
    SOAName = repRateName(2:end);

    combName = ["Low";"High";"Low_o1st";"High_o1st"];
    nComb = size(combName,1); % number of different stimulus combination to be processed

    dataAv_cell = cell(1,nComb); % preallocate
    dataGrAv_cell = cell(1,nComb); % preallocate
    dataSe_cell = cell(1,nComb); % preallocate
    filenames_stim_cell = cell(1,nComb); % preallocate
    recID_cell = cell(1,nComb); % preallocate
    timetrBlck_cell = cell(1,nComb); % preallocate
    timetrUnit_cell  = cell(1,nComb); % preallocate
    pts2begin_cell = cell(1,nComb); % preallocate
    nBlcks_cell = cell(1,nComb); % preallocate
    nCh_s_cell = cell(1,nComb); % preallocate
    nFiles_cell = cell(1,nComb); % preallocate
    nFilt_cell = cell(1,nComb); % preallocate
    nResp_cell = cell(1,nComb); % preallocate
    SOA_cell = cell(1,nComb); % preallocate
    stimDur_cell = cell(1,nComb); % preallocate
    stimDelay_cell = cell(1,nComb); % preallocate
    nBlcksAcc_cell = cell(1,nComb); % preallocate
    fs_stim_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
    pSiResp_cell = cell(1,nComb); % preallocate
    stimCat_cell = cell(1,nComb); % preallocate
    
    for c = 1:nComb % run analysis once for each stimulation-combination
        switch c % switch dataset depending on loop iteration
            case 1
                data = Data_Low;
            case 2
                data = Data_High;
            case 3
                data = Data_Low_o1st;
            case 4
                data = Data_High_o1st;
        end
        
        % Organise dataset in multiple variables (in separate script)
        RS_OrgData              
        
        % extract data of channel to be processed NOT FINISHED YET
        ch_idx = zeros(nCh_s(1),nFiles,'single'); % preallocate
        rec_raw_ch = cell(nFiles,1); % preallocate
        recRawCut = zeros(pUnit,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) = cutblockOnset(rec_raw_ch{f},trigs_all(f,:),pUnit,nBlcks); % cut out buffers and overlengths from raw data and split data into step-blocks; ATTENTION: different routine to DD processing (cutblockOnset vs. cutblockTrial) 
%             if repDaSe==2
%                 figure; plot(recRawCut(:,:,f)) % plot raw trace to identify artefacts
%             end
        end
        
        % artefact rejection
        nBlcksAcc = zeros(1,nFiles); % preallocate
        switch artRejI
            case 0 % if artefact rejection is turned off
                nBlcksAcc(:) = nBlcks;
            case 1 % run if artefact rejection is turned on
                % preallocate
                maxV = zeros(1,nBlcks,nFiles);
                rejIdx = zeros(1,nBlcks,nFiles);
                for f = 1:nFiles
                    for t = 1:nBlcks
                        maxV(1,t,f) = max(abs(recRawCut(:,t,f))); % find (absolute) max value in each trial
                        rejIdx(1,t,f) = maxV(:,t,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
                    nBlcksAcc(f) = sum(rejIdx(:,:,f)); % count number of accepted trials for current file and step
                    recRawCut(:,1:nBlcksAcc(f),f) = recRawCut(:,logical(rejIdx(:,:,f)),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
                end
        end

        % remove trigger artifact
        if remTrArtI==1
            for f = 1:nFiles
                switch artDateI(f)
                    case 0 % if no Stimtrak was used
                        nArt = 1; % fix only first artifact
                    case 1 % if Stimtrak was used
                        nArt = 2; % fix first and second artifact
                end
                for art = 1:nArt % once for each artifact per stimulus
                    for r = 1:nResp % once for each stimulus in the block
                        switch art
                            case 1
                                artStart = 0+SOA*(r-1); % start of trigger artifact in s
                            case 2
                                artStart = 0.0107+SOA*(r-1); % start of trigger artifact in s
                        end
                        if art==1&&r==1 % for the very first artifact per block, add 1 point to not start at 0
                            artStartPts = round(artStart*fs)+1; % convert start of trigger artifact  to points
                        else
                            artStartPts = round(artStart*fs); % convert start of trigger artifact  to points
                        end
                        artDur = 0.0015; % duration of trigger artifact in s
                        artDurPts = round(artDur*fs); % convert duration of trigger artifact to points
                        artEndPts = artStartPts+artDurPts; % calculate end of trigger artifact in points
                        for t = 1:nBlcks
                            recRawCut(artStartPts:artEndPts,t,f) = linspace(recRawCut(artStartPts,t,f),recRawCut(artEndPts,t,f),artDurPts+1); % replace artifact by linear vector from start to end
                        end
                    end
                end
            end
        end
                
        % process data
        recPro = recRawCut; % rename
        recFilt = zeros(pUnit*nBlcks,nFiles,nFilt-1); % preallocate
        for f = 1:nFiles
            % detrending
            if detrendI==1
                recPro(:,1:nBlcksAcc(f),f) = detrend(recPro(:,1:nBlcksAcc(f),f));
            end
            % z score normalisation
            if zsNormI==1
                recPro(:,1:nBlcksAcc(f),f) = zscoreJo(recPro(:,1:nBlcksAcc(f),f),pUnit);
            end
            % demeaning
            if demeanI==1
                recPro(:,1:nBlcksAcc(f),f) = demean(recPro(:,1:nBlcksAcc(f),f));
            end
            % filtering
            recProRS = reshape(recPro,nBlcks*pUnit,nFiles);
            % recProRS(1:nBlcksAcc(f)*pUnit,f) = Filter_Butter(recProRS(1:nBlcksAcc(f)*pUnit,f),20,45,55,"stop",fs); % apply 50 Hz notch filter to the data to get rid of electrical noise
            for ft = 1:nFilt-1 % "-1" because "nFilt" includes raw-array which must be excluded here
                recFilt(1:nBlcksAcc(f)*pUnit,f,ft) = Filter_Butter(recProRS(1:nBlcksAcc(f)*pUnit,f),order,lowc(ft),hic,"bandpass",fs);
            end
        end
        recProFilt = cat(3,recProRS,recFilt); % concatenate unfiltered and filtered data in 4th dimension
        recProFilt = reshape(recProFilt,pUnit,nBlcks,nFiles,nFilt); % reshape to split raw-data into separate trials (2nd dimension)
        

        % perform baseline correction: depends on processed stimulus (different to bat analysis)
        blcStartPtsAll = round(pts2begin(1)+(blcEnd-blcWin)*fs); % transform time window for baseline correction into points
        blcEndPtsAll = round(pts2begin(1)+blcEnd*fs);
        switch c
            case {1,3}
                blcStartPts = blcStartPtsAll(1);
                blcEndPts = blcEndPtsAll(1);
            case {2,4}
                blcStartPts = blcStartPtsAll(2);
                blcEndPts = blcEndPtsAll(2);
        end
        for f = 1:nFiles
            for ft = 1:nFilt
                recProFilt(:,1:nBlcksAcc(f),f,ft) = blCorr(recProFilt(:,1:nBlcksAcc(f),f,ft),blcStartPts,blcEndPts);
            end
        end
        
        
        % calculate average & grand-average of sorted responses
        dataAv = zeros(pUnit,nFiles,nFilt);
        for f = 1:nFiles
            dataAv(:,f,:) = mean(recProFilt(:,1:nBlcksAcc(f),f,:),2); % averaging is performed only on those vectors that contain an actual response and not only zeros (2nd dim)
        end
        % calculate grand average
        dataGrAv = reshape(mean(dataAv,2),pUnit,nFilt);
        % calculate standard deviation
        dataStdD = reshape(std(dataAv,0,2),pUnit,nFilt);
        % calculate standard error
        dataSe = dataStdD/sqrt(nFiles);
        
        
        % concatenate important data in another cell array with one cell
        % for each stimulation-condition
        dataAv_cell{c} = dataAv;
        dataGrAv_cell{c} = dataGrAv;
        dataSe_cell{c} = dataSe;
        filenames_stim_cell{c} = filenames_stim;
        recID_cell{c} = recID;
        timetrBlck_cell{c} = timetrBlck;
        timetrUnit_cell{c} = timetrUnit;
        pts2begin_cell{c} = pts2begin;
        nCh_s_cell{c} = nCh_s;
        nFiles_cell{c} = nFiles;
        nFilt_cell{c} = nFilt;
        nResp_cell{c} = nResp;
        stimDur_cell{c} = stimDur;
        stimDelay_cell{c} = stimDelay;
        fs_stim_cell{c} = fs_stim;
        fs_cell{c} = fs;
        ch_idx_cell{c} = ch_idx;
        nameCh_cell{c} = nameCh;
        chLbl_cell{c} = chLbl;
        nBlcksAcc_cell{c} = nBlcksAcc;
        pSiResp_cell{c} = pSiResp;

        %% calculate stimulus template to be plotted next to response

        % load stim data
        switch c
            case {1,3} % Low
                [stim] = audioread('Broadband Low.wav');
            case {2,4} % High
                [stim] = audioread('Broadband High.wav');
        end
        % add ISI after the stimulus
        stimLen = length(stim); % calculate length of stim
        stimLong = zeros(round(SOA*fs_stim),1); % create array with the length of stim + ISI (= SOA)
        stimLong(1:stimLen) = stim; % replace the first points by the actual stimulus
        stimLongLen = length(stimLong); % calculate length of stim + ISI
        stimCat_cell{c} = zeros(round((stimLongLen*nResp)+(IBI*fs_stim)),1); % preallocate
        % create concatenated stimulus of block-stimulation
        for r = 1:nResp
            stimCat_cell{c}(1+(r-1)*stimLongLen:r*stimLongLen) = stimLong;
            stimLong = -stimLong; % inverse stimulus for the next concatenation
        end
    end

    save(fileName,'dataAv_cell','dataGrAv_cell','dataSe_cell','filenames_stim_cell','recID_cell',...
        'timetrBlck_cell','timetrUnit_cell','pts2begin_cell','combName',...
        'nComb','nCh_s_cell','nFiles_cell','nFilt_cell','nResp_cell','nBlcksAcc_cell',...
        'stimDur_cell','stimDelay_cell','fs_cell','ch_idx_cell','nameCh_cell',...
        'chLbl_cell','artRejI','zsNormI','detrendI','demeanI','blcEnd','SOAName','stimCat_cell',...
        'lowc','fs_stim_cell','pSiResp_cell')
end
beep