123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249 |
- function a6_automatic_artifact_correction(id, rootpath)
- if ismac % run if function is not pre-compiled
- id = '1'; % 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]);
- 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
- %% 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
|