|
@@ -0,0 +1,253 @@
|
|
|
+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
|