18. Suppression of artifacts in EEG collected inside MRI
Electroencephalography (EEG) concurrently collected with functional magnetic resonance imaging (fMRI) is heavily distorted by the repetitive gradient coil switching during the fMRI acquisition. EEG signals are also degraded by the ballistocardiogram (BCG), the induced electric potentials by the head motion originating from heartbeats. Here are methods to process EEG collected inside MRI.
Our related works on this topic include:
We explained part codes, which use Average Artifact Subtract (AAS) and Optimal Basis Sets (OBS) or AAS and Dynamic Modeling of Heartbeats to suppress MRI gradient coil and pulse artifacts, receptively.
1. Gradient artifact (GA) suppression
Onsets of each MRI volume acquisition (represented by trigger_token) are identified. We code each onset by a number "1000". Then, ECG and EEG data are respectively processed by [Average Artifact Suppression]. The ECG signal after AAS is temporally smoothed by Savitzky-Golay filtering to smooth high-frequency oscillations while maintaining the shape. This process is meant to improve the identification of cardiac peaks in the following pulse artifact suppression.
gradient_trigger{f_idx}=zeros(size(eeg{f_idx},2),1);
iidx=find(trigger{f_idx}.event==trigger_token);
gradient_trigger{f_idx}(trigger{f_idx}.time(iidx))=1e3;
ecg{f_idx}=eeg_ga(ecg{f_idx},gradient_trigger{f_idx},TR,fs(f_idx),'flag_display',0,'flag_ma_aas',1,'flag_aas_svd',0,'flag_anchor_bnd',0,'n_ma_aas',7); %AAS on ECG
ecg{f_idx}=sgolayfilt(ecg{f_idx},6,301); %smoothing out residual GA in ECG
eeg{f_idx}=eeg_ga(eeg{f_idx},gradient_trigger{f_idx},TR,fs(f_idx),'flag_display',0,'flag_ma_aas',1,'flag_aas_svd',0,'flag_anchor_bnd',0,'n_ma_aas',7); %AAS on EEG
2. Pulse artifact (PA) suppression
Here we present two versions of PA suppression: [Optimal Basis Sets (OBS)](https://pubmed.ncbi.nlm.nih.gov/16150610/) and Dynamic Modeling of Heartbeats (DMH).We first use the Pan-Tompkins method to identify QRS complexes in the ECG. Then, OBS is used to estimate and suppress BCG artifacts across EEG data aligned at each QRS complex with 20% and 80% of the minimum of inter-QRS interval, respectively.
[qrs_amp_raw,qrs{f_idx}]=pan_tompkin(ecg{f_idx},fs(f_idx),0);
%%plot ECG and QRS detection
%tt=[1:length(ecg{f_idx})]./fs(f_idx);
%figure; plot(tt,ecg{f_idx}); hold on;
%line(repmat(tt(qrs{f_idx}),[2 1]),repmat([min(ecg{f_idx}-mean(ecg{f_idx}))/2; max(ecg{f_idx}-mean(ecg{f_idx}))/2],size(qrs{f_idx})),'LineWidth',2.5,'LineStyle','-.','Color','r');
if(flag_bcg(f_idx))
fprintf('\tBCG...\n');
iqi=diff(qrs{f_idx})./fs(f_idx);
iqi=sort(iqi);
iqi_min=iqi(round(length(iqi).*0.1));
tpre=0.2.*iqi_min;
tpost=0.8.*iqi_min;
[eeg{f_idx}, qrs_i_raw, bcg_all, ecg_all, bad_trials]=eeg_bcg(eeg{f_idx},ecg{f_idx},fs(f_idx),'flag_display',0,'bcg_nsvd',2,'n_ma_bcg',21,'flag_anchor_ends',0,'flag_post_ssp',0,'flag_badrejection',0,'BCG_tPre',tpre,'BCG_tPost',tpost); %BCG correction on EEG
end;
3. Examples of artifact suppression
Here are examples of EEG artifacts caused by MRI. Data were from subject s007 in the emotion experiment, where subjects were watching movie clips of various genres during concurrent EEG-MRI collection. EEG were collected by a 31-ch cap at 5,000 Hz. MRI were collected in the 3T scanner with a 64-ch head coil array. We used a temporally sparse MRI (2-s TR with each MRI volume taking 0.1 s) to minimize the gradient artifacts and their error propagation. Data were 30-s segements
- Before gradient artifact suppression by Average Artifact Subtraction, prominent intermittent (once every 2 s) gradient artifacts were observed.
- After gradient artifact suppression by Average Artifact Subtraction and low-pass filtering (100-Hz), gradient artifacts were much reduced. Now pulse artifacts around 1-Hz became more clear.
- After gradient artifact suppression by Average Artifact Subtraction, low-pass filtering (100-Hz), and pulse artifact suppression by Dynamic Modeling of Heartbeats, EEG became more like what it should be outside the MRI.
Detailed steps of processing EEG data collected inside an MRI scanner are explained below. The complete script is [here]. This is based on the processing of EEG-MRI data using simultaneous multi-slice inverse imaging (sparse fMRI) using music as the stimuli in 3T MRI with 32-channel EEG on October 6, 2023. Example scripts and data are within Sunnybrook server at /space_lin1/fhlin/eegmri_music/s008/eeg_analysis and /space_lin1/fhlin/eegmri_music/s008/eeg_raw, respectively.
1. Set input data files
It reads one set of EEG data (run 0001 for sub-008). One set of file including .eeg, .vhdr, and .vmrk three files.
headerFile = {
'../eeg_raw/sub-008_eeg-0001.vhdr';
};
markerFile={
'../eeg_raw/sub-008_eeg-0001.vmrk';
};
2. Set parameters
It reads 32 channels. Names of selected channels are noted in select_channel. We also indicate that EEG channels are #1 to #31 and ECG channel is #32.
n_chan=32; %32-channel EEG;
select_channel={ 'Fp1' 'Fp2' 'F3' 'F4' 'C3' 'C4' 'P3' 'P4' 'O1' 'O2' 'F7' 'F8' 'T7' 'T8' 'P7' 'P8' 'Fz' 'Cz' 'Pz' 'Oz' 'FC1' 'FC2' 'CP1' 'CP2' 'FC5' 'FC6' 'CP5' 'CP6' 'TP9' 'TP10' 'POz' 'ECG'};
eeg_channel=[1:31];
ecg_channel=[32];
Then we specified that MRI TR was 1.5 s. This parameter is important for gradient artifact suppression.
TR=[1.5];
The followings are markers in EEG data. During MRI acquisition, the MRI system sent a trigger named "R128" at the beginning of each TR. We used code 1000 to indicate this event. At the same time, there is "SYNC" trigger, which is numerically coded by 100 in our analysis.
%these two tokens are required (but values are arbitrary) for data
%collected inside MRI
trigger_token=1e3;
trigger_token_str='R128';
sync_token=1e2;
Finally, there are flags to indicate:
- re-referencing: subtract the average signals from each channel (flag_reref).
- high-pass filtering: default cut-off frequency at 0.1 Hz (flag_hp).
- low-pass filtering: default cut-off frequency at 70 Hz (flag_lp).
- gradient artifact suppression using Average Artifact Suppression (AAS) method (flag_aas).
- pulse artifact suppression using Optimal Basis Sets (OBS) method (flag_bcg).
- pulse artifact suppression using Dynamic Modeling of Heartbeats (DMH) method (flag_bcg_ccm).
A variable time_trim is used to indicate the beginning of time series to be trimmed off (in seconds; time_trim).
flag_reref= [1 1 1 1 1 1 1 1]; %referencing: remove the average time course across all electrodes from the time course at each individual electrode
flag_hp= [1 1 1 1 1 1 1 1]; %high-pass filtering at 0.1 Hz
flag_lp= [0 0 0 0 0 0 0 0]; %low-pass filtering at 70 Hz
flag_aas= [1 1 1 1 1 1 1 1]; %AAS for gradient artifact suppression
flag_bcg= [1 1 1 1 1 1 1 1]; %BCG artifact suppression
flag_bcg_ccm = [0 0 0 0 0 0 0 0]; %BCG artiface suppression using DMH
time_trim=[0]; %trim off the first few seconds...
3. Read data and basic info
Read data and related info (EEG channel name in `info`; sampling rate in `fs`). [dummy,fstem]=fileparts(headerFile{f_idx});
fprintf('reading [%s]...\n',fstem);
fstem=sprintf('%s%s',fstem,file_suffix);
% first get the continuous data as a matlab array
eeg{f_idx} = double(bva_loadeeg(headerFile{f_idx}));
% meta information such as samplingRate (fs), labels, etc
[fs(f_idx) label meta] = bva_readheader(headerFile{f_idx});
4. Select EEG channels from the raw data
Find the EEG channel specified in `select_channel`. found_channel={};
for s_idx=1:length(select_channel)
IndexC = strcmp(lower(label),lower(select_channel{s_idx})); %change all labels into lower case
Index = find(IndexC);
if(~isempty(Index))
fprintf('\tChannel [%s] found:: index=%03d \r',select_channel{s_idx},Index);
if(strcmp(lower(select_channel{s_idx}),'ecg'))
ecg_channel=Index;
else
eeg_channel(s_idx)=Index;
end;
found_channel{end+1}=select_channel{s_idx};
else
fprintf('\tChannel [%s] not found! \r',select_channel{s_idx});
end;
end;
5. Create EEG and ECG data variables
EEG and ECG data are set to eeg_orig and ecg_orig, respectively. They are also set to eeg and ecg variables. Note that these two variables are updated along the script.
ecg_orig=eeg{f_idx}(ecg_channel,:)';
eeg_orig=eeg{f_idx}(eeg_channel,:);
ecg{f_idx}=eeg{f_idx}(ecg_channel,:)';
eeg{f_idx}=eeg{f_idx}(eeg_channel,:);
You can now examine the data visually. Note that ECG has much larger amplitude than EEG.
etc_trace([eeg{1}; ecg{1}.'],'fs',fs(1),'ch_names',label);
6. Read EEG triggers/events
Events and triggers are now read from .vmrk files. .
%read maker file
fprintf('\treading triggers...\n');
trigger{f_idx}=etc_read_vmrk(markerFile{f_idx});
EEG_before_aas=eeg{f_idx};
7. Gradient artifact suppression
First find the triggers for each onset of MRI gradient switching (`gradient_trigger`). Then call `eeg_ga.m` to suppress gradient artifacts by AAS. if(flag_aas(f_idx))
fprintf('\tAAS...\n');
gradient_trigger{f_idx}=zeros(size(eeg{f_idx},2),1);
IndexC = strfind(trigger{f_idx}.event_str,trigger_token_str);
iidx = find(not(cellfun('isempty',IndexC)));
%iidx=find(trigger{f_idx}.event==trigger_token);
gradient_trigger{f_idx}(trigger{f_idx}.time(iidx))=trigger_token;
ecg{f_idx}=eeg_ga(ecg{f_idx},gradient_trigger{f_idx},TR,fs(f_idx),'flag_display',0,'flag_ma_aas',1,'flag_aas_svd',0,'flag_anchor_bnd',0,'n_ma_aas',7); %AAS on ECG
ecg{f_idx}=sgolayfilt(ecg{f_idx},6,301); %smoothing out residual GA in ECG
eeg{f_idx}=eeg_ga(eeg{f_idx},gradient_trigger{f_idx},TR,fs(f_idx),'flag_display',0,'flag_ma_aas',1,'flag_aas_svd',0,'flag_anchor_bnd',0,'n_ma_aas',7); %AAS on EEG
end;
Now we can examine the EEG after GA. Note that a larger scale in y-axis [-200 200] is now specified. The default y-axis range is automatically set to [-50 50]. For visualization, we only show EEG traces here.
etc_trace(eeg{1},'fs',fs(1),'ch_names',label(1:31),'ylim',[-200 200]);
To compare the EEG before AAS, we can load the variable EEG_before_aas in the viewer (by clicking "l" button at the far left of the command window and select the "data" list). Pressing and keys can switch between variables.
This is before AAS.
This is after AAS.
8. High-pass filtering
Perform high-pass filtering if `flag_hp` is set. %high-pass filtering
if(flag_hp(f_idx))
fprintf('\tHP...\n');
nn=round(fs(f_idx).*1); %1 Hz
for ch_idx=1:size(eeg{f_idx},1)
fprintf('*');
lp_data=filtfilt(ones(1,nn)./nn,1,eeg{f_idx}(ch_idx,:));
eeg{f_idx}(ch_idx,:)=eeg{f_idx}(ch_idx,:)-lp_data;
end;
fprintf('~');
lp_data=filtfilt(ones(1,nn)./nn,1,ecg{f_idx}(:)');
ecg{f_idx}=ecg{f_idx}(:)'-lp_data;
fprintf('\n');
else
ecg{f_idx}=ecg{f_idx}(:)';
end;
9. Trim data
Trim the first few seconds of EEG/ECG data if `time_trim` is specified. EEG_before_bcg=eeg{f_idx};
%remove first few seconds (if needed....)
if(~isempty(time_trim))
time_trim_idx=round(time_trim*fs(f_idx));
fprintf('trimming [%1.1f] s data {(%d) samples}....\n',time_trim, time_trim_idx);
ecg{f_idx}=ecg{f_idx}(time_trim_idx+1:end);
eeg{f_idx}=eeg{f_idx}(:,time_trim_idx+1:end);
%trigger{f_idx}.time=trigger{f_idx}.time-time_trim_idx;
end;
10. Pulse artifact suppression
First use ECG to detect QRS complexes. This is very critical in the pulse artifact suppression! It might be good to visualize the ECG waveform and the detected QRS complexes before proceeding next steps.NOTE: The required pan_tompkin2.m is [here] .
[qrs_amp_raw,qrs{f_idx}]=pan_tompkin2(ecg{f_idx},fs(f_idx));
Then use the OBS method to suppress pulse artifacts. Here we first down sample the data (with low-pass filtering to avoid aliasing) to improve the computational efficiency.
NOTE In the low-pass filtering, we assumed here that the sampling rate is 5,000 Hz and we down-sampled data to 50 Hz. This 100-fold down sampling may be different in other settings.
iqi=diff(qrs{f_idx})./fs(f_idx);
iqi=sort(iqi);
iqi_min=iqi(round(length(iqi).*0.5));
tpre=0.2.*iqi_min;
tpost=0.8.*iqi_min;
%down-sampling
tmp=[];
for ch_idx=1:size(eeg{f_idx},1)
tmp0=filtfilt(ones(100,1)./100,1,eeg{f_idx}(ch_idx,:));
tmp(ch_idx,:)=tmp0(1:100:end);
end;
eeg_dec=tmp;
tmp=filtfilt(ones(100,1)./100,1,ecg{f_idx});
ecg_dec=tmp(1:100:end);
fs_dec=fs(f_idx)/100;
eeg_dec_bcg_pred=[];
if(~flag_bcg_ccm(f_idx))
[eeg_dec_bcg, qrs_i_raw, bcg_all, ecg_all, bad_trials]=eeg_bcg(eeg_dec,ecg_dec,fs_dec,'flag_display',0,'bcg_nsvd',3,'n_ma_bcg',21,'flag_anchor_ends',0,'flag_post_ssp',0,'flag_badrejection',0,'BCG_tPre',tpre,'BCG_tPost',tpost,'flag_pan_tompkin2',1); %BCG correction on EEG
else
e=8;
tau=1;
nn=10;
[eeg_dec_bcg, qrs_i_raw, eeg_dec_bcg_pred, eeg_manifold_neighbor, eeg_manifold_now,eeg_manifold_now_approx, ecg_manifold_neighbor, ecg_manifold_now,ecg_manifold_now_approx, ccm_D, ccm_IDX]=eeg_bcg_ccm2(eeg_dec,ecg_dec,fs_dec,'e',e,'tau',tau,'nn',nn,'flag_display',0,'n_ecg',30,'flag_pan_tompkin2',1);
%[eeg_dec_bcg, eeg_dec_bcg_pred]=eeg_bcg_ccm5(eeg_dec,ecg_dec,fs_dec,'flag_display',0,'flag_auto_hp',1,'flag_pan_tompkin2',1);
end;
%up-samping
for ch_idx=1:size(eeg_dec_bcg,1)
tmp=interp(eeg_dec_bcg(ch_idx,:),100);
%eeg{f_idx}(ch_idx,:)=sgolayfilt(tmp,4,301);
eeg{f_idx}(ch_idx,:)=tmp(1:size(eeg{f_idx},2));
if(~isempty(eeg_dec_bcg_pred))
tmp=interp(eeg_dec_bcg_pred(ch_idx,:),100);
eeg_bcg_pred{f_idx}(ch_idx,:)=tmp(1:size(eeg{f_idx},2));
else
eeg_bcg_pred{f_idx}=[];
end;
end;
EEG_after_bcg=eeg{f_idx};
To visualize EEG traces after suppressing BCG artifacts, use the following script. The data before suppressing BCG artifacts should be low-pass filtered (50 Hz cut-off) for fair comparison, because BCG suppression in the script above automaticallyd smooth and down/up sample data at between 5,000 Hz and 50 Hz.
for ch_idx=1:size(EEG_before_bcg,1)
tmp(ch_idx,:)=filtfilt(ones(100,1)./100,1,EEG_before_bcg(ch_idx,:));
end;
data_before_bcg=cat(1,tmp,ecg{1});
data_after_bcg=cat(1,EEG_after_bcg,ecg{1});
etc_trace(data_after_bcg,'fs',fs(1));
To compare the EEG before BCG suppression, we can load the variable data_before_bcg in the viewer (by clicking "l" button at the far left of the command window and select the "data" list). Pressing and keys can switch between variables.
This is before BCG suppression.
This is after BCG suppression.
11. Re-reference
Perform re-referencing (subtracting the average of all channels from each channel) if `flag_reref` is set. %re-referencing
if(flag_reref(f_idx))
fprintf('\treferencing...\n');
eeg_ref{f_idx}=mean(eeg{f_idx},1);
for ch_idx=1:size(eeg{f_idx},1)
eeg{f_idx}(ch_idx,:)=eeg{f_idx}(ch_idx,:)-eeg_ref{f_idx};
end;
else
eeg_ref{f_idx}=[];
end;
EEG_after_reref=eeg{f_idx};
12. Update ECG trigger
Add the ECG QRS peaks as new triggers.%assign a trigger to R-peaks of ECG
trigger_ecg{f_idx}.event=33.*ones(1,length(qrs{f_idx}(:)));
trigger_ecg{f_idx}.time=qrs{f_idx}(:)';
%if data has been trimmed, ECG trigger should be adjusted.
if(~isempty(time_trim))
trigger_ecg{f_idx}.time=trigger_ecg{f_idx}.time+round(fs(f_idx).*time_trim);
end;
%merge ECG triggers
trigger{end}=etc_trigger_append(trigger{end},trigger_ecg{f_idx});
%update trigger information by including ECG triggers
TRIGGER_ECG=trigger_ecg{end};
TRIGGER=trigger{end};
EEG=eeg{end}; %<--------
ECG=ecg{end};
13. Perform low-pass filtering
Perform high-pass filtering if `flag_lp` is set. If BCG has been suppressed, this may not be necessary, because low-pass filtering has been implemented there.%low-poass filtering EEG (cut-off = 70 Hz)
if(flag_lp(f_idx))
a=1;
b=ones(round(fs(f_idx)./70),1)./round(fs(f_idx)./70);
for idx=1:size(EEG,1)
EEG(idx,:)=filtfilt(b,a,EEG(idx,:));
end;
end;
14. Finalize procesing
- Compenset trimmed time series, if any, by zeros.
- Save results into an MAT file.
%compensate the EEC/ECG data with zeros for the truncated time, if any
EEG=cat(2,eeg_orig(:,1:fs(f_idx).*time_trim),EEG);
ECG=cat(2,ecg_orig(1:fs(f_idx).*time_trim)',ECG);
HeaderFile=headerFile{f_idx};
MarkerFile=markerFile{f_idx};
sfreq=fs(f_idx);
save(sprintf('%s_nolp.mat',fstem),'label','sfreq','HeaderFile','MarkerFile','EEG','ECG','time_trim','TRIGGER_ECG','TRIGGER','-v7.3','EEG_before_aas','EEG_before_bcg','EEG_bcg_pred','EEG_after_reref','EEG_after_bcg');
This is the final EEG traces after all processing steps.
There are still noticeable artifacts, mostly related to heartbeats. So there is always room to improve!!