|
@@ -1,307 +0,0 @@
|
|
|
-function H_STSW_automatic_artifact_correction_180123_noTardis()
|
|
|
-
|
|
|
-%% H_SS_automatic_artifact_correction_20170922
|
|
|
-
|
|
|
-% 180326 | instead of equalizing to smallest amount of samples, interpolate
|
|
|
-% for the same timing across subjects. Note that this does not
|
|
|
-% counteract any natural jitter in the paradigm and should only
|
|
|
-% equalize the size of matrices (fill with zeros! at ends).
|
|
|
-% Artefact detection does NOT currently work with potential NaN
|
|
|
-% values that would naturally arise from the interpolation.
|
|
|
-
|
|
|
-%% initialize
|
|
|
-
|
|
|
-%restoredefaultpath;
|
|
|
-%clear all; close all; pack; clc;
|
|
|
-
|
|
|
-%% pathdef
|
|
|
-
|
|
|
-if ismac
|
|
|
- pn.study = '/Volumes/LNDG/Projects/StateSwitch/';
|
|
|
- pn.eeg_root = [pn.study, 'dynamic/data/eeg/task/A_preproc/SA_preproc_study/'];
|
|
|
- pn.dynamic_In = [pn.eeg_root, 'B_data/B_EEG_ET_ByRun/'];
|
|
|
- pn.triggerTiming= [pn.eeg_root, 'C_figures/D_TriggerTiming/'];
|
|
|
- pn.EEG = [pn.eeg_root, 'B_data/C_EEG_FT/']; mkdir(pn.EEG);
|
|
|
- pn.History = [pn.eeg_root, 'B_data/D_History/']; mkdir(pn.History);
|
|
|
- pn.logs = [pn.eeg_root, 'Y_logs/H_ArtCorr/'];
|
|
|
- % add ConMemEEG tools
|
|
|
- pn.MWBtools = [pn.eeg_root, 'T_tools/fnct_MWB/']; addpath(genpath(pn.MWBtools));
|
|
|
- pn.THGtools = [pn.eeg_root, 'T_tools/fnct_THG/']; addpath(genpath(pn.THGtools));
|
|
|
- pn.commontools = [pn.eeg_root, 'T_tools/fnct_common/']; addpath(genpath(pn.commontools));
|
|
|
- pn.fnct_JQK = [pn.eeg_root, 'T_tools/fnct_JQK/']; addpath(genpath(pn.fnct_JQK));
|
|
|
- pn.FT = [pn.eeg_root, 'T_tools/fieldtrip-20170904/']; addpath(pn.FT); ft_defaults;
|
|
|
- pn.helper = [pn.eeg_root, 'A_scripts/helper/']; addpath(pn.helper);
|
|
|
-else
|
|
|
- pn.root = '/home/mpib/LNDG/StateSwitch/WIP_eeg/SA_preproc_study/';
|
|
|
- pn.EEG = [pn.root, 'B_data/C_EEG_FT/'];
|
|
|
- pn.History = [pn.root, 'B_data/D_History/'];
|
|
|
- pn.triggerTiming= [pn.root, 'C_figures/D_TriggerTiming/'];
|
|
|
- pn.logs = [pn.root, 'Y_logs/H_ArtCorr/'];
|
|
|
-end
|
|
|
-
|
|
|
-%% define Condition & IDs for preprocessing
|
|
|
-
|
|
|
-condEEG = 'dynamic';
|
|
|
-
|
|
|
-%% define IDs for segmentation
|
|
|
-
|
|
|
-% N = 47 YAs + 53 OAs;
|
|
|
-IDs = {'1117';'1118';'1120';'1124';'1126';'1131';'1132';'1135';'1136';'1138';...
|
|
|
- '1144';'1151';'1158';'1160';'1163';'1164';'1167';'1169';'1172';'1173';...
|
|
|
- '1178';'1182';'1215';'1216';'1219';'1221';'1223';'1227';'1228';'1233';...
|
|
|
- '1234';'1237';'1239';'1240';'1243';'1245';'1247';'1250';'1252';'1257';...
|
|
|
- '1261';'1265';'1266';'1268';'1270';'1276';'1281';...
|
|
|
- '2104';'2107';'2108';'2112';'2118';'2120';'2121';'2123';'2125';'2129';...
|
|
|
- '2130';'2131';'2132';'2133';'2134';'2135';'2139';'2140';'2145';'2147';...
|
|
|
- '2149';'2157';'2160';'2201';'2202';'2203';'2205';'2206';'2209';'2210';...
|
|
|
- '2211';'2213';'2214';'2215';'2216';'2217';'2219';'2222';'2224';'2226';...
|
|
|
- '2227';'2236';'2237';'2238';'2241';'2244';'2246';'2248';'2250';'2251';...
|
|
|
- '2252';'2258';'2261'};
|
|
|
-
|
|
|
-%% loop IDs
|
|
|
-
|
|
|
-ID_unavailable = cell(length(IDs),1);
|
|
|
-for id = 1:length(IDs)
|
|
|
- try
|
|
|
- display(['processing ID ' num2str(IDs{id})]);
|
|
|
- for iRun = 1:4
|
|
|
-
|
|
|
- % load config
|
|
|
- load([pn.History, IDs{id}, '_r',num2str(iRun), '_', condEEG, '_config'],'config')
|
|
|
-
|
|
|
- if ~isfield(config, 'ArtDect')
|
|
|
-
|
|
|
- % load data
|
|
|
- load([pn.EEG, IDs{id}, '_r',num2str(iRun), '_', condEEG, '_EEG_Rlm_Fhl_rdSeg'],'data_EEG');
|
|
|
-
|
|
|
- dataICA = load([pn.EEG, IDs{id}, '_', condEEG, '_EEG_Rlm_Fhl_Ica'],'data');
|
|
|
-
|
|
|
- data = data_EEG; clear data_EEG;
|
|
|
- data.elec = dataICA.data.elec;
|
|
|
- data.chanlocs = dataICA.data.chanlocs; clear dataICA;
|
|
|
-
|
|
|
- %% ------------------ ARTIFACT DETECTION - PREPARATION ----------------- %%
|
|
|
-
|
|
|
- config.elec = data.elec;
|
|
|
-
|
|
|
- %% ICA (from weights)
|
|
|
-
|
|
|
- % ica config
|
|
|
- cfg.method = 'runica';
|
|
|
- cfg.channel = {'all','-ECG','-A2'};
|
|
|
- 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','-IOR','-LHEOG','-RHEOG','-A1'};
|
|
|
- cfg.demean = 'yes';
|
|
|
-
|
|
|
- % remove channels
|
|
|
- tmpdat = ft_preprocessing(cfg,data);
|
|
|
-
|
|
|
- % clear cfg & data variable
|
|
|
- clear cfg data
|
|
|
-
|
|
|
-
|
|
|
- %% ------------------------- ARTIFACT DETECTION ------------------------ %%
|
|
|
-
|
|
|
- % open log file
|
|
|
- fid = fopen([pn.logs, 'log_' IDs{id}, '_r',num2str(iRun), '_' condEEG '_ArtCorr.txt'],'a');
|
|
|
-
|
|
|
- % write log
|
|
|
- % fprintf(fid,['********************* ' BLOCK ' *********************\n']); % Undefined function or variable 'BLOCK'. but see: %% H_prep_data_for_analysis_20141218
|
|
|
-
|
|
|
- fprintf(fid,['********************* BLOCK *********************\n']);
|
|
|
- fprintf(fid,['function: H_SS_automatic_artifact_correction_20170922.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
|
|
|
-
|
|
|
- %% interpolate data to requested time (necessary due to acquisition jitter)
|
|
|
-
|
|
|
- % 2160 has only 62 trials, i.e. the first two trials were
|
|
|
- % not recorded.
|
|
|
-
|
|
|
- timeRequested = -1.5:1/500:9.566;
|
|
|
-
|
|
|
- for indTrial = 1:size(tmpdat.trial,2)
|
|
|
- recTime = tmpdat.time{indTrial};
|
|
|
- recData = tmpdat.trial{indTrial};
|
|
|
- for indChannel = 1:size(recData,1)
|
|
|
- interpData(indChannel,:) = interp1(recTime,recData(indChannel,:),timeRequested,'linear');
|
|
|
- end
|
|
|
- interpData(isnan(interpData)) = 0;
|
|
|
- tmpdat.trial{indTrial} = interpData;
|
|
|
- tmpdat.time{indTrial} = timeRequested;
|
|
|
- end
|
|
|
-
|
|
|
- %figure; plot(timeRequested,interpData, 'k'); hold on; plot(recTime,recData(1,:), 'r')
|
|
|
-
|
|
|
- %% 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([pn.History, IDs{id}, '_r',num2str(iRun), '_', condEEG, '_config'],'config')
|
|
|
-
|
|
|
- %% finalize log
|
|
|
-
|
|
|
- % log
|
|
|
- fprintf(fid,'Artifact detection completed.\n');
|
|
|
- fprintf(fid,['Information saved to: ' IDs{id} '_config.mat\n\n']);
|
|
|
- fclose(fid);
|
|
|
-
|
|
|
-
|
|
|
- else
|
|
|
- disp([IDs{id}, ' Run ',num2str(iRun), ' has already run through automatic artifact detection']);
|
|
|
- end % if data was not processed yet
|
|
|
- end % run
|
|
|
- catch ME
|
|
|
- warning('Error occured. Please check.');
|
|
|
- rethrow(ME)
|
|
|
- fprintf('\n ID not availble! Skip! \n')
|
|
|
- ID_unavailable{id,1} = (IDs{id});
|
|
|
- end
|
|
|
-
|
|
|
-end; clear id
|
|
|
-
|
|
|
-end
|