diff --git a/doc/changes/latest.inc b/doc/changes/latest.inc index 173f1b11679..c72ce384acd 100644 --- a/doc/changes/latest.inc +++ b/doc/changes/latest.inc @@ -297,6 +297,8 @@ Bugs - Fix bug where ``backend='notebook'`` could not be used in :meth:`mne.SourceEstimate.plot` (:gh:`9305` by `Jean-Remi King`_) +- `mne.preprocessing.find_eog_events`, `mne.preprocessing.find_ecg_events`, `mne.preprocessing.create_eog_epochs`, and `mne.preprocessing.create_ecg_epochs` now have better handling for situations where the signal is flat or the provided threshold value is too small due to floating point imprecisions (:gh:`xxx` by `Richard Höchenberger`_) + - `mne.preprocessing.compute_proj_eog` and `mne.preprocessing.compute_proj_ecg` now return empty lists if no EOG or ECG events, respectively, could be found. Previously, we'd return ``None`` in these situations, which does not match the documented behavior of returning a list of projectors (:gh:`9277` by `Richard Höchenberger`_) API changes diff --git a/doc/conf.py b/doc/conf.py index de86ec90dfe..18422dbcf04 100644 --- a/doc/conf.py +++ b/doc/conf.py @@ -793,6 +793,7 @@ def reset_warnings(gallery_conf, fname): 'ignore', message='.*mne-realtime.*', category=DeprecationWarning) + reset_warnings(None, None) diff --git a/mne/preprocessing/_peak_finder.py b/mne/preprocessing/_peak_finder.py index 9bac177988a..3335b075a9f 100644 --- a/mne/preprocessing/_peak_finder.py +++ b/mne/preprocessing/_peak_finder.py @@ -50,11 +50,26 @@ def peak_finder(x0, thresh=None, extrema=1, verbose=None): s = x0.size if x0.ndim >= 2 or s == 0: - raise ValueError('The input data must be a non empty 1D vector') + raise ValueError('The input data must be a non-empty 1D vector') + + if x0.dtype.kind in 'uib': # integer dtype + x0 = x0.astype(float) + + eps = 1e2 * np.finfo(x0.dtype).eps + + if thresh is not None: + assert thresh >= 0. if thresh is None: thresh = (np.max(x0) - np.min(x0)) / 4 + if thresh < eps: + thresh = eps logger.debug('Peak finder automatic threshold: %0.2g' % (thresh,)) + elif thresh < eps: + raise ValueError( + f'The selected threshold value is too small for peak detection ' + f'due to inevitable floating-point imprecisions.\nPlease pick a ' + f'value >= {eps:g} (You passed thresh={thresh})') assert extrema in [-1, 1] @@ -160,10 +175,6 @@ def peak_finder(x0, thresh=None, extrema=1, verbose=None): peak_mags = [] peak_inds = [] - # Change sign of data if was finding minima - if extrema < 0: - peak_mags *= -1.0 - # ensure output type array if not isinstance(peak_inds, np.ndarray): peak_inds = np.atleast_1d(peak_inds).astype('int64') @@ -171,8 +182,11 @@ def peak_finder(x0, thresh=None, extrema=1, verbose=None): if not isinstance(peak_mags, np.ndarray): peak_mags = np.atleast_1d(peak_mags).astype('float64') - # Plot if no output desired - if len(peak_inds) == 0: + # Change sign of data if was finding minima + if extrema < 0 and peak_mags.size > 0: + peak_mags *= -1.0 + + if peak_inds.size == 0: logger.info('No significant peaks found') else: logger.info('Found %d significant peak%s' diff --git a/mne/preprocessing/ecg.py b/mne/preprocessing/ecg.py index 0ed1e0e52e7..3559c990f72 100644 --- a/mne/preprocessing/ecg.py +++ b/mne/preprocessing/ecg.py @@ -31,7 +31,7 @@ def qrs_detector(sfreq, ecg, thresh_value=0.6, levels=2.5, n_thresh=3, Sampling rate ecg : array ECG signal - thresh_value : float | str + thresh_value : float | 'auto' qrs detection threshold. Can also be "auto" for automatic selection of threshold. levels : float @@ -65,18 +65,32 @@ def qrs_detector(sfreq, ecg, thresh_value=0.6, levels=2.5, n_thresh=3, n_points = len(ecg_abs) + # get peak amplitudes in the 1st, 2nd, and 3rd second of the recording maxpt = np.empty(3) maxpt[0] = np.max(ecg_abs[:init]) maxpt[1] = np.max(ecg_abs[init:init * 2]) maxpt[2] = np.max(ecg_abs[init * 2:init * 3]) init_max = np.mean(maxpt) + eps = 1e2 * np.finfo(ecg_abs.dtype).eps + + if init_max < eps: + warn('The first 3 seconds of your ECG data appear to be flat. Please ' + 'inspect your data to ensure this is intentional.') + init_max = eps if thresh_value == 'auto': thresh_runs = np.arange(0.3, 1.1, 0.05) elif isinstance(thresh_value, str): raise ValueError('threshold value must be "auto" or a float') else: + if init_max * thresh_value < eps: + raise ValueError( + f'The selected threshold value is too small for peak ' + f'detection with the provided data, due to inevitable ' + f'floating-point imprecisions.\n' + f'Please pick a value >= {eps/init_max:g} ' + f'(You passed thresh={thresh_value})') thresh_runs = [thresh_value] # Try a few thresholds (or just one) @@ -100,7 +114,7 @@ def qrs_detector(sfreq, ecg, thresh_value=0.6, levels=2.5, n_thresh=3, else: ii += 1 - if len(rms) == 0: + if not rms: rms.append(0.0) time.append(0.0) time = np.array(time) @@ -151,7 +165,7 @@ def find_ecg_events(raw, event_id=999, ch_name=None, tstart=0.0, %(ecg_ch_name)s %(ecg_tstart)s %(ecg_filter_freqs)s - qrs_threshold : float | str + qrs_threshold : float | 'auto' Between 0 and 1. qrs detection threshold. Can also be "auto" to automatically choose the threshold that generates a reasonable number of heartbeats (40-160 beats / min). diff --git a/mne/preprocessing/ica.py b/mne/preprocessing/ica.py index 4d8241ad79d..29e22473140 100644 --- a/mne/preprocessing/ica.py +++ b/mne/preprocessing/ica.py @@ -1246,7 +1246,7 @@ def find_bads_ecg(self, inst, ch_name=None, threshold='auto', start=None, The name of the channel to use for ECG peak detection. The argument is mandatory if the dataset contains no ECG channels. - threshold : float | str + threshold : float | 'auto' The value above which a feature is classified as outlier. If 'auto' and method is 'ctps', automatically compute the threshold. If 'auto' and method is 'correlation', defaults to 3.0. The default diff --git a/mne/preprocessing/tests/test_ecg.py b/mne/preprocessing/tests/test_ecg.py index 75148098b16..fa18e66f0fc 100644 --- a/mne/preprocessing/tests/test_ecg.py +++ b/mne/preprocessing/tests/test_ecg.py @@ -86,12 +86,20 @@ def test_find_ecg(): assert 'ECG-SYN' not in raw.ch_names assert 'ECG-SYN' not in ecg_epochs.ch_names - # Test behavior if no peaks can be found -> achieve this by providing - # all-zero'd data - raw._data[ecg_idx] = 0. + # test behavior if no peaks can be found -> achieve this by setting a high + # threshold ecg_events, _, average_pulse, ecg = find_ecg_events( - raw, ch_name=raw.ch_names[ecg_idx], return_ecg=True + raw, ch_name=raw.ch_names[ecg_idx], return_ecg=True, qrs_threshold=1e5 ) assert ecg_events.size == 0 assert average_pulse == 0 + + # test behavior for an entirely flat channel + raw._data[ecg_idx] = 0. + with pytest.warns(RuntimeWarning, match='first 3 seconds.*flat'): + ecg_events, _, average_pulse, ecg = find_ecg_events( + raw, ch_name=raw.ch_names[ecg_idx], return_ecg=True + ) + assert ecg_events.size == 0 + assert average_pulse == 0 assert np.allclose(ecg, np.zeros_like(ecg)) diff --git a/mne/preprocessing/tests/test_eog.py b/mne/preprocessing/tests/test_eog.py index f0435c42878..8eb05ae4300 100644 --- a/mne/preprocessing/tests/test_eog.py +++ b/mne/preprocessing/tests/test_eog.py @@ -1,5 +1,6 @@ import os.path as op +import numpy as np import pytest from mne import Annotations @@ -39,3 +40,14 @@ def test_find_eog(): events = find_eog_events(raw, ch_name=['EEG 060', 'EOG 061']) assert len(events) == 4 + + # test too low threshold + eps = np.finfo(raw.get_data().dtype).eps + with pytest.raises(ValueError, match='threshold value is too small'): + find_eog_events(raw, ch_name=None, thresh=eps) + + # test with an entirely flat signal + eog = raw.copy().pick_types(meg=False, eog=True).load_data() + eog._data[:] = 0. + events = find_eog_events(eog, ch_name=None) + assert events.size == 0 diff --git a/mne/preprocessing/tests/test_ssp.py b/mne/preprocessing/tests/test_ssp.py index 4307dc9f6ea..934084717a4 100644 --- a/mne/preprocessing/tests/test_ssp.py +++ b/mne/preprocessing/tests/test_ssp.py @@ -63,11 +63,7 @@ def test_compute_proj_ecg(short_raw, average): projs, events, drop_log = compute_proj_ecg( raw, n_mag=2, n_grad=2, n_eeg=2, ch_name='MEG 1531', bads=[], average=average, avg_ref=True, no_proj=True, l_freq=None, - h_freq=None, tmax=dur_use, return_drop_log=True, - # XXX can be removed once - # XXX https://github.com/mne-tools/mne-python/issues/9273 - # XXX has been resolved: - qrs_threshold=1e-15) + h_freq=None, tmax=dur_use, return_drop_log=True) assert projs == [] assert len(events) == len(drop_log)