ENH: muscle artifact detection

doc/changes/latest.inc

@@ -145,6 +145,8 @@ Changelog
 
 - Improve repr of :class:`mne.Info` (empty items are not shown anymore) by `Clemens Brunner`_
 
+- Add function :func:`mne.preprocessing.annotate_muscle_zscore` to annotate periods with muscle artifacts. by `Adonay Nunes`_
+
 Bug
 ~~~
 

doc/references.bib

@@ -1037,6 +1037,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},

examples/preprocessing/plot_muscle_detection.py

@@ -5,21 +5,21 @@
 
 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 head muscle. Muscle artifacts are most noticeable
-in the range of 110-140 Hz.
+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, we then low-pass
-filter 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.
+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]>
@@ -29,7 +29,6 @@
 import os.path as op
 import matplotlib.pyplot as plt
 import numpy as np
-from mne import pick_types
 from mne.datasets.brainstorm import bst_auditory
 from mne.io import read_raw_ctf
 from mne.preprocessing import annotate_muscle_zscore
@@ -45,14 +44,7 @@
 raw.resample(300, npad="auto")
 
 ###############################################################################
-# The inputted raw data has to be from a single channel type as there might be
-# different channel type magnitudes.
-# If the MEG has axial gradiometers and magnetometers, select magnetometers as
-# they are more sensitive to muscle activity.
-picks = pick_types(raw.info, meg=True, ref_meg=False)
-
-###############################################################################
-# Notch filter the data.
+# Notch filter the data:
 #
 # .. note::
 #     If line noise is present, you should perform notch-filtering *before*
@@ -61,10 +53,14 @@
 
 raw.notch_filter([50, 100])
 
+###############################################################################
+
 # The threshold is data dependent, check the optimal threshold by plotting
 # ``scores_muscle``.
 threshold_muscle = 5  # z-score
-out = annotate_muscle_zscore(raw, picks=picks, threshold=threshold_muscle,
+# Choose one channel type, if there are axial gradiometers and magnetometers,
+# select magnetometers as they are more sensitive to muscle activity.
+out = annotate_muscle_zscore(raw, ch_type="mag", threshold=threshold_muscle,
                              min_length_good=0.2, filter_freq=[110, 140])
 annot_muscle, scores_muscle = out
 

mne/preprocessing/artifact_detection.py

@@ -13,21 +13,19 @@
 
 
 @verbose
-def annotate_muscle_zscore(raw, threshold=4, picks=None, min_length_good=0.1,
+def annotate_muscle_zscore(raw, threshold=4, ch_type=None, min_length_good=0.1,
                            filter_freq=(110, 140), n_jobs=1, verbose=None):
-    """Detect segments that likely contain muscle artifacts.
+    """Create annotations for segments that likely contain muscle artifacts.
 
     Detects data segments containing activity in the frequency range given by
-    ``filter_freq``, that exceeds the specified z-score threshold.
-
-    Raw data is band-pass filtered between ``filter_freq`` especified
-    frequencies (default is 110 - 140 Hz), the signal envelope computed,
-    z-scored across samples, channel summation and division by the square root
-    of the channel number, and low-pass filtered to better capture beginning
-    and end of muscle activity and false positive transient peaks.
-
-    .. note::
-        Use a single channel type.
+    `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
     ----------
@@ -36,7 +34,9 @@ def annotate_muscle_zscore(raw, threshold=4, picks=None, min_length_good=0.1,
     threshold : float
         The threshold in z-scores for marking segments as containing muscle
         activity artifacts.
-    %(picks_all)s
+    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
@@ -54,31 +54,39 @@ def annotate_muscle_zscore(raw, threshold=4, picks=None, min_length_good=0.1,
         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()
-    raw_copy.pick(picks)
-    ch_type = raw_copy.get_channel_types()
 
-    # Remove ref chans if MEG data just in case
-    meg = [True for e in ['mag', 'grad'] if (e in ch_type)]
-    if bool(meg) is True:
-        raw_copy.pick_types(ref_meg=False)
-        ch_type = raw_copy.get_channel_types()
-
-    # Only one type of channel, otherwise z-score will be biased
-    assert(len(np.unique(ch_type)) == 1), \
-        'Different channel types, pick one type'
+    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'
+        logger.info(u'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)
+    elif ch_type == 'eeg':
+        raw_copy.pick_types(meg=False, eeg=True)
 
     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)
-    sfreq = raw_copy.info['sfreq']
 
     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])

mne/preprocessing/tests/test_artifact_detection.py

@@ -54,10 +54,10 @@ def test_movement_annotation_head_correction():
 def test_muscle_annotation():
     """Test correct detection muscle artifacts."""
     raw = read_raw_fif(raw_fname, allow_maxshield='yes').load_data()
-    raw.pick_types(meg='mag', ref_meg=False)
     raw.notch_filter([50, 110, 150])
     # Check 2 muscle segments are detected
-    annot_muscle, scores = annotate_muscle_zscore(raw, threshold=10)
+    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,