基于眼动数据在ICA中识别并移除EEG中的眼动伪迹

开篇之前，先在小黑板上划重点：

欲用此法，在数据采集时必须同时记录眼动数和EEG数据！

欲用此法，在数据采集时必须同时记录眼动数和EEG数据！！

欲用此法，在数据采集时必须同时记录眼动数和EEG数据！！！

如果眼动数据和EEG数据没有同时记录，那接下来的操作可能不完全适用于你的数据。尤其是眼动数据不可或缺。已经和工具包（或示例代码）的开发者 Olaf 联系了，如果缺少眼动数据，那目前使用此方法仍有较大难度，但未来可期。

鉴于此部分内容较多，我们将会推出的是系列组稿，具体时间安排未知，因为小编也不知道自己什么时候有时间来更新，毕竟组稿不易，且行且珍惜，且看且转推，点赞、在看、转发是我的源泉和动力。

以上闲篇结束，步入正题。

1.  准备工作

1.1 环境配置

文中的方法是基于Matlab、EEGlab和FiledTrip的数据分析环境，并辅之EYE-EEG toolbox(EEGLAB插件，其安装方法见 >> TFA：一个基于EEGLAB的时频分析工具包)、Unflod(EEG分析集成工具包，详见 参考资料 [3]，后续会有推文对其的使用进行详细介绍) ，相关工具包的下载地址：

• EEGLab

https://eeglab.org/

• FiledTrip

https://www.fieldtriptoolbox.org/

• EYE-EEG toolbox

https://www.eyetracking-eeg.org/

• Unflod

https://www.unfoldtoolbox.org/

图源：https://www.eyetracking-eeg.org/

为了方便使用和下载，我们对上述工具包进行了下载，示例中的Matlab代码也一并上传到了百度云，链接和提取密码如下：

链接: https://pan.baidu.com/s/1F7hB_n0-s8f413s1lEqf5Q 

提取码: 8sv3 

*如遇链接失效，请后台联系小编处理。

1.2 眼动数据

1.2.1 眼动数据导出为文本格式

将不同设备采集的不同数据格式(如：*.edf or *.idf)的眼动数据转换为text文本数据（ASCII）。（*此部分教程我们会在后续的教程中推出，等不及的小伙伴也可以先去EYE-EEG工具包官网查阅。）

1.2.2 眼动数据导入Matlab并另存为.mat文件

将导出的后的眼动数据转存为Matlab数据格式并另存后备用，具体操作为：

• 点击 Eyetracker > Parse eyetracker raw data 并选择相应的眼动文本数据（如下图所示）
• 选择相应的文件夹，保存.mat格式的眼动数据。

2. 实例代码解读

*本实例代码来源于Federica Degno等人公开共享的数据分析代码(见参考资料[5] [6])。数据分析流程示意图如下图所示。鉴于篇幅关系，仅对核心过程进行注释和解读。鉴于个人能力有限，若有不当之处，恳请批评指正。

图源：https://doi.org/10.1525/collabra.18032

如图所示，基于眼动数据在ICA中识别并移除EEG中的眼动伪迹的核心步骤主要有以下11个环节：

Step 01:对眼动数据进行重新编码

Step 02: 对EEG数据进行重新编码

Step 03: 同步眼动数据和EEG数据

Step 04: （如果必要）对坏导进行检测并删除。

Step 05: 对连续的EEG数据进行滤波。

Step 06: 对记录不完整或未记录到等的眼动数据以及任务相关EEG记录期间的眼动进行检测。

Step 07: 优化ICA中的EEG训练数据集。

Step 08: 在EEG训练数据集结果的基础上进行ICA运算。

Step 09: 在EEG训练数据集和ICA计算结果的基础上进行眼动伪迹等伪迹成分的检测。

Step 10: ICA转换和眼动伪迹IC成分标记。

Step 11: 在连续的EEG原始数据中对标记为眼动伪迹等IC成分进行移除。

接下来，我们将对每一步的代码进行详细的展示和注解。

---

Step 01:对眼动数据进行重新编码

%% 1. RECODING OF EM DATA
%Here we replace all trial onset triggers (i.e., 'TTLon') with the same code (i.e., 601 regardless of the trial index) and all trial offset triggers (i.e., 'TTLoff') with the same code (i.e., 602 regardless of the trial index),
%we remove practice and filler trials from the dataset,
%and we covert the raw EM data from ASCII to MATLAB format (.mat).

clear
eeglab

rawPATH = 'C:\Users\Federica\Documents\DegnoLobergLiversedge\EMdata\ASCIIfiles\'; %Folder with data to process.
cd(rawPATH); %Go to folder.
allraw = dir('*.asc'); %Detect all files with this extension in the folder.

tempPATH = 'C:\Users\Federica\Documents\DegnoLobergLiversedge\EMdata\TemporaryFiles\'; %Folder where to save the parsed data (EYE1.mat and EYE2.mat).
savePATH = 'C:\Users\Federica\Documents\DegnoLobergLiversedge\EMdata\EMCONVERTED\'; %Folder where to save the processed data.

for master_iterator =1:length(allraw) %Start the loop.

    %Define strings to load/save files
    sourcefile = strcat(rawPATH,allraw(master_iterator).name); %Create the string to load eye-data.
    fileEYE1 = strcat(tempPATH,'EYE1.mat'); %Create the string to save the EYE1 data.
    fileEYE2 = strcat(tempPATH,'EYE2.mat'); %Create the string to save the EYE2 data.

    %Parse raw eye tracking data from Eyelink
    ET = parseeyelink(sourcefile,fileEYE1,'TTLon'); %Keyword string for trial onset: 'TTLon'. Saved as EYE1.mat
    ET = parseeyelink(sourcefile,fileEYE2,'TTLoff'); %Keyword string for trial offset: 'TTLoff'. Saved as EYE2.mat

    %Recode trial onset and offset messages
    cd(tempPATH); %Move to folder where the parsed data are saved.
    load('EYE1.mat') %Load the temporary file with the trial onset message (i.e., TTLon).
    event2=event; %TTLon events to matrix (122 events in the current experiment).
    load('EYE2.mat') %Load the temporary file of the trial offset message (i.e., TTLoff).
    event2=[event2; event]; %Combine TTLon (event2) and TTLoff (event) messages into same structure (244 events in the current experiment: the first 122 are TTLon, from 123 are TTLoff).
    event2(1:122,2)=601;% Now create synchronization code 601 for TTLon.
    event2(123:end,2)=602;% Now create synchronization code 602 for TTLoff.
    event2=sortrows(event2);% Sort events by time (first column of the matrix event2). Now 601 and 602 are alternated.

    %Remove practice and filler trials
    %In this experiment, practice: first 10 trials; fillers: TTLon/off11, TTLon/off39, TTLon/off67, TTLon/off95.
    %Because each trial has 2 events (i.e., TTLon and TTLoff), we need to remove
    %practice trials: events 1-20,
    %filler trial 11 (TTLon/off11): events 21-22, filler trial 39 (TTLon/off39): events 77-78,
    %filler trial 67 (TTLon/off67): events 133-134, filler trial 95 (TTLon/off95): events 189-190.
    b = [1:20 21 22 77 78 133 134 189 190]; %Create vector of practice and filler trial event indexes.
    event2(b,:) = []; %Remove practice and filler trials. Now you should obtain 216 rows, as we only have 108 experimental trials in the current experiment.
    event=event2; %Replace the original event structure with our new synchronization events (216 rows).
    clear event2

    %Save the processed data
    newStr=regexprep(allraw(master_iterator).name,'.asc',''); %Remove .asc from filename.
    savename=strcat(newStr,'.mat'); %Add .mat to the filename.
    cd(savePATH); %%Go to folder.
    save(savename) %Save the synchronizable eye-file SYNCHRONIZABLE? OR RECODED EM DATA? CAN WE SAVE WITH savePATH?
end

clear

Step 02: 对EEG数据进行重新编码

%% 2. RECODING OF EEG DATA
%Here we add recoded duplicates of original events (replacing all trial onset triggers (i.e., 'TTLon') with the same code (i.e., 601 regardless of the trial index) and all trial offset triggers (i.e., 'TTLoff') with the same code (i.e., 602 regardless of the trial index)),
%we remove practice and filler trials from the dataset,
%and we covert the raw EEG data from NEUROSCAN (.cnt) to MATLAB EEGLAB format (.set).

clear
eeglab

rawPATH = 'C:\Users\Federica\Documents\DegnoLobergLiversedge\EEGdata\'; %Folder with data to process.
cd(rawPATH); %Go to folder.
allraw = dir('*.cnt'); %Detect all files with this extension in the folder.

savePATH = 'C:\Users\Federica\Documents\DegnoLobergLiversedge\EEGdata\EEGSET\'; %Folder where to save the processed data.

for master_iterator =1:length(allraw) %Start the loop.

    %Load each dataset
    EEG = pop_loadcnt(allraw(master_iterator).name , 'dataformat', 'auto', 'memmapfile', ''); %Load each Neuroscan .cnt file. If you encounter problems here, check that you have Neuroscan import plugin in your EEGLAB build
    [ALLEEG EEG CURRENTSET] = pop_newset(ALLEEG, EEG, 0,'gui','off'); %Save EEG dataset structure information.
    EEG = eeg_checkset( EEG ); %Check the consistency of the fields of the EEG dataset.

    %Remove practice and filler trials
    allraw(master_iterator).events =length(EEG.event); %Save the number of events.
    b = [1:10 11 66 121 176]; %Create vector of practice (i.e., the first 10 rows) and filler trials (trial 11: event 11, trial 39: event 66, trial 67: event 121, trial 95: event 176).
    EEG.event(b)=[]; %Remove practice and filler trials. Now you should obtain 216 rows, as we only have 108 experimental trials in the current experiment.
    event2=EEG.event; %Copy events to secondary structure.

    %Recode trial onset and offset triggers
    %In the current experiment, the EEG.event.type (now 'event2.type') comprised
    %the values 1, 2, 3, 4, 5, 6 (which represent the marker sent at the
    %beginning of each trial and corresponding to the condition number of that
    %trial), and the values 101, 102, 103, 104, 105, 106 (which represent the marker sent at the
    %end of each trial and corresponding to the condition number of that
    %trial). 
    for i = 1:length(event2) %Loop for the length of the structure.
        if event2(i).type < 50    %If the event type is less than 50 (i.e., it is the marker representing the trial onset),
           event2(i).type = 601;  %then replace the value with the new event2.type equal to 601 (i.e., syncronization start)
        else                      %otherwise,
            event2(i).type = 602; %replace the value with 602 (i.e., syncronization end).
        end
    end
    EEG.event = [EEG.event event2]; %Here we merge new and old EEG events.

    %Save the processed data
    newStr = regexprep(allraw(master_iterator).name,'.cnt',''); %Remove '.cnt' from filename.
    savename =strcat(newStr,'_raw.set'); %Add '_raw.set' to the filename.
    EEG = pop_saveset( EEG, 'filename',savename,'filepath',savePATH); %Save as a new dataset.
    STUDY = []; CURRENTSTUDY = 0; ALLEEG = []; EEG=[]; CURRENTSET=[]; %Initialize EEGLAB/STUDY variables. This clears the memory of the data after saving, so only one file is kept in memory, for fluid computations.
end

Step 03: 同步眼动数据和EEG数据

%% 3. SYNCHRONISATION OF EM AND EEG DATA
%Here we import the eye tracking (ET) data as additional ET data channels in the EEG,
%we synchronise EM and EEG data based on common events present in both recordings,
%and we obtain an estimate of the quality of synchronization between the two recordings.

clear
eeglab

rawPATH1 = 'C:\Users\Federica\Documents\DegnoLobergLiversedge\EEGdata\EEGSET\'; %Folder with data to process.
cd(rawPATH1); %Go to folder.
allset = dir('*.set'); %Detect all files with this extension in the folder.

rawPATH2 = 'C:\Users\Federica\Documents\DegnoLobergLiversedge\EMdata\EMCONVERTED\'; %Folder with data to process.
cd(rawPATH2); %Go to folder.
alleye = dir('*.mat'); %Detect all files with this extension in the folder.

savePATH = 'C:\Users\Federica\Documents\DegnoLobergLiversedge\EEGdata\SYNC\'; %%Folder where to save the processed data.
savePATHfig = 'C:\Users\Federica\Documents\DegnoLobergLiversedge\EEGdata\SYNC\SYNCFIGURES\'; %Folder where to save the figures of synchronisation quality.

%Load the channel information as a variable.
%If your .set data already contains channel locations, this line of code is not necessary.
load('C:\Users\Federica\Documents\DegnoLobergLiversedge\chanlocs.mat')  

for master_iterator =1:length(allset) %Start the loop.

    %Load each dataset
    EEG = pop_loadset('filename',allset(master_iterator).name,'filepath',rawPATH1); %Load each existing dataset.
    [ALLEEG, EEG, CURRENTSET] = eeg_store( ALLEEG, EEG, 0 ); %Save EEG dataset structure information.

    %Add the channel location information to the EEG structure
    EEG.chanlocs=chanlocs;

    %Define strings to load EM files
    loadname = strcat(rawPATH2,alleye(master_iterator).name); %EM data file name.

    %Import and synchronise ET data
    %Here 601 (STARTevent) and 602 (ENDevent) are the common events between EM and EEG data.
    EEG = pop_importeyetracker(EEG,loadname,[601 602] ,[2:5] ,{'R_GAZE_X' 'R_GAZE_Y' 'R_AREA' 'INPUT'},1,1,0,1,5); %Here we set the searchRadius to identify matching triggers to 5 samples.
    [ALLEEG EEG CURRENTSET] = pop_newset(ALLEEG, EEG, 1,'overwrite','on','gui','off'); %Save EEG dataset structure information.
    EEG = eeg_checkset( EEG ); %Check the consistency of the fields of the EEG dataset.

    %Save the processed data
    name = regexprep(allset(master_iterator).name,'_raw.set','');  %Remove '_raw.set'.
    savename = strcat(name,'_sync.set'); %Generate the new file name for saving.
    EEG = pop_saveset( EEG, 'filename',savename,'filepath',savePATH); %Save.

    %Save the synchronization picture
    cd(savePATHfig); %Go to folder.
    name = regexprep(allset(master_iterator).name,'.set',''); %Remove .set. 
    savename = strcat(name,'_syncresult.fig'); %Generate the new file name for saving.
    saveas(gcf,savename); %Save the current figure as FIG-file, i.e., Matlab figure.
    close gcf 
end

Step 04: （如果必要）对坏导进行检测并删除。注：进行此步骤必须将 Unfold工具包添加到路径。（关于Unfold工具包的安装和使用细节我们会在近期推出，尽情期待！）

%% 4. DETECTION AND REMOVAL OF BAD EEG CHANNELS, IF ANY
%Here we check whether there is any EEG channel with a high percentage of 'bad' data (i.e., with extreme values).
%If so, we identify those channels and remove them from the dataset.
%To identify the bad EEG channels, we use the C.R.A.P algorithm channel by
%channel, and see whether any channel has considerably larger proportion of bad data than the others.,
%here greater than 3 standard deviations from mean amount of bad data is considered to be indicative of a bad channel.  

%Unfold needs to be in the path for this part of the script to work.
%If an error occurs when running uf_continuousArtifactDetect, please ensure
%that erplab plugin is removed from the eeglab plugin folder and path.

clear
eeglab

rawPATH = 'C:\Users\Federica\Documents\DegnoLobergLiversedge\EEGdata\SYNC\'; %Folder with data to process.
cd(rawPATH); %Go to folder.
myDataSet = dir('*_sync.set'); %Detect all files with this extension in the folder.

savepath = 'C:\Users\Federica\Documents\DegnoLobergLiversedge\EEGdata\BAD_CH\'; %Folder where to save the processed data.

for master_iterator = 1:length(myDataSet) %Start the loop.

    %Load each dataset
    EEG = pop_loadset('filename',myDataSet(master_iterator).name,'filepath',rawPATH); %Load each dataset.
    [ALLEEG, EEG, CURRENTSET] = eeg_store( ALLEEG, EEG, 0 ); %Save EEG dataset structure information.
    EEG = eeg_checkset( EEG );

    %Save original channel locations
    EEG.urchanlocs = EEG.chanlocs; %This will allow to keep track of the removed channels.

    %Look up the EEG channel indexes
    EEGch = find(strcmp({EEG.chanlocs.type},'EEG')); 
    badsamples_per_channel = []; %Create a variable for holding the amount of bad data in each channel.

    %Use the C.R.A.P. algorithm, as implemented in the Unfold toolbox, channel by channel to detect intervals with bad EEG data
    %Default moving window width of 2000 ms.
    for channel_iterator= EEGch %Start loop, to go through the channel indexes.
        bad_intervals = uf_continuousArtifactDetect(EEG,'amplitudeThreshold',250,'channels',channel_iterator); %Get intervals with bad data in each channel (i.e., the intervals where the amplitude differences are greater than 250 microVolts).
        bad_samples=sum(bad_intervals(:,2)-bad_intervals(:,1)); %Calculate the sum of bad data in each channel (in samples).
        badsamples_per_channel(channel_iterator)=bad_samples; %Store sum of bad data in a channel to a corresponding cell.
    end

    %Identify and mark the 'bad' EEG channels
    %That is, any channel with a considerably larger proportion of bad data than the others - here greater than 3 standard deviations from mean amount of bad data.  
    badsample_percentage_per_channel = (badsamples_per_channel *100)/ (EEG.pnts); %Calculate the percentage of bad data in each channel.
    badsample_percentage_zscored = zscore(badsample_percentage_per_channel); %Zscore the percentages for standardizing the amount of bad data.
    badchannel_index=find(badsample_percentage_zscored>3); %Get indexes of channels with Zscored bad data percentage higher than 3.

    %Remove the 'bad' EEG channels     
    EEG = pop_select( EEG, 'nochannel',badchannel_index); %Remove the 'bad' EEG channels. If badchannel_index_is_empty, then none are removed.

    %Mark the 'bad' EEG channels down to the datafile
    %and note down which EEG electrodes have been removed, and their percentage of 'bad' data
    EEG.etc.removed_bad_channels_index=badchannel_index; %Add the indexes of channels to be rejected.   
    EEG.etc.percentage_of_bad_data_per_channel= badsample_percentage_per_channel; %Save the bad sample percentages per channel.

    %Save the processed data
    name = regexprep(myDataSet(master_iterator).name,'_sync.set','');
    savename = strcat(name,'_sync_badCH.set');
    EEG = pop_saveset( EEG, 'filename',savename,'filepath',savepath); %Save the file.
    [ALLEEG EEG] = eeg_store(ALLEEG, EEG, CURRENTSET);
    STUDY = []; CURRENTSTUDY = 0; ALLEEG = []; EEG=[]; CURRENTSET=[];
end

Step 05: 对连续的EEG数据进行滤波。在此处我们进行的是0.1Hz 和 100Hz。

%% 5. FILTER ORIGINAL EEG CONTINUOUS DATA
%Here we filter the original EEG data with 0.1Hz and 100Hz.
%with the EYE-EEG extension, to avoid filtering the new ET channels.

clear
eeglab

rawPATH = 'C:\Users\Federica\Documents\DegnoLobergLiversedge\EEGdata\BAD_CH\'; %Folder with data to process.
cd(rawPATH); %Go to folder.
myDataSet = dir('*_sync_badCH.set'); %Detect all files with this extension in the folder.

savePATH = 'C:\Users\Federica\Documents\DegnoLobergLiversedge\EEGdata\FILT1\'; %Folder where to save the processed data.

%Parameters for filtering:
HIGHPASS = 0.1;
%LOWPASS = 100; %Our data were already low-pass filtered at 100 Hz, so we do not apply any low-pass filter here.

for master_iterator =1:length(myDataSet) %Start the loop.

    %Load each dataset
    EEG = pop_loadset('filename',myDataSet(master_iterator).name,'filepath',rawPATH); %Load each dataset.
    [ALLEEG, EEG, CURRENTSET] = eeg_store( ALLEEG, EEG, 0 ); %Save EEG dataset structure information.

    %Identify the indexes of channels to filter (EEG and EOG, not ET)
    %Note that this is necessary because we might have removed some 'bad' EEG channels.
    EEG_EOG_ch = [find(strcmp({EEG.chanlocs.type},'EEG')) find(strcmp({EEG.chanlocs.type},'EOG'))]; %Identify the indexes of EEG and EOG channels.

    %High pass filter
    %Note that applytochannels is implemented as a wrapper function, thus,
    %we need to build an adapting string that calls the filter function from within applytochannels  
    funcvar = strcat('pop_eegfiltnew(EEG, [],', num2str(HIGHPASS), ', [], 1, [], 1);'); %Using New Basic FIR filter
    EEG = applytochannels(EEG, EEG_EOG_ch ,funcvar);  %Apply high-pass filter to the EEG and EOG channels

    [ALLEEG EEG CURRENTSET] = pop_newset(ALLEEG, EEG, 1,'overwrite','on','gui','off'); %Save EEG dataset structure information.
    EEG = eeg_checkset( EEG ); %Check the consistency of the fields of the EEG dataset.
    close gcf %Close current picture

    %Low pass filter
    %Our data were already low-pass filtered at 100 Hz, so we do not apply any low-pass filter here
    %funcvar = strcat('pop_eegfiltnew(EEG,', num2str(LOWPASS), ',[], [], 1, [], 1);');  %Using New Basic FIR filter
    %EEG = applytochannels(EEG, EEG_EOG_ch ,funcvar); %Apply low-pass filter to the EEG and EOG channels
    %[ALLEEG EEG CURRENTSET] = pop_newset(ALLEEG, EEG, 1,'overwrite','on','gui','off'); %Save EEG dataset structure information.

    %EEG = eeg_checkset( EEG ); %Check the consistency of the fields of the EEG dataset.
    %close gcf %Close current picture

    %Save the processed data
    name = regexprep(myDataSet(master_iterator).name,'_sync_badCH.set','');  %Remove '_sync.set'.
    savename = strcat(name,'_filt.set'); %Generate the new file name for saving.
    EEG = pop_saveset( EEG, 'filename',savename,'filepath',savePATH); %Save.
    [ALLEEG EEG] = eeg_store(ALLEEG, EEG, CURRENTSET); %Store the current dataset in ALLEEG structure.
end

Step 06: 对记录不完整或未记录到等的眼动数据以及任务相关EEG记录期间的眼动进行检测。

%% 6. DETECTION OF BAD ET AND EEG INTERVALS
%Here we detect the intervals of continuous data with extreme (out-of-range) values as recorded by the ET system.
%A marker "bad_ET" will be added to EEG.event for each bad data interval.
%It is important to only detect, and not remove, these segments if one is going to conduct deconvolution later on.

clear
eeglab

rawPATH = 'C:\Users\Federica\Documents\DegnoLobergLiversedge\EEGdata\FILT1\'; %Folder with data to process.
cd(rawPATH); %Go to folder.
myDataSet = dir('*_filt.set'); %Detect all files with this extension in the folder.

savePATH = 'C:\Users\Federica\Documents\DegnoLobergLiversedge\EEGdata\BAD_ET1\'; %Folder where to save the processed data.

for master_iterator =1:length(myDataSet) %Start the loop.

    %Load each dataset
    EEG = pop_loadset('filename',myDataSet(master_iterator).name,'filepath',rawPATH); %Load each dataset.
    [ALLEEG, EEG, CURRENTSET] = eeg_store( ALLEEG, EEG, 0 ); %Save EEG dataset structure information.

    %Select channels
    chanlocs = EEG.chanlocs;
    idx_X = find(strcmp({chanlocs.labels},'R-GAZE-X')); %This corresponds to the imported ET channel 'R_GAZE_X'.
    idx_Y = find(strcmp({chanlocs.labels},'R-GAZE-Y')); %This corresponds to the imported ET channel 'R_GAZE_Y'.
    EEGch = find(strcmp({chanlocs.type},'EEG')); %This corresponds to the EEG channels.

    %Identify the bad ET intervals from the ET record
    EEG = pop_rej_eyecontin(EEG,[(idx_X) (idx_Y)] ,[1 1] ,[1024 768] ,100,2); %Identify (method '2') the bad ET intervals.

    %Detect intervals with bad data with C.R.A.P. as modified by Benedikt Ehinger in the Unfold toolbox
    %Use the algorithm on EEG channels, with default moving window width of 2000 ms.
    craprej = uf_continuousArtifactDetect(EEG,'amplitudeThreshold',250, 'channels',EEGch); %Detect intervals with amplitude differences in the EEG signal greater than 250 microvolts.
    EEG.craprej = craprej;

    %Estimate the ET-EEG synchronisation lag via the EYE-EEG cross-correlation 
    try
        substitute_HEOG_left = find(strcmp({chanlocs.labels},'f7')); %Look up index of f7.
        substitute_HEOG_right = find(strcmp({chanlocs.labels},'f8')); %Look up index of f8.
        EEG = pop_checksync(EEG,idx_X,substitute_HEOG_left,substitute_HEOG_right,0); %Check synchronisation lag.
        clear substitute_HEOG_left substitute_HEOG_right
        [M,Index]=max(EEG.etc.xcorr_eyeeeg(1,:)) ;
        EEG.etc.estimate_of_lag=EEG.etc.xcorr_eyeeeg(2,Index); % 
        clear M Index
    catch %if error is encountered in lag estimation, put a note to the EEG.etc.xcorr_eyeeeg cell
        EEG.etc.xcorr_eyeeeg = 'Estimation of the ET-EEG failed, probably due either f7 or f8 -channel missing';
    end
    %If f7 and/or f8 channels are missing, the researcher will need to use different electrodes.

    %Save the processed data
    name = regexprep(myDataSet(master_iterator).name,'_filt.set','');  %Remove '_sync.set'.
    savename = strcat(name,'_badET1.set'); %Generate the new file name for saving.
    EEG = pop_saveset( EEG, 'filename',savename,'filepath',savePATH); %Save.
    [ALLEEG, EEG] = eeg_store(ALLEEG, EEG, CURRENTSET);
    STUDY = []; CURRENTSTUDY = 0; ALLEEG = []; EEG=[]; CURRENTSET=[]; %Reset the eeglab data-structures.
end

Step 07: 优化ICA中的EEG训练数据集。

%% 7. OPTIMISED ICA TRAINING DATA
%This part of code is based on Dimigen, O. (2020). Optimizing the ICA-based removal of ocular EEG artifacts during from viewing experiments. NeuroImage 
%Here we create the training data to run an optimised ICA.

clear
eeglab

rawPATH = 'C:\Users\Federica\Documents\DegnoLobergLiversedge\EEGdata\BAD_ET1\'; %Folder with data to process.
cd(rawPATH); %Go to the folder.
myDataSet  = dir('*_badET1.set'); %Detect all files with this extension in the folder.

savePATH = 'C:\Users\Federica\Documents\DegnoLobergLiversedge\EEGdata\TRAINING_DATA\'; %Folder where to save the processed data.

%Set up the constants
HIPASS = 3;                  %Filter's passband edge (in Hz).
OW_FACTOR = 0.3;             %Value for overweighting of SPs (0.3 = add spike potentials corresponding to 30% of original data length)
REMOVE_EPOCHMEAN = true;     %Mean-center the appended peri-saccadic epochs.

for master_iterator =1:length(myDataSet) %Start the loop.

    %Load each dataset    
    EEG = pop_loadset('filename',myDataSet(master_iterator).name,'filepath',rawPATH); %Load each dataset.
    [ALLEEG, EEG, CURRENTSET] = eeg_store( ALLEEG, EEG, 0 ); %Save EEG dataset structure information.

    %Indices of all EEG and EOGchannels
    EEG_CHANNELS = [find(strcmp({EEG.chanlocs.type},'EEG')) find(strcmp({EEG.chanlocs.type},'EOG'))];         % indices of all EEG and EOGchannels

    %Create copy of data used as training data and high pass filter it
    fprintf('\nCreating optimized ICA training data...')
    EEG_training = EEG; %Make a copy of the data.

    %High pass filter
    %Note that applytochannels is implemented as a wrapper function, thus,
    %we need to build an adapting string that calls the filter function from within applytochannels  
    funcvar = strcat('pop_eegfiltnew(EEG, [],', num2str(HIPASS), ', [], 1, [], 1);'); %Using New Basic FIR filter
    EEG_training = applytochannels(EEG_training, EEG_CHANNELS ,funcvar); %Apply high-pass filter to the EEG and EOG channels.
    close gcf % close current picture

    %Remove the previously defined segments of bad data
    EEG_training = pop_select(EEG_training,'nopoint',EEG_training.craprej); 

    %Cut training data into epochs (e.g., around stimulus onsets)
    %Here stimulus onset triggers are "601", and we cut from stimulus onset to +2 seconds.
    EEG_training = pop_epoch(EEG_training,{'601'},[0 2],'newname', 'ICA Training FILE', 'epochinfo', 'yes');
    [ALLEEG EEG_training CURRENTSET] = pop_newset(ALLEEG, EEG_training,1,'overwrite','on','gui','off'); %v Overwrite the existing file.

    %Overweight spike potentials
    %Repeatedly append intervals around saccade onsets (-20 to +10 ms) to training data.
    EEG_training = pop_overweightevents(EEG_training,'R_saccade',[-0.02 0.01],OW_FACTOR,REMOVE_EPOCHMEAN);  %This has a known bug related to the remove mean and new versions of EEGLAB, load most current version of EYE-EEG from github to have a patched version.

    %Save the processed data
    name = regexprep(myDataSet(master_iterator).name,'_badET1.set',''); %Remove '_sync_filt.set'.
    savename = strcat(name,'_training30.set'); %Create the savename: training30 represent: (0 2) from stimulus onset and 30% overweight.
    EEG_training = pop_saveset(EEG_training, 'filename',savename,'filepath',savePATH); %Save in the folder ready for ICA-decomposition.
    [ALLEEG EEG] = eeg_store(ALLEEG, EEG, CURRENTSET); %Save EEG dataset structure information.
    STUDY = []; CURRENTSTUDY = 0; ALLEEG = []; EEG=[];EEG_training=[]; CURRENTSET=[]; %Reset the eeglab data-structures.
end

Step 08: 在EEG训练数据集结果的基础上进行ICA运算。

%% 8. RUN ICA ON TRAINING DATA
%Here we run the extended Infomax ICA on the optimised ICA training data.

clear
eeglab

rawPATH = 'C:\Users\Federica\Documents\DegnoLobergLiversedge\EEGdata\TRAINING_DATA\'; %Folder with data to process.
cd(rawPATH); %Go to folder.
myDataSet  = dir('*_training30.set'); %Detect all files with this extension in the folder.

savePATH = 'C:\Users\Federica\Documents\DegnoLobergLiversedge\EEGdata\INFOMAX_ICA\'; %Folder where to save the processed data.

for master_iterator =1:length(myDataSet ) %Start the loop.

    %Load each dataset
    EEG = pop_loadset('filename',myDataSet (master_iterator).name,'filepath',rawPATH); %Load each existing dataset.
    [ALLEEG, EEG, CURRENTSET] = eeg_store( ALLEEG, EEG, 0 ); %Save EEG dataset structure information.

    %Select channels   
    MyCHANNELS = [find(strcmp({EEG.chanlocs.type},'EEG')) find(strcmp({EEG.chanlocs.type},'EOG'))]; %Indices of all EEG and EOG channels.

    %Downsample the data
    EEG_training2 = pop_resample(EEG, 500); %To reduce computational resources and time. 
    [ALLEEG EEG_training2 CURRENTSET] = pop_newset(ALLEEG, EEG_training2, 1,'overwrite','on','gui','off'); 

    %Run ICA
    fprintf('\nRunning ICA on optimized training data...')
    EEG_training2 = pop_runica(EEG_training2,'chanind',MyCHANNELS,'extended',1,'interupt','on','chanind',MyCHANNELS); %Run the extended Infomax ICA.

    %Transfer the ICA information from training2 data to training data
    EEG.icasphere = EEG_training2.icasphere;
    EEG.icaweights = EEG_training2.icaweights;  
    EEG.icachansind = EEG_training2.icachansind; 
    EEG = eeg_checkset(EEG);  
    clear('EEG_training2')

    %Save the processed data
    name = regexprep(myDataSet(master_iterator).name,'_training30.set',''); %Remove previous extension name.
    savename = strcat(name,'_training30_icasegments.set'); %Create the savename. _training30_icasegments (30%).
    EEG = pop_saveset(EEG, 'filename',savename,'filepath',savePATH); %Save in the folder ready for ICA-decomposition.
    [ALLEEG EEG] = eeg_store(ALLEEG, EEG, CURRENTSET); %Save EEG dataset structure information.
    STUDY = []; CURRENTSTUDY = 0; ALLEEG = []; EEG=[]; CURRENTSET=[]; %Reset the eeglab data-structures.
end

Step 09: 在EEG训练数据集和ICA计算结果的基础上进行眼动伪迹等伪迹成分的检测。在这个环节需要使用到EYE-EEG工具包。

%% 9. DETECTION OF OCULAR ICs ON TRAINING DATA
%Here we detect ocular-artifacts with the temporal covariance criterion from the EYE-EEG extension.

clear
eeglab

rawPATH= 'C:\Users\Federica\Documents\DegnoLobergLiversedge\EEGdata\INFOMAX_ICA\'; %Folder with data to process.
cd(rawPATH); %Go to folder.
myDataSet = dir('*_training30_icasegments.set'); %Detect all files with this extension in the folder.

rawpathFIG = 'C:\Users\Federica\Documents\DegnoLobergLiversedge\EEGdata\INFOMAX_ICA\ICAFIGURES\';

savePATH = 'C:\Users\Federica\Documents\DegnoLobergLiversedge\EEGdata\IC_DETECTION\'; %Folder where to save the processed data.

for master_iterator = 1:length(myDataSet) %Start the loop.

    %Load each dataset
    EEG = pop_loadset('filename',myDataSet(master_iterator).name,'filepath',rawPATH); %Load EEG training dataset with the ICs weights.
    [ALLEEG, EEG, CURRENTSET] = eeg_store( ALLEEG, EEG, 0 ); 
    close gcf %To close the eeglab GUI window.

    %Detect ocular ICs
    %Please note that there is a known bug with pop+eyetrackerica.
    %Do not use this function on resampled EEG data. 
    [EEG vartable] = pop_eyetrackerica(EEG,'R_saccade','R_fixation',[10 10] ,1.1,3,1,4); %Temporal covariance threshold detection at 1.1.
    [ALLEEG EEG] = eeg_store(ALLEEG, EEG, CURRENTSET);
    reject = find(EEG.reject.gcompreject==1); %Find the flagged ICs to reject
    EEG.ocularICs = reject; %Add the ICs indexes to the dataset, to be used for removing the ocular ICs.

    %Save the ICA pictures
    cd(rawpathFIG);
    name = regexprep(myDataSet(master_iterator).name,'_training30_icasegments.set',''); %Remove .set. 
    savename1 = strcat(name,'_ica.fig'); %Generate the new file name for saving.
    saveas(figure(1),savename1); %To save the current figure as FIG-file, i.e., Matlab figure.
    close gcf 

    %Save the processed data
    savename2 = strcat(name,'_detectedICs.set'); %Create the savename.
    EEG = pop_saveset( EEG, 'filename',savename2,'filepath',savePATH); %Save in the folder ready for ICA-decomposition.
    %[ALLEEG EEG] = eeg_store(ALLEEG, EEG, CURRENTSET);
    STUDY = []; CURRENTSTUDY = 0; ALLEEG = []; EEG=[]; CURRENTSET=[]; %Clear eeglab memory.
end

%It is recommended to check the flagged ICs to evaluate whether you agree or not with their removal.
%To check this, first, load each dataset from each single participant from the folder INFOMAX_ICA (i.e., Load existing dataset --> C:\Users\Federica\Documents\DegnoLobergLiversedge\EEGdata\INFOMAX_ICA\ppt20_training30_icasegments.set,
%then use "Eyetracker --> Reject ICA components with eye tracker" tool to check, and select as "Saccade event": R_saccade, and as "Fixation event": R_fixation.

%If you agree with the detected components to be rejected, leave everything as it is. If you don't agree, and you would like to remove additional components, then
%add the index number (associated with that component to reject) in the
%EEG.ocularICs to keep a record of what you have rejected.

Step 10: ICA转换和眼动伪迹IC成分标记。将优化ICA训练EEG数据中的IC成分数量和检测到的眼动伪迹等IC成分传递到原始的连续EEG数据中并进行标记。

%% 10. ICA TRANSFER
%Here we transfer the IC weights and the detected ocular ICs from the optimised ICA training EEG data to the original continuous EEG data.
%This part of code is based on Dimigen, O. (2020). Optimizing the ICA-based removal of ocular EEG artifacts during from viewing experiments. NeuroImage.

clear
eeglab

rawPATH1 = 'C:\Users\Federica\Documents\DegnoLobergLiversedge\EEGdata\IC_DETECTION\'; %Folder with data to process.
cd(rawPATH1); %Go to folder.
myDataSet = dir('*_detectedICs.set'); %Detect all files with this extension in the folder.

rawPATH2 = 'C:\Users\Federica\Documents\DegnoLobergLiversedge\EEGdata\BAD_ET1\'; %Folder with data to process.
cd(rawPATH2); %Go to folder.
mySyncSet = dir('*_badET1.set'); %Detect all files with this extension in the folder.

%Before running further, make sure that the order of the subjects in
%myDataSet and mySyncSet is  exactly the same.

savePATH = 'C:\Users\Federica\Documents\DegnoLobergLiversedge\EEGdata\IC_TRANSFER\'; %Folder where to save the processed data.

for master_iterator = 1:length(myDataSet) %Start the loop.

    %Load each dataset
    EEG = pop_loadset('filename',myDataSet(master_iterator).name,'filepath',rawPATH1); %Load each ICA training dataset with the detected ICs and weights.
    [ALLEEG, EEG, CURRENTSET] = eeg_store( ALLEEG, EEG, 0 ); 

    %Select channels
    MyCHANNELS = [find(strcmp({EEG.chanlocs.type},'EEG')) find(strcmp({EEG.chanlocs.type},'EOG'))]; %Indices of all EEG and EOG channels.

    %Store ICA information and flagged ocular ICs of training data
    %The true weight ICA matrix is equal to EEG.icaweights*EEG.icasphere.
    Sphere=EEG.icasphere;
    Weights=EEG.icaweights;
    ocularICs=EEG.ocularICs;

    %%Clear eeglab memory from the training data
    STUDY = []; CURRENTSTUDY = 0; ALLEEG = []; EEG=[]; CURRENTSET=[]; 

    %Load original EEG data
    EEG = pop_loadset('filename',mySyncSet(master_iterator).name,'filepath',rawPATH2); %Load the continuous original EEG data.
    [ALLEEG, EEG, CURRENTSET] = eeg_store( ALLEEG, EEG, 0 ); 

    %Remove any existing ICA solutions from the original dataset
    EEG.icaact      = [];
    EEG.icasphere   = [];
    EEG.icaweights  = [];
    EEG.icachansind = [];
    EEG.icawinv     = [];

    %Transfer unmixing weights (ICA information) from training data to original data 
    fprintf('\nTransfering ICA weights from training data to original data...')
    EEG.icasphere = Sphere;
    EEG.icaweights = Weights;  
    EEG.icachansind = MyCHANNELS; 
    EEG = eeg_checkset(EEG); %Let EEGLAB re-compute EEG.icaact & EEG.icawinv

    %Transfer detected ocular ICs indexes
    EEG.ocularICs = ocularICs;

    %Save the processed data
    name = regexprep(mySyncSet(master_iterator).name,'_badET1.set',''); %Remove '_badET1.set'.
    savename = strcat(name,'_ICweights.set'); %Create the savename.
    EEG = pop_saveset( EEG, 'filename',savename,'filepath',savePATH); %Save.
    [ALLEEG EEG] = eeg_store(ALLEEG, EEG, CURRENTSET);
    STUDY = []; CURRENTSTUDY = 0; ALLEEG = []; EEG=[]; CURRENTSET=[]; %clear eeglab memory from datafiles
end

Step 11: 在连续的EEG原始数据中对标记为眼动伪迹等IC成分进行移除。

%% 11. PRUNING OF THE OCULAR ICs
%Here we remove the ocular ICs.

clear
eeglab

rawPATH = 'C:\Users\Federica\Documents\DegnoLobergLiversedge\EEGdata\IC_TRANSFER\'; %Folder with data to process.
cd(rawPATH); %Go to folder.
myDataSet = dir('*_ICweights.set'); %Detect all files with this extension in the folder.

savePATH = 'C:\Users\Federica\Documents\DegnoLobergLiversedge\EEGdata\PRUNED\'; %Folder where to save the processed data.

for master_iterator = 1:length(myDataSet) %Start the loop.

    %Load each dataset
    EEG = pop_loadset('filename',myDataSet(master_iterator).name,'filepath',rawPATH); %Load each dataset.
    [ALLEEG, EEG, CURRENTSET] = eeg_store( ALLEEG, EEG, 0 ); 

    %Remove ocular ICs
    EEG = pop_subcomp( EEG, EEG.ocularICs, 0); %Remove the components that had their index numbers stored in 'ocularICs'.

    %Save the processed data
    name = regexprep(myDataSet(master_iterator).name,'_ICweights.set',''); 
    savename = strcat(name,'_pruned.set');
    EEG = pop_saveset( EEG, 'filename',savename,'filepath',savePATH); %Save the file.
    [ALLEEG EEG] = eeg_store(ALLEEG, EEG, CURRENTSET);
    STUDY = []; CURRENTSTUDY = 0; ALLEEG = []; EEG=[]; CURRENTSET=[];%Clear eeglab memory from datafiles.
end

以上11步就是基于眼动数据在ICA中识别并移除EEG中的眼动伪迹的全部步骤，关于基于ICA伪迹拒绝的后续分析步骤，我们会在后续的推文中进行详细的讲解。希望以上内容能对你有一定的帮助。下期见！！！

参考资料

[1]   Dimigen, O. (2020). Optimizing the ICA-based removal of ocular artifacts from free viewing EEG. 116117. NeuroImage, https://doi.org/10.1016/j.neuroimage.2019.116117

[2] Ehinger BV, Dimigen O. 2019. Unfold: an integrated toolbox for overlap correction, non-linear modeling, and regression-based EEG analysis. PeerJ 7:e7838 https://doi.org/10.7717/peerj.7838

[3] Dimigen, O. & Ehinger, B.V. (2021). Regression-based analysis of combined EEG and eye-tracking data: Theory and Applications. Journal of Vision, 21, 3

[4] OPTICAT: https://github.com/olafdimigen/OPTICAT

[5] Federica Degno, Otto Loberg, Simon P. Liversedge; Co-Registration of Eye Movements and Fixation—Related Potentials in Natural Reading: Practical Issues of Experimental Design and Data Analysis. Collabra: Psychology 4 January 2021; 7 (1): 18032. doi: https://doi.org/10.1525/collabra.18032

[6] preprocessing_FRPdata_v2.m:https://osf.io/k9ujr/?view_only=d28689a07d7a43bc9f8060bdf2a506a6

不用于商业行为，转载请联系后台

若有侵权，请后台留言，管理员即时删侵！

更多阅读

投稿通道

如何让你的工作让更多人知晓和受益？

脑机接口社区就是这样一个连接学界、

企业界和爱好者的平台渠道。

社区鼓励高校实验室、企业或个人在我们平台上分享优质内容。

稿件要求

稿件系个人原创作品，若已在其他平台发表，请明确标注。

稿件一经录取，便提供稿费！

投稿通道

微信扫码，备注：投稿+姓名+单位

微信交流群，请扫码上方微信

(备注：姓名+单位+专业/领域/行业)

QQ交流群：913607986

你的每一次在看，我都很在意！