Use independent component analysis (ICA) to remove ECG artifacts

This script demonstrates how you can use ICA for cleaning the ECG artifacts from your MEG data. It consists of three steps:

1. decomposition of the MEG data
2. identifying the components that reflect heart artifacts
3. removing those components and backprojecting the data

You can run the code below on your own data. Alternatively, try with the example MEG dataset ArtifactMEG.ds (available from ftp://ftp.fcdonders.nl/pub/fieldtrip/tutorial/ArtifactMEG.zip). This dataset was acquired continuously with trials of 10 seconds. All figures in this example script are based on these data.

To load this dataset into matlab and preprocess with FieldTrip, use:

% ft_preprocessing of example dataset
cfg = [];
cfg.dataset            = 'ArtifactMEG.ds';
cfg.trialdef.eventtype = 'trial';
cfg = ft_definetrial(cfg);

cfg.channel            = 'MEG';
cfg.continuous         = 'yes';
data = ft_preprocessing(cfg);

This script demonstrates how you can use ICA for cleaning the ECG artifacts from your MEG data. It starts by first doing a decomposition of the MEG data in the data segments of interest (i.e. the real trials in your experiment). Subsequently goes back to the original raw datafile and it reads the data segments around the QRS peaks that can easily be detected in the ECG channel. It uses the decomposition from the original data to estimate the timecourse of the components around the ECG artifacts. By looking at the component timecourses (averaged), the coherence between the components and the ECG channel, and the spatial topographies, it is possible to determine which components are responsible for the ECG artifact in the MEG channels. Those components can then be removed from the original data.

% read the already preprocessed MEG data from a Matlab file
load datfile.mat data

% keep the original data prior to the downsampling, we will use it at the end
data_orig = data;

% resample to speed up the ICA decomposition, especially usefull for 1200Hz MEG data
cfg            = [];
cfg.resamplefs = 300;
cfg.detrend    = 'no';
data           = ft_resampledata(cfg, data);

% perform the independent component analysis, i.e. decompose the downsampled original data
% this results in a mixing matrix and in the component timecourses
cfg            = [];
comp           = ft_componentanalysis(cfg, data);

% go back to the raw data on disk and detect the peaks in the ECG channel, i.e. the QRS-complex
cfg                      = [];
cfg.trl                  = data.cfg.previous.trl;
cfg.dataset              = data.cfg.previous.dataset;
cfg.datatype             = 'continuous';
cfg.artfctdef.ecg.pretim = 0.25;
cfg.artfctdef.ecg.psttim = 0.50-1/1200;
[cfg, artifact]          = ft_artifact_ecg(cfg);

% preproces the data around the QRS-complex, i.e. read the segments of raw data containing the ECG artifact
cfg           = [];
cfg.dataset   = data.cfg.previous.dataset;
cfg.datatype  = 'continuous';
cfg.padding   = 10;
cfg.dftfilter = 'yes';
cfg.blc       = 'yes';
cfg.trl       = [artifact zeros(size(artifact,1),1)];
cfg.channel   = {'MEG' 'ECG'};
data_ecg      = ft_preprocessing(cfg);

% resample to speed up the decomposition and frequency analysis, especially usefull for 1200Hz MEG data
cfg            = [];
cfg.resamplefs = 300;
cfg.detrend    = 'no';
data_ecg       = ft_resampledata(cfg, data_ecg);

% decompose the ECG-locked datasegments into components, using the previously found (un)mixing matrix
cfg           = [];
cfg.topo      = comp.topo;
cfg.topolabel = comp.topolabel;
comp_ecg      = ft_componentanalysis(cfg, data_ecg);

% average the components timelocked to the QRS-complex
cfg      = [];
timelock = ft_timelockanalysis(cfg, comp_ecg);

% look at the timelocked/averaged components and compare them with the ECG
figure
subplot(2,1,1); plot(timelock.time, timelock.avg(1,:))
subplot(2,1,2); plot(timelock.time, timelock.avg(2:end,:))
figure
subplot(2,1,1); plot(timelock.time, timelock.avg(1,:))
subplot(2,1,2); imagesc(timelock.avg(2:end,:));

% compute a frequency decomposition of all components and the ECG
cfg            = [];
cfg.method     = 'mtmfft';
cfg.output     = 'fourier';
cfg.foilim     = [0 100];
cfg.taper      = 'hanning';
cfg.pad        = 'maxperlen';
freq           = ft_freqanalysis(cfg, comp_ecg);
% compute coherence between all components and the ECG
cfg            = [];
cfg.channelcmb = {'all' 'ECG'};
cfg.jackknife  = 'no';
fdcomp         = ft_freqdescriptives(cfg, freq);

% look at the coherence spectrum between all components and the ECG
figure;
subplot(2,1,1); plot(fdcomp.freq, abs(fdcomp.cohspctrm));
subplot(2,1,2); imagesc(abs(fdcomp.cohspctrm));

% look at the spatial topography of the components
cfg = [];
cfg.layout = 'CTF151.lay';
lay = ft_prepare_layout(cfg);
figure;
subplot(2,3,1);topoplot([], lay.pos(1:151,1), lay.pos(1:151,2), comp.topo(:,1));title(num2str(1));
subplot(2,3,2);topoplot([], lay.pos(1:151,1), lay.pos(1:151,2), comp.topo(:,2));title(num2str(2));
subplot(2,3,3);topoplot([], lay.pos(1:151,1), lay.pos(1:151,2), comp.topo(:,3));title(num2str(3));
subplot(2,3,4);topoplot([], lay.pos(1:151,1), lay.pos(1:151,2), comp.topo(:,4));title(num2str(4));
subplot(2,3,5);topoplot([], lay.pos(1:151,1), lay.pos(1:151,2), comp.topo(:,5));title(num2str(5));
subplot(2,3,6);topoplot([], lay.pos(1:151,1), lay.pos(1:151,2), comp.topo(:,6));title(num2str(6));
% ... you may want to plot more topographies ...

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% based on the figures, you now should select the compontents that explain the ECG artifact
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% decompose the original data as it was prior to downsampling to 300Hz
cfg = [];
cfg.topo      = comp.topo;
cfg.topolabel = comp.topolabel;
comp_orig     = ft_componentanalysis(cfg, data_orig);

% the original data can now be reconstructed, excluding those components
cfg = [];
cfg.component = [2 3 6];
data_clean = ft_rejectcomponent(cfg, comp_orig);