% This script uses EEG data in which all trials (segments) of all
% participants is kept. Baseline correction, artifact rejection on
% eye-movements (EOG) and on individual channel mode have been performed.
% The data was also down-sampeled to 500Hz. All channels are also included.
% The script requires a data set in which trial markers are unique. Eye-
% movement contaminated trials are excluded. The rest of the artefacts are
% marked per channel. The script first checks what segments are removed by
% looking for the marker 'Bad Interval', gives all segments a name
% corresponding to the item, and finally replaces the values in bad
% segments with NA values. Segments are then averaged (101:600, 0:100 is
% baseline) to do further trial-based analysis

% Matlab v.R2019a
%% Processing Object data

% Specify folder for trial data
workdirparent = ('L:\projects\jvbproject\Nicola Project CBN\CBNtrialdata_all_pp');
cd('L:\projects\jvbproject\Nicola Project CBN\CBNtrialdata_all_pp');

% Collect files and lists
files = dir('*trialdata.mat');
numberofdatasets = size(files,1);

% Create the brain regions
% Anterior (P2): Fp1, AF3, AF7, F1, F3, F5, F7, FC1, FC3, FC5, FT7, Fp2,
% AF4, AF8, F2, F4, F6, F8, FC2,  FC4, FC6, FT8, FPz, AFz, Fz, FCz
% Central Posterior (N400 en LPP): C1, C3, C5, T7, CP1, CP3, CP5, TP7, P1,
% P3, P5, P7, PO3, PO7, O1, C2, C4, C6, T8, CP2, CP4, CP6, TP8, P2, P4, P6,
% P8, PO4, PO8, O2, Cz, CPz, Pz, POz, Oz

Anterior = ["Fp1","AF3","AF7","F1","F3","F5","F7","FC1","FC3","FC5",...
    "FT7","Fp2","AF4","AF8","F2","F4","F6","F8","FC2","F4","FC6","FT8",...
    "Fpz","AFz","Fz","FCz"];
CPosterior = ["C1","C3","C5","T7","CP1","CP3","CP5","TP7","P1","P3",...
    "P5","P7","PO3","PO7","O1","C2","C4","C6","T8","CP2","CP4","CP6",...
    "TP8","P2","P4","P6","P8","PO4","PO8","O2","Cz","CPz","Pz","POz","Oz"];
% size of regions
regionAnt = size(Anterior,2);
regionCPos = size(CPosterior,2);
anteriority = ["Anterior","Central Posterior"]';

% Predefine variabeles
%newstimfield = "stimuluslist";
endbaseline = 100;
numseg = zeros(1,1);
channels = zeros(64,1);
brainregions = regionAnt+regionCPos;
pattern = '(CBN(\d+)([A-Z]+)(\d+))'; % necessary to isolate pp name

% file for output
sjsFileTrials = ['CBN_ERP_trial_data_', date, '.txt']; %text file for analysis
fidTrials = fopen(sjsFileTrials, 'wt+');
all_erptrialdata_labels = ["Participant","P2_200-250","N400_350-450","LPP_500-650","block","item"];
all_sjs_erptrialdata = [];


% Loop through files
for f = 45 %1:length(files) %files' %
    file = files(f);
    % Clear variables
    clearvars -except files numberofdatasets Anterior CPosterior regionAnt ...
        regionCPos anteriority endbaseline numseg channels brainregions ...
        pattern all_erptrialdata_labels all_sjs_erptrialdata ...
        sjsFileTrials fidTrials f file
    
    % Load in one file. This is from 1 participant, from 1 word, in 1
    % condition.
    data = load(file.name, '-mat');
    
    % Count number of segments. This differs from 540 for some participants
    % since eye movement trials are excluded
    numseg = data.SegmentCount;
    ERP = zeros(numseg*numberofdatasets*brainregions,1);
    
    % Make datavariables
    segdata_Ant = zeros(numseg,1);
    segdata_CPos = zeros(numseg,1);
    segchannel = {data.Markers.ChannelNumber};
    
    % ARTIFACT REMOVAL OF SEGMENTS
    % The individual channel mode artifact rejection needs to be
    % done. In the IC AR, segments that were removed are labeled
    % by Markers in the field Description. These can vary, but the most
    % common one is 'Bad Amplitude'. EOG removed segments do not seem to
    % have these descriptions, so the script will gather the removed
    % channels in a different variable.
    
    % First, the script gathers the marker types and the marker
    % descriptions.
    segtype = {data.Markers.Type};
    segdes = {data.Markers.Description};
    
    % Now we prepare variabeles in which we will later write down
    % codes for what segments and channels will be removed.
    segmenttype = zeros(size(segtype,2),1);
    segmentdescription = zeros(size(segdes,2),1);
    
    % For IC AR, the script goes over the segment types and whenever there
    % is a 'Bad Interval', it will note down 99 in the segmenttype variabele.
    % If there is no 'Bad Interval' it notes down a 0.
    for marker = 1:size(segtype,2)
        if strcmp(segtype(marker),'Bad Interval')
            segmenttype(marker) = 99;
        end
    end
    
    % Reshape variabele segmenttype to go back to the structures of
    % segments in the original data structure.
    badchannels = reshape(segmenttype,numseg,size(data.Markers,2));
    
    % The data from the structure is isolated in a separate variabele
    % because I could not find a way to integrate it in the structure.
    dataperchannel = struct2cell(data);
    % Only the first 64 rows are kept because these are the 64 scalp electrodes.
    dataperchannel = dataperchannel(1:64);
    
    % The script goes over the data per channel and replaces all values with
    % NaN if there was a 99 noted down in the variabele bad intervals.
    for elec = 1:64
        for seg = 1:numseg
            % IC artifact rejection
            for samplepoint=1:size(badchannels,2)
                if (badchannels(seg,samplepoint) == 99)
                    dataperchannel{data.Markers(seg,samplepoint).ChannelNumber}(seg,:)=NaN;
                end
            end
        end
    end
    
    % AVERAGE WITHIN SEGMENTS OVER TIME-WINDOWS
    % average results per time window per channel. Pay attention t=0 is the
    % start of the baseline, Time 0 of the ERP is at t=101, so 200ms afer
    % onset is at t=101+200ms*0.5Hz=201
    for x = 1:64
        dataperchannel_200_250{x} = nanmean(dataperchannel{x}(:,201:225),2);
        dataperchannel_350_450{x} = nanmean(dataperchannel{x}(:,276:326),2);
        dataperchannel_500_650{x} = nanmean(dataperchannel{x}(:,351:426),2);
    end
    
    % AVERAGE ACROSS THE CHANNELS FOR THE 2 BRAIN REGIONS
    % Isolates the 64 scalp channels
    channels = {data.Channels(1:64).Name}';
    
    % Isolates the data of channels and puts it in their respective
    % variables that are related to the brain regions.
    for channel_index = 1:64
        for group_id = 1:regionAnt
            if Anterior(group_id) == channels(channel_index)
                dataperchan_Ant(:,group_id) = dataperchannel_200_250{channel_index};
            end
        end
        for group_id = 1:regionCPos
            if CPosterior(group_id) == channels(channel_index)
                dataperchan_CPos_350_450(:,group_id) = dataperchannel_350_450{channel_index};
                dataperchan_CPos_500_650(:,group_id) = dataperchannel_500_650{channel_index};
            end
        end
    end
    
    % Save segmented data in the 2 brain regions, but not NaN channel data.
    for seg = 1:numseg
        segdata_Ant(seg,:)=nanmean(dataperchan_Ant(seg,:));
        segdata_CPos_350_450(seg,:)=nanmean(dataperchan_CPos_350_450(seg,:));
        segdata_CPos_550_650(seg,:)=nanmean(dataperchan_CPos_500_650(seg,:));
    end
    
    % MAKE R-FRIENDLY FORMAT
    % Save data in R-friendly format
    Pp(:,1)=repmat(regexp(file.name, pattern, 'tokens'),numseg,1);
    ERP_pp_200_250 = num2cell(segdata_Ant);
    ERP_pp_350_450 = num2cell(segdata_CPos_350_450);
    ERP_pp_500_650 = num2cell(segdata_CPos_550_650);
    Pp_ERP = horzcat(Pp,ERP_pp_200_250,ERP_pp_350_450, ERP_pp_500_650);
    
    % Add item number and block to R-friendly format
    % Isolate the relevant markers to check what the item nr and block was.
    for t = 1:numseg
        blockmarkers{t} = NaN;
        trialmarkers{t} = NaN;
        for m = 1:size(data.Markers,2)
            data.Markers(t,m).Description
            if ~isempty(data.Markers(t,m).Description)
                switch data.Markers(t,m).Description(1)
                    case 'S'
                    case '1'
                        if size(data.Markers(t,m).Description,2)==3
                            blockmarkers{t} = str2num(data.Markers(t,m).Description);
                        else
                            trialmarkers{t} = str2num(data.Markers(t,m).Description);
                        end
                    case 'B'
                    otherwise
                        trialmarkers{t} = str2num(data.Markers(t,m).Description);
                end
            end
        end
    end

    Pp_ERP_item_block = horzcat(Pp_ERP,blockmarkers',trialmarkers');
    all_sjs_erptrialdata = vertcat(all_sjs_erptrialdata,Pp_ERP_item_block);
    
    
    %f = f + 1;
end
all_erptrialdata = vertcat(all_erptrialdata_labels,all_sjs_erptrialdata);
writematrix(all_erptrialdata, sjsFileTrials);
fclose(fidTrials);
