ENH: muscle artifact detection (#7407)

* muscle artifact detection * docstrings * docstrings * Alex suggestions Co-Authored-By: Alexandre Gramfort <[email protected]> * example * example * for eeg * skip_by_annotation * new func added * docstring * cross ref * cross ref2 * Alex suggestions Co-Authored-By: Alexandre Gramfort <[email protected]> * new parameters * docstring * API name change, zscore * test update * sqrt nchans * handle nans zscore * example note * Apply suggestions from code review drammock Co-Authored-By: Daniel McCloy <[email protected]> * revision changes * docstring * docstring * ch_type param * logger * docstrings * ch type * test no meeg data * test no meeg data * fix docstring Co-authored-by: Alexandre Gramfort <[email protected]> Co-authored-by: Daniel McCloy <[email protected]>
mne-tools · Apr 7, 2020 · e3cc4b1 · e3cc4b1
1 parent ab50b23
commit e3cc4b1
Show file tree

Hide file tree

Showing 8 changed files with 232 additions and 4 deletions.
diff --git a/doc/changes/latest.inc b/doc/changes/latest.inc
@@ -25,6 +25,8 @@ Changelog
 
 - Add functionality to interpolate bad NIRS channels by `Robert Luke`_
 
+- Add function :func:`mne.preprocessing.annotate_muscle_zscore` to annotate periods with muscle artifacts. by `Adonay Nunes`_
+
 Bug
 ~~~
 

diff --git a/doc/python_reference.rst b/doc/python_reference.rst
@@ -346,6 +346,7 @@ Projections:
    ICA
    Xdawn
    annotate_movement
+   annotate_muscle_zscore
    compute_average_dev_head_t
    compute_current_source_density
    compute_proj_ecg

diff --git a/doc/references.bib b/doc/references.bib
@@ -1048,6 +1048,16 @@ @article{MourtazaevEtAl1995
   year = {1995}
 }
 
+@article{Muthukumaraswamy2013,
+  author = {Muthukumaraswamy, Suresh},
+  doi = {10.3389/fnhum.2013.00138},
+  journal = {Frontiers in Human Neuroscience},
+  pages = {138},
+  title = {High-frequency brain activity and muscle artifacts in {{MEG}}/{{EEG}}: A review and recommendations},
+  volume = {7},
+  year = {2013}
+}
+
 @inproceedings{NdiayeEtAl2016,
   author = {Ndiaye, Eugene and Fercoq, Olivier and Gramfort, Alexandre and Salmon, Joseph},
   booktitle = {Advances in Neural Information Processing Systems 29},

diff --git a/examples/preprocessing/plot_muscle_detection.py b/examples/preprocessing/plot_muscle_detection.py
@@ -0,0 +1,80 @@
+"""
+===========================
+Annotate muscle artifacts
+===========================
+
+Muscle contractions produce high frequency activity that can mask brain signal
+of interest. Muscle artifacts can be produced when clenching the jaw,
+swallowing, or twitching a cranial muscle. Muscle artifacts are most
+noticeable in the range of 110-140 Hz.
+
+This example uses :func:`~mne.preprocessing.annotate_muscle_zscore` to annotate
+segments where muscle activity is likely present. This is done by band-pass
+filtering the data in the 110-140 Hz range. Then, the envelope is taken using
+the hilbert analytical signal to only consider the absolute amplitude and not
+the phase of the high frequency signal. The envelope is z-scored and summed
+across channels and divided by the square root of the number of channels.
+Because muscle artifacts last several hundred milliseconds, a low-pass filter
+is applied on the averaged z-scores at 4 Hz, to remove transient peaks.
+Segments above a set threshold are annotated as ``BAD_muscle``. In addition,
+the ``min_length_good`` parameter determines the cutoff for whether short
+spans of "good data" in between muscle artifacts are included in the
+surrounding "BAD" annotation.
+
+"""
+# Authors: Adonay Nunes <[email protected]>
+#          Luke Bloy <[email protected]>
+# License: BSD (3-clause)
+
+import os.path as op
+import matplotlib.pyplot as plt
+import numpy as np
+from mne.datasets.brainstorm import bst_auditory
+from mne.io import read_raw_ctf
+from mne.preprocessing import annotate_muscle_zscore
+
+
+# Load data
+data_path = bst_auditory.data_path()
+raw_fname = op.join(data_path, 'MEG', 'bst_auditory', 'S01_AEF_20131218_01.ds')
+
+raw = read_raw_ctf(raw_fname, preload=False)
+
+raw.crop(130, 160).load_data()  # just use a fraction of data for speed here
+raw.resample(300, npad="auto")
+
+###############################################################################
+# Notch filter the data:
+#
+# .. note::
+#     If line noise is present, you should perform notch-filtering *before*
+#     detecting muscle artifacts. See :ref:`tut-section-line-noise` for an
+#     example.
+
+raw.notch_filter([50, 100])
+
+###############################################################################
+
+# The threshold is data dependent, check the optimal threshold by plotting
+# ``scores_muscle``.
+threshold_muscle = 5  # z-score
+# Choose one channel type, if there are axial gradiometers and magnetometers,
+# select magnetometers as they are more sensitive to muscle activity.
+annot_muscle, scores_muscle = annotate_muscle_zscore(
+    raw, ch_type="mag", threshold=threshold_muscle, min_length_good=0.2,
+    filter_freq=[110, 140])
+
+###############################################################################
+# Plot muscle z-scores across recording
+# --------------------------------------------------------------------------
+
+fig, ax = plt.subplots()
+ax.plot(raw.times, scores_muscle)
+ax.axhline(y=threshold_muscle, color='r')
+ax.set(xlabel='time, (s)', ylabel='zscore', title='Muscle activity')
+###############################################################################
+# View the annotations
+# --------------------------------------------------------------------------
+order = np.arange(144, 164)
+raw.set_annotations(annot_muscle)
+raw.plot(start=5, duration=20, order=order)
diff --git a/mne/preprocessing/__init__.py b/mne/preprocessing/__init__.py
@@ -22,4 +22,5 @@
 from .xdawn import Xdawn
 from ._csd import compute_current_source_density
 from . import nirs
-from .artifact_detection import (annotate_movement, compute_average_dev_head_t)
+from .artifact_detection import (annotate_movement, compute_average_dev_head_t,
+                                 annotate_muscle_zscore)
diff --git a/mne/preprocessing/artifact_detection.py b/mne/preprocessing/artifact_detection.py
@@ -4,12 +4,117 @@
 
 
 import numpy as np
-
 from ..annotations import (Annotations, _annotations_starts_stops)
 from ..transforms import (quat_to_rot, _average_quats, _angle_between_quats,
                           apply_trans, _quat_to_affine)
+from ..filter import filter_data
 from .. import Transform
-from ..utils import (_mask_to_onsets_offsets, logger)
+from ..utils import (_mask_to_onsets_offsets, logger, verbose)
+
+
+@verbose
+def annotate_muscle_zscore(raw, threshold=4, ch_type=None, min_length_good=0.1,
+                           filter_freq=(110, 140), n_jobs=1, verbose=None):
+    """Create annotations for segments that likely contain muscle artifacts.
+
+    Detects data segments containing activity in the frequency range given by
+    ``filter_freq`` whose envelope magnitude exceeds the specified z-score
+    threshold, when summed across channels and divided by ``sqrt(n_channels)``.
+    False-positive transient peaks are prevented by low-pass filtering the
+    resulting z-score time series at 4 Hz. Only operates on a single channel
+    type, if ``ch_type`` is ``None`` it will select the first type in the list
+    ``mag``, ``grad``, ``eeg``.
+    See :footcite:`Muthukumaraswamy2013` for background on choosing
+    ``filter_freq`` and ``threshold``.
+
+    Parameters
+    ----------
+    raw : instance of Raw
+        Data to estimate segments with muscle artifacts.
+    threshold : float
+        The threshold in z-scores for marking segments as containing muscle
+        activity artifacts.
+    ch_type : 'mag' | 'grad' | 'eeg' | None
+        The type of sensors to use. If ``None`` it will take the first type in
+        ``mag``, ``grad``, ``eeg``.
+    min_length_good : float | None
+        The shortest allowed duration of "good data" (in seconds) between
+        adjacent annotations; shorter segments will be incorporated into the
+        surrounding annotations.``None`` is equivalent to ``0``.
+        Default is ``0.1``.
+    filter_freq : array-like, shape (2,)
+        The lower and upper frequencies of the band-pass filter.
+        Default is ``(110, 140)``.
+    %(n_jobs)s
+    %(verbose)s
+
+    Returns
+    -------
+    annot : mne.Annotations
+        Periods with muscle artifacts annotated as BAD_muscle.
+    scores_muscle : array
+        Z-score values averaged across channels for each sample.
+
+    References
+    ----------
+    .. footbibliography::
+    """
+    from scipy.stats import zscore
+    from scipy.ndimage.measurements import label
+
+    raw_copy = raw.copy()
+
+    if ch_type is None:
+        raw_ch_type = raw_copy.get_channel_types()
+        if 'mag' in raw_ch_type:
+            ch_type = 'mag'
+        elif 'grad' in raw_ch_type:
+            ch_type = 'grad'
+        elif 'eeg' in raw_ch_type:
+            ch_type = 'eeg'
+        else:
+            raise ValueError('No M/EEG channel types found, please specify a'
+                             ' ch_type or provide M/EEG sensor data')
+        logger.info('Using %s sensors for muscle artifact detection'
+                    % (ch_type))
+
+    if ch_type in ('mag', 'grad'):
+        raw_copy.pick_types(meg=ch_type, ref_meg=False)
+    else:
+        ch_type = {'meg': False, ch_type: True}
+        raw_copy.pick_types(**ch_type)
+
+    raw_copy.filter(filter_freq[0], filter_freq[1], fir_design='firwin',
+                    pad="reflect_limited", n_jobs=n_jobs)
+    raw_copy.apply_hilbert(envelope=True, n_jobs=n_jobs)
+
+    data = raw_copy.get_data(reject_by_annotation="NaN")
+    nan_mask = ~np.isnan(data[0])
+    sfreq = raw_copy.info['sfreq']
+
+    art_scores = zscore(data[:, nan_mask], axis=1)
+    art_scores = art_scores.sum(axis=0) / np.sqrt(art_scores.shape[0])
+    art_scores = filter_data(art_scores, sfreq, None, 4)
+
+    scores_muscle = np.zeros(data.shape[1])
+    scores_muscle[nan_mask] = art_scores
+
+    art_mask = scores_muscle > threshold
+    # return muscle scores with NaNs
+    scores_muscle[~nan_mask] = np.nan
+
+    # remove artifact free periods shorter than min_length_good
+    min_length_good = 0 if min_length_good is None else min_length_good
+    min_samps = min_length_good * sfreq
+    comps, num_comps = label(art_mask == 0)
+    for com in range(1, num_comps + 1):
+        l_idx = np.nonzero(comps == com)[0]
+        if len(l_idx) < min_samps:
+            art_mask[l_idx] = True
+
+    annot = _annotations_from_mask(raw_copy.times, art_mask, 'BAD_muscle')
+
+    return annot, scores_muscle
 
 
 def annotate_movement(raw, pos, rotation_velocity_limit=None,

diff --git a/mne/preprocessing/tests/test_artifact_detection.py b/mne/preprocessing/tests/test_artifact_detection.py
@@ -5,11 +5,13 @@
 
 import os.path as op
 import numpy as np
+import pytest
 from numpy.testing import assert_allclose
 from mne.chpi import read_head_pos
 from mne.datasets import testing
 from mne.io import read_raw_fif
-from mne.preprocessing import (annotate_movement, compute_average_dev_head_t)
+from mne.preprocessing import (annotate_movement, compute_average_dev_head_t,
+                               annotate_muscle_zscore)
 
 data_path = testing.data_path(download=False)
 sss_path = op.join(data_path, 'SSS')
@@ -47,3 +49,28 @@ def test_movement_annotation_head_correction():
                               [0., 0., 0., 1.]])
 
     assert_allclose(dev_head_t_ori, dev_head_t['trans'], rtol=1e-5, atol=0)
+
+
+@testing.requires_testing_data
+def test_muscle_annotation():
+    """Test correct detection muscle artifacts."""
+    raw = read_raw_fif(raw_fname, allow_maxshield='yes').load_data()
+    raw.notch_filter([50, 110, 150])
+    # Check 2 muscle segments are detected
+    annot_muscle, scores = annotate_muscle_zscore(raw, ch_type='mag',
+                                                  threshold=10)
+    onset = annot_muscle.onset * raw.info['sfreq']
+    onset = onset.astype(int)
+    np.testing.assert_array_equal(scores[onset].astype(int), np.array([23,
+                                                                       10]))
+    assert(annot_muscle.duration.size == 2)
+
+
+@testing.requires_testing_data
+def test_muscle_annotation_without_meeg_data():
+    """Call annotate_muscle_zscore with data without meg or eeg."""
+    raw = read_raw_fif(raw_fname, allow_maxshield='yes')
+    raw.crop(0, .1).load_data()
+    raw.pick_types(meg=False, stim=True)
+    with pytest.raises(ValueError, match="No M/EEG channel types found"):
+        annot_muscle, scores = annotate_muscle_zscore(raw, threshold=10)
diff --git a/tutorials/preprocessing/plot_30_filtering_resampling.py b/tutorials/preprocessing/plot_30_filtering_resampling.py
@@ -101,6 +101,8 @@
 mne.viz.plot_filter(filter_params, raw.info['sfreq'], flim=(0.01, 5))
 
 ###############################################################################
+# .. _tut-section-line-noise:
+#
 # Power line noise
 # ~~~~~~~~~~~~~~~~
 #
-Original file line number
+Diff line change
@@ Expand Up / @@ -25,6 +25,8 @@ Changelog @@
     - Add functionality to interpolate bad NIRS channels by `Robert Luke`_
+    - Add function :func:`mne.preprocessing.annotate_muscle_zscore` to annotate periods with muscle artifacts. by `Adonay Nunes`_
     Bug
     ~~~
@@ Expand Down @@