From 7c5e680d15177b872af3190e148bc435dc4c652b Mon Sep 17 00:00:00 2001 From: AdoNunes Date: Sun, 8 Mar 2020 19:24:28 -0700 Subject: [PATCH 01/31] muscle artifact detection --- doc/python_reference.rst | 1 + mne/preprocessing/__init__.py | 3 +- mne/preprocessing/artifact_detection.py | 62 ++++++++++++++++++- .../tests/test_artifact_detection.py | 14 ++++- 4 files changed, 77 insertions(+), 3 deletions(-) diff --git a/doc/python_reference.rst b/doc/python_reference.rst index f74b6f22cdd..fd3d0bc2d64 100644 --- a/doc/python_reference.rst +++ b/doc/python_reference.rst @@ -346,6 +346,7 @@ Projections: ICA Xdawn annotate_movement + annotate_muscle compute_average_dev_head_t compute_current_source_density compute_proj_ecg diff --git a/mne/preprocessing/__init__.py b/mne/preprocessing/__init__.py index 5c439aa3250..3941ba9b153 100644 --- a/mne/preprocessing/__init__.py +++ b/mne/preprocessing/__init__.py @@ -24,4 +24,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) diff --git a/mne/preprocessing/artifact_detection.py b/mne/preprocessing/artifact_detection.py index 2aa12fab89a..c439b2f257e 100644 --- a/mne/preprocessing/artifact_detection.py +++ b/mne/preprocessing/artifact_detection.py @@ -4,14 +4,74 @@ import numpy as np - +from scipy.stats import zscore +from scipy.ndimage.measurements import label 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) +def annotate_muscle(raw, threshold=1.5, picks=None, min_length_good=.1): + """Detect segments with muscle artifacts. + + Detects segments periods that contains high frequency activity beyond the + specified threshold. Muscle artifacts are most notable in the range of 110- + 140Hz. + + Raw data is band pass filtered between 110 and 140 Hz, the signal envelope + computed, z-scored across samples, channel averaged and low-pass + filtered to smooth transient peaks. + + Parameters + ---------- + raw : instance of Raw + Data to compute head position. + threshold : float + The threshod for selecting segments with muscle activity artifacts. + picks: + Channels to use for artifact detection. + min_length_good : int | float | None, in seconds + The minimal good segment length between annotations, smaller segments will + be included in the movement annotation. + + Returns + ------- + annot : mne.Annotations + Periods with muscle artifacts. + scores_muscle : array shape (1 x n samples) + Z-score values averaged accros channels for each sample. + """ + raw_copy = raw.copy() + raw_copy.pick(picks) + raw_copy.pick_types(ref_meg=False) # Remove ref chans just in case + # Only one type of channel, otherwise z-score will be biased + assert(len(set(raw_copy.get_channel_types())) == 1), 'Different channel ' \ + 'types, pick one type' + + raw_copy.filter(110, 140, fir_design='firwin') + raw_copy.apply_hilbert(envelope=True) + sfreq = raw_copy.info['sfreq'] + + art_scores = zscore(raw_copy._data, axis=1) + scores_muscle = filter_data(art_scores.mean(axis=0), sfreq, None, 4) + art_mask = scores_muscle > threshold + + # remove artifact free periods shorter than min_length_good + idx_min = min_length_good * sfreq + comps, num_comps = label(art_mask == 0) + for l in range(1, num_comps + 1): + l_idx = np.nonzero(comps == l)[0] + if len(l_idx) < idx_min: + 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, translation_velocity_limit=None, mean_distance_limit=None): diff --git a/mne/preprocessing/tests/test_artifact_detection.py b/mne/preprocessing/tests/test_artifact_detection.py index 6fbfa627d57..b30eca6c2c1 100755 --- a/mne/preprocessing/tests/test_artifact_detection.py +++ b/mne/preprocessing/tests/test_artifact_detection.py @@ -9,7 +9,8 @@ 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) data_path = testing.data_path(download=False) sss_path = op.join(data_path, 'SSS') @@ -47,3 +48,14 @@ 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.pick_types(meg='grad', ref_meg=False) + raw.notch_filter([50, 110, 150]) + # Check 2 muscle segments are detected + annot_muscle, scores = annotate_muscle(raw, threshold=1) + assert(annot_muscle.duration.size == 2) From e42ed5cc54aa6ce174076f793c8a6ff897117058 Mon Sep 17 00:00:00 2001 From: AdoNunes Date: Sun, 8 Mar 2020 19:49:42 -0700 Subject: [PATCH 02/31] docstrings --- mne/preprocessing/artifact_detection.py | 8 ++++---- mne/preprocessing/tests/test_artifact_detection.py | 2 +- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/mne/preprocessing/artifact_detection.py b/mne/preprocessing/artifact_detection.py index c439b2f257e..5265c082c7b 100644 --- a/mne/preprocessing/artifact_detection.py +++ b/mne/preprocessing/artifact_detection.py @@ -33,15 +33,15 @@ def annotate_muscle(raw, threshold=1.5, picks=None, min_length_good=.1): The threshod for selecting segments with muscle activity artifacts. picks: Channels to use for artifact detection. - min_length_good : int | float | None, in seconds - The minimal good segment length between annotations, smaller segments will - be included in the movement annotation. + min_length_good : int | float | None + The minimal good segment length between annotations, smaller segments + will be included in the movement annotation. Returns ------- annot : mne.Annotations Periods with muscle artifacts. - scores_muscle : array shape (1 x n samples) + scores_muscle : array Z-score values averaged accros channels for each sample. """ raw_copy = raw.copy() diff --git a/mne/preprocessing/tests/test_artifact_detection.py b/mne/preprocessing/tests/test_artifact_detection.py index b30eca6c2c1..1612c811a43 100755 --- a/mne/preprocessing/tests/test_artifact_detection.py +++ b/mne/preprocessing/tests/test_artifact_detection.py @@ -52,7 +52,7 @@ def test_movement_annotation_head_correction(): @testing.requires_testing_data def test_muscle_annotation(): - """Test correct detection muscle artifacts""" + """Test correct detection muscle artifacts.""" raw = read_raw_fif(raw_fname, allow_maxshield='yes').load_data() raw.pick_types(meg='grad', ref_meg=False) raw.notch_filter([50, 110, 150]) From b705ae708ed4575bee7c312a043629f25d67703e Mon Sep 17 00:00:00 2001 From: AdoNunes Date: Sun, 8 Mar 2020 20:49:13 -0700 Subject: [PATCH 03/31] docstrings --- mne/preprocessing/artifact_detection.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/mne/preprocessing/artifact_detection.py b/mne/preprocessing/artifact_detection.py index 5265c082c7b..54db302747a 100644 --- a/mne/preprocessing/artifact_detection.py +++ b/mne/preprocessing/artifact_detection.py @@ -31,7 +31,7 @@ def annotate_muscle(raw, threshold=1.5, picks=None, min_length_good=.1): Data to compute head position. threshold : float The threshod for selecting segments with muscle activity artifacts. - picks: + picks : array Channels to use for artifact detection. min_length_good : int | float | None The minimal good segment length between annotations, smaller segments From 73351dbbf50f1d8bac69420d2e0e4a4fa266f346 Mon Sep 17 00:00:00 2001 From: Adonay Nunes Date: Mon, 9 Mar 2020 16:20:42 -0700 Subject: [PATCH 04/31] Alex suggestions Co-Authored-By: Alexandre Gramfort --- mne/preprocessing/artifact_detection.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/mne/preprocessing/artifact_detection.py b/mne/preprocessing/artifact_detection.py index 54db302747a..05c6cb2d124 100644 --- a/mne/preprocessing/artifact_detection.py +++ b/mne/preprocessing/artifact_detection.py @@ -21,7 +21,7 @@ def annotate_muscle(raw, threshold=1.5, picks=None, min_length_good=.1): specified threshold. Muscle artifacts are most notable in the range of 110- 140Hz. - Raw data is band pass filtered between 110 and 140 Hz, the signal envelope + Raw data is band-pass filtered between 110 and 140 Hz, the signal envelope computed, z-scored across samples, channel averaged and low-pass filtered to smooth transient peaks. @@ -30,7 +30,7 @@ def annotate_muscle(raw, threshold=1.5, picks=None, min_length_good=.1): raw : instance of Raw Data to compute head position. threshold : float - The threshod for selecting segments with muscle activity artifacts. + The threshold for selecting segments with muscle activity artifacts. picks : array Channels to use for artifact detection. min_length_good : int | float | None @@ -42,7 +42,7 @@ def annotate_muscle(raw, threshold=1.5, picks=None, min_length_good=.1): annot : mne.Annotations Periods with muscle artifacts. scores_muscle : array - Z-score values averaged accros channels for each sample. + Z-score values averaged across channels for each sample. """ raw_copy = raw.copy() raw_copy.pick(picks) From 24da9f697e5033cf2ca721a8688ad07d659d6a94 Mon Sep 17 00:00:00 2001 From: AdoNunes Date: Mon, 9 Mar 2020 16:48:17 -0700 Subject: [PATCH 05/31] example --- .../preprocessing/plot_muscle_detection.py | 82 +++++++++++++++++++ 1 file changed, 82 insertions(+) create mode 100644 examples/preprocessing/plot_muscle_detection.py diff --git a/examples/preprocessing/plot_muscle_detection.py b/examples/preprocessing/plot_muscle_detection.py new file mode 100644 index 00000000000..048c851aa80 --- /dev/null +++ b/examples/preprocessing/plot_muscle_detection.py @@ -0,0 +1,82 @@ +""" +=========================== +Annotate muscle artifacts +=========================== + +Muscle contractions produce high frequency activity that can contaminate the +brain signal of interest. Muscle artifacts can be produced when clinching the +jaw, swalloing, or twitching a head muscle. Muscle artifacts are most notable +in teh range of 110-140Hz. + +This example uses `annotate_muscle` to annotate segments where muscle activity +likely occurred. This is done by band-pass filtering the data in the 110-140 Hz +range. Then, the envelope is taken to account for peacks and troughs. The +envelope is z-scored and averaged across channels. To remove noisy tansient +peak, the z-scored average is low-pass filtered to 4 Hz. Segments above a set +threshold are annotated as BAD_motion. In addition, `min_length_good` allows to +discard god segments of data between bad segments that are to transient. + +Note: +The raw data needs to be notched to remove line activity, otherwise, changes of +AC power could bias the estimation. + +The inputted raw data has to be from a single channel type as there might be +different channel type magnitudes. + +""" +# Authors: Adonay Nunes +# Luke Bloy +# License: BSD (3-clause) + +import os.path as op +import matplotlib.pyplot as plt +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 +from mne.io import concatenate_raws + +# Load data + +data_path = bst_auditory.data_path() +data_path_MEG = op.join(data_path, 'MEG') +subject = 'bst_auditory' +subjects_dir = op.join(data_path, 'subjects') +trans_fname = op.join(data_path, 'MEG', 'bst_auditory', + 'bst_auditory-trans.fif') +raw_fname1 = op.join(data_path_MEG, 'bst_auditory', 'S01_AEF_20131218_01.ds') +raw_fname2 = op.join(data_path_MEG, 'bst_auditory', 'S01_AEF_20131218_02.ds') + +raw = read_raw_ctf(raw_fname1, preload=False) + +raw = concatenate_raws([raw, read_raw_ctf(raw_fname2, preload=False)]) +raw.crop(350, 410).load_data() +raw.resample(300, npad="auto") +raw.notch_filter([50, 110, 140]) + +picks = pick_types(raw.info, meg=True, ref_meg=False) + +# detect muscle artifacts +threshold_muscle = 1.5 # z-score +annotation_muscle, scores_muscle = annotate_muscle(raw, picks=picks, + threshold=threshold_muscle, + min_length_good=0.2) + +############################################################################### +# Plot movement z-scores across recording +# -------------------------------------------------------------------------- + +plt.figure() +plt.plot(raw.times, scores_muscle) +plt.axhline(y=threshold_muscle, color='r') +plt.show(block=False) +plt.title('High frequency ') +plt.xlabel('time, s.') +plt.ylabel('zscore') + +############################################################################### +# Plot raw with annotated muscles +# -------------------------------------------------------------------------- + +raw.set_annotations(annotation_muscle) +raw.plot(n_channels=100, duration=20) From 72c8c4d2b08cf05f2a90d03da9b04fb047fa0016 Mon Sep 17 00:00:00 2001 From: AdoNunes Date: Mon, 9 Mar 2020 16:54:57 -0700 Subject: [PATCH 06/31] example --- examples/preprocessing/plot_muscle_detection.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/examples/preprocessing/plot_muscle_detection.py b/examples/preprocessing/plot_muscle_detection.py index 048c851aa80..3eb296dfe76 100644 --- a/examples/preprocessing/plot_muscle_detection.py +++ b/examples/preprocessing/plot_muscle_detection.py @@ -6,13 +6,13 @@ Muscle contractions produce high frequency activity that can contaminate the brain signal of interest. Muscle artifacts can be produced when clinching the jaw, swalloing, or twitching a head muscle. Muscle artifacts are most notable -in teh range of 110-140Hz. +in the range of 110-140Hz. This example uses `annotate_muscle` to annotate segments where muscle activity likely occurred. This is done by band-pass filtering the data in the 110-140 Hz -range. Then, the envelope is taken to account for peacks and troughs. The +range. Then, the envelope is taken to account for peaks and troughs. The envelope is z-scored and averaged across channels. To remove noisy tansient -peak, the z-scored average is low-pass filtered to 4 Hz. Segments above a set +peaks, the z-scored average is low-pass filtered to 4 Hz. Segments above a set threshold are annotated as BAD_motion. In addition, `min_length_good` allows to discard god segments of data between bad segments that are to transient. From 39be3d33f61ab82ed2ffa5ab959f73405b68442a Mon Sep 17 00:00:00 2001 From: Adonay Nunes Date: Mon, 9 Mar 2020 19:25:33 -0700 Subject: [PATCH 07/31] for eeg --- .../preprocessing/plot_muscle_detection.py | 38 ++++++++++--------- mne/preprocessing/artifact_detection.py | 33 ++++++++++------ 2 files changed, 41 insertions(+), 30 deletions(-) diff --git a/examples/preprocessing/plot_muscle_detection.py b/examples/preprocessing/plot_muscle_detection.py index 3eb296dfe76..b7ce8218891 100644 --- a/examples/preprocessing/plot_muscle_detection.py +++ b/examples/preprocessing/plot_muscle_detection.py @@ -32,29 +32,31 @@ import matplotlib.pyplot as plt from mne import pick_types from mne.datasets.brainstorm import bst_auditory -from mne.io import read_raw_ctf +from mne.datasets import eegbci +from mne.io import read_raw_ctf, read_raw_edf from mne.preprocessing import annotate_muscle from mne.io import concatenate_raws # Load data - -data_path = bst_auditory.data_path() -data_path_MEG = op.join(data_path, 'MEG') -subject = 'bst_auditory' -subjects_dir = op.join(data_path, 'subjects') -trans_fname = op.join(data_path, 'MEG', 'bst_auditory', - 'bst_auditory-trans.fif') -raw_fname1 = op.join(data_path_MEG, 'bst_auditory', 'S01_AEF_20131218_01.ds') -raw_fname2 = op.join(data_path_MEG, 'bst_auditory', 'S01_AEF_20131218_02.ds') - -raw = read_raw_ctf(raw_fname1, preload=False) - -raw = concatenate_raws([raw, read_raw_ctf(raw_fname2, preload=False)]) -raw.crop(350, 410).load_data() +dataset = 'brainstorm' + +if dataset == 'eegbci': + fname = eegbci.load_data(2, runs=[3])[0] + raw = read_raw_edf(fname).load_data() + picks = pick_types(raw.info, eeg=True) +elif dataset == 'brainstorm': + data_path = bst_auditory.data_path() + data_path_MEG = op.join(data_path, 'MEG') + raw_fname1 = op.join(data_path_MEG, 'bst_auditory', 'S01_AEF_20131218_01.ds') + raw_fname2 = op.join(data_path_MEG, 'bst_auditory', 'S01_AEF_20131218_02.ds') + + raw = read_raw_ctf(raw_fname1, preload=False) + raw = concatenate_raws([raw, read_raw_ctf(raw_fname2, preload=False)]) + raw.crop(350, 410).load_data() + picks = pick_types(raw.info, meg=True, ref_meg=False) + raw.resample(300, npad="auto") -raw.notch_filter([50, 110, 140]) - -picks = pick_types(raw.info, meg=True, ref_meg=False) +raw.notch_filter([50, 100]) # detect muscle artifacts threshold_muscle = 1.5 # z-score diff --git a/mne/preprocessing/artifact_detection.py b/mne/preprocessing/artifact_detection.py index 05c6cb2d124..439908b001b 100644 --- a/mne/preprocessing/artifact_detection.py +++ b/mne/preprocessing/artifact_detection.py @@ -4,17 +4,17 @@ import numpy as np -from scipy.stats import zscore -from scipy.ndimage.measurements import label 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) -def annotate_muscle(raw, threshold=1.5, picks=None, min_length_good=.1): +@verbose +def annotate_muscle(raw, threshold=1.5, picks=None, min_length_good=.1, + verbose=None): """Detect segments with muscle artifacts. Detects segments periods that contains high frequency activity beyond the @@ -30,32 +30,41 @@ def annotate_muscle(raw, threshold=1.5, picks=None, min_length_good=.1): raw : instance of Raw Data to compute head position. threshold : float - The threshold for selecting segments with muscle activity artifacts. - picks : array - Channels to use for artifact detection. + The threshold in z-scores for selecting segments with muscle activity + artifacts. + %(picks_all)s min_length_good : int | float | None The minimal good segment length between annotations, smaller segments will be included in the movement annotation. + %(verbose)s Returns ------- annot : mne.Annotations - Periods with muscle artifacts. + Periods with muscle artifacts. If the z-score of scores_muscle : array Z-score values averaged across channels for each sample. """ + from scipy.stats import zscore + from scipy.ndimage.measurements import label + raw_copy = raw.copy() raw_copy.pick(picks) - raw_copy.pick_types(ref_meg=False) # Remove ref chans just in case + chan_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 chan_type)] + if meg is True: + raw_copy.pick_types(ref_meg=False) + # Only one type of channel, otherwise z-score will be biased - assert(len(set(raw_copy.get_channel_types())) == 1), 'Different channel ' \ - 'types, pick one type' + assert(len(set(chan_type)) == 1), 'Different channel types, pick one type' raw_copy.filter(110, 140, fir_design='firwin') raw_copy.apply_hilbert(envelope=True) sfreq = raw_copy.info['sfreq'] - art_scores = zscore(raw_copy._data, axis=1) + art_scores = zscore(raw_copy.get_data(), axis=1) scores_muscle = filter_data(art_scores.mean(axis=0), sfreq, None, 4) art_mask = scores_muscle > threshold From ab04ad736f94af883c3becf5b1b11c1d6439aca2 Mon Sep 17 00:00:00 2001 From: AdoNunes Date: Tue, 10 Mar 2020 21:13:28 -0700 Subject: [PATCH 08/31] skip_by_annotation --- examples/preprocessing/plot_muscle_detection.py | 11 ++++++----- mne/preprocessing/artifact_detection.py | 9 +++++---- 2 files changed, 11 insertions(+), 9 deletions(-) diff --git a/examples/preprocessing/plot_muscle_detection.py b/examples/preprocessing/plot_muscle_detection.py index b7ce8218891..a03733a1dec 100644 --- a/examples/preprocessing/plot_muscle_detection.py +++ b/examples/preprocessing/plot_muscle_detection.py @@ -11,7 +11,7 @@ This example uses `annotate_muscle` to annotate segments where muscle activity likely occurred. This is done by band-pass filtering the data in the 110-140 Hz range. Then, the envelope is taken to account for peaks and troughs. The -envelope is z-scored and averaged across channels. To remove noisy tansient +envelope is z-scored and averaged across channels. To remove noisy transient peaks, the z-scored average is low-pass filtered to 4 Hz. Segments above a set threshold are annotated as BAD_motion. In addition, `min_length_good` allows to discard god segments of data between bad segments that are to transient. @@ -47,14 +47,15 @@ elif dataset == 'brainstorm': data_path = bst_auditory.data_path() data_path_MEG = op.join(data_path, 'MEG') - raw_fname1 = op.join(data_path_MEG, 'bst_auditory', 'S01_AEF_20131218_01.ds') - raw_fname2 = op.join(data_path_MEG, 'bst_auditory', 'S01_AEF_20131218_02.ds') - + subject = 'bst_auditory' + raw_fname1 = op.join(data_path_MEG, subject, 'S01_AEF_20131218_01.ds') + raw_fname2 = op.join(data_path_MEG, subject, 'S01_AEF_20131218_02.ds') + raw = read_raw_ctf(raw_fname1, preload=False) raw = concatenate_raws([raw, read_raw_ctf(raw_fname2, preload=False)]) raw.crop(350, 410).load_data() picks = pick_types(raw.info, meg=True, ref_meg=False) - + raw.resample(300, npad="auto") raw.notch_filter([50, 100]) diff --git a/mne/preprocessing/artifact_detection.py b/mne/preprocessing/artifact_detection.py index 439908b001b..a0ee979aa77 100644 --- a/mne/preprocessing/artifact_detection.py +++ b/mne/preprocessing/artifact_detection.py @@ -41,13 +41,13 @@ def annotate_muscle(raw, threshold=1.5, picks=None, min_length_good=.1, Returns ------- annot : mne.Annotations - Periods with muscle artifacts. If the z-score of + Periods with muscle artifacts annotated as BAD_muscle. scores_muscle : array Z-score values averaged across channels for each sample. """ from scipy.stats import zscore from scipy.ndimage.measurements import label - + raw_copy = raw.copy() raw_copy.pick(picks) chan_type = raw_copy.get_channel_types() @@ -55,12 +55,13 @@ def annotate_muscle(raw, threshold=1.5, picks=None, min_length_good=.1, # Remove ref chans if MEG data just in case meg = [True for e in ['mag', 'grad'] if (e in chan_type)] if meg is True: - raw_copy.pick_types(ref_meg=False) + raw_copy.pick_types(ref_meg=False) # Only one type of channel, otherwise z-score will be biased assert(len(set(chan_type)) == 1), 'Different channel types, pick one type' - raw_copy.filter(110, 140, fir_design='firwin') + raw_copy.filter(110, 140, fir_design='firwin', + skip_by_annotation=['edge', 'BAD_']) raw_copy.apply_hilbert(envelope=True) sfreq = raw_copy.info['sfreq'] From 2a6691080f240fdd83ab13d8da73acab674f45ae Mon Sep 17 00:00:00 2001 From: AdoNunes Date: Mon, 16 Mar 2020 17:59:47 -0700 Subject: [PATCH 09/31] new func added --- .../preprocessing/plot_muscle_detection.py | 97 ++++++++++--------- mne/preprocessing/artifact_detection.py | 33 ++++--- .../plot_30_filtering_resampling.py | 1 + 3 files changed, 69 insertions(+), 62 deletions(-) diff --git a/examples/preprocessing/plot_muscle_detection.py b/examples/preprocessing/plot_muscle_detection.py index a03733a1dec..9cf05724fcf 100644 --- a/examples/preprocessing/plot_muscle_detection.py +++ b/examples/preprocessing/plot_muscle_detection.py @@ -3,25 +3,28 @@ Annotate muscle artifacts =========================== -Muscle contractions produce high frequency activity that can contaminate the -brain signal of interest. Muscle artifacts can be produced when clinching the -jaw, swalloing, or twitching a head muscle. Muscle artifacts are most notable -in the range of 110-140Hz. - -This example uses `annotate_muscle` to annotate segments where muscle activity -likely occurred. This is done by band-pass filtering the data in the 110-140 Hz -range. Then, the envelope is taken to account for peaks and troughs. The -envelope is z-scored and averaged across channels. To remove noisy transient -peaks, the z-scored average is low-pass filtered to 4 Hz. Segments above a set -threshold are annotated as BAD_motion. In addition, `min_length_good` allows to -discard god segments of data between bad segments that are to transient. - -Note: -The raw data needs to be notched to remove line activity, otherwise, changes of -AC power could bias the estimation. - -The inputted raw data has to be from a single channel type as there might be -different channel type magnitudes. +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 notable +in the range of 110-140 Hz. + +This example uses :func:`~mne.preprocessing.annotate_muscle` 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 averaged +across channels. To remove noisy transient peaks (muscle artifacts last +several hundred milliseconds), the channel averaged z-scored is low-pass +filtered to 4 Hz. Segments above a set threshold are annotated as +``BAD_muscle``. In addition, ``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. + + +.. note:: + If line noise is present, you should perform notch-filtering *before* + detecting muscle artifacts. See :ref:`tut-section-line-noise` for an + example. """ # Authors: Adonay Nunes @@ -32,34 +35,33 @@ import matplotlib.pyplot as plt from mne import pick_types from mne.datasets.brainstorm import bst_auditory -from mne.datasets import eegbci -from mne.io import read_raw_ctf, read_raw_edf +from mne.io import read_raw_ctf from mne.preprocessing import annotate_muscle -from mne.io import concatenate_raws + # Load data -dataset = 'brainstorm' - -if dataset == 'eegbci': - fname = eegbci.load_data(2, runs=[3])[0] - raw = read_raw_edf(fname).load_data() - picks = pick_types(raw.info, eeg=True) -elif dataset == 'brainstorm': - data_path = bst_auditory.data_path() - data_path_MEG = op.join(data_path, 'MEG') - subject = 'bst_auditory' - raw_fname1 = op.join(data_path_MEG, subject, 'S01_AEF_20131218_01.ds') - raw_fname2 = op.join(data_path_MEG, subject, 'S01_AEF_20131218_02.ds') - - raw = read_raw_ctf(raw_fname1, preload=False) - raw = concatenate_raws([raw, read_raw_ctf(raw_fname2, preload=False)]) - raw.crop(350, 410).load_data() - picks = pick_types(raw.info, meg=True, ref_meg=False) +data_path = bst_auditory.data_path() +data_path_MEG = op.join(data_path, 'MEG') +subject = 'bst_auditory' +raw_fname = op.join(data_path_MEG, subject, '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.filter(1, 150) 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) + +# Remove line noise as line activity changes could bias the artifact detection. raw.notch_filter([50, 100]) -# detect muscle artifacts +# The threshold of 1.5 is generally well suited for magnetometer, axial +# gradiometer and eeg data. Planar gradiometers will need a lower threshold. threshold_muscle = 1.5 # z-score annotation_muscle, scores_muscle = annotate_muscle(raw, picks=picks, threshold=threshold_muscle, @@ -69,17 +71,16 @@ # Plot movement z-scores across recording # -------------------------------------------------------------------------- -plt.figure() -plt.plot(raw.times, scores_muscle) -plt.axhline(y=threshold_muscle, color='r') -plt.show(block=False) -plt.title('High frequency ') -plt.xlabel('time, s.') -plt.ylabel('zscore') +fig, ax = plt.subplots() +ax.plot(raw.times, scores_muscle) +ax.axhline(y=threshold_muscle, color='r') +ax.set_title('Muscle activity') +ax.set_xlabel('time, s.') +ax.set_ylabel('zscore') ############################################################################### # Plot raw with annotated muscles # -------------------------------------------------------------------------- raw.set_annotations(annotation_muscle) -raw.plot(n_channels=100, duration=20) +raw.plot(n_channels=150, duration=60) diff --git a/mne/preprocessing/artifact_detection.py b/mne/preprocessing/artifact_detection.py index a0ee979aa77..0a07775e365 100644 --- a/mne/preprocessing/artifact_detection.py +++ b/mne/preprocessing/artifact_detection.py @@ -14,16 +14,18 @@ @verbose def annotate_muscle(raw, threshold=1.5, picks=None, min_length_good=.1, - verbose=None): + filter_freq=[110, 140], verbose=None): """Detect segments with muscle artifacts. Detects segments periods that contains high frequency activity beyond the - specified threshold. Muscle artifacts are most notable in the range of 110- - 140Hz. + specified threshold. Muscle artifacts are most notable in the range of + 110-140 Hz. - Raw data is band-pass filtered between 110 and 140 Hz, the signal envelope - computed, z-scored across samples, channel averaged and low-pass - filtered to smooth transient peaks. + 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 averaged and low-pass filtered to better + capture beggining and end of muscle activity and false positive transient + peaks. Parameters ---------- @@ -36,6 +38,9 @@ def annotate_muscle(raw, threshold=1.5, picks=None, min_length_good=.1, min_length_good : int | float | None The minimal good segment length between annotations, smaller segments will be included in the movement annotation. + filter_freq : list + The lower and upper high frequency to filter the signal for muscle + detection. %(verbose)s Returns @@ -50,30 +55,30 @@ def annotate_muscle(raw, threshold=1.5, picks=None, min_length_good=.1, raw_copy = raw.copy() raw_copy.pick(picks) - chan_type = raw_copy.get_channel_types() + 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 chan_type)] + meg = [True for e in ['mag', 'grad'] if (e in ch_type)] if meg is True: raw_copy.pick_types(ref_meg=False) # Only one type of channel, otherwise z-score will be biased - assert(len(set(chan_type)) == 1), 'Different channel types, pick one type' + assert(len(set(ch_type)) == 1), 'Different channel types, pick one type' - raw_copy.filter(110, 140, fir_design='firwin', - skip_by_annotation=['edge', 'BAD_']) + raw_copy.filter(filter_freq[0], filter_freq[1], fir_design='firwin', + pad="reflect_limited") raw_copy.apply_hilbert(envelope=True) sfreq = raw_copy.info['sfreq'] - art_scores = zscore(raw_copy.get_data(), axis=1) + art_scores = zscore(raw_copy.get_data(reject_by_annotation="NaN"), axis=1) scores_muscle = filter_data(art_scores.mean(axis=0), sfreq, None, 4) art_mask = scores_muscle > threshold # remove artifact free periods shorter than min_length_good idx_min = min_length_good * sfreq comps, num_comps = label(art_mask == 0) - for l in range(1, num_comps + 1): - l_idx = np.nonzero(comps == l)[0] + for ith in range(1, num_comps + 1): + l_idx = np.nonzero(comps == ith)[0] if len(l_idx) < idx_min: art_mask[l_idx] = True diff --git a/tutorials/preprocessing/plot_30_filtering_resampling.py b/tutorials/preprocessing/plot_30_filtering_resampling.py index 0ead1e601af..628b4951624 100644 --- a/tutorials/preprocessing/plot_30_filtering_resampling.py +++ b/tutorials/preprocessing/plot_30_filtering_resampling.py @@ -100,6 +100,7 @@ mne.viz.plot_filter(filter_params, raw.info['sfreq'], flim=(0.01, 5)) +# .. _tut-section-line-noise ############################################################################### # Power line noise # ~~~~~~~~~~~~~~~~ From d477c6097b188a8e1aaed03909c2aa797db61aa8 Mon Sep 17 00:00:00 2001 From: AdoNunes Date: Mon, 16 Mar 2020 18:08:13 -0700 Subject: [PATCH 10/31] docstring --- mne/preprocessing/artifact_detection.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/mne/preprocessing/artifact_detection.py b/mne/preprocessing/artifact_detection.py index 0a07775e365..4e6e5dc1b41 100644 --- a/mne/preprocessing/artifact_detection.py +++ b/mne/preprocessing/artifact_detection.py @@ -24,7 +24,7 @@ def annotate_muscle(raw, threshold=1.5, picks=None, min_length_good=.1, 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 averaged and low-pass filtered to better - capture beggining and end of muscle activity and false positive transient + capture beginning and end of muscle activity and false positive transient peaks. Parameters @@ -77,8 +77,8 @@ def annotate_muscle(raw, threshold=1.5, picks=None, min_length_good=.1, # remove artifact free periods shorter than min_length_good idx_min = min_length_good * sfreq comps, num_comps = label(art_mask == 0) - for ith in range(1, num_comps + 1): - l_idx = np.nonzero(comps == ith)[0] + for com in range(1, num_comps + 1): + l_idx = np.nonzero(comps == com)[0] if len(l_idx) < idx_min: art_mask[l_idx] = True From 62a00eb7c2e85436b9931cd6b35e63e505bf0b65 Mon Sep 17 00:00:00 2001 From: AdoNunes Date: Mon, 16 Mar 2020 21:27:07 -0700 Subject: [PATCH 11/31] cross ref --- tutorials/preprocessing/plot_30_filtering_resampling.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tutorials/preprocessing/plot_30_filtering_resampling.py b/tutorials/preprocessing/plot_30_filtering_resampling.py index 628b4951624..39db54459f4 100644 --- a/tutorials/preprocessing/plot_30_filtering_resampling.py +++ b/tutorials/preprocessing/plot_30_filtering_resampling.py @@ -100,7 +100,7 @@ mne.viz.plot_filter(filter_params, raw.info['sfreq'], flim=(0.01, 5)) -# .. _tut-section-line-noise +# .. _tut-section-line-noise: ############################################################################### # Power line noise # ~~~~~~~~~~~~~~~~ From 8f289c5e8834da0f8e1a1a9a5db3d2bb6fce82be Mon Sep 17 00:00:00 2001 From: AdoNunes Date: Mon, 16 Mar 2020 21:46:27 -0700 Subject: [PATCH 12/31] cross ref2 --- tutorials/preprocessing/plot_30_filtering_resampling.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tutorials/preprocessing/plot_30_filtering_resampling.py b/tutorials/preprocessing/plot_30_filtering_resampling.py index 39db54459f4..a7e35d16fc4 100644 --- a/tutorials/preprocessing/plot_30_filtering_resampling.py +++ b/tutorials/preprocessing/plot_30_filtering_resampling.py @@ -100,8 +100,8 @@ mne.viz.plot_filter(filter_params, raw.info['sfreq'], flim=(0.01, 5)) -# .. _tut-section-line-noise: ############################################################################### +# .. _tut-section-line-noise: # Power line noise # ~~~~~~~~~~~~~~~~ # From a0722fcc9723af186aae7ad5064bd57e60414fa4 Mon Sep 17 00:00:00 2001 From: Adonay Nunes Date: Tue, 17 Mar 2020 07:56:30 -0700 Subject: [PATCH 13/31] Alex suggestions Co-Authored-By: Alexandre Gramfort --- mne/preprocessing/artifact_detection.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/mne/preprocessing/artifact_detection.py b/mne/preprocessing/artifact_detection.py index 4e6e5dc1b41..6df90748081 100644 --- a/mne/preprocessing/artifact_detection.py +++ b/mne/preprocessing/artifact_detection.py @@ -14,10 +14,10 @@ @verbose def annotate_muscle(raw, threshold=1.5, picks=None, min_length_good=.1, - filter_freq=[110, 140], verbose=None): + filter_freq=(110, 140), verbose=None): """Detect segments with muscle artifacts. - Detects segments periods that contains high frequency activity beyond the + Detects segments periods that contain high frequency activity beyond the specified threshold. Muscle artifacts are most notable in the range of 110-140 Hz. @@ -38,7 +38,7 @@ def annotate_muscle(raw, threshold=1.5, picks=None, min_length_good=.1, min_length_good : int | float | None The minimal good segment length between annotations, smaller segments will be included in the movement annotation. - filter_freq : list + filter_freq : list | tuple, default (110, 140) The lower and upper high frequency to filter the signal for muscle detection. %(verbose)s From 1eaee829580d390d490a5c00315d8ecafeaf8205 Mon Sep 17 00:00:00 2001 From: AdoNunes Date: Tue, 17 Mar 2020 10:31:28 -0700 Subject: [PATCH 14/31] new parameters --- examples/preprocessing/plot_muscle_detection.py | 5 +++-- mne/preprocessing/artifact_detection.py | 10 ++++++---- .../preprocessing/plot_30_filtering_resampling.py | 1 + 3 files changed, 10 insertions(+), 6 deletions(-) diff --git a/examples/preprocessing/plot_muscle_detection.py b/examples/preprocessing/plot_muscle_detection.py index 9cf05724fcf..d0f23fdcccc 100644 --- a/examples/preprocessing/plot_muscle_detection.py +++ b/examples/preprocessing/plot_muscle_detection.py @@ -33,6 +33,7 @@ import os.path as op import matplotlib.pyplot as plt +from numpy import arange from mne import pick_types from mne.datasets.brainstorm import bst_auditory from mne.io import read_raw_ctf @@ -81,6 +82,6 @@ ############################################################################### # Plot raw with annotated muscles # -------------------------------------------------------------------------- - +order = arange(220, 240) raw.set_annotations(annotation_muscle) -raw.plot(n_channels=150, duration=60) +raw.plot( duration=30, order=order) diff --git a/mne/preprocessing/artifact_detection.py b/mne/preprocessing/artifact_detection.py index 4e6e5dc1b41..2e042f561f6 100644 --- a/mne/preprocessing/artifact_detection.py +++ b/mne/preprocessing/artifact_detection.py @@ -14,7 +14,7 @@ @verbose def annotate_muscle(raw, threshold=1.5, picks=None, min_length_good=.1, - filter_freq=[110, 140], verbose=None): + filter_freq=[110, 140], n_jobs=1, verbose=None): """Detect segments with muscle artifacts. Detects segments periods that contains high frequency activity beyond the @@ -33,7 +33,8 @@ def annotate_muscle(raw, threshold=1.5, picks=None, min_length_good=.1, Data to compute head position. threshold : float The threshold in z-scores for selecting segments with muscle activity - artifacts. + artifacts. Check ``scores_muscle`` to see optimal thesholding for the + data. %(picks_all)s min_length_good : int | float | None The minimal good segment length between annotations, smaller segments @@ -41,6 +42,7 @@ def annotate_muscle(raw, threshold=1.5, picks=None, min_length_good=.1, filter_freq : list The lower and upper high frequency to filter the signal for muscle detection. + %(n_jobs)s %(verbose)s Returns @@ -66,8 +68,8 @@ def annotate_muscle(raw, threshold=1.5, picks=None, min_length_good=.1, assert(len(set(ch_type)) == 1), 'Different channel types, pick one type' raw_copy.filter(filter_freq[0], filter_freq[1], fir_design='firwin', - pad="reflect_limited") - raw_copy.apply_hilbert(envelope=True) + pad="reflect_limited", n_jobs=n_jobs) + raw_copy.apply_hilbert(envelope=True, n_jobs=n_jobs) sfreq = raw_copy.info['sfreq'] art_scores = zscore(raw_copy.get_data(reject_by_annotation="NaN"), axis=1) diff --git a/tutorials/preprocessing/plot_30_filtering_resampling.py b/tutorials/preprocessing/plot_30_filtering_resampling.py index a7e35d16fc4..ccb0170423d 100644 --- a/tutorials/preprocessing/plot_30_filtering_resampling.py +++ b/tutorials/preprocessing/plot_30_filtering_resampling.py @@ -102,6 +102,7 @@ ############################################################################### # .. _tut-section-line-noise: +# # Power line noise # ~~~~~~~~~~~~~~~~ # From a9d5d4d60dd2546a82d253320b442dcb34e80ecb Mon Sep 17 00:00:00 2001 From: AdoNunes Date: Tue, 17 Mar 2020 10:51:58 -0700 Subject: [PATCH 15/31] docstring --- examples/preprocessing/plot_muscle_detection.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/preprocessing/plot_muscle_detection.py b/examples/preprocessing/plot_muscle_detection.py index d0f23fdcccc..3f90f7c645e 100644 --- a/examples/preprocessing/plot_muscle_detection.py +++ b/examples/preprocessing/plot_muscle_detection.py @@ -84,4 +84,4 @@ # -------------------------------------------------------------------------- order = arange(220, 240) raw.set_annotations(annotation_muscle) -raw.plot( duration=30, order=order) +raw.plot(duration=30, order=order) From 38a328b7039761223a8a036c4a9add6b619a5fbe Mon Sep 17 00:00:00 2001 From: AdoNunes Date: Thu, 19 Mar 2020 22:04:52 -0700 Subject: [PATCH 16/31] API name change, zscore --- doc/python_reference.rst | 2 +- .../preprocessing/plot_muscle_detection.py | 21 ++++++++++++------- mne/preprocessing/__init__.py | 2 +- mne/preprocessing/artifact_detection.py | 17 +++++++++------ 4 files changed, 26 insertions(+), 16 deletions(-) diff --git a/doc/python_reference.rst b/doc/python_reference.rst index fdce1f388f7..6b0e8e445d9 100644 --- a/doc/python_reference.rst +++ b/doc/python_reference.rst @@ -346,7 +346,7 @@ Projections: ICA Xdawn annotate_movement - annotate_muscle + annotate_muscle_zscore compute_average_dev_head_t compute_current_source_density compute_proj_ecg diff --git a/examples/preprocessing/plot_muscle_detection.py b/examples/preprocessing/plot_muscle_detection.py index 3f90f7c645e..32c324278bc 100644 --- a/examples/preprocessing/plot_muscle_detection.py +++ b/examples/preprocessing/plot_muscle_detection.py @@ -37,7 +37,7 @@ 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 +from mne.preprocessing import annotate_muscle_zscore # Load data @@ -58,15 +58,20 @@ # they are more sensitive to muscle activity picks = pick_types(raw.info, meg=True, ref_meg=False) -# Remove line noise as line activity changes could bias the artifact detection. +""" +.. 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 of 1.5 is generally well suited for magnetometer, axial +# The threshold of 4 is generally well suited for magnetometer, axial # gradiometer and eeg data. Planar gradiometers will need a lower threshold. -threshold_muscle = 1.5 # z-score -annotation_muscle, scores_muscle = annotate_muscle(raw, picks=picks, - threshold=threshold_muscle, - min_length_good=0.2) +threshold_muscle = 4 # z-score +out = annotate_muscle_zscore(raw, picks=picks, threshold=threshold_muscle, + min_length_good=0.2, filter_freq=[90, 120]) +annot_muscle, scores_muscle = out ############################################################################### # Plot movement z-scores across recording @@ -83,5 +88,5 @@ # Plot raw with annotated muscles # -------------------------------------------------------------------------- order = arange(220, 240) -raw.set_annotations(annotation_muscle) +raw.set_annotations(annot_muscle) raw.plot(duration=30, order=order) diff --git a/mne/preprocessing/__init__.py b/mne/preprocessing/__init__.py index 6210ad8d871..387341c8eac 100644 --- a/mne/preprocessing/__init__.py +++ b/mne/preprocessing/__init__.py @@ -25,4 +25,4 @@ from ._csd import compute_current_source_density from . import nirs from .artifact_detection import (annotate_movement, compute_average_dev_head_t, - annotate_muscle) + annotate_muscle_zscore) diff --git a/mne/preprocessing/artifact_detection.py b/mne/preprocessing/artifact_detection.py index a726b5282c2..6c52e0ca515 100644 --- a/mne/preprocessing/artifact_detection.py +++ b/mne/preprocessing/artifact_detection.py @@ -13,12 +13,12 @@ @verbose -def annotate_muscle(raw, threshold=1.5, picks=None, min_length_good=.1, - filter_freq=[110, 140], n_jobs=1, verbose=None): +def annotate_muscle_zscore(raw, threshold=4, picks=None, min_length_good=.1, + filter_freq=[110, 140], n_jobs=1, verbose=None): """Detect segments with muscle artifacts. Detects segments periods that contain high frequency activity beyond the - specified threshold. Muscle artifacts are most notable in the range of + specified threshold. Muscle artifacts are most detectable in the range of 110-140 Hz. Raw data is band-pass filtered between ``filter_freq`` especified @@ -52,7 +52,7 @@ def annotate_muscle(raw, threshold=1.5, picks=None, min_length_good=.1, scores_muscle : array Z-score values averaged across channels for each sample. """ - from scipy.stats import zscore + from scipy.stats.mstats import zscore from scipy.ndimage.measurements import label raw_copy = raw.copy() @@ -72,8 +72,13 @@ def annotate_muscle(raw, threshold=1.5, picks=None, min_length_good=.1, raw_copy.apply_hilbert(envelope=True, n_jobs=n_jobs) sfreq = raw_copy.info['sfreq'] - art_scores = zscore(raw_copy.get_data(reject_by_annotation="NaN"), axis=1) - scores_muscle = filter_data(art_scores.mean(axis=0), sfreq, None, 4) + art_scores = zscore(raw_copy.get_data(reject_by_annotation="NaN"), axis=1, + nan_policy='omit') + + art_scores = art_scores.sum(axis=0) / np.sqrt(art_scores.shape[0]) + + scores_muscle = filter_data(art_scores, sfreq, None, 4) + scores_muscle[np.isnan(scores_muscle)] = threshold art_mask = scores_muscle > threshold # remove artifact free periods shorter than min_length_good From b232e4f0d3c03911a5f4e50fe42c33e768b29c3c Mon Sep 17 00:00:00 2001 From: AdoNunes Date: Thu, 19 Mar 2020 22:25:25 -0700 Subject: [PATCH 17/31] test update --- mne/preprocessing/tests/test_artifact_detection.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/mne/preprocessing/tests/test_artifact_detection.py b/mne/preprocessing/tests/test_artifact_detection.py index 1612c811a43..a9e99cc7af9 100755 --- a/mne/preprocessing/tests/test_artifact_detection.py +++ b/mne/preprocessing/tests/test_artifact_detection.py @@ -10,7 +10,7 @@ from mne.datasets import testing from mne.io import read_raw_fif from mne.preprocessing import (annotate_movement, compute_average_dev_head_t, - annotate_muscle) + annotate_muscle_zscore) data_path = testing.data_path(download=False) sss_path = op.join(data_path, 'SSS') @@ -54,8 +54,8 @@ 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='grad', ref_meg=False) + 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(raw, threshold=1) + annot_muscle, scores = annotate_muscle_zscore(raw, threshold=10) assert(annot_muscle.duration.size == 2) From d9412dbb0b4a10869fd04950901751ebc2a046c7 Mon Sep 17 00:00:00 2001 From: AdoNunes Date: Sun, 22 Mar 2020 17:14:18 -0700 Subject: [PATCH 18/31] sqrt nchans --- .../preprocessing/plot_muscle_detection.py | 41 ++++++++----------- mne/preprocessing/artifact_detection.py | 8 ++-- .../tests/test_artifact_detection.py | 3 ++ 3 files changed, 23 insertions(+), 29 deletions(-) diff --git a/examples/preprocessing/plot_muscle_detection.py b/examples/preprocessing/plot_muscle_detection.py index 32c324278bc..02d847a531b 100644 --- a/examples/preprocessing/plot_muscle_detection.py +++ b/examples/preprocessing/plot_muscle_detection.py @@ -8,23 +8,17 @@ swallowing, or twitching a head muscle. Muscle artifacts are most notable in the range of 110-140 Hz. -This example uses :func:`~mne.preprocessing.annotate_muscle` to annotate +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 averaged -across channels. To remove noisy transient peaks (muscle artifacts last -several hundred milliseconds), the channel averaged z-scored is low-pass -filtered to 4 Hz. Segments above a set threshold are annotated as -``BAD_muscle``. In addition, ``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. - - -.. note:: - If line noise is present, you should perform notch-filtering *before* - detecting muscle artifacts. See :ref:`tut-section-line-noise` for an - example. +the phase of the high frequency signal. The envelope is z-scored and summed +across channels and divedid by the square root of the channel numer. To remove +noisy transient peaks (muscle artifacts last several hundred milliseconds), +the channel averaged z-scored is low-pass filtered to 4 Hz. Segments above a +set threshold are annotated as ``BAD_muscle``. In addition, ``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 @@ -33,7 +27,7 @@ import os.path as op import matplotlib.pyplot as plt -from numpy import arange +import numpy as np from mne import pick_types from mne.datasets.brainstorm import bst_auditory from mne.io import read_raw_ctf @@ -49,7 +43,6 @@ raw = read_raw_ctf(raw_fname, preload=False) raw.crop(130, 160).load_data() # just use a fraction of data for speed here -raw.filter(1, 150) raw.resample(300, npad="auto") # The inputted raw data has to be from a single channel type as there might be @@ -58,16 +51,17 @@ # they are more sensitive to muscle activity picks = pick_types(raw.info, meg=True, ref_meg=False) +############################################################################### +raw.notch_filter([50, 100]) """ .. 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 of 4 is generally well suited for magnetometer, axial -# gradiometer and eeg data. Planar gradiometers will need a lower threshold. +# The threshold is data dependent, check the optimal threshold by plotting +# ``threshold_muscle`` threshold_muscle = 4 # z-score out = annotate_muscle_zscore(raw, picks=picks, threshold=threshold_muscle, min_length_good=0.2, filter_freq=[90, 120]) @@ -80,13 +74,10 @@ fig, ax = plt.subplots() ax.plot(raw.times, scores_muscle) ax.axhline(y=threshold_muscle, color='r') -ax.set_title('Muscle activity') -ax.set_xlabel('time, s.') -ax.set_ylabel('zscore') - +ax.set(xlabel='time, s.', ylabel='zscore', title='Muscle activity') ############################################################################### # Plot raw with annotated muscles # -------------------------------------------------------------------------- -order = arange(220, 240) +order = np.arange(220, 240) raw.set_annotations(annot_muscle) -raw.plot(duration=30, order=order) +raw.plot(duration=20, order=order) diff --git a/mne/preprocessing/artifact_detection.py b/mne/preprocessing/artifact_detection.py index 6c52e0ca515..fbc60661f46 100644 --- a/mne/preprocessing/artifact_detection.py +++ b/mne/preprocessing/artifact_detection.py @@ -14,7 +14,7 @@ @verbose def annotate_muscle_zscore(raw, threshold=4, picks=None, min_length_good=.1, - filter_freq=[110, 140], n_jobs=1, verbose=None): + filter_freq=(110, 140), n_jobs=1, verbose=None): """Detect segments with muscle artifacts. Detects segments periods that contain high frequency activity beyond the @@ -23,9 +23,9 @@ def annotate_muscle_zscore(raw, threshold=4, picks=None, min_length_good=.1, 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 averaged and low-pass filtered to better - capture beginning and end of muscle activity and false positive transient - peaks. + 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. Parameters ---------- diff --git a/mne/preprocessing/tests/test_artifact_detection.py b/mne/preprocessing/tests/test_artifact_detection.py index a9e99cc7af9..8f9231b9e14 100755 --- a/mne/preprocessing/tests/test_artifact_detection.py +++ b/mne/preprocessing/tests/test_artifact_detection.py @@ -58,4 +58,7 @@ def test_muscle_annotation(): raw.notch_filter([50, 110, 150]) # Check 2 muscle segments are detected annot_muscle, scores = annotate_muscle_zscore(raw, threshold=10) + onset = annot_muscle.onset * raw.info['sfreq'] + onset = onset.astype(int) + assert(np.all(scores[onset].astype(int) == np.array([23, 10]))) assert(annot_muscle.duration.size == 2) From b5e378e388fe589d5c06d3f3edaac0199e12bdae Mon Sep 17 00:00:00 2001 From: AdoNunes Date: Mon, 23 Mar 2020 09:05:38 -0700 Subject: [PATCH 19/31] handle nans zscore --- examples/preprocessing/plot_muscle_detection.py | 16 ++++++++-------- mne/preprocessing/artifact_detection.py | 15 ++++++++++----- 2 files changed, 18 insertions(+), 13 deletions(-) diff --git a/examples/preprocessing/plot_muscle_detection.py b/examples/preprocessing/plot_muscle_detection.py index 02d847a531b..1d704c28505 100644 --- a/examples/preprocessing/plot_muscle_detection.py +++ b/examples/preprocessing/plot_muscle_detection.py @@ -52,13 +52,13 @@ picks = pick_types(raw.info, meg=True, ref_meg=False) ############################################################################### +# 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]) -""" -.. note:: - If line noise is present, you should perform notch-filtering *before* - detecting muscle artifacts. See :ref:`tut-section-line-noise` for an - example. -""" # The threshold is data dependent, check the optimal threshold by plotting # ``threshold_muscle`` @@ -78,6 +78,6 @@ ############################################################################### # Plot raw with annotated muscles # -------------------------------------------------------------------------- -order = np.arange(220, 240) +order = np.arange(144, 164) raw.set_annotations(annot_muscle) -raw.plot(duration=20, order=order) +raw.plot(start=5, duration=20, order=order) diff --git a/mne/preprocessing/artifact_detection.py b/mne/preprocessing/artifact_detection.py index fbc60661f46..9eadd409a71 100644 --- a/mne/preprocessing/artifact_detection.py +++ b/mne/preprocessing/artifact_detection.py @@ -52,7 +52,7 @@ def annotate_muscle_zscore(raw, threshold=4, picks=None, min_length_good=.1, scores_muscle : array Z-score values averaged across channels for each sample. """ - from scipy.stats.mstats import zscore + from scipy.stats import zscore from scipy.ndimage.measurements import label raw_copy = raw.copy() @@ -72,14 +72,19 @@ def annotate_muscle_zscore(raw, threshold=4, picks=None, min_length_good=.1, raw_copy.apply_hilbert(envelope=True, n_jobs=n_jobs) sfreq = raw_copy.info['sfreq'] - art_scores = zscore(raw_copy.get_data(reject_by_annotation="NaN"), axis=1, - nan_policy='omit') + data = raw_copy.get_data(reject_by_annotation="NaN") + nan_mask = ~np.isnan(data[0, :]) + 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 - scores_muscle = filter_data(art_scores, sfreq, None, 4) - scores_muscle[np.isnan(scores_muscle)] = threshold 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 idx_min = min_length_good * sfreq From b055b53d745ab1ce62b03b47d2287638953be408 Mon Sep 17 00:00:00 2001 From: AdoNunes Date: Mon, 23 Mar 2020 09:42:49 -0700 Subject: [PATCH 20/31] example note --- examples/preprocessing/plot_muscle_detection.py | 1 + 1 file changed, 1 insertion(+) diff --git a/examples/preprocessing/plot_muscle_detection.py b/examples/preprocessing/plot_muscle_detection.py index 1d704c28505..41d862b2aeb 100644 --- a/examples/preprocessing/plot_muscle_detection.py +++ b/examples/preprocessing/plot_muscle_detection.py @@ -53,6 +53,7 @@ ############################################################################### # 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 From 5c36c037d9a15396f703b04dc93e4aedacfe166a Mon Sep 17 00:00:00 2001 From: Adonay Nunes Date: Mon, 23 Mar 2020 19:55:28 -0700 Subject: [PATCH 21/31] Apply suggestions from code review drammock Co-Authored-By: Daniel McCloy --- .../preprocessing/plot_muscle_detection.py | 17 +++++++------ mne/preprocessing/artifact_detection.py | 25 +++++++++---------- 2 files changed, 21 insertions(+), 21 deletions(-) diff --git a/examples/preprocessing/plot_muscle_detection.py b/examples/preprocessing/plot_muscle_detection.py index 41d862b2aeb..5c0c7a757b1 100644 --- a/examples/preprocessing/plot_muscle_detection.py +++ b/examples/preprocessing/plot_muscle_detection.py @@ -13,10 +13,12 @@ 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 divedid by the square root of the channel numer. To remove -noisy transient peaks (muscle artifacts last several hundred milliseconds), -the channel averaged z-scored is low-pass filtered to 4 Hz. Segments above a -set threshold are annotated as ``BAD_muscle``. In addition, ``min_length_good`` +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. @@ -36,15 +38,14 @@ # Load data data_path = bst_auditory.data_path() -data_path_MEG = op.join(data_path, 'MEG') -subject = 'bst_auditory' -raw_fname = op.join(data_path_MEG, subject, 'S01_AEF_20131218_01.ds') +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") +############################################################################### # 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 @@ -77,7 +78,7 @@ ax.axhline(y=threshold_muscle, color='r') ax.set(xlabel='time, s.', ylabel='zscore', title='Muscle activity') ############################################################################### -# Plot raw with annotated muscles +# View the annotations # -------------------------------------------------------------------------- order = np.arange(144, 164) raw.set_annotations(annot_muscle) diff --git a/mne/preprocessing/artifact_detection.py b/mne/preprocessing/artifact_detection.py index 9eadd409a71..3671c82f6f7 100644 --- a/mne/preprocessing/artifact_detection.py +++ b/mne/preprocessing/artifact_detection.py @@ -13,13 +13,12 @@ @verbose -def annotate_muscle_zscore(raw, threshold=4, picks=None, min_length_good=.1, +def annotate_muscle_zscore(raw, threshold=4, picks=None, min_length_good=0.1, filter_freq=(110, 140), n_jobs=1, verbose=None): - """Detect segments with muscle artifacts. + """Detect segments that likely contain muscle artifacts. - Detects segments periods that contain high frequency activity beyond the - specified threshold. Muscle artifacts are most detectable in the range of - 110-140 Hz. + 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, @@ -32,16 +31,16 @@ def annotate_muscle_zscore(raw, threshold=4, picks=None, min_length_good=.1, raw : instance of Raw Data to compute head position. threshold : float - The threshold in z-scores for selecting segments with muscle activity - artifacts. Check ``scores_muscle`` to see optimal thesholding for the - data. + The threshold in z-scores for marking segments as containg muscle + activity artifacts. %(picks_all)s min_length_good : int | float | None - The minimal good segment length between annotations, smaller segments - will be included in the movement annotation. - filter_freq : list | tuple, default (110, 140) - The lower and upper high frequency to filter the signal for muscle - detection. + The shortest allowed duration of "good data" (in seconds) between + adjacent annotations; shorter segments will be incorporated into the + surrounding annotations. + 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 From ee8e871734d8e103329ec646fb830c099cca9c05 Mon Sep 17 00:00:00 2001 From: AdoNunes Date: Mon, 23 Mar 2020 21:19:33 -0700 Subject: [PATCH 22/31] revision changes --- .../preprocessing/plot_muscle_detection.py | 23 +++++++++---------- mne/preprocessing/artifact_detection.py | 23 ++++++++++++------- .../tests/test_artifact_detection.py | 3 ++- 3 files changed, 28 insertions(+), 21 deletions(-) diff --git a/examples/preprocessing/plot_muscle_detection.py b/examples/preprocessing/plot_muscle_detection.py index 5c0c7a757b1..88b5ea9c0f8 100644 --- a/examples/preprocessing/plot_muscle_detection.py +++ b/examples/preprocessing/plot_muscle_detection.py @@ -5,7 +5,7 @@ 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 notable +swallowing, or twitching a head muscle. Muscle artifacts are most noticeable in the range of 110-140 Hz. This example uses :func:`~mne.preprocessing.annotate_muscle_zscore` to annotate @@ -16,11 +16,10 @@ 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. +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 @@ -49,7 +48,7 @@ # 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 +# they are more sensitive to muscle activity. picks = pick_types(raw.info, meg=True, ref_meg=False) ############################################################################### @@ -63,20 +62,20 @@ raw.notch_filter([50, 100]) # The threshold is data dependent, check the optimal threshold by plotting -# ``threshold_muscle`` -threshold_muscle = 4 # z-score +# ``scores_muscle``. +threshold_muscle = 5 # z-score out = annotate_muscle_zscore(raw, picks=picks, threshold=threshold_muscle, - min_length_good=0.2, filter_freq=[90, 120]) + min_length_good=0.2, filter_freq=[110, 140]) annot_muscle, scores_muscle = out ############################################################################### -# Plot movement z-scores across recording +# 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') +ax.set(xlabel='time, (s)', ylabel='zscore', title='Muscle activity') ############################################################################### # View the annotations # -------------------------------------------------------------------------- diff --git a/mne/preprocessing/artifact_detection.py b/mne/preprocessing/artifact_detection.py index 3671c82f6f7..5c2c5192c9d 100644 --- a/mne/preprocessing/artifact_detection.py +++ b/mne/preprocessing/artifact_detection.py @@ -26,18 +26,22 @@ def annotate_muscle_zscore(raw, threshold=4, picks=None, min_length_good=0.1, 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. + Parameters ---------- raw : instance of Raw - Data to compute head position. + Data to estimate segments with muscle artifacts. threshold : float The threshold in z-scores for marking segments as containg muscle activity artifacts. %(picks_all)s - min_length_good : int | float | None + 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. + 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)``. @@ -60,11 +64,13 @@ def annotate_muscle_zscore(raw, threshold=4, picks=None, min_length_good=0.1, # Remove ref chans if MEG data just in case meg = [True for e in ['mag', 'grad'] if (e in ch_type)] - if meg is True: + 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(set(ch_type)) == 1), 'Different channel types, pick one type' + assert(len(np.unique(ch_type)) == 1), \ + 'Different channel types, pick one type' raw_copy.filter(filter_freq[0], filter_freq[1], fir_design='firwin', pad="reflect_limited", n_jobs=n_jobs) @@ -72,7 +78,7 @@ def annotate_muscle_zscore(raw, threshold=4, picks=None, min_length_good=0.1, sfreq = raw_copy.info['sfreq'] data = raw_copy.get_data(reject_by_annotation="NaN") - nan_mask = ~np.isnan(data[0, :]) + nan_mask = ~np.isnan(data[0]) art_scores = zscore(data[:, nan_mask], axis=1) art_scores = art_scores.sum(axis=0) / np.sqrt(art_scores.shape[0]) @@ -86,11 +92,12 @@ def annotate_muscle_zscore(raw, threshold=4, picks=None, min_length_good=0.1, scores_muscle[~nan_mask] = np.nan # remove artifact free periods shorter than min_length_good - idx_min = min_length_good * sfreq + 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) < idx_min: + if len(l_idx) < min_samps: art_mask[l_idx] = True annot = _annotations_from_mask(raw_copy.times, art_mask, 'BAD_muscle') diff --git a/mne/preprocessing/tests/test_artifact_detection.py b/mne/preprocessing/tests/test_artifact_detection.py index 8f9231b9e14..35742fc7f5c 100755 --- a/mne/preprocessing/tests/test_artifact_detection.py +++ b/mne/preprocessing/tests/test_artifact_detection.py @@ -60,5 +60,6 @@ def test_muscle_annotation(): annot_muscle, scores = annotate_muscle_zscore(raw, threshold=10) onset = annot_muscle.onset * raw.info['sfreq'] onset = onset.astype(int) - assert(np.all(scores[onset].astype(int) == np.array([23, 10]))) + np.testing.assert_array_equal(scores[onset].astype(int), np.array([23, + 10])) assert(annot_muscle.duration.size == 2) From 0953d9a8a013b60af74590c246a67b37bdbd73dd Mon Sep 17 00:00:00 2001 From: AdoNunes Date: Mon, 23 Mar 2020 21:26:21 -0700 Subject: [PATCH 23/31] docstring --- mne/preprocessing/artifact_detection.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/mne/preprocessing/artifact_detection.py b/mne/preprocessing/artifact_detection.py index 5c2c5192c9d..1a35ce6bc4c 100644 --- a/mne/preprocessing/artifact_detection.py +++ b/mne/preprocessing/artifact_detection.py @@ -26,8 +26,8 @@ def annotate_muscle_zscore(raw, threshold=4, picks=None, min_length_good=0.1, 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. + .. note:: + Use a single channel type. Parameters ---------- From 647b9abd67ccf2caa78bb0c470f0508f57f4e106 Mon Sep 17 00:00:00 2001 From: AdoNunes Date: Mon, 23 Mar 2020 21:27:35 -0700 Subject: [PATCH 24/31] docstring --- mne/preprocessing/artifact_detection.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/mne/preprocessing/artifact_detection.py b/mne/preprocessing/artifact_detection.py index 1a35ce6bc4c..b1700fb9595 100644 --- a/mne/preprocessing/artifact_detection.py +++ b/mne/preprocessing/artifact_detection.py @@ -34,7 +34,7 @@ def annotate_muscle_zscore(raw, threshold=4, picks=None, min_length_good=0.1, raw : instance of Raw Data to estimate segments with muscle artifacts. threshold : float - The threshold in z-scores for marking segments as containg muscle + The threshold in z-scores for marking segments as containing muscle activity artifacts. %(picks_all)s min_length_good : float | None From 973c1d7adcfd60f1269a2952f452ab46c1790cec Mon Sep 17 00:00:00 2001 From: AdoNunes Date: Wed, 25 Mar 2020 16:38:46 -0700 Subject: [PATCH 25/31] ch_type param --- doc/changes/latest.inc | 3 + doc/references.bib | 10 ++++ .../preprocessing/plot_muscle_detection.py | 32 +++++------ mne/preprocessing/artifact_detection.py | 57 +++++++++++-------- .../tests/test_artifact_detection.py | 4 +- 5 files changed, 61 insertions(+), 45 deletions(-) diff --git a/doc/changes/latest.inc b/doc/changes/latest.inc index 98797cdb4ef..0b3b4491bb8 100644 --- a/doc/changes/latest.inc +++ b/doc/changes/latest.inc @@ -131,6 +131,9 @@ Changelog - Expose the number of ICA iterations during the fitting procedure via the ``n_iter_`` attribute of :class:`mne.preprocessing.ICA` by `Richard Höchenberger`_ +- Add function :func:`mne.preprocessing.annotate_muscle_zscore` to annotate periods with muscle artifacts. by `Adonay Nunes`_ + + Bug ~~~ diff --git a/doc/references.bib b/doc/references.bib index 091bcd00df5..f28ccc3a8f7 100644 --- a/doc/references.bib +++ b/doc/references.bib @@ -1026,6 +1026,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 index 88b5ea9c0f8..c62f11c52db 100644 --- a/examples/preprocessing/plot_muscle_detection.py +++ b/examples/preprocessing/plot_muscle_detection.py @@ -5,8 +5,8 @@ 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 @@ -14,12 +14,12 @@ 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 @@ -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 diff --git a/mne/preprocessing/artifact_detection.py b/mne/preprocessing/artifact_detection.py index b1700fb9595..57ce5e4aa6d 100644 --- a/mne/preprocessing/artifact_detection.py +++ b/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,38 @@ 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(f"Using {ch_type} sensors for muscle artifact detection.") + + 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]) diff --git a/mne/preprocessing/tests/test_artifact_detection.py b/mne/preprocessing/tests/test_artifact_detection.py index 35742fc7f5c..92bd03dba7a 100755 --- a/mne/preprocessing/tests/test_artifact_detection.py +++ b/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, From 5701cb018bfb36a84ff58b5f2a767f66b06a3be7 Mon Sep 17 00:00:00 2001 From: Adonay Nunes Date: Wed, 25 Mar 2020 21:25:03 -0700 Subject: [PATCH 26/31] logger --- mne/preprocessing/artifact_detection.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/mne/preprocessing/artifact_detection.py b/mne/preprocessing/artifact_detection.py index 57ce5e4aa6d..d7f27da49c7 100644 --- a/mne/preprocessing/artifact_detection.py +++ b/mne/preprocessing/artifact_detection.py @@ -72,7 +72,8 @@ def annotate_muscle_zscore(raw, threshold=4, ch_type=None, min_length_good=0.1, ch_type = 'grad' elif 'eeg' in raw_ch_type: ch_type = 'eeg' - logger.info(f"Using {ch_type} sensors for muscle artifact detection.") + 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) From 416ab4704149b5f4b0710cb80476c8e21917fab0 Mon Sep 17 00:00:00 2001 From: Adonay Nunes Date: Thu, 26 Mar 2020 12:43:07 -0700 Subject: [PATCH 27/31] docstrings --- mne/preprocessing/artifact_detection.py | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/mne/preprocessing/artifact_detection.py b/mne/preprocessing/artifact_detection.py index d7f27da49c7..02e25b8e7b2 100644 --- a/mne/preprocessing/artifact_detection.py +++ b/mne/preprocessing/artifact_detection.py @@ -18,14 +18,14 @@ def annotate_muscle_zscore(raw, threshold=4, ch_type=None, min_length_good=0.1, """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)`). + ``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``. + 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`. + ``filter_freq`` and ``threshold``. Parameters ---------- @@ -36,7 +36,7 @@ def annotate_muscle_zscore(raw, threshold=4, ch_type=None, min_length_good=0.1, 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``. + ``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 @@ -72,7 +72,7 @@ def annotate_muscle_zscore(raw, threshold=4, ch_type=None, min_length_good=0.1, ch_type = 'grad' elif 'eeg' in raw_ch_type: ch_type = 'eeg' - logger.info(u'Using %s sensors for muscle artifact detection' + logger.info('Using %s sensors for muscle artifact detection' % (ch_type)) if ch_type in ('mag', 'grad'): From 269277e689e8c3d7fd7217aa78891e29d7b6bf41 Mon Sep 17 00:00:00 2001 From: AdoNunes Date: Mon, 30 Mar 2020 14:48:21 -0700 Subject: [PATCH 28/31] ch type --- examples/preprocessing/plot_muscle_detection.py | 6 +++--- mne/preprocessing/artifact_detection.py | 8 ++++++-- 2 files changed, 9 insertions(+), 5 deletions(-) diff --git a/examples/preprocessing/plot_muscle_detection.py b/examples/preprocessing/plot_muscle_detection.py index c62f11c52db..7befc482440 100644 --- a/examples/preprocessing/plot_muscle_detection.py +++ b/examples/preprocessing/plot_muscle_detection.py @@ -60,9 +60,9 @@ 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. -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 +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 diff --git a/mne/preprocessing/artifact_detection.py b/mne/preprocessing/artifact_detection.py index 02e25b8e7b2..c3bebcdb827 100644 --- a/mne/preprocessing/artifact_detection.py +++ b/mne/preprocessing/artifact_detection.py @@ -72,13 +72,17 @@ def annotate_muscle_zscore(raw, threshold=4, ch_type=None, min_length_good=0.1, 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) - elif ch_type == 'eeg': - raw_copy.pick_types(meg=False, eeg=True) + 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) From d82c6147067b7e28e6faec2348775b19ff3cc4db Mon Sep 17 00:00:00 2001 From: AdoNunes Date: Tue, 7 Apr 2020 10:08:47 -0700 Subject: [PATCH 29/31] test no meeg data --- mne/preprocessing/tests/test_artifact_detection.py | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/mne/preprocessing/tests/test_artifact_detection.py b/mne/preprocessing/tests/test_artifact_detection.py index 92bd03dba7a..7329a9b7890 100755 --- a/mne/preprocessing/tests/test_artifact_detection.py +++ b/mne/preprocessing/tests/test_artifact_detection.py @@ -5,6 +5,7 @@ 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 @@ -63,3 +64,14 @@ def test_muscle_annotation(): 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) + + \ No newline at end of file From 2e1531d4cb768f8301aee3e637ea4cd63256e16f Mon Sep 17 00:00:00 2001 From: AdoNunes Date: Tue, 7 Apr 2020 10:11:31 -0700 Subject: [PATCH 30/31] test no meeg data --- mne/preprocessing/tests/test_artifact_detection.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/mne/preprocessing/tests/test_artifact_detection.py b/mne/preprocessing/tests/test_artifact_detection.py index 7329a9b7890..db716e461fa 100755 --- a/mne/preprocessing/tests/test_artifact_detection.py +++ b/mne/preprocessing/tests/test_artifact_detection.py @@ -65,6 +65,7 @@ def test_muscle_annotation(): 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.""" @@ -73,5 +74,3 @@ def test_muscle_annotation_without_meeg_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) - - \ No newline at end of file From d1d448d7373431caeaf1d8636178a60cb3e9a9ee Mon Sep 17 00:00:00 2001 From: AdoNunes Date: Tue, 7 Apr 2020 10:16:10 -0700 Subject: [PATCH 31/31] fix docstring --- mne/preprocessing/tests/test_artifact_detection.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/mne/preprocessing/tests/test_artifact_detection.py b/mne/preprocessing/tests/test_artifact_detection.py index db716e461fa..4a39830454a 100755 --- a/mne/preprocessing/tests/test_artifact_detection.py +++ b/mne/preprocessing/tests/test_artifact_detection.py @@ -68,7 +68,7 @@ def test_muscle_annotation(): @testing.requires_testing_data def test_muscle_annotation_without_meeg_data(): - """ Call annotate_muscle_zscore with data without meg or eeg.""" + """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)