eegmp_preproc/code/a6_automatic_artifact_correction/a6_automatic_artifact_correction.m

function a6_automatic_artifact_correction(id, rootpath)

if ismac % run if function is not pre-compiled
    id = '21'; % test for example subject
    currentFile = mfilename('fullpath');
    [pathstr,~,~] = fileparts(currentFile);
    cd(fullfile(pathstr,'..', '..'))
    rootpath = pwd;
end

% inputs
pn.eeg_BIDS = fullfile(rootpath, 'data', 'inputs', 'rawdata', 'eeg_BIDS');
pn.channel_locations = fullfile(rootpath, 'data', 'inputs', 'rawdata', 'channel_locations');
pn.events = fullfile(rootpath, 'data', 'inputs', 'rawdata', 'events');
pn.tools = fullfile(rootpath, 'tools');
% outputs
pn.eeg_ft = fullfile(rootpath, 'data', 'outputs', 'eeg');
pn.history = fullfile(rootpath, 'data', 'outputs', 'history');

if ismac % run if function is not pre-compiled
     addpath(fullfile(pn.tools, 'fieldtrip')); ft_defaults;
     addpath(fullfile(pn.tools, 'helpers'));
end

% N = 33
IDs = tdfread(fullfile(pn.eeg_BIDS, 'participants.tsv'));
IDs = cellstr(IDs.participant_id);

id = str2num(id);

display(['processing ID ' num2str(IDs{id})]);

%% load data

load(fullfile(pn.eeg_ft, [IDs{id}, '_task-xxxx_eeg_seg.mat']),'data_eeg', 'events');

% remove M2 (collinear with M1 due to avg. mastoid REF: mean(M1,M2)=0)
cfg_preproc               = [];
cfg_preproc.channel       = {'all', '-M2'};
data_eeg = ft_preprocessing(cfg_preproc, data_eeg);

data_ica = load(fullfile(pn.eeg_ft, [IDs{id}, '_task-xxxx_eeg_ica.mat']),'data');

data = data_eeg; clear data_eeg;
data.elec = data_ica.data.elec;
data.chanlocs = data_ica.data.chanlocs; clear data_ica;

%% ------------------ ARTIFACT DETECTION - PREPARATION ----------------- %%

% load config

load(fullfile(pn.history, [IDs{id}, '_task-xxxx_config.mat']),'config');

config.elec = data.elec;

%% ICA (from weights)

% ica config
cfg.method           = 'runica';
cfg.channel          = {'all'}; % additional channels should be excluded already...
cfg.trials           = 'all';
cfg.numcomponent     = 'all';
cfg.demean           = 'yes';
cfg.runica.extended  = 1;

% ICA solution
cfg.unmixing     = config.ica1.unmixing;
cfg.topolabel    = config.ica1.topolabel;

% components
comp = ft_componentanalysis(cfg,data);

% clear cfg
clear cfg data

%% remove components

% get IC labels
iclabels = config.ica1.iclabels.manual;

% cfg for rejecting components (reject: blinks, eye movements,
% ecg, ref); JQK: added artifacts and emg
cfg.component = sortrows([iclabels.bli(:); iclabels.mov(:); iclabels.hrt(:); iclabels.ref(:); iclabels.art(:); iclabels.emg(:)])';
cfg.demean    = 'yes';

% reject components
data = ft_rejectcomponent(cfg,comp);

% clear cfg
clear cfg comp

%% remove eye & reference channels

cfg.channel     = {'all','-IO1','-IO2','-M1','-VEOG','-HEOG'};
cfg.demean      = 'yes';

% remove channels
tmpdat = ft_preprocessing(cfg,data);

% clear cfg & data variable
clear cfg data

%% ------------------------- ARTIFACT DETECTION ------------------------ %%

% open log file
fid = fopen(fullfile(pn.history, [IDs{id}, '_task-xxxx_ArtCorr.txt']),'a');

fprintf(fid,['*********************   BLOCK   *********************\n']);
fprintf(fid,['function: a6_automatic_artifact_correction.m \n\n']);                
fprintf(fid,['eeg file = ' config.data_file '\n\n']);

n_trl = length(tmpdat.trial);

%%  get artifact contaminated channels by kurtosis, low & high frequency artifacts

cfg.criterion = 3;
cfg.recursive = 'no';

[index0 parm0 zval0] = cm_MWB_channel_x_epoch_artifacts_20170922(cfg,tmpdat);

% write log
tmp_log = '';
for j = 1:length(index0.c)
    tmp_log = [tmp_log num2str(index0.c(j)) ' '];
end; clear j
tmp_log = [tmp_log(1:end-1) '\n'];
fprintf(fid,'(1) automatic bad channel detection:\n');
fprintf(fid,['MWB:          channel(s) ' tmp_log]);

% clear cfg
clear cfg tmp_log

%%  get artifact contaminated channels by FASTER

cfg.criterion = 3;
cfg.recursive = 'no';

[index1 parm1 zval1] = THG_FASTER_1_channel_artifacts_20140302(cfg,tmpdat);

% write log
tmp_log = '';
for j = 1:length(index1)
    tmp_log = [tmp_log num2str(index1(j)) ' '];
end; clear j
tmp_log = [tmp_log(1:end-1) '\n'];
fprintf(fid,['FASTER:      channel(s) ' tmp_log]);

% clear cfg
clear cfg tmp_log

%%  interpolate artifact contaminated channels

% collect bad channels
badchan = unique([index0.c; index1]);

if ~isempty(badchan)

    fprintf(fid,['--> ' num2str(length(badchan)) ' channels interpolated\n\n']);

    cfg.method     = 'spline';
    cfg.badchannel = tmpdat.label(badchan);
    cfg.trials     = 'all';
    cfg.lambda     = 1e-5; 
    cfg.order      = 4; 
    cfg.elec       = config.elec;

    % interpolate
    tmpdat = ft_channelrepair(cfg,tmpdat);

    % clear cfg
    clear cfg
    
end

%% equalize duration of trials to the trial with fewest samples

for indTrial = 1:numel(tmpdat.trial)
    TrialLength(indTrial) = size(tmpdat.trial{indTrial},2);
end
[val, idx] = min(TrialLength);
for indTrial = 1:numel(tmpdat.trial)
    tmpdat.trial{indTrial} = tmpdat.trial{indTrial}(:,1:val);
    tmpdat.time{indTrial} = tmpdat.time{idx};
end

%%  get artifact contaminated epochs & exclude epochs
% includes: - THG_MWB_channel_x_epoch_artifacts_20140311
%           - THG_FASTER_2_epoch_artifacts_20140302
% recursive epoch exclusion!

[tmpdat index] = THG_automatic_artifact_correction_trials_20170922(tmpdat);

% write log
fprintf(fid,'(2) automatic recursive bad epoch detection:\n');
fprintf(fid,['MWB & FASTER: ' num2str(n_trl-length(index)) ' bad epoch(s) detected\n\n']);

%%  get channel x epoch artifacts

% cfg
cfg.criterion = 3;
cfg.recursive = 'yes';

[index3 parm3 zval3] = THG_FASTER_4_channel_x_epoch_artifacts_20140302(cfg,tmpdat);

% write log
fprintf(fid,'(3) automatic single epoch/channel detection:\n');
fprintf(fid,['FASTER:       ' num2str(sum(sum(index3))) ' channel(s) x trial(s) detected\n\n']);

% clear cfg
clear cfg

%%  collect and save detected artifacts & artifact correction infos

% include ArtDect.parameters

% bad channels
ArtDect.channels = badchan;

% bad trials
tmp  = zeros(n_trl,1); tmp(index,1) = 1;
badtrl = find(tmp==0);
ArtDect.trials  = badtrl; 
clear tmp

% bad single epoch(s)/channel(s) - after exclusion of bad epochs
ArtDect.channels_x_trials = index3;

% overview
ind = [1:n_trl];
ind(badtrl) = [];
tmp = ones(length(tmpdat.label),n_trl);
tmp(:,ind) = index3;
tmp(badchan,:) = 1;
ArtDect.channels_x_trials_all = tmp;
clear tmp

% save version
ArtDect.version = 'JQK 20180123';

% add to config
config.ArtDect = ArtDect;

% save config
save(fullfile(pn.history, [IDs{id}, '_task-xxxx_config.mat']),'config');

%%  finalize log

% log
fprintf(fid,'Artifact detection completed.\n');
fprintf(fid,['Information saved to: ' IDs{id} '_config.mat\n\n']);
fclose(fid);

end