ENH: Automatically compute threshold for CTPS cardiac artifacts detection

[yh-luo]

Describe the problem

Current default threshold is 0.25 for mne.preprocessing.ICA.find_bads_ecg with method='ctps'. However, the value seems arbitrary. Would it make sense to implement automatic computation for CTPS threshold?

0.25 is actually the threshold for data with 400 Hz sampling rate. I found that 0.25 is sometimes not sensitive enough for clean data with sampling rate higher than 400 Hz.

For example, the sampling rate of sample dataset is 600 Hz. After filtering, mne.preprocessing.ICA.find_bads_ecg does not detect any ECG related IC using default (0.25) threshold. Using 0.21 (computed from the formula), it is able to detect one ECG related IC.

The ECG Evoked

The ECG scores of ICs

The detected IC (ICA001) with threshold=0.21

Describe your solution

In the reference, the threshold was set to Pk >= 10^-20. Pk can be computed from formula (13) with λ and λ will change according to the number of data points. Thus, the threshold will also change according to the sampling rate.

The threshold of the Kuiper index can be determined easily:

import numpy as np

def _get_ctps_threshold(N=1000, Pk_threshold=1e-20):
    vs = [x / 100 for x in range(1, 100)]
    Pks = list()
    C = np.sqrt(N) + 0.155 + 0.24 / np.sqrt(N)
    # because when k gets large, only k=1 matters for the summation
    # k*v*C thus becomes v*C
    for v in vs:
        Pk = 2 * (4 * (v * C)**2 - 1) * (np.exp(-2 * (v * C)**2))
        Pks.append(abs(Pk - Pk_threshold))
    return vs[Pks.index(min(Pks))]

In [2]: _get_ctps_threshold(400)                                                                                     
Out[2]: 0.25

In [3]: _get_ctps_threshold(600)                                                                                     
Out[3]: 0.21

In [4]: _get_ctps_threshold(1000)                                                                                    
Out[4]: 0.16

If that's reasonable, I can make a PR for this implementation.

Additional context

The test script:

import os.path as op

from IPython import get_ipython

import mne
from autoreject import get_rejection_threshold

mne.set_log_level('INFO')
get_ipython().run_line_magic('matplotlib', 'qt')
get_ipython().run_line_magic('gui', 'qt')

h_freq = None
l_freq = 40

# Calibration files for Maxwell filter
sample_dir = mne.datasets.sample.data_path()
ctc = op.join(sample_dir, 'SSS', 'ct_sparse_mgh.fif')
cal = op.join(sample_dir, 'SSS', 'sss_cal_mgh.dat')

raw_fname = op.join(sample_dir, 'MEG', 'sample', 'sample_audvis_raw.fif')
sss_fname = op.join(sample_dir, 'MEG', 'sample', 'sample_audvis_raw_sss.fif')

raw = mne.io.read_raw_fif(raw_fname, preload=True)
# remove bad channels
raw.info['bads'].extend(['MEG 1032', 'MEG 2313'])
raw_sss = mne.preprocessing.maxwell_filter(raw,
                                           calibration=cal,
                                           cross_talk=ctc)
# filter
raw_sss.filter(None, 40)

# ICA
n_components = 0.999
ica = mne.preprocessing.ICA(n_components=n_components, random_state=99)

picks_eog = mne.pick_types(raw.info,
                           meg=True,
                           eeg=False,
                           eog=False,
                           stim=False,
                           exclude='bads')
# high-pass for ICA
raw_sss.filter(1, None, picks=picks_eog, fir_window='hann')

# Only remove ECG artifacts for now
picks = mne.pick_types(raw.info,
                       meg=True,
                       eeg=False,
                       eog=False,
                       stim=False,
                       exclude='bads')
# use autoreject to find the rejection threshold to get better ICA results
tstep = 1.0
events = mne.make_fixed_length_events(raw_sss, duration=tstep)
even_epochs = mne.Epochs(raw_sss,
                         events,
                         baseline=(0, 0),
                         tmin=0.0,
                         tmax=tstep,
                         reject_by_annotation=True)
reject = get_rejection_threshold(even_epochs,
                                 ch_types=['mag', 'grad'])
ica.fit(raw_sss, picks=picks, reject=reject, tstep=tstep)

# ECG
ecg_epochs = mne.preprocessing.create_ecg_epochs(raw_sss)
ecg_epochs.apply_baseline((None, None))

# threshold = 0.25
# found 0 components
ecg_inds, ecg_scores = ica.find_bads_ecg(ecg_epochs,
                                         method='ctps',
                                         threshold=0.25)
ecg_evoked = ecg_epochs.average()
ecg_evoked.plot()
# barplot of ICA component "ECG match" scores
ica.plot_scores(ecg_scores)
print(ecg_inds)

# threshold = automatically for 600Hz data
# Using threshold: 0.21 for CTPS ECG detection
# found 1 components [1]
ecg_inds, ecg_scores = ica.find_bads_ecg(ecg_epochs, method='ctps')
print(ecg_inds)

# plot diagnostics
ica.plot_properties(raw_sss, picks=ecg_inds)

# plot ICs applied to raw data, with ECG matches highlighted
ica.plot_sources(raw_sss)

# plot ICs applied to the averaged ECG epochs, with ECG matches highlighted
ica.plot_sources(ecg_evoked)

ica.plot_overlay(raw_sss, exclude=ecg_inds, picks='mag')

[agramfort]

can you send us a PR?

[yh-luo]

sure, I can make a PR for this.

[yh-luo]

Closed by #7819