%% Subscript: Organises data from datasets in multiple variables

% subdivide dataset into its components
if remI==1
    data(remRec{repDaSe},:) = []; % remove certain recordings from analysis
end
filenames_stim = cat(1,data{:,1});
bsPos = cell2mat(strfind(filenames_stim,'\')); % find positions of back slashes in filename
% bsPos = strfind(filenames_stim,'\'); % find positions of back slashes in filename
startPos = bsPos(:,end); % use last back slash position as starting point
ptPos = cell2mat(strfind(filenames_stim,'.')); % find positions of points in filename (here only one)
% ptPos = strfind(filenames_stim,'.'); % find positions of points in filename (here only one)
endPos = ptPos; % use point position as ending point
recID = extractBetween(filenames_stim,startPos,endPos,'Boundaries','Exclusive'); % extract between Positions
switch runType
    case {1,5,6}
        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,4}
        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_stim = filenames_stim(fileI_Ctrl); % use only those control filenames that have a matching oddball response
end
param_stim = cat(1,data{:,2});
date = cat(1,data{:,8});
artDateI = date<datetime(trArtDate,'Format','yyyyMMdd'); % logical array that identifies recordings with artefact
seq = cat(1,data{:,3});
filenames_rec = cat(1,data{:,4});
param_rec = struct2cell(cat(1,(data{:,5}))); % concatenate struct arrays, then convert to cell arrays
rec_raw = cat(1,data(:,7)); % concatenate recorded raw data
rec_raw = cellfun(@(x) x/100,rec_raw,'un',0); % divide data by 100 because data from EP-Preamp are amplified by factor 100
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_stim = param_stim(1,4);
fdown_stim = param_stim(1,17);
fIndx = param_stim(1,2);
Alvl = param_stim(1,9);
Afrq = param_stim(1,10);
Blvl = param_stim(1,26);
Bfrq = param_stim(1,25);
AfrqI = 999; % frequency index in MR control --> redundant but required in other runTypes
BfrqI = 999;
nDev = param_stim(1,23);
devProb = param_stim(1,24);
nTrials = param_stim(1,12);
nStd = nTrials-nDev;
nBlcks = 2; % number of blocks (AB and BA)
pts2begin_stim = param_stim(1,41);
stimDur = round(param_stim(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_stim(1,13)+trDel,4); % stimulus delay in s (corrected for sound travelling time)
repper = round(param_stim(1,8),4);
recdur = repper;
recdur = recdur*recDurMult; % use double recdur

% recording-parameters
fs = cat(2,param_rec{1,1})';
chLbl = string(cat(2,param_rec{3,:}))';

% down sample data if activated
switch dwnSmpleI
    case 0
        dwnSmplF = 1; % down sampling factor is 1 if no down sampling is performed
    case 1
        dwnSmplF = dwnSmplRate/fs(1); % down sampling factor
        for f = 1:nFiles
            rec_raw{f} = resample(rec_raw{f},dwnSmplRate/1000,fs(1)/1000); % down sample data
        end
        fs = dwnSmplRate; % change sampling rate to down sampling rate
end

% additional parameters
nCh = zeros(1,nFiles); % preallocate
rawLen = zeros(1,nFiles); % preallocate
for f = 1:nFiles
    nCh(f) = size(rec_raw{f},1); % number of channes in raw file
    rawLen(f) = size(rec_raw{f},2); % number of samples in each channel of raw file
end
blockLen = round(recdur.*nTrials.*fs); % length of one block (AB or BA) in samples
pSmpl = blockLen(1)/nTrials(1); % number of samples of each single response
if channels=="all"
    nameCh = chLbl;
    nCh_s = nCh;
else
    nameCh = channels;
    nCh_s = size(channels,2);
end
pts2begin = round(stimDelay.*fs); % samples before stimulus onset
timetr = (0:pSmpl-1)/fs(1); % time trace of one trial

% triggers
nTrigsE = 20000; % number of expected triggers
trigs_all = zeros(nFiles,nBlcks,nTrials); % preallocate
for f = 1:nFiles % extract trigger events for each file individually
    eventsTemp = struct2cell(cat(1,data{f,6})); % extract all events from struct array "data" and convert it to cell array
    eventsTemp = eventsTemp(:,:,eventsTemp(1,1,:)=="Stimulus"); % find all triggers "Stimulus" within events array
    trigsIdx = eventsTemp(2,1,:)=="S  1"; % find all triggers "S  1" within events array
    trigsTemp = eventsTemp(3,1,trigsIdx(1,1,:)); % apply logical array on events array to extract triggers
    trigsTemp = cell2mat(trigsTemp); % convert trigger array to matrix
    trigsDiff = diff(trigsTemp,1,3); % calculate differences between adjacent triggers to find too small gaps (additional triggers that must be removed)
    corTrIdx = trigsDiff>(pSmpl/recDurMult*0.8); % find indices of correct triggers; those with a distance of at least 80 % of the expected distance
    corTrIdx = logical(cat(3,1,corTrIdx)); % add "1" at the beginning of the indexing array (since diff has one value less than original array) and make it logical
    trigsTemp = trigsTemp(corTrIdx); % select only those triggers that have the correct distance to the adjacent ones
    nTrigs = size(trigsTemp,3); % determine number of recorded triggers
    trigsTempRs = reshape(trigsTemp,1,nTrigs); % reshape trigger array to 2 dimensional array

    if nTrigs>nTrigsE % if more triggers were recorded than expected, stop and print error
        error("At least one file contains too many trigger events. Please check and remove those files.")
    end
    % if less triggers were recorded than expected, go on with the
    % following routine to add missing triggers
    while nTrigs<nTrigsE % repeat loop as long as number of recorded triggers is lower than number of expected triggers (only if 2 or more consecutive triggers are missing)
        trigsDiff = diff(trigsTempRs); % calculate difference-array of trigger array
        mTrigIdx = find(trigsDiff>1000&trigsDiff<100000)+1; % find indices of missing trigger events ("<100000" because gap between sequences must be excluded)
        nMTrigs = numel(mTrigIdx);
        if isempty(mTrigIdx) % if no indices were found, either first or last trigger is missing. In this case it is not possible to reconstruct trigger array. File must be excluded
            error("Either the first or last trigger is missing in at least one file. Please check and remove those files.")
        end
        valPreTrig = trigsTempRs(mTrigIdx-1); % check if this works for more than one missing trigger; determine trigger value preceding the missing trigger
        valTrig = valPreTrig+pSmpl; % calculate value for missing trigger
        for t = 1:nMTrigs % insert missing triggers in trigger array
            trigsTempRs = [trigsTempRs(1:mTrigIdx(t)+(t-2)),valTrig(t),trigsTempRs(mTrigIdx(t)+(t-1):end)]; % "(t-2)" and "(t-1)" because every added trigger shifts next trigger index one the the right
        end
        nTrigs = size(trigsTempRs,2); % count number of triggers again
    end % repeat loop if number of triggers is still less than number of expected triggers
    for b = 0:nBlcks-1
        trigs_all(f,b+1,:) = trigsTempRs(b*nTrials+1:b*nTrials+nTrials); % after all missing triggers are added, concatenate triggers of all files
    end
end
trigs_all = round(trigs_all*dwnSmplF); % multiply triggers with down sampling factor (which is 1 if no down sampling is performed)