diff --git a/.codeclimate.yml b/.codeclimate.yml index a60dd55b42..c4946e57d8 100644 --- a/.codeclimate.yml +++ b/.codeclimate.yml @@ -1,43 +1,43 @@ engines: - pep8: - enabled: true - checks: - E501: + pep8: + enabled: true + checks: + E501: + enabled: false + E203: + enabled: false + E302: + enabled: false + E303: + enabled: false + +checks: + argument-count: enabled: false - E203: + complex-logic: enabled: false - E302: + file-lines: enabled: false - E303: + method-complexity: + enabled: false + method-count: + enabled: false + method-lines: + enabled: false + nested-control-flow: + enabled: true + return-statements: + enabled: true + similar-code: + enabled: false + identical-code: enabled: false - -checks: - argument-count: - enabled: false - complex-logic: - enabled: false - file-lines: - enabled: true - method-complexity: - enabled: false - method-count: - enabled: true - method-lines: - enabled: false - nested-control-flow: - enabled: true - return-statements: - enabled: true - similar-code: - enabled: false - identical-code: - enabled: false exclude_patterns: -- "tests/" -- "scripts/" -- "benchmarks/" -- "studies/" -- "data/" -- "docs/img/" -- "docs/readme/" + - "tests/" + - "scripts/" + - "benchmarks/" + - "studies/" + - "data/" + - "docs/img/" + - "docs/readme/" diff --git a/.github/workflows/docs-build.yml b/.github/workflows/docs-build.yml index 96d26c165c..99b4f636fc 100644 --- a/.github/workflows/docs-build.yml +++ b/.github/workflows/docs-build.yml @@ -27,10 +27,11 @@ jobs: pip install sphinx-book-theme pip install sphinxemoji pip install sphinx-copybutton + pip install ipykernel==6.23.2 pip install ipython pip install myst-parser pip install myst-nb - pip install numpy + pip install numpy==1.24.3 pip install pandas pip install scipy pip install scikit-learn diff --git a/.github/workflows/docs-check.yml b/.github/workflows/docs-check.yml index db2d3c4782..b7026c7033 100644 --- a/.github/workflows/docs-check.yml +++ b/.github/workflows/docs-check.yml @@ -22,10 +22,11 @@ jobs: pip install sphinx-book-theme pip install sphinxemoji pip install sphinx-copybutton + pip install ipykernel==6.23.2 pip install ipython pip install myst-parser pip install myst-nb - pip install numpy + pip install numpy==1.24.3 pip install pandas pip install scipy pip install scikit-learn diff --git a/docs/conf.py b/docs/conf.py index af75667886..99d45d4142 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -12,6 +12,8 @@ import os import re import sys +import asyncio +import platform # -- Path setup -------------------------------------------------------------- @@ -105,6 +107,9 @@ def find_version(): # googleanalytics_id = "G-DVXSEGN5M9" +# Address asyncio warning +if platform.system() == "Windows": + asyncio.set_event_loop_policy(asyncio.WindowsSelectorEventLoopPolicy()) # NumPyDoc configuration ----------------------------------------------------- diff --git a/docs/functions/signal.rst b/docs/functions/signal.rst index 6be83fa5c1..49848be29b 100644 --- a/docs/functions/signal.rst +++ b/docs/functions/signal.rst @@ -19,6 +19,10 @@ Preprocessing """"""""""""""""""""" .. autofunction:: neurokit2.signal_resample +*signal_fillmissing()* +"""""""""""""""""""" +.. autofunction:: neurokit2.signal_fillmissing + Transformation ^^^^^^^^^^^^^^^^ diff --git a/neurokit2/__init__.py b/neurokit2/__init__.py index 7fa6113c27..2aff8433e9 100644 --- a/neurokit2/__init__.py +++ b/neurokit2/__init__.py @@ -32,7 +32,7 @@ from .video import * # Info -__version__ = "0.2.4" +__version__ = "0.2.5" # Maintainer info diff --git a/neurokit2/complexity/complexity.py b/neurokit2/complexity/complexity.py index 9aae8ad880..649cbc5d68 100644 --- a/neurokit2/complexity/complexity.py +++ b/neurokit2/complexity/complexity.py @@ -37,10 +37,10 @@ def complexity(signal, which="makowski2022", delay=1, dimension=2, tolerance="sd depending on what the literature suggests. We recommend using this function only for quick exploratory analyses, but then replacing it by the calls to the individual functions. - Check-out our open-access `study `_ + Check-out our `open-access study `_ explaining the selection of indices. - The categorization by "computation time" is based on our `study + The categorization by "computation time" is based on `our study `_ results: .. figure:: https://raw.githubusercontent.com/DominiqueMakowski/ComplexityStructure/main/figures/time1-1.png diff --git a/neurokit2/complexity/fractal_tmf.py b/neurokit2/complexity/fractal_tmf.py index 46e4118e4f..17a868ecde 100644 --- a/neurokit2/complexity/fractal_tmf.py +++ b/neurokit2/complexity/fractal_tmf.py @@ -109,7 +109,7 @@ def fractal_tmf(signal, n=40, show=False, **kwargs): w = np.zeros(n) for i in range(n): surro = signal_surrogate(signal, method="IAAFT") - w[i] = fractal_dfa(surro, multifractal=True, show=False)[0]["Width"] + w[i] = float(fractal_dfa(surro, multifractal=True, show=False)[0]["Width"].iloc[0]) # Run t-test # TODO: adjust in the future diff --git a/neurokit2/complexity/optim_complexity_tolerance.py b/neurokit2/complexity/optim_complexity_tolerance.py index a6a232b856..aa10d4607d 100644 --- a/neurokit2/complexity/optim_complexity_tolerance.py +++ b/neurokit2/complexity/optim_complexity_tolerance.py @@ -157,7 +157,7 @@ def complexity_tolerance( .. ipython:: python - @savefig p_complexity_tolerance4.png scale=100% + @savefig p_complexity_tolerance5.png scale=100% r, info = nk.complexity_tolerance(signal, delay=1, dimension=3, method = 'bin', show=True) @suppress @@ -168,7 +168,7 @@ def complexity_tolerance( .. ipython:: python # Slow method - @savefig p_complexity_tolerance5.png scale=100% + @savefig p_complexity_tolerance6.png scale=100% r, info = nk.complexity_tolerance(signal, delay=8, dimension=6, method = 'maxApEn', show=True) @suppress @@ -184,7 +184,7 @@ def complexity_tolerance( .. ipython:: python # Narrower range - @savefig p_complexity_tolerance6.png scale=100% + @savefig p_complexity_tolerance7.png scale=100% r, info = nk.complexity_tolerance(signal, delay=8, dimension=6, method = 'maxApEn', r_range=np.linspace(0.002, 0.8, 30), show=True) @suppress diff --git a/neurokit2/complexity/utils.py b/neurokit2/complexity/utils.py index c6af29d79a..f3801bfb09 100644 --- a/neurokit2/complexity/utils.py +++ b/neurokit2/complexity/utils.py @@ -2,6 +2,7 @@ import numpy as np import sklearn.metrics import sklearn.neighbors +from packaging import version from .utils_complexity_embedding import complexity_embedding @@ -109,11 +110,16 @@ def _get_count( # Get neighbors count # ------------------- # Sanity checks - if distance not in sklearn.neighbors.KDTree.valid_metrics + ["range"]: + sklearn_version = version.parse(sklearn.__version__) + if sklearn_version >= version.parse("1.3.0"): + valid_metrics = sklearn.neighbors.KDTree.valid_metrics() + ["range"] + else: + valid_metrics = sklearn.neighbors.KDTree.valid_metrics + ["range"] + if distance not in valid_metrics: raise ValueError( "The given metric (%s) is not valid." "The valid metric names are: %s" - % (distance, sklearn.neighbors.KDTree.valid_metrics + ["range"]) + % (distance, valid_metrics) ) if fuzzy is True: diff --git a/neurokit2/complexity/utils_fractal_mandelbrot.py b/neurokit2/complexity/utils_fractal_mandelbrot.py index f83bfffd72..6893dc59de 100644 --- a/neurokit2/complexity/utils_fractal_mandelbrot.py +++ b/neurokit2/complexity/utils_fractal_mandelbrot.py @@ -260,7 +260,7 @@ def _buddhabrot_initialize(size=1000, iterations=100, real_range=(-2, 2), imagin def _mandelbrot_optimize(c): # Optimizations: most of the mset points lie within the # within the cardioid or in the period-2 bulb. (The two most - # prominant shapes in the mandelbrot set. We can eliminate these + # prominent shapes in the mandelbrot set. We can eliminate these # from our search straight away and save alot of time. # see: http://en.wikipedia.org/wiki/Mandelbrot_set#Optimizations diff --git a/neurokit2/ecg/ecg_clean.py b/neurokit2/ecg/ecg_clean.py index c42acd4452..e27f80deed 100644 --- a/neurokit2/ecg/ecg_clean.py +++ b/neurokit2/ecg/ecg_clean.py @@ -17,7 +17,9 @@ def ecg_clean(ecg_signal, sampling_rate=1000, method="neurokit", **kwargs): * ``'neurokit'`` (default): 0.5 Hz high-pass butterworth filter (order = 5), followed by powerline filtering (see ``signal_filter()``). By default, ``powerline = 50``. - * ``'biosppy'``: Same as in the biosppy package. **Please help providing a better description!** + * ``'biosppy'``: Method used in the BioSPPy package. A FIR filter ([0.67, 45] Hz; order = 1.5 * + SR). The 0.67 Hz cutoff value was selected based on the fact that there are no morphological + features below the heartrate (assuming a minimum heart rate of 40 bpm). * ``'pantompkins1985'``: Method used in Pan & Tompkins (1985). **Please help providing a better description!** * ``'hamilton2002'``: Method used in Hamilton (2002). **Please help providing a better @@ -121,7 +123,7 @@ def ecg_clean(ecg_signal, sampling_rate=1000, method="neurokit", **kwargs): "swt", "kalidas", "kalidastamil", - "kalidastamil2017" + "kalidastamil2017", ]: clean = ecg_signal else: @@ -137,7 +139,6 @@ def ecg_clean(ecg_signal, sampling_rate=1000, method="neurokit", **kwargs): # Handle missing data # ============================================================================= def _ecg_clean_missing(ecg_signal): - ecg_signal = pd.DataFrame.pad(pd.Series(ecg_signal)) return ecg_signal @@ -147,11 +148,18 @@ def _ecg_clean_missing(ecg_signal): # NeuroKit # ============================================================================= def _ecg_clean_nk(ecg_signal, sampling_rate=1000, **kwargs): - # Remove slow drift and dc offset with highpass Butterworth. - clean = signal_filter(signal=ecg_signal, sampling_rate=sampling_rate, lowcut=0.5, method="butterworth", order=5) + clean = signal_filter( + signal=ecg_signal, + sampling_rate=sampling_rate, + lowcut=0.5, + method="butterworth", + order=5, + ) - clean = signal_filter(signal=clean, sampling_rate=sampling_rate, method="powerline", **kwargs) + clean = signal_filter( + signal=clean, sampling_rate=sampling_rate, method="powerline", **kwargs + ) return clean @@ -160,18 +168,24 @@ def _ecg_clean_nk(ecg_signal, sampling_rate=1000, **kwargs): # ============================================================================= def _ecg_clean_biosppy(ecg_signal, sampling_rate=1000): """Adapted from https://github.com/PIA- - Group/BioSPPy/blob/e65da30f6379852ecb98f8e2e0c9b4b5175416c3/biosppy/signals/ecg.py#L69.""" + Group/BioSPPy/blob/e65da30f6379852ecb98f8e2e0c9b4b5175416c3/biosppy/signals/ecg.py#L69. + """ - order = int(0.3 * sampling_rate) + # The order and frequency was recently changed + # (see https://github.com/scientisst/BioSPPy/pull/12) + + order = int(1.5 * sampling_rate) if order % 2 == 0: order += 1 # Enforce odd number # -> filter_signal() - frequency = [3, 45] + frequency = [0.67, 45] # -> get_filter() # -> _norm_freq() - frequency = 2 * np.array(frequency) / sampling_rate # Normalize frequency to Nyquist Frequency (Fs/2). + frequency = ( + 2 * np.array(frequency) / sampling_rate + ) # Normalize frequency to Nyquist Frequency (Fs/2). # -> get coeffs a = np.array([1]) @@ -180,6 +194,9 @@ def _ecg_clean_biosppy(ecg_signal, sampling_rate=1000): # _filter_signal() filtered = scipy.signal.filtfilt(b, a, ecg_signal) + # DC offset + filtered -= np.mean(filtered) + return filtered @@ -188,11 +205,17 @@ def _ecg_clean_biosppy(ecg_signal, sampling_rate=1000): # ============================================================================= def _ecg_clean_pantompkins(ecg_signal, sampling_rate=1000): """Adapted from https://github.com/PIA- - Group/BioSPPy/blob/e65da30f6379852ecb98f8e2e0c9b4b5175416c3/biosppy/signals/ecg.py#L69.""" + Group/BioSPPy/blob/e65da30f6379852ecb98f8e2e0c9b4b5175416c3/biosppy/signals/ecg.py#L69. + """ order = 1 clean = signal_filter( - signal=ecg_signal, sampling_rate=sampling_rate, lowcut=5, highcut=15, method="butterworth_zi", order=order + signal=ecg_signal, + sampling_rate=sampling_rate, + lowcut=5, + highcut=15, + method="butterworth_zi", + order=order, ) return clean # Return filtered @@ -212,7 +235,12 @@ def _ecg_clean_elgendi(ecg_signal, sampling_rate=1000): order = 2 clean = signal_filter( - signal=ecg_signal, sampling_rate=sampling_rate, lowcut=8, highcut=20, method="butterworth_zi", order=order + signal=ecg_signal, + sampling_rate=sampling_rate, + lowcut=8, + highcut=20, + method="butterworth_zi", + order=order, ) return clean # Return filtered @@ -223,11 +251,17 @@ def _ecg_clean_elgendi(ecg_signal, sampling_rate=1000): # ============================================================================= def _ecg_clean_hamilton(ecg_signal, sampling_rate=1000): """Adapted from https://github.com/PIA- - Group/BioSPPy/blob/e65da30f6379852ecb98f8e2e0c9b4b5175416c3/biosppy/signals/ecg.py#L69.""" + Group/BioSPPy/blob/e65da30f6379852ecb98f8e2e0c9b4b5175416c3/biosppy/signals/ecg.py#L69. + """ order = 1 clean = signal_filter( - signal=ecg_signal, sampling_rate=sampling_rate, lowcut=8, highcut=16, method="butterworth_zi", order=order + signal=ecg_signal, + sampling_rate=sampling_rate, + lowcut=8, + highcut=16, + method="butterworth_zi", + order=order, ) return clean # Return filtered @@ -249,7 +283,12 @@ def _ecg_clean_engzee(ecg_signal, sampling_rate=1000): order = 4 clean = signal_filter( - signal=ecg_signal, sampling_rate=sampling_rate, lowcut=52, highcut=48, method="butterworth_zi", order=order + signal=ecg_signal, + sampling_rate=sampling_rate, + lowcut=52, + highcut=48, + method="butterworth_zi", + order=order, ) return clean # Return filtered @@ -270,6 +309,12 @@ def _ecg_clean_vgraph(ecg_signal, sampling_rate=1000): """ order = 2 - clean = signal_filter(signal=ecg_signal, sampling_rate=sampling_rate, lowcut=4, method="butterworth", order=order) + clean = signal_filter( + signal=ecg_signal, + sampling_rate=sampling_rate, + lowcut=4, + method="butterworth", + order=order, + ) return clean # Return filtered diff --git a/neurokit2/ecg/ecg_delineate.py b/neurokit2/ecg/ecg_delineate.py index aee288c5be..24218f82a6 100644 --- a/neurokit2/ecg/ecg_delineate.py +++ b/neurokit2/ecg/ecg_delineate.py @@ -153,15 +153,20 @@ def ecg_delineate( method = method.lower() # remove capitalised letters if method in ["peak", "peaks", "derivative", "gradient"]: - waves = _ecg_delineator_peak(ecg_cleaned, rpeaks=rpeaks, sampling_rate=sampling_rate) + waves = _ecg_delineator_peak( + ecg_cleaned, rpeaks=rpeaks, sampling_rate=sampling_rate + ) elif method in ["cwt", "continuous wavelet transform"]: - waves = _ecg_delineator_cwt(ecg_cleaned, rpeaks=rpeaks, sampling_rate=sampling_rate) + waves = _ecg_delineator_cwt( + ecg_cleaned, rpeaks=rpeaks, sampling_rate=sampling_rate + ) elif method in ["dwt", "discrete wavelet transform"]: waves = _dwt_ecg_delineator(ecg_cleaned, rpeaks, sampling_rate=sampling_rate) else: raise ValueError( - "NeuroKit error: ecg_delineate(): 'method' should be one of 'peak'," "'cwt' or 'dwt'." + "NeuroKit error: ecg_delineate(): 'method' should be one of 'peak'," + "'cwt' or 'dwt'." ) # Ensure that all indices are not larger than ECG signal indices @@ -172,7 +177,9 @@ def ecg_delineate( # Remove NaN in Peaks, Onsets, and Offsets waves_noNA = waves.copy() for feature in waves_noNA.keys(): - waves_noNA[feature] = [int(x) for x in waves_noNA[feature] if ~np.isnan(x) and x > 0] + waves_noNA[feature] = [ + int(x) for x in waves_noNA[feature] if ~np.isnan(x) and x > 0 + ] instant_peaks = signal_formatpeaks(waves_noNA, desired_length=len(ecg_cleaned)) signals = instant_peaks @@ -202,14 +209,18 @@ def ecg_delineate( # ============================================================================= def _dwt_resample_points(peaks, sampling_rate, desired_sampling_rate): """Resample given points to a different sampling rate.""" - if isinstance(peaks, np.ndarray): # peaks are passed in from previous processing steps + if isinstance( + peaks, np.ndarray + ): # peaks are passed in from previous processing steps # Prevent overflow by converting to np.int64 (peaks might be passed in containing np.int32). peaks = peaks.astype(dtype=np.int64) elif isinstance(peaks, list): # peaks returned from internal functions # Cannot be converted to int since list might contain np.nan. Automatically cast to np.float64 if list contains np.nan. peaks = np.array(peaks) peaks_resample = peaks * desired_sampling_rate / sampling_rate - peaks_resample = [np.nan if np.isnan(x) else int(x) for x in peaks_resample.tolist()] + peaks_resample = [ + np.nan if np.isnan(x) else int(x) for x in peaks_resample.tolist() + ] return peaks_resample @@ -242,7 +253,9 @@ def _dwt_ecg_delineator(ecg, rpeaks, sampling_rate, analysis_sampling_rate=2000) for i, rpeak in enumerate(rpeaks): heartbeat = heartbeats[str(i + 1)] # Get index of R peaks - R = heartbeat.index.get_loc(np.min(heartbeat.index.values[heartbeat.index.values > 0])) + R = heartbeat.index.get_loc( + np.min(heartbeat.index.values[heartbeat.index.values > 0]) + ) # Q wave Q_index, Q = _ecg_delineator_peak_Q(rpeak, heartbeat, R) qpeaks.append(Q_index) @@ -263,13 +276,23 @@ def _dwt_ecg_delineator(ecg, rpeaks, sampling_rate, analysis_sampling_rate=2000) # plt.legend() # plt.grid(True) # plt.show() - rpeaks_resampled = _dwt_resample_points(rpeaks, sampling_rate, analysis_sampling_rate) + rpeaks_resampled = _dwt_resample_points( + rpeaks, sampling_rate, analysis_sampling_rate + ) + qpeaks_resampled = _dwt_resample_points( + qpeaks, sampling_rate, analysis_sampling_rate + ) tpeaks, ppeaks = _dwt_delineate_tp_peaks( ecg, rpeaks_resampled, dwtmatr, sampling_rate=analysis_sampling_rate ) qrs_onsets, qrs_offsets = _dwt_delineate_qrs_bounds( - rpeaks_resampled, dwtmatr, ppeaks, tpeaks, sampling_rate=analysis_sampling_rate + rpeaks_resampled, + dwtmatr, + ppeaks, + tpeaks, + qpeaks_resampled, + sampling_rate=analysis_sampling_rate, ) ponsets, poffsets = _dwt_delineate_tp_onsets_offsets( ppeaks, rpeaks_resampled, dwtmatr, sampling_rate=analysis_sampling_rate @@ -383,12 +406,11 @@ def _dwt_delineate_tp_peaks( srch_idx_start = rpeak_ + srch_bndry srch_idx_end = rpeak_ + 2 * int(rt_duration * sampling_rate) dwt_local = dwtmatr[degree_tpeak + degree_add, srch_idx_start:srch_idx_end] - height = epsilon_T_weight * np.sqrt(np.mean(np.square(dwt_local))) - if len(dwt_local) == 0: tpeaks.append(np.nan) continue + height = epsilon_T_weight * np.sqrt(np.mean(np.square(dwt_local))) ecg_local = ecg[srch_idx_start:srch_idx_end] peaks, __ = scipy.signal.find_peaks(np.abs(dwt_local), height=height) peaks = list( @@ -406,7 +428,8 @@ def _dwt_delineate_tp_peaks( ) # pylint: disable=R1716 if correct_sign: idx_zero = ( - signal_zerocrossings(dwt_local[idx_peak : idx_peak_nxt + 1])[0] + idx_peak + signal_zerocrossings(dwt_local[idx_peak : idx_peak_nxt + 1])[0] + + idx_peak ) # This is the score assigned to each peak. The peak with the highest score will be # selected. @@ -420,7 +443,9 @@ def _dwt_delineate_tp_peaks( tpeaks.append(np.nan) continue - tpeaks.append(candidate_peaks[np.argmax(candidate_peaks_scores)] + srch_idx_start) + tpeaks.append( + candidate_peaks[np.argmax(candidate_peaks_scores)] + srch_idx_start + ) ppeaks = [] for rpeak in rpeaks: @@ -432,15 +457,16 @@ def _dwt_delineate_tp_peaks( srch_idx_start = rpeak - 2 * int(p2r_duration * sampling_rate) srch_idx_end = rpeak - srch_bndry dwt_local = dwtmatr[degree_ppeak + degree_add, srch_idx_start:srch_idx_end] - height = epsilon_P_weight * np.sqrt(np.mean(np.square(dwt_local))) - if len(dwt_local) == 0: ppeaks.append(np.nan) continue + height = epsilon_P_weight * np.sqrt(np.mean(np.square(dwt_local))) ecg_local = ecg[srch_idx_start:srch_idx_end] peaks, __ = scipy.signal.find_peaks(np.abs(dwt_local), height=height) - peaks = list(filter(lambda p: np.abs(dwt_local[p]) > 0.025 * max(dwt_local), peaks)) + peaks = list( + filter(lambda p: np.abs(dwt_local[p]) > 0.025 * max(dwt_local), peaks) + ) if dwt_local[0] > 0: # just append peaks = [0] + peaks @@ -453,7 +479,8 @@ def _dwt_delineate_tp_peaks( ) # pylint: disable=R1716 if correct_sign: idx_zero = ( - signal_zerocrossings(dwt_local[idx_peak : idx_peak_nxt + 1])[0] + idx_peak + signal_zerocrossings(dwt_local[idx_peak : idx_peak_nxt + 1])[0] + + idx_peak ) # This is the score assigned to each peak. The peak with the highest score will be # selected. @@ -467,7 +494,9 @@ def _dwt_delineate_tp_peaks( ppeaks.append(np.nan) continue - ppeaks.append(candidate_peaks[np.argmax(candidate_peaks_scores)] + srch_idx_start) + ppeaks.append( + candidate_peaks[np.argmax(candidate_peaks_scores)] + srch_idx_start + ) return tpeaks, ppeaks @@ -510,7 +539,9 @@ def _dwt_delineate_tp_onsets_offsets( if not (dwt_local[: onset_slope_peaks[-1]] < epsilon_onset).any(): onsets.append(np.nan) continue - candidate_onsets = np.where(dwt_local[: onset_slope_peaks[-1]] < epsilon_onset)[0] + candidate_onsets = np.where(dwt_local[: onset_slope_peaks[-1]] < epsilon_onset)[ + 0 + ] onsets.append(candidate_onsets[-1] + srch_idx_start) # # only for debugging @@ -548,13 +579,15 @@ def _dwt_delineate_tp_onsets_offsets( return onsets, offsets -def _dwt_delineate_qrs_bounds(rpeaks, dwtmatr, ppeaks, tpeaks, sampling_rate=250): +def _dwt_delineate_qrs_bounds( + rpeaks, dwtmatr, ppeaks, tpeaks, qpeaks, sampling_rate=250 +): degree = _dwt_adjust_parameters(rpeaks, sampling_rate, target="degree") onsets = [] - for i in range(len(rpeaks)): # pylint: disable=C0200 + for i in range(len(qpeaks)): # pylint: disable=C0200 # look for onsets srch_idx_start = ppeaks[i] - srch_idx_end = rpeaks[i] + srch_idx_end = qpeaks[i] if srch_idx_start is np.nan or srch_idx_end is np.nan: onsets.append(np.nan) continue @@ -567,7 +600,9 @@ def _dwt_delineate_qrs_bounds(rpeaks, dwtmatr, ppeaks, tpeaks, sampling_rate=250 if not (-dwt_local[: onset_slope_peaks[-1]] < epsilon_onset).any(): onsets.append(np.nan) continue - candidate_onsets = np.where(-dwt_local[: onset_slope_peaks[-1]] < epsilon_onset)[0] + candidate_onsets = np.where( + -dwt_local[: onset_slope_peaks[-1]] < epsilon_onset + )[0] onsets.append(candidate_onsets[-1] + srch_idx_start) # only for debugging @@ -595,7 +630,8 @@ def _dwt_delineate_qrs_bounds(rpeaks, dwtmatr, ppeaks, tpeaks, sampling_rate=250 offsets.append(np.nan) continue candidate_offsets = ( - np.where(dwt_local[onset_slope_peaks[0] :] < epsilon_offset)[0] + onset_slope_peaks[0] + np.where(dwt_local[onset_slope_peaks[0] :] < epsilon_offset)[0] + + onset_slope_peaks[0] ) offsets.append(candidate_offsets[0] + srch_idx_start) @@ -642,7 +678,9 @@ def _apply_G_filter(signal_i, power=0): T_deg = _apply_H_filter(intermediate_ret, power=deg) dwtmatr.append(S_deg) intermediate_ret = np.array(T_deg) - dwtmatr = [arr[: len(ecg)] for arr in dwtmatr] # rescale transforms to the same length + dwtmatr = [ + arr[: len(ecg)] for arr in dwtmatr + ] # rescale transforms to the same length return np.array(dwtmatr) @@ -650,7 +688,6 @@ def _apply_G_filter(signal_i, power=0): # WAVELET METHOD (CWT) # ============================================================================= def _ecg_delineator_cwt(ecg, rpeaks=None, sampling_rate=1000): - # P-Peaks and T-Peaks tpeaks, ppeaks = _peaks_delineator(ecg, rpeaks, sampling_rate=sampling_rate) @@ -677,7 +714,9 @@ def _ecg_delineator_cwt(ecg, rpeaks=None, sampling_rate=1000): for i, rpeak in enumerate(rpeaks): heartbeat = heartbeats[str(i + 1)] # Get index of R peaks - R = heartbeat.index.get_loc(np.min(heartbeat.index.values[heartbeat.index.values > 0])) + R = heartbeat.index.get_loc( + np.min(heartbeat.index.values[heartbeat.index.values > 0]) + ) # Q wave Q_index, Q = _ecg_delineator_peak_Q(rpeak, heartbeat, R) q_peaks.append(Q_index) @@ -760,11 +799,15 @@ def _onset_offset_delineator(ecg, peaks, peak_type="rpeaks", sampling_rate=1000) leftbase = wt_peaks_data["left_bases"][-1] + index_peak - half_wave_width if peak_type == "rpeaks": candidate_onsets = ( - np.where(cwtmatr[2, nfirst - 100 : nfirst] < epsilon_onset)[0] + nfirst - 100 + np.where(cwtmatr[2, nfirst - 100 : nfirst] < epsilon_onset)[0] + + nfirst + - 100 ) elif peak_type in ["tpeaks", "ppeaks"]: candidate_onsets = ( - np.where(-cwtmatr[4, nfirst - 100 : nfirst] < epsilon_onset)[0] + nfirst - 100 + np.where(-cwtmatr[4, nfirst - 100 : nfirst] < epsilon_onset)[0] + + nfirst + - 100 ) candidate_onsets = candidate_onsets.tolist() + [leftbase] @@ -804,11 +847,13 @@ def _onset_offset_delineator(ecg, peaks, peak_type="rpeaks", sampling_rate=1000) rightbase = wt_peaks_data["right_bases"][0] + index_peak if peak_type == "rpeaks": candidate_offsets = ( - np.where((-cwtmatr[2, nlast : nlast + 100]) < epsilon_offset)[0] + nlast + np.where((-cwtmatr[2, nlast : nlast + 100]) < epsilon_offset)[0] + + nlast ) elif peak_type in ["tpeaks", "ppeaks"]: candidate_offsets = ( - np.where((cwtmatr[4, nlast : nlast + 100]) < epsilon_offset)[0] + nlast + np.where((cwtmatr[4, nlast : nlast + 100]) < epsilon_offset)[0] + + nlast ) candidate_offsets = candidate_offsets.tolist() + [rightbase] @@ -845,13 +890,17 @@ def _peaks_delineator(ecg, rpeaks, sampling_rate=1000): end = rpeaks[i + 1] - search_boundary search_window = cwtmatr[4, start:end] height = 0.25 * np.sqrt(np.mean(np.square(search_window))) - peaks_tp, heights_tp = scipy.signal.find_peaks(np.abs(search_window), height=height) + peaks_tp, heights_tp = scipy.signal.find_peaks( + np.abs(search_window), height=height + ) peaks_tp = peaks_tp + rpeaks[i] + search_boundary # set threshold for heights of peaks to find significant peaks in wavelet threshold = 0.125 * max(search_window) significant_peaks_tp = [] significant_peaks_tp = [ - peaks_tp[j] for j in range(len(peaks_tp)) if heights_tp["peak_heights"][j] > threshold + peaks_tp[j] + for j in range(len(peaks_tp)) + if heights_tp["peak_heights"][j] > threshold ] significant_peaks_groups.append( @@ -888,12 +937,13 @@ def _find_tppeaks(ecg, keep_tp, sampling_rate=1000): # if near and correct_sign: if correct_sign: index_zero_cr = ( - signal_zerocrossings(cwtmatr[4, :][index_cur : index_next + 1])[0] + index_cur + signal_zerocrossings(cwtmatr[4, :][index_cur : index_next + 1])[0] + + index_cur ) nb_idx = int(max_search_duration * sampling_rate) - index_max = np.argmax(ecg[index_zero_cr - nb_idx : index_zero_cr + nb_idx]) + ( - index_zero_cr - nb_idx - ) + index_max = np.argmax( + ecg[index_zero_cr - nb_idx : index_zero_cr + nb_idx] + ) + (index_zero_cr - nb_idx) tppeaks.append(index_max) if len(tppeaks) == 0: tppeaks = [np.nan] @@ -904,7 +954,6 @@ def _find_tppeaks(ecg, keep_tp, sampling_rate=1000): # PEAK METHOD # ============================================================================= def _ecg_delineator_peak(ecg, rpeaks=None, sampling_rate=1000): - # Initialize heartbeats = ecg_segment(ecg, rpeaks, sampling_rate) @@ -920,7 +969,9 @@ def _ecg_delineator_peak(ecg, rpeaks=None, sampling_rate=1000): heartbeat = heartbeats[str(i + 1)] # Get index of heartbeat - R = heartbeat.index.get_loc(np.min(heartbeat.index.values[heartbeat.index.values > 0])) + R = heartbeat.index.get_loc( + np.min(heartbeat.index.values[heartbeat.index.values > 0]) + ) # Peaks ------ # Q wave @@ -980,7 +1031,8 @@ def _ecg_delineator_peak_P(rpeak, heartbeat, R, Q): segment = heartbeat.iloc[:Q] # Select left of Q wave P = signal_findpeaks( - segment["Signal"], height_min=0.05 * (segment["Signal"].max() - segment["Signal"].min()) + segment["Signal"], + height_min=0.05 * (segment["Signal"].max() - segment["Signal"].min()), ) if len(P["Peaks"]) == 0: @@ -993,7 +1045,8 @@ def _ecg_delineator_peak_P(rpeak, heartbeat, R, Q): def _ecg_delineator_peak_S(rpeak, heartbeat): segment = heartbeat[0:] # Select right hand side S = signal_findpeaks( - -segment["Signal"], height_min=0.05 * (segment["Signal"].max() - segment["Signal"].min()) + -segment["Signal"], + height_min=0.05 * (segment["Signal"].max() - segment["Signal"].min()), ) if len(S["Peaks"]) == 0: @@ -1008,7 +1061,8 @@ def _ecg_delineator_peak_T(rpeak, heartbeat, R, S): segment = heartbeat.iloc[R + S :] # Select right of S wave T = signal_findpeaks( - segment["Signal"], height_min=0.05 * (segment["Signal"].max() - segment["Signal"].min()) + segment["Signal"], + height_min=0.05 * (segment["Signal"].max() - segment["Signal"].min()), ) if len(T["Peaks"]) == 0: @@ -1070,7 +1124,6 @@ def _ecg_delineate_plot( window_start=-0.35, window_end=0.55, ): - """ import neurokit2 as nk import numpy as np @@ -1155,7 +1208,9 @@ def _ecg_delineate_plot( ax.plot(epoch_data.Time, epoch_data.Signal, color="grey", alpha=0.2) for i, feature_type in enumerate(features.columns.values): # pylint: disable=W0612 event_data = data[data[feature_type] == 1.0] - ax.scatter(event_data.Time, event_data.Signal, label=feature_type, alpha=0.5, s=200) + ax.scatter( + event_data.Time, event_data.Signal, label=feature_type, alpha=0.5, s=200 + ) ax.legend() return fig @@ -1188,5 +1243,7 @@ def _calculate_abs_z(df, columns): """This function helps to calculate the absolute standardized distance between R-peaks and other delineated waves features by `ecg_delineate()`""" for column in columns: - df["Dist_R_" + column] = np.abs(standardize(df[column].sub(df["ECG_R_Peaks"], axis=0))) + df["Dist_R_" + column] = np.abs( + standardize(df[column].sub(df["ECG_R_Peaks"], axis=0)) + ) return df diff --git a/neurokit2/ecg/ecg_intervalrelated.py b/neurokit2/ecg/ecg_intervalrelated.py index 81880e9fd0..d735b2869e 100644 --- a/neurokit2/ecg/ecg_intervalrelated.py +++ b/neurokit2/ecg/ecg_intervalrelated.py @@ -127,6 +127,7 @@ def _ecg_intervalrelated_hrv(data, sampling_rate, output={}): results = hrv(rpeaks, sampling_rate=sampling_rate) for column in results.columns: - output[column] = float(results[column]) + # Add and convert to float + output[column] = results[[column]].values return output diff --git a/neurokit2/ecg/ecg_process.py b/neurokit2/ecg/ecg_process.py index 247f436265..488f287da8 100644 --- a/neurokit2/ecg/ecg_process.py +++ b/neurokit2/ecg/ecg_process.py @@ -60,7 +60,7 @@ def ecg_process(ecg_signal, sampling_rate=1000, method="neurokit"): * **This list is not up-to-date. Help us improve the documentation!** info : dict A dictionary containing the samples at which the R-peaks occur, accessible with the key - ``"ECG_Peaks"``, as well as the signals' sampling rate. + ``"ECG_R_Peaks"``, as well as the signals' sampling rate. See Also -------- diff --git a/neurokit2/ecg/ecg_quality.py b/neurokit2/ecg/ecg_quality.py index 3c50058483..41594d35d5 100644 --- a/neurokit2/ecg/ecg_quality.py +++ b/neurokit2/ecg/ecg_quality.py @@ -406,4 +406,4 @@ def _ecg_quality_basSQI( num_power = psd.iloc[0][0] dem_power = psd.iloc[0][1] - return 1 - num_power / dem_power + return (1 - num_power) / dem_power diff --git a/neurokit2/eda/eda_process.py b/neurokit2/eda/eda_process.py index e978dbb676..abff6d7e28 100644 --- a/neurokit2/eda/eda_process.py +++ b/neurokit2/eda/eda_process.py @@ -20,7 +20,7 @@ def eda_process(eda_signal, sampling_rate=1000, method="neurokit", report=None, eda_signal : Union[list, np.array, pd.Series] The raw EDA signal. sampling_rate : int - The sampling frequency of ``"rsp_signal"`` (in Hz, i.e., samples/second). + The sampling frequency of ``"eda_signal"`` (in Hz, i.e., samples/second). method : str The processing pipeline to apply. Can be one of ``"biosppy"`` or ``"neurokit"`` (default). report : str diff --git a/neurokit2/eeg/mne_channel_extract.py b/neurokit2/eeg/mne_channel_extract.py index 8e93a18cab..f403913abd 100644 --- a/neurokit2/eeg/mne_channel_extract.py +++ b/neurokit2/eeg/mne_channel_extract.py @@ -64,7 +64,7 @@ def mne_channel_extract(raw, what, name=None, add_firstsamples=False): "check channel names in raw.info['ch_names']. " ) - channels, __ = raw.copy().pick_channels(what)[:] + channels, __ = raw.copy().pick_channels(what, ordered=False)[:] if len(what) > 1: channels = pd.DataFrame(channels.T, columns=what) if name is not None: diff --git a/neurokit2/emg/emg_clean.py b/neurokit2/emg/emg_clean.py index 3293679862..f5c2c33d66 100644 --- a/neurokit2/emg/emg_clean.py +++ b/neurokit2/emg/emg_clean.py @@ -9,7 +9,7 @@ from ..signal import signal_detrend -def emg_clean(emg_signal, sampling_rate=1000): +def emg_clean(emg_signal, sampling_rate=1000, method="biosppy"): """**Preprocess an electromyography (emg) signal** Clean an EMG signal using a set of parameters. Only one method is available at the moment. @@ -24,6 +24,10 @@ def emg_clean(emg_signal, sampling_rate=1000): sampling_rate : int The sampling frequency of ``emg_signal`` (in Hz, i.e., samples/second). Defaults to 1000. + method : str + The processing pipeline to apply. Can be one of ``"biosppy"`` or ``"none"``. + Defaults to ``"biosppy"``. If ``"none"`` is passed, the raw signal will be returned without + any cleaning. Returns ------- @@ -61,6 +65,32 @@ def emg_clean(emg_signal, sampling_rate=1000): ) emg_signal = _emg_clean_missing(emg_signal) + method = str(method).lower() + if method in ["none"]: + clean = emg_signal + elif method in ["biosppy"]: + clean = _emg_clean_biosppy(emg_signal, sampling_rate=sampling_rate) + else: + raise ValueError( + "NeuroKit error: emg_clean(): 'method' should be one of 'biosppy' or 'none'." + ) + return clean + + +# ============================================================================= +# Handle missing data +# ============================================================================= +def _emg_clean_missing(emg_signal): + + emg_signal = pd.DataFrame.pad(pd.Series(emg_signal)) + + return emg_signal + + +# ============================================================================= +# BioSPPy +# ============================================================================= +def _emg_clean_biosppy(emg_signal, sampling_rate=1000): # Parameters order = 4 frequency = 100 @@ -76,13 +106,3 @@ def emg_clean(emg_signal, sampling_rate=1000): clean = signal_detrend(filtered, order=0) return clean - - -# ============================================================================= -# Handle missing data -# ============================================================================= -def _emg_clean_missing(emg_signal): - - emg_signal = pd.DataFrame.pad(pd.Series(emg_signal)) - - return emg_signal diff --git a/neurokit2/emg/emg_methods.py b/neurokit2/emg/emg_methods.py new file mode 100644 index 0000000000..9bee832f4c --- /dev/null +++ b/neurokit2/emg/emg_methods.py @@ -0,0 +1,135 @@ +# -*- coding: utf-8 -*- +import numpy as np + +from ..misc.report import get_kwargs +from .emg_activation import emg_activation + + +def emg_methods( + sampling_rate=1000, + method_cleaning="biosppy", + method_activation="threshold", + **kwargs, +): + """**EMG Preprocessing Methods** + + This function analyzes and specifies the methods used in the preprocessing, and create a + textual description of the methods used. It is used by :func:`eda_process()` to dispatch the + correct methods to each subroutine of the pipeline and to create a + preprocessing report. + + Parameters + ---------- + sampling_rate : int + The sampling frequency of the raw EMG signal (in Hz, i.e., samples/second). + method_cleaning : str + The method used for cleaning the raw EMG signal. Can be one of ``"biosppy"`` or ``"none"``. + Defaults to ``"biosppy"``. If ``"none"`` is passed, the raw signal will be used without + any cleaning. + method_activation: str + The method used for locating EMG activity. Defaults to ``"threshold"``. + For more information, see the ``"method"`` argument + of :func:`.emg_activation`. + **kwargs + Other arguments to be passed to :func:`.emg_activation`, + + Returns + ------- + report_info : dict + A dictionary containing the keyword arguments passed to the cleaning and activation + functions, text describing the methods, and the corresponding references. + + See Also + -------- + + """ + # Sanitize inputs + method_cleaning = str(method_cleaning).lower() + method_activation = str(method_activation).lower() + + # Create dictionary with all inputs + report_info = { + "sampling_rate": sampling_rate, + "method_cleaning": method_cleaning, + "method_activation": method_activation, + **kwargs, + } + + # Get arguments to be passed to activation function + kwargs_activation, report_info = get_kwargs(report_info, emg_activation) + + # Save keyword arguments in dictionary + report_info["kwargs_activation"] = kwargs_activation + + # Initialize refs list with NeuroKit2 reference + refs = [ + """Makowski, D., Pham, T., Lau, Z. J., Brammer, J. C., Lespinasse, F., Pham, H., + Schölzel, C., & Chen, S. A. (2021). NeuroKit2: A Python toolbox for neurophysiological signal processing. + Behavior Research Methods, 53(4), 1689–1696. https://doi.org/10.3758/s13428-020-01516-y + """ + ] + + # 1. Cleaning + # ------------ + # If no cleaning + report_info["text_cleaning"] = f"The raw signal, sampled at {sampling_rate} Hz," + if method_cleaning in ["none"]: + report_info["text_cleaning"] += " was directly used without any cleaning." + else: + report_info["text_cleaning"] += ( + " was cleaned using the " + method_cleaning + " method." + ) + + # 2. Activation + # ------------- + report_info["text_activation"] = ( + "EMG activity was detected using the " + method_activation + " method. " + ) + if method_activation in ["silva"]: + if str(report_info["threshold"]) == "default": + threshold_str = "0.05" + else: + threshold_str = str(report_info["threshold"]) + report_info["text_activation"] += f"""The threshold was {threshold_str}. """ + + refs.append( + """Silva H, Scherer R, Sousa J, Londral A , "Towards improving the ssability of + electromyographic interfacess", Journal of Oral Rehabilitation, pp. 1-2, 2012.""" + ) + if method_activation in ["mixture"]: + report_info[ + "text_activation" + ] += """A Gaussian mixture model was used to discriminate between activity and baseline. """ + if str(report_info["threshold"]) == "default": + threshold_str = "0.33" + else: + threshold_str = str(report_info["threshold"]) + report_info[ + "text_activation" + ] += f"""The minimum probability required to + be considered as activated was {threshold_str}. """ + elif method_activation in ["threshold"]: + report_info[ + "text_activation" + ] += """The signal was considered as activated when the amplitude exceeded a threshold. """ + if str(report_info["threshold"]) == "default": + threshold_str = "one tenth of the standard deviation of emg_amplitude" + else: + threshold_str = str(report_info["threshold"]) + report_info[ + "text_activation" + ] += f"""The minimum amplitude to detect as onset was set to {threshold_str}.""" + elif method_activation in ["biosppy"]: + if str(report_info["threshold"]) == "default": + threshold_str = "1.2 times of the mean of the absolute of the smoothed, full-wave-rectified signal" + else: + threshold_str = str(report_info["threshold"]) + report_info[ + "text_activation" + ] += f"""The threshold was set to {threshold_str}.""" + + # 3. References + # ------------- + report_info["references"] = list(np.unique(refs)) + + return report_info diff --git a/neurokit2/emg/emg_plot.py b/neurokit2/emg/emg_plot.py index d897eee90d..9b493202d0 100644 --- a/neurokit2/emg/emg_plot.py +++ b/neurokit2/emg/emg_plot.py @@ -4,7 +4,7 @@ import pandas as pd -def emg_plot(emg_signals, sampling_rate=None): +def emg_plot(emg_signals, sampling_rate=None, static=True): """**EMG Graph** Visualize electromyography (EMG) data. @@ -17,6 +17,10 @@ def emg_plot(emg_signals, sampling_rate=None): The sampling frequency of the EMG (in Hz, i.e., samples/second). Needs to be supplied if the data should be plotted over time in seconds. Otherwise the data is plotted over samples. Defaults to ``None``. + static : bool + If True, a static plot will be generated with matplotlib. + If False, an interactive plot will be generated with plotly. + Defaults to True. See Also -------- @@ -70,6 +74,31 @@ def emg_plot(emg_signals, sampling_rate=None): else: x_axis = np.arange(0, emg_signals.shape[0]) + if static is True: + return _emg_plot_static(emg_signals, x_axis, onsets, offsets, sampling_rate) + else: + return _emg_plot_interactive(emg_signals, x_axis, onsets, offsets, sampling_rate) + + +# ============================================================================= +# Internals +# ============================================================================= +def _emg_plot_activity(emg_signals, onsets, offsets): + + activity_signal = pd.Series(np.full(len(emg_signals), np.nan)) + activity_signal[onsets] = emg_signals["EMG_Amplitude"][onsets].values + activity_signal[offsets] = emg_signals["EMG_Amplitude"][offsets].values + activity_signal = activity_signal.fillna(method="backfill") + + if np.any(activity_signal.isna()): + index = np.min(np.where(activity_signal.isna())) - 1 + value_to_fill = activity_signal[index] + activity_signal = activity_signal.fillna(value_to_fill) + + return activity_signal + + +def _emg_plot_static(emg_signals, x_axis, onsets, offsets, sampling_rate): # Prepare figure. fig, (ax0, ax1) = plt.subplots(nrows=2, ncols=1, sharex=True) if sampling_rate is not None: @@ -126,21 +155,114 @@ def emg_plot(emg_signals, sampling_rate=None): ax1.axvline(i, color="#4a4a4a", linestyle="--", label=None, zorder=2) ax1.axvline(j, color="#4a4a4a", linestyle="--", label=None, zorder=2) ax1.legend(loc="upper right") + plt.close() + return fig -# ============================================================================= -# Internals -# ============================================================================= -def _emg_plot_activity(emg_signals, onsets, offsets): +def _emg_plot_interactive(emg_signals, x_axis, onsets, offsets, sampling_rate): + try: + import plotly.graph_objects as go + from plotly.subplots import make_subplots + except ImportError: + raise ImportError( + "NeuroKit error: emg_plot(): the 'plotly' " + "module is required for this feature." + "Please install it first (`pip install plotly`)." + ) - activity_signal = pd.Series(np.full(len(emg_signals), np.nan)) - activity_signal[onsets] = emg_signals["EMG_Amplitude"][onsets].values - activity_signal[offsets] = emg_signals["EMG_Amplitude"][offsets].values - activity_signal = activity_signal.fillna(method="backfill") + # Prepare figure. + fig = make_subplots(rows=2, cols=1, shared_xaxes=True) + fig.update_layout(title="Electromyography (EMG)", font=dict(size=18), height=600) - if np.any(activity_signal.isna()): - index = np.min(np.where(activity_signal.isna())) - 1 - value_to_fill = activity_signal[index] - activity_signal = activity_signal.fillna(value_to_fill) + # Plot cleaned and raw EMG. + fig.add_trace( + go.Scatter( + x=x_axis, + y=emg_signals["EMG_Raw"], + mode="lines", + name="Raw", + line=dict(color="#B0BEC5"), + ), + row=1, + col=1, + ) + fig.add_trace( + go.Scatter( + x=x_axis, + y=emg_signals["EMG_Clean"], + mode="lines", + name="Cleaned", + line=dict(color="#FFC107"), + ), + row=1, + col=1, + ) - return activity_signal + # Plot Amplitude. + fig.add_trace( + go.Scatter( + x=x_axis, + y=emg_signals["EMG_Amplitude"], + mode="lines", + name="Amplitude", + line=dict(color="#FF9800"), + ), + row=2, + col=1, + ) + + # Mark onsets and offsets. + fig.add_trace( + go.Scatter( + x=x_axis[onsets], + y=emg_signals["EMG_Amplitude"][onsets], + mode="markers", + name="Onsets", + marker=dict(color="#f03e65", size=10), + ), + row=2, + col=1, + ) + fig.add_trace( + go.Scatter( + x=x_axis[offsets], + y=emg_signals["EMG_Amplitude"][offsets], + mode="markers", + name="Offsets", + marker=dict(color="#f03e65", size=10), + ), + row=2, + col=1, + ) + + if sampling_rate is not None: + onsets = onsets / sampling_rate + offsets = offsets / sampling_rate + fig.update_xaxes(title_text="Time (seconds)", row=2, col=1) + elif sampling_rate is None: + fig.update_xaxes(title_text="Samples", row=2, col=1) + + for i, j in zip(list(onsets), list(offsets)): + fig.add_shape( + type="line", + x0=i, + y0=0, + x1=i, + y1=1, + line=dict(color="#4a4a4a", width=2, dash="dash"), + row=2, + col=1, + ) + fig.add_shape( + type="line", + x0=j, + y0=0, + x1=j, + y1=1, + line=dict(color="#4a4a4a", width=2, dash="dash"), + row=2, + col=1, + ) + + fig.update_yaxes(title_text="Amplitude", row=2, col=1) + return fig diff --git a/neurokit2/emg/emg_process.py b/neurokit2/emg/emg_process.py index dccbc07c30..4742770216 100644 --- a/neurokit2/emg/emg_process.py +++ b/neurokit2/emg/emg_process.py @@ -1,13 +1,16 @@ # -*- coding: utf-8 -*- import pandas as pd +from ..misc.report import create_report from ..signal import signal_sanitize from .emg_activation import emg_activation from .emg_amplitude import emg_amplitude from .emg_clean import emg_clean +from .emg_methods import emg_methods +from .emg_plot import emg_plot -def emg_process(emg_signal, sampling_rate=1000): +def emg_process(emg_signal, sampling_rate=1000, report=None, **kwargs): """**Process a electromyography (EMG) signal** Convenience function that automatically processes an electromyography signal. @@ -18,6 +21,14 @@ def emg_process(emg_signal, sampling_rate=1000): The raw electromyography channel. sampling_rate : int The sampling frequency of ``emg_signal`` (in Hz, i.e., samples/second). + report : str + The filename of a report containing description and figures of processing + (e.g. ``"myreport.html"``). Needs to be supplied if a report file + should be generated. Defaults to ``None``. Can also be ``"text"`` to + just print the text in the console without saving anything. + **kwargs + Other arguments to be passed to specific methods. For more information, + see :func:`.emg_methods`. Returns ------- @@ -57,16 +68,23 @@ def emg_process(emg_signal, sampling_rate=1000): """ # Sanitize input emg_signal = signal_sanitize(emg_signal) + methods = emg_methods(sampling_rate=sampling_rate, **kwargs) # Clean signal - emg_cleaned = emg_clean(emg_signal, sampling_rate=sampling_rate) + emg_cleaned = emg_clean( + emg_signal, sampling_rate=sampling_rate, method=methods["method_cleaning"] + ) # Get amplitude amplitude = emg_amplitude(emg_cleaned) # Get onsets, offsets, and periods of activity activity_signal, info = emg_activation( - amplitude, sampling_rate=sampling_rate, threshold="default" + emg_amplitude=amplitude, + emg_cleaned=emg_cleaned, + sampling_rate=sampling_rate, + method=methods["method_activation"], + **methods["kwargs_activation"] ) info["sampling_rate"] = sampling_rate # Add sampling rate in dict info @@ -77,4 +95,12 @@ def emg_process(emg_signal, sampling_rate=1000): signals = pd.concat([signals, activity_signal], axis=1) + if report is not None: + # Generate report containing description and figures of processing + if ".html" in str(report): + fig = emg_plot(signals, sampling_rate=sampling_rate, static=False) + else: + fig = None + create_report(file=report, signals=signals, info=methods, fig=fig) + return signals, info diff --git a/neurokit2/epochs/eventrelated_utils.py b/neurokit2/epochs/eventrelated_utils.py index eea4c8d104..1d5b743ba8 100644 --- a/neurokit2/epochs/eventrelated_utils.py +++ b/neurokit2/epochs/eventrelated_utils.py @@ -12,7 +12,14 @@ def _eventrelated_sanitizeinput(epochs, what="ecg", silent=False): # Sanity checks if isinstance(epochs, pd.DataFrame): - epochs = _df_to_epochs(epochs) # Convert df to dict + if "Time" in epochs.columns.values: + # Assumpe it is a dataframe of epochs created by `epochs_to_df` + epochs = _df_to_epochs(epochs) # Convert df back to dict + else: + raise ValueError( + "It seems that you are trying to pass a single epoch to an eventrelated function." + + 'If this is what you want, please wrap it in a dict: `{"0": epoch}` ' + ) if not isinstance(epochs, dict): raise ValueError( @@ -37,7 +44,6 @@ def _eventrelated_sanitizeinput(epochs, what="ecg", silent=False): def _eventrelated_addinfo(epoch, output={}): - # Add label if "Index" in epoch.columns: output["Event_Onset"] = epoch.loc[np.min(np.abs(epoch.index))]["Index"] @@ -58,7 +64,6 @@ def _eventrelated_addinfo(epoch, output={}): def _eventrelated_sanitizeoutput(data): - df = pd.DataFrame.from_dict(data, orient="index") # Convert to a dataframe colnames = df.columns.values @@ -76,12 +81,12 @@ def _eventrelated_sanitizeoutput(data): def _eventrelated_rate(epoch, output={}, var="ECG_Rate"): - # Sanitize input colnames = epoch.columns.values if len([i for i in colnames if var in i]) == 0: warn( - "Input does not have an `" + var + "` column." " Will skip all rate-related features.", + "Input does not have an `" + var + "` column." + " Will skip all rate-related features.", category=NeuroKitWarning, ) return output diff --git a/neurokit2/events/events_find.py b/neurokit2/events/events_find.py index e4fea5d4bf..da94c98060 100644 --- a/neurokit2/events/events_find.py +++ b/neurokit2/events/events_find.py @@ -36,8 +36,8 @@ def events_find( threshold_keep : str ``"above"`` or ``"below"``, define the events as above or under the threshold. For photosensors, a white screen corresponds usually to higher values. Therefore, if your - events are signaled by a black colour, events values are the lower ones, and you should set - the cut to ``"below"``. + events are signaled by a black colour, events values are the lower ones (i.e., the signal + "drops" when the events onset), and you should set the cut to ``"below"``. start_at : int Keep events which onset is after a particular time point. end_at : int @@ -109,7 +109,9 @@ def events_find( """ - events = _events_find(event_channel, threshold=threshold, threshold_keep=threshold_keep) + events = _events_find( + event_channel, threshold=threshold, threshold_keep=threshold_keep + ) # Warning when no events detected if len(events["onset"]) == 0: @@ -218,7 +220,12 @@ def _events_find_label( def _events_find(event_channel, threshold="auto", threshold_keep="above"): binary = signal_binarize(event_channel, threshold=threshold) - if threshold_keep.lower() != "above": + if threshold_keep not in ["above", "below"]: + raise ValueError( + "In events_find(), 'threshold_keep' must be one of 'above' or 'below'." + ) + + if threshold_keep != "above": binary = np.abs(binary - 1) # Reverse if events are below # Initialize data diff --git a/neurokit2/hrv/hrv_frequency.py b/neurokit2/hrv/hrv_frequency.py index d56667bc2d..5b61ad6286 100644 --- a/neurokit2/hrv/hrv_frequency.py +++ b/neurokit2/hrv/hrv_frequency.py @@ -39,6 +39,7 @@ def hrv_frequency( * **LF**: The spectral power of low frequencies (by default, .04 to .15 Hz). * **HF**: The spectral power of high frequencies (by default, .15 to .4 Hz). * **VHF**: The spectral power of very high frequencies (by default, .4 to .5 Hz). + * **TP**: The total spectral power. * **LFHF**: The ratio obtained by dividing the low frequency power by the high frequency power. * **LFn**: The normalized low frequency, obtained by dividing the low frequency power by the total power. @@ -212,6 +213,7 @@ def hrv_frequency( # Normalized total_power = np.nansum(power.values) + out["TP"] = total_power out["LFHF"] = out["LF"] / out["HF"] out["LFn"] = out["LF"] / total_power out["HFn"] = out["HF"] / total_power diff --git a/neurokit2/hrv/hrv_nonlinear.py b/neurokit2/hrv/hrv_nonlinear.py index b79b8aaade..177a08954d 100644 --- a/neurokit2/hrv/hrv_nonlinear.py +++ b/neurokit2/hrv/hrv_nonlinear.py @@ -500,9 +500,9 @@ def _hrv_nonlinear_show(rri, rri_time=None, rri_missing=False, out={}, ax=None, sd1 = out["SD1"] sd2 = out["SD2"] if isinstance(sd1, pd.Series): - sd1 = float(sd1) + sd1 = float(sd1.iloc[0]) if isinstance(sd2, pd.Series): - sd2 = float(sd2) + sd2 = float(sd2.iloc[0]) # Poincare values ax1 = rri[:-1] diff --git a/neurokit2/hrv/hrv_time.py b/neurokit2/hrv/hrv_time.py index c4fc76bfde..aa977dbbcc 100644 --- a/neurokit2/hrv/hrv_time.py +++ b/neurokit2/hrv/hrv_time.py @@ -70,6 +70,8 @@ def hrv_time(peaks, sampling_rate=1000, show=False, **kwargs): * **MCVNN**: The median absolute deviation of the RR intervals (**MadNN**) divided by the median of the RR intervals (**MedianNN**). * **IQRNN**: The interquartile range (**IQR**) of the RR intervals. + * **SDRMSSD**: SDNN / RMSSD, a time-domain equivalent for the low Frequency-to-High + Frequency (LF/HF) Ratio (Sollers et al., 2007). * **Prc20NN**: The 20th percentile of the RR intervals (Han, 2017; Hovsepian, 2015). * **Prc80NN**: The 80th percentile of the RR intervals (Han, 2017; Hovsepian, 2015). * **pNN50**: The proportion of RR intervals greater than 50ms, out of the total number of @@ -132,6 +134,10 @@ def hrv_time(peaks, sampling_rate=1000, show=False, **kwargs): * Subramaniam, S. D., & Dass, B. (2022). An Efficient Convolutional Neural Network for Acute Pain Recognition Using HRV Features. In Proceedings of the International e-Conference on Intelligent Systems and Signal Processing (pp. 119-132). Springer, Singapore. + * Sollers, J. J., Buchanan, T. W., Mowrer, S. M., Hill, L. K., & Thayer, J. F. (2007). + Comparison of the ratio of the standard deviation of the RR interval and the root mean + squared successive differences (SD/rMSSD) to the low frequency-to-high frequency (LF/HF) + ratio in a patient population and normal healthy controls. Biomed Sci Instrum, 43, 158-163. """ # Sanitize input @@ -166,6 +172,7 @@ def hrv_time(peaks, sampling_rate=1000, show=False, **kwargs): out["MadNN"] = mad(rri) out["MCVNN"] = out["MadNN"] / out["MedianNN"] # Normalized out["IQRNN"] = scipy.stats.iqr(rri) + out["SDRMSSD"] = out["SDNN"] / out["RMSSD"] # Sollers (2007) out["Prc20NN"] = np.nanpercentile(rri, q=20) out["Prc80NN"] = np.nanpercentile(rri, q=80) @@ -200,7 +207,6 @@ def hrv_time(peaks, sampling_rate=1000, show=False, **kwargs): def _hrv_time_show(rri, **kwargs): - fig = summary_plot(rri, **kwargs) plt.xlabel("R-R intervals (ms)") fig.suptitle("Distribution of R-R intervals") @@ -209,7 +215,6 @@ def _hrv_time_show(rri, **kwargs): def _sdann(rri, rri_time=None, window=1): - window_size = window * 60 * 1000 # Convert window in min to ms if rri_time is None: # Compute the timestamps of the R-R intervals in seconds @@ -230,7 +235,6 @@ def _sdann(rri, rri_time=None, window=1): def _sdnni(rri, rri_time=None, window=1): - window_size = window * 60 * 1000 # Convert window in min to ms if rri_time is None: # Compute the timestamps of the R-R intervals in seconds @@ -267,10 +271,14 @@ def _hrv_TINN(rri, bar_x, bar_y, binsize): while m < np.max(rri): n_start = np.where(bar_x == n)[0][0] n_end = np.where(bar_x == X)[0][0] - qn = np.polyval(np.polyfit([n, X], [0, Y], deg=1), bar_x[n_start : n_end + 1]) + qn = np.polyval( + np.polyfit([n, X], [0, Y], deg=1), bar_x[n_start : n_end + 1] + ) m_start = np.where(bar_x == X)[0][0] m_end = np.where(bar_x == m)[0][0] - qm = np.polyval(np.polyfit([X, m], [Y, 0], deg=1), bar_x[m_start : m_end + 1]) + qm = np.polyval( + np.polyfit([X, m], [Y, 0], deg=1), bar_x[m_start : m_end + 1] + ) q = np.zeros(len(bar_x)) q[n_start : n_end + 1] = qn q[m_start : m_end + 1] = qm diff --git a/neurokit2/hrv/intervals_process.py b/neurokit2/hrv/intervals_process.py index 5e89bd3704..c5131e953d 100644 --- a/neurokit2/hrv/intervals_process.py +++ b/neurokit2/hrv/intervals_process.py @@ -51,6 +51,8 @@ def intervals_process( Preprocessed intervals, in milliseconds. np.ndarray Preprocessed timestamps corresponding to intervals, in seconds. + int + Sampling rate (Hz) of the interpolated interbeat intervals. Examples -------- diff --git a/neurokit2/hrv/intervals_to_peaks.py b/neurokit2/hrv/intervals_to_peaks.py index fddf3ad728..f07e6e0308 100644 --- a/neurokit2/hrv/intervals_to_peaks.py +++ b/neurokit2/hrv/intervals_to_peaks.py @@ -13,9 +13,10 @@ def intervals_to_peaks(intervals, intervals_time=None, sampling_rate=1000): Parameters ---------- intervals : list or array - List or numpy array of intervals, in milliseconds. + List of intervals (by default in milliseconds). intervals_time : list or array, optional - List or numpy array of timestamps corresponding to intervals, in seconds. + Optional list of timestamps corresponding to intervals, in seconds. If None (default), the + cumulative sum of the intervals is used. sampling_rate : int, optional Sampling rate (Hz) of the continuous signal in which the peaks occur. @@ -31,9 +32,11 @@ def intervals_to_peaks(intervals, intervals_time=None, sampling_rate=1000): import neurokit2 as nk + # Suppose we have a vector of RRi from data sampled at 1000 Hz ibi = [500, 400, 700, 500, 300, 800, 500] - peaks = nk.intervals_to_peaks(ibi) + peaks = nk.intervals_to_peaks(ibi, sampling_rate=1000) + # We can then use NeuroKit's functionalities to compute HRV indices @savefig p_intervals_to_peaks.png scale=100% hrv_indices = nk.hrv_time(peaks, sampling_rate=100, show=True) @suppress @@ -41,6 +44,13 @@ def intervals_to_peaks(intervals, intervals_time=None, sampling_rate=1000): hrv_indices + .. ipython:: python + + # We can also use the timestamps of the intervals + rri = [400, 500, 700, 800, 900] + rri_idx = [0.7, 1.2, 2.5, 3.3, 4.2] + nk.intervals_to_peaks(rri, rri_idx, sampling_rate=1000) + """ if intervals is None: return None diff --git a/neurokit2/hrv/intervals_utils.py b/neurokit2/hrv/intervals_utils.py index 160b5ed413..4f41d9b947 100644 --- a/neurokit2/hrv/intervals_utils.py +++ b/neurokit2/hrv/intervals_utils.py @@ -53,13 +53,15 @@ def _intervals_successive(intervals, intervals_time=None, thresh_unequal=10, n_d if intervals_time is None: intervals_time = np.nancumsum(intervals / 1000) else: - intervals_time = np.array(intervals_time) + intervals_time = np.array(intervals_time).astype(float) intervals_time[np.isnan(intervals)] = np.nan diff_intervals_time_ms = np.diff(intervals_time, n=n_diff) * 1000 - abs_error_intervals_ref_time = abs(diff_intervals_time_ms - np.diff(intervals[1:], n=n_diff - 1)) + abs_error_intervals_ref_time = abs( + diff_intervals_time_ms - np.diff(intervals[1:], n=n_diff - 1) + ) successive_intervals = abs_error_intervals_ref_time <= thresh_unequal @@ -125,7 +127,9 @@ def _intervals_sanitize(intervals, intervals_time=None, remove_missing=True): intervals = np.array(intervals) if intervals_time is None: # Impute intervals with median in case of missing values to calculate timestamps - imputed_intervals = np.where(np.isnan(intervals), np.nanmedian(intervals, axis=0), intervals) + imputed_intervals = np.where( + np.isnan(intervals), np.nanmedian(intervals, axis=0), intervals + ) # Compute the timestamps of the intervals in seconds intervals_time = np.nancumsum(imputed_intervals / 1000) else: @@ -141,7 +145,9 @@ def _intervals_sanitize(intervals, intervals_time=None, remove_missing=True): # If none of the differences between timestamps match # the length of the R-R intervals in seconds, # try converting milliseconds to seconds - converted_successive_intervals = _intervals_successive(intervals, intervals_time=intervals_time / 1000) + converted_successive_intervals = _intervals_successive( + intervals, intervals_time=intervals_time / 1000 + ) # Check if converting to seconds increased the number of differences # between timestamps that match the length of the R-R intervals in seconds diff --git a/neurokit2/misc/report.py b/neurokit2/misc/report.py index 444a92cc47..2da99ea3cf 100644 --- a/neurokit2/misc/report.py +++ b/neurokit2/misc/report.py @@ -19,13 +19,14 @@ def create_report(file="myreport.html", signals=None, info={"sampling_rate": 100 Name of the file to save the report to. Can also be ``"text"`` to simply print the text in the console. signals : pd.DataFrame - A DataFrame of signals. Usually obtained from :func:`.rsp_process` or :func:`.ppg_process` + A DataFrame of signals. Usually obtained from :func:`.rsp_process`, :func:`.ppg_process`, or + :func:`.emg_process`. info : dict A dictionary containing the information of peaks and the signals' sampling rate. Usually obtained from :func:`.rsp_process` or :func:`.ppg_process`. fig : matplotlib.figure.Figure or plotly.graph_objects.Figure - A figure containing the processed signals. Usually obtained from :func:`.rsp_plot` or - :func:`.ppg_plot`. + A figure containing the processed signals. Usually obtained from :func:`.rsp_plot`, + :func:`.ppg_plot`, or :func:`.emg_plot`. Returns ------- @@ -34,7 +35,7 @@ def create_report(file="myreport.html", signals=None, info={"sampling_rate": 100 See Also -------- - rsp_process, ppg_process + rsp_process, ppg_process, emg_process Examples -------- @@ -71,7 +72,7 @@ def create_report(file="myreport.html", signals=None, info={"sampling_rate": 100 def summarize_table(signals): - """Create table to summarize statistics of a RSP signal.""" + """Create table to summarize statistics of a signal.""" # TODO: add more features summary = {} diff --git a/neurokit2/ppg/ppg_clean.py b/neurokit2/ppg/ppg_clean.py index 0d9f14eb24..81742780a0 100644 --- a/neurokit2/ppg/ppg_clean.py +++ b/neurokit2/ppg/ppg_clean.py @@ -2,10 +2,9 @@ from warnings import warn import numpy as np -import pandas as pd from ..misc import NeuroKitWarning, as_vector -from ..signal import signal_filter +from ..signal import signal_fillmissing, signal_filter def ppg_clean(ppg_signal, sampling_rate=1000, heart_rate=None, method="elgendi"): @@ -23,8 +22,9 @@ def ppg_clean(ppg_signal, sampling_rate=1000, heart_rate=None, method="elgendi") sampling_rate : int The sampling frequency of the PPG (in Hz, i.e., samples/second). The default is 1000. method : str - The processing pipeline to apply. Can be one of ``"elgendi"`` or ``"nabian2018"``. The - default is ``"elgendi"``. + The processing pipeline to apply. Can be one of ``"elgendi"``, ``"nabian2018"``, or ``"none"``. + The default is ``"elgendi"``. If ``"none"`` is passed, the raw signal will be returned without + any cleaning. Returns ------- @@ -72,41 +72,32 @@ def ppg_clean(ppg_signal, sampling_rate=1000, heart_rate=None, method="elgendi") if n_missing > 0: warn( "There are " + str(n_missing) + " missing data points in your signal." - " Filling missing values by using the forward filling method.", + " Filling missing values using `signal_fillmissing`.", category=NeuroKitWarning, ) - ppg_signal = _ppg_clean_missing(ppg_signal) + ppg_signal = signal_fillmissing(ppg_signal, method="both") - method = method.lower() + method = str(method).lower() if method in ["elgendi"]: clean = _ppg_clean_elgendi(ppg_signal, sampling_rate) elif method in ["nabian2018"]: clean = _ppg_clean_nabian2018(ppg_signal, sampling_rate, heart_rate=heart_rate) - elif method is None or method == "none": + elif method in ["none"]: clean = ppg_signal else: - raise ValueError("`method` not found. Must be one of 'elgendi' or 'nabian2018'.") + raise ValueError( + "`method` not found. Must be one of 'elgendi', 'nabian2018', or 'none'." + ) return clean -# ============================================================================= -# Handle missing data -# ============================================================================= -def _ppg_clean_missing(ppg_signal): - - ppg_signal = pd.DataFrame.pad(pd.Series(ppg_signal)) - - return ppg_signal - - # ============================================================================= # Methods # ============================================================================= def _ppg_clean_elgendi(ppg_signal, sampling_rate): - filtered = signal_filter( ppg_signal, sampling_rate=sampling_rate, diff --git a/neurokit2/ppg/ppg_findpeaks.py b/neurokit2/ppg/ppg_findpeaks.py index ca636659be..0746681f63 100644 --- a/neurokit2/ppg/ppg_findpeaks.py +++ b/neurokit2/ppg/ppg_findpeaks.py @@ -220,8 +220,8 @@ def _ppg_findpeaks_bishop( # - column-wise summation m_max_sum = np.sum(m_max == False, axis=0) m_min_sum = np.sum(m_min == False, axis=0) - peaks = np.asarray(np.where(m_max_sum == 0)).astype(int) - onsets = np.asarray(np.where(m_min_sum == 0)).astype(int) + peaks = np.where(m_max_sum == 0)[0].astype(int) + onsets = np.where(m_min_sum == 0)[0].astype(int) if show: _, ax0 = plt.subplots(nrows=1, ncols=1, sharex=True) diff --git a/neurokit2/ppg/ppg_intervalrelated.py b/neurokit2/ppg/ppg_intervalrelated.py index 31d31a024b..7f91483140 100644 --- a/neurokit2/ppg/ppg_intervalrelated.py +++ b/neurokit2/ppg/ppg_intervalrelated.py @@ -21,11 +21,11 @@ def ppg_intervalrelated(data, sampling_rate=1000): Returns ------- DataFrame - A dataframe containing the analyzed ECG features. The analyzed features consist of the following: + A dataframe containing the analyzed PPG features. The analyzed features consist of the following: * ``"PPG_Rate_Mean"``: the mean heart rate. - * ``"ECG_HRV"``: the different heart rate variability metrices. + * ``"HRV"``: the different heart rate variability metrices. See :func:`.hrv` docstrings for details. @@ -61,7 +61,7 @@ def ppg_intervalrelated(data, sampling_rate=1000): rate_cols = [col for col in data.columns if "PPG_Rate" in col] if len(rate_cols) == 1: intervals.update(_ppg_intervalrelated_formatinput(data)) - intervals.update(_ppg_intervalrelated_hrv(data, sampling_rate)) + intervals.update(*_ppg_intervalrelated_hrv(data, sampling_rate)) else: raise ValueError( "NeuroKit error: ppg_intervalrelated(): Wrong input," @@ -71,20 +71,18 @@ def ppg_intervalrelated(data, sampling_rate=1000): ppg_intervals = pd.DataFrame.from_dict(intervals, orient="index").T elif isinstance(data, dict): - for index in data: - intervals[index] = {} # Initialize empty container + for epoch in data: + intervals[epoch] = {} # Initialize empty container # Add label info - intervals[index]["Label"] = data[index]["Label"].iloc[0] + intervals[epoch]["Label"] = data[epoch]["Label"].iloc[0] # Rate - intervals[index] = _ppg_intervalrelated_formatinput( - data[index], intervals[index] - ) + intervals[epoch].update(_ppg_intervalrelated_formatinput(data[epoch])) # HRV - intervals[index] = _ppg_intervalrelated_hrv( - data[index], sampling_rate, intervals[index] + intervals[epoch].update( + *_ppg_intervalrelated_hrv(data[epoch], sampling_rate) ) ppg_intervals = pd.DataFrame.from_dict(intervals, orient="index") @@ -97,31 +95,32 @@ def ppg_intervalrelated(data, sampling_rate=1000): # ============================================================================= -def _ppg_intervalrelated_formatinput(data, output={}): - +def _ppg_intervalrelated_formatinput(data): # Sanitize input colnames = data.columns.values - if len([i for i in colnames if "PPG_Rate" in i]) == 0: + + if "PPG_Rate" not in colnames: raise ValueError( "NeuroKit error: ppg_intervalrelated(): Wrong input," "we couldn't extract heart rate. Please make sure" - "your DataFrame contains an `PPG_Rate` column." + "your DataFrame contains a `PPG_Rate` column." ) signal = data["PPG_Rate"].values - output["PPG_Rate_Mean"] = np.mean(signal) - return output + PPG_Rate_Mean = np.mean(signal) + return {"PPG_Rate_Mean": PPG_Rate_Mean} -def _ppg_intervalrelated_hrv(data, sampling_rate, output={}): +def _ppg_intervalrelated_hrv(data, sampling_rate): # Sanitize input colnames = data.columns.values - if len([i for i in colnames if "PPG_Peaks" in i]) == 0: + + if "PPG_Peaks" not in colnames: raise ValueError( "NeuroKit error: ppg_intervalrelated(): Wrong input," "we couldn't extract peaks. Please make sure" - "your DataFrame contains an `PPG_Peaks` column." + "your DataFrame contains a `PPG_Peaks` column." ) # Transform rpeaks from "signal" format to "info" format. @@ -129,7 +128,5 @@ def _ppg_intervalrelated_hrv(data, sampling_rate, output={}): peaks = {"PPG_Peaks": peaks} results = hrv(peaks, sampling_rate=sampling_rate) - for column in results.columns: - output[column] = float(results[column]) - return output + return results.astype("float").to_dict("records") diff --git a/neurokit2/ppg/ppg_methods.py b/neurokit2/ppg/ppg_methods.py index 8b4e778516..57554c7cd3 100644 --- a/neurokit2/ppg/ppg_methods.py +++ b/neurokit2/ppg/ppg_methods.py @@ -65,8 +65,14 @@ def ppg_methods( """ # Sanitize inputs - method_cleaning = str(method).lower() if method_cleaning == "default" else str(method_cleaning).lower() - method_peaks = str(method).lower() if method_peaks == "default" else str(method_peaks).lower() + method_cleaning = ( + str(method).lower() + if method_cleaning == "default" + else str(method_cleaning).lower() + ) + method_peaks = ( + str(method).lower() if method_peaks == "default" else str(method_peaks).lower() + ) # Create dictionary with all inputs report_info = { @@ -85,8 +91,11 @@ def ppg_methods( report_info["kwargs_cleaning"] = kwargs_cleaning report_info["kwargs_peaks"] = kwargs_peaks - # Initialize refs list - refs = [] + # Initialize refs list with NeuroKit2 reference + refs = ["""Makowski, D., Pham, T., Lau, Z. J., Brammer, J. C., Lespinasse, F., Pham, H., + Schölzel, C., & Chen, S. A. (2021). NeuroKit2: A Python toolbox for neurophysiological signal processing. + Behavior Research Methods, 53(4), 1689–1696. https://doi.org/10.3758/s13428-020-01516-y + """] # 1. Cleaning # ------------ @@ -121,10 +130,14 @@ def ppg_methods( IEEE Journal of Translational Engineering in Health and Medicine, 6, 1-11.""" ) elif method_cleaning in ["none"]: - report_info["text_cleaning"] += " was directly used for peak detection without preprocessing." + report_info[ + "text_cleaning" + ] += " was directly used for peak detection without preprocessing." else: # just in case more methods are added - report_info["text_cleaning"] = "was cleaned following the " + method + " method." + report_info["text_cleaning"] = ( + "was cleaned following the " + method + " method." + ) # 2. Peaks # ---------- @@ -141,7 +154,9 @@ def ppg_methods( elif method_peaks in ["none"]: report_info["text_peaks"] = "There was no peak detection carried out." else: - report_info["text_peaks"] = f"The peak detection was carried out using the method {method_peaks}." + report_info[ + "text_peaks" + ] = f"The peak detection was carried out using the method {method_peaks}." report_info["references"] = list(np.unique(refs)) return report_info diff --git a/neurokit2/ppg/ppg_process.py b/neurokit2/ppg/ppg_process.py index a7ffee2c42..a361a3d8d2 100644 --- a/neurokit2/ppg/ppg_process.py +++ b/neurokit2/ppg/ppg_process.py @@ -10,7 +10,10 @@ from .ppg_methods import ppg_methods from .ppg_plot import ppg_plot -def ppg_process(ppg_signal, sampling_rate=1000, method="elgendi", report=None, **kwargs): + +def ppg_process( + ppg_signal, sampling_rate=1000, method="elgendi", report=None, **kwargs +): """**Process a photoplethysmogram (PPG) signal** Convenience function that automatically processes a photoplethysmogram signal. @@ -87,7 +90,9 @@ def ppg_process(ppg_signal, sampling_rate=1000, method="elgendi", report=None, * info["sampling_rate"] = sampling_rate # Add sampling rate in dict info # Mark peaks - peaks_signal = _signal_from_indices(info["PPG_Peaks"], desired_length=len(ppg_cleaned)) + peaks_signal = _signal_from_indices( + info["PPG_Peaks"], desired_length=len(ppg_cleaned) + ) # Rate computation rate = signal_rate( diff --git a/neurokit2/ppg/ppg_simulate.py b/neurokit2/ppg/ppg_simulate.py index b30aeaad27..cb928633ac 100644 --- a/neurokit2/ppg/ppg_simulate.py +++ b/neurokit2/ppg/ppg_simulate.py @@ -2,10 +2,9 @@ import matplotlib.pyplot as plt import numpy as np -import scipy.interpolate from ..misc import check_random_state, check_random_state_children -from ..signal import signal_distort +from ..signal import signal_distort, signal_interpolate def ppg_simulate( @@ -149,9 +148,8 @@ def ppg_simulate( # Interpolate a continuous signal between the landmarks (i.e., Cartesian # coordinates). - f = scipy.interpolate.Akima1DInterpolator(x_all, y_all) samples = np.arange(int(np.ceil(duration * sampling_rate))) - ppg = f(samples) + ppg = signal_interpolate(x_values=x_all, y_values=y_all, x_new=samples, method="akima") # Remove NAN (values outside interpolation range, i.e., after last sample). ppg[np.isnan(ppg)] = np.nanmean(ppg) diff --git a/neurokit2/rsp/rsp_intervalrelated.py b/neurokit2/rsp/rsp_intervalrelated.py index 87bd2c3e45..bb061d6554 100644 --- a/neurokit2/rsp/rsp_intervalrelated.py +++ b/neurokit2/rsp/rsp_intervalrelated.py @@ -93,7 +93,7 @@ def _rsp_intervalrelated_features(data, sampling_rate, output={}): output["RSP_Rate_Mean"] = np.nanmean(data["RSP_Rate"].values) rrv = rsp_rrv(data, sampling_rate=sampling_rate) for column in rrv.columns: - output[column] = float(rrv[column]) + output[column] = rrv[column].values.astype("float") if "RSP_Amplitude" in colnames: output["RSP_Amplitude_Mean"] = np.nanmean(data["RSP_Amplitude"].values) diff --git a/neurokit2/rsp/rsp_methods.py b/neurokit2/rsp/rsp_methods.py index 169dd35c22..81f4a48972 100644 --- a/neurokit2/rsp/rsp_methods.py +++ b/neurokit2/rsp/rsp_methods.py @@ -95,8 +95,11 @@ def rsp_methods( report_info["kwargs_peaks"] = kwargs_peaks report_info["kwargs_rvt"] = kwargs_rvt - # Initialize refs list - refs = [] + # Initialize refs list with NeuroKit2 reference + refs = ["""Makowski, D., Pham, T., Lau, Z. J., Brammer, J. C., Lespinasse, F., Pham, H., + Schölzel, C., & Chen, S. A. (2021). NeuroKit2: A Python toolbox for neurophysiological signal processing. + Behavior Research Methods, 53(4), 1689–1696. https://doi.org/10.3758/s13428-020-01516-y + """] # 1. Cleaning # ------------ diff --git a/neurokit2/rsp/rsp_process.py b/neurokit2/rsp/rsp_process.py index 3de69c8854..2751218a14 100644 --- a/neurokit2/rsp/rsp_process.py +++ b/neurokit2/rsp/rsp_process.py @@ -66,7 +66,7 @@ def rsp_process( * ``"RSP_Rate"``: breathing rate interpolated between inhalation peaks. * ``"RSP_Amplitude"``: breathing amplitude interpolated between inhalation peaks. * ``"RSP_Phase"``: breathing phase, marked by "1" for inspiration and "0" for expiration. - * ``"RSP_PhaseCompletion"``: breathing phase completion, expressed in percentage (from 0 to + * ``"RSP_Phase_Completion"``: breathing phase completion, expressed in percentage (from 0 to 1), representing the stage of the current respiratory phase. * ``"RSP_RVT"``: respiratory volume per time (RVT). info : dict diff --git a/neurokit2/signal/__init__.py b/neurokit2/signal/__init__.py index cfefdf8fcf..68d56cde6d 100644 --- a/neurokit2/signal/__init__.py +++ b/neurokit2/signal/__init__.py @@ -5,6 +5,7 @@ from .signal_decompose import signal_decompose from .signal_detrend import signal_detrend from .signal_distort import signal_distort +from .signal_fillmissing import signal_fillmissing from .signal_filter import signal_filter from .signal_findpeaks import signal_findpeaks from .signal_fixpeaks import signal_fixpeaks @@ -59,4 +60,5 @@ "signal_timefrequency", "signal_sanitize", "signal_flatline", + "signal_fillmissing", ] diff --git a/neurokit2/signal/signal_fillmissing.py b/neurokit2/signal/signal_fillmissing.py new file mode 100644 index 0000000000..fdf1289fd2 --- /dev/null +++ b/neurokit2/signal/signal_fillmissing.py @@ -0,0 +1,38 @@ +import pandas as pd + + +def signal_fillmissing(signal, method="both"): + """**Handle missing values** + + Fill missing values in a signal using specific methods. + + Parameters + ---------- + signal : Union[list, np.array, pd.Series] + The signal (i.e., a time series) in the form of a vector of values. + method : str + The method to use to fill missing values. Can be one of ``"forward"``, ``"backward"``, + or ``"both"``. The default is ``"both"``. + + Returns + ------- + signal + + Examples + -------- + .. ipython:: python + + import neurokit2 as nk + + signal = [np.nan, 1, 2, 3, np.nan, np.nan, 6, 7, np.nan] + nk.signal_fillmissing(signal, method="forward") + nk.signal_fillmissing(signal, method="backward") + nk.signal_fillmissing(signal, method="both") + + """ + if method in ["forward", "forwards", "ffill", "both"]: + signal = pd.Series(signal).fillna(method="ffill").values + + if method in ["backward", "backwards", "back", "bfill", "both"]: + signal = pd.Series(signal).fillna(method="bfill").values + return signal diff --git a/neurokit2/signal/signal_filter.py b/neurokit2/signal/signal_filter.py index 37ada44b7f..94065c9d6f 100644 --- a/neurokit2/signal/signal_filter.py +++ b/neurokit2/signal/signal_filter.py @@ -6,6 +6,7 @@ import scipy.signal from ..misc import NeuroKitWarning +from .signal_interpolate import signal_interpolate def signal_filter( @@ -143,10 +144,12 @@ def signal_filter( """ method = method.lower() + signal_sanitized, missing = _signal_filter_missing(signal) + if method in ["sg", "savgol", "savitzky-golay"]: - filtered = _signal_filter_savgol(signal, sampling_rate, order, window_size=window_size) + filtered = _signal_filter_savgol(signal_sanitized, sampling_rate, order, window_size=window_size) elif method in ["powerline"]: - filtered = _signal_filter_powerline(signal, sampling_rate, powerline) + filtered = _signal_filter_powerline(signal_sanitized, sampling_rate, powerline) else: # Sanity checks @@ -154,15 +157,15 @@ def signal_filter( raise ValueError("NeuroKit error: signal_filter(): you need to specify a 'lowcut' or a 'highcut'.") if method in ["butter", "butterworth"]: - filtered = _signal_filter_butterworth(signal, sampling_rate, lowcut, highcut, order) + filtered = _signal_filter_butterworth(signal_sanitized, sampling_rate, lowcut, highcut, order) elif method in ["butter_ba", "butterworth_ba"]: - filtered = _signal_filter_butterworth_ba(signal, sampling_rate, lowcut, highcut, order) + filtered = _signal_filter_butterworth_ba(signal_sanitized, sampling_rate, lowcut, highcut, order) elif method in ["butter_zi", "butterworth_zi"]: - filtered = _signal_filter_butterworth_zi(signal, sampling_rate, lowcut, highcut, order) + filtered = _signal_filter_butterworth_zi(signal_sanitized, sampling_rate, lowcut, highcut, order) elif method in ["bessel"]: - filtered = _signal_filter_bessel(signal, sampling_rate, lowcut, highcut, order) + filtered = _signal_filter_bessel(signal_sanitized, sampling_rate, lowcut, highcut, order) elif method in ["fir"]: - filtered = _signal_filter_fir(signal, sampling_rate, lowcut, highcut, window_size=window_size) + filtered = _signal_filter_fir(signal_sanitized, sampling_rate, lowcut, highcut, window_size=window_size) else: raise ValueError( "NeuroKit error: signal_filter(): 'method' should be", @@ -170,6 +173,8 @@ def signal_filter( " 'savgol' or 'fir'.", ) + filtered[missing] = np.nan + if show is True: plt.plot(signal, color="lightgrey") plt.plot(filtered, color="red", alpha=0.9) @@ -353,3 +358,12 @@ def _signal_filter_windowsize(window_size="default", sampling_rate=1000): if (window_size % 2) == 0: window_size + 1 # pylint: disable=W0104 return window_size + + +def _signal_filter_missing(signal): + """Interpolate missing data and save the indices of the missing data.""" + missing = np.where(np.isnan(signal))[0] + if len(missing) > 0: + return signal_interpolate(signal, method="linear"), missing + else: + return signal, missing diff --git a/neurokit2/signal/signal_interpolate.py b/neurokit2/signal/signal_interpolate.py index 291b259ff5..8110cc2b6e 100644 --- a/neurokit2/signal/signal_interpolate.py +++ b/neurokit2/signal/signal_interpolate.py @@ -1,10 +1,16 @@ # -*- coding: utf-8 -*- +from warnings import warn + import numpy as np import pandas as pd import scipy.interpolate +from ..misc import NeuroKitWarning + -def signal_interpolate(x_values, y_values=None, x_new=None, method="quadratic", fill_value=None): +def signal_interpolate( + x_values, y_values=None, x_new=None, method="quadratic", fill_value=None +): """**Interpolate a signal** Interpolate a signal using different methods. @@ -25,13 +31,15 @@ def signal_interpolate(x_values, y_values=None, x_new=None, method="quadratic", first and the last values of x_values. method : str Method of interpolation. Can be ``"linear"``, ``"nearest"``, ``"zero"``, ``"slinear"``, - ``"quadratic"``, ``"cubic"``, ``"previous"``, ``"next"`` or ``"monotone_cubic"``. The - methods ``"zero"``, ``"slinear"``,``"quadratic"`` and ``"cubic"`` refer to a spline + ``"quadratic"``, ``"cubic"``, ``"previous"``, ``"next"``, ``"monotone_cubic"``, or ``"akima"``. + The methods ``"zero"``, ``"slinear"``,``"quadratic"`` and ``"cubic"`` refer to a spline interpolation of zeroth, first, second or third order; whereas ``"previous"`` and ``"next"`` simply return the previous or next value of the point. An integer specifying the order of the spline interpolator to use. See `here `_ for details on the ``"monotone_cubic"`` method. + See `here `_ for details on the ``"akima"`` method. fill_value : float or tuple or str If a ndarray (or float), this value will be used to fill in for requested points outside of the data range. @@ -82,7 +90,9 @@ def signal_interpolate(x_values, y_values=None, x_new=None, method="quadratic", """ # Sanity checks if x_values is None: - raise ValueError("NeuroKit error: signal_interpolate(): x_values must be provided.") + raise ValueError( + "NeuroKit error: signal_interpolate(): x_values must be provided." + ) if y_values is None: # for interpolating NaNs return _signal_interpolate_nan(x_values, method=method, fill_value=fill_value) @@ -100,8 +110,16 @@ def signal_interpolate(x_values, y_values=None, x_new=None, method="quadratic", x_new = np.linspace(x_values[0], x_values[-1], x_new) else: # if x_values is identical to x_new, no need for interpolation - if np.all(x_values == x_new): + if np.array_equal(x_values, x_new): return y_values + elif np.any(x_values[1:] == x_values[:-1]): + warn( + "Duplicate x values detected. Averaging their corresponding y values.", + category=NeuroKitWarning, + ) + x_values, y_values = _signal_interpolate_average_duplicates( + x_values, y_values + ) # If only one value, return a constant signal if len(x_values) == 1: @@ -111,6 +129,10 @@ def signal_interpolate(x_values, y_values=None, x_new=None, method="quadratic", interpolation_function = scipy.interpolate.PchipInterpolator( x_values, y_values, extrapolate=True ) + elif method == "akima": + interpolation_function = scipy.interpolate.Akima1DInterpolator( + x_values, y_values + ) else: if fill_value is None: fill_value = ([y_values[0]], [y_values[-1]]) @@ -160,8 +182,18 @@ def _signal_interpolate_nan(values, method="quadratic", fill_value=None): # interpolate to get the values at the indices where they are missing return signal_interpolate( - x_values=x_values, y_values=y_values, x_new=x_new, method=method, fill_value=fill_value + x_values=x_values, + y_values=y_values, + x_new=x_new, + method=method, + fill_value=fill_value, ) else: # if there are no missing values, return original values return values + + +def _signal_interpolate_average_duplicates(x_values, y_values): + unique_x, indices = np.unique(x_values, return_inverse=True) + mean_y = np.bincount(indices, weights=y_values) / np.bincount(indices) + return unique_x, mean_y diff --git a/neurokit2/signal/signal_period.py b/neurokit2/signal/signal_period.py index 78071c2750..8d02422748 100644 --- a/neurokit2/signal/signal_period.py +++ b/neurokit2/signal/signal_period.py @@ -8,9 +8,17 @@ from .signal_interpolate import signal_interpolate -def signal_period(peaks, sampling_rate=1000, desired_length=None, interpolation_method="monotone_cubic"): +def signal_period( + peaks, + sampling_rate=1000, + desired_length=None, + interpolation_method="monotone_cubic", +): """**Calculate signal period from a series of peaks** + Calculate the period of a signal from a series of peaks. The period is defined as the time + in seconds between two consecutive peaks. + Parameters ---------- peaks : Union[list, np.array, pd.DataFrame, pd.Series, dict] @@ -51,12 +59,20 @@ def signal_period(peaks, sampling_rate=1000, desired_length=None, interpolation_ import neurokit2 as nk - signal = nk.signal_simulate(duration=10, sampling_rate=1000, frequency=1) - info = nk.signal_findpeaks(signal) + # Generate 2 signals (with fixed and variable period) + sig1 = nk.signal_simulate(duration=20, sampling_rate=200, frequency=1) + sig2 = nk.ecg_simulate(duration=20, sampling_rate=200, heart_rate=60) + + # Find peaks + info1 = nk.signal_findpeaks(sig1) + info2 = nk.ecg_findpeaks(sig2, sampling_rate=200) + + # Compute period + period1 = nk.signal_period(peaks=info1["Peaks"], desired_length=len(sig1), sampling_rate=200) + period2 = nk.signal_period(peaks=info2["ECG_R_Peaks"], desired_length=len(sig2), sampling_rate=200) - period = nk.signal_period(peaks=info["Peaks"], desired_length=len(signal)) @savefig p_signal_period.png scale=100% - nk.signal_plot(period) + nk.signal_plot([period1, period2], subplots=True) @suppress plt.close() @@ -67,13 +83,15 @@ def signal_period(peaks, sampling_rate=1000, desired_length=None, interpolation_ if np.size(peaks) <= 3: warn( "Too few peaks detected to compute the rate. Returning empty vector.", - category=NeuroKitWarning + category=NeuroKitWarning, ) return np.full(desired_length, np.nan) if isinstance(desired_length, (int, float)): if desired_length <= peaks[-1]: - raise ValueError("NeuroKit error: desired_length must be None or larger than the index of the last peak.") + raise ValueError( + "NeuroKit error: desired_length must be None or larger than the index of the last peak." + ) # Calculate period in sec, based on peak to peak difference and make sure # that rate has the same number of elements as peaks (important for @@ -83,6 +101,8 @@ def signal_period(peaks, sampling_rate=1000, desired_length=None, interpolation_ # Interpolate all statistics to desired length. if desired_length is not None: - period = signal_interpolate(peaks, period, x_new=np.arange(desired_length), method=interpolation_method) + period = signal_interpolate( + peaks, period, x_new=np.arange(desired_length), method=interpolation_method + ) return period diff --git a/neurokit2/signal/signal_plot.py b/neurokit2/signal/signal_plot.py index d9057d8523..12b7810f6c 100644 --- a/neurokit2/signal/signal_plot.py +++ b/neurokit2/signal/signal_plot.py @@ -138,7 +138,7 @@ def signal_plot( vector = signal[col] events.append(np.where(vector == np.max(vector.unique()))[0]) events_plot(events, signal=signal[continuous_columns]) - if sampling_rate is None and signal.index.is_integer(): + if sampling_rate is None and pd.api.types.is_integer_dtype(signal.index): plt.gca().set_xlabel("Samples") else: plt.gca().set_xlabel(title_x) @@ -172,7 +172,7 @@ def signal_plot( else: _ = signal[continuous_columns].plot(subplots=False, sharex=True, **kwargs) - if sampling_rate is None and signal.index.is_integer(): + if sampling_rate is None and pd.api.types.is_integer_dtype(signal.index): plt.xlabel("Samples") else: plt.xlabel(title_x) diff --git a/setup.py b/setup.py index e65aaf0dbc..c943f1d5fc 100644 --- a/setup.py +++ b/setup.py @@ -20,7 +20,8 @@ def find_version(): result = re.search( - r'{}\s*=\s*[\'"]([^\'"]*)[\'"]'.format("__version__"), open("neurokit2/__init__.py").read() + r'{}\s*=\s*[\'"]([^\'"]*)[\'"]'.format("__version__"), + open("neurokit2/__init__.py").read(), ) return result.group(1) @@ -77,9 +78,8 @@ def find_version(): "Intended Audience :: Developers", "License :: OSI Approved :: MIT License", "Natural Language :: English", - "Programming Language :: Python :: 3.7", - "Programming Language :: Python :: 3.8", "Programming Language :: Python :: 3.9", "Programming Language :: Python :: 3.10", + "Programming Language :: Python :: 3.11", ], ) diff --git a/tests/tests_complexity.py b/tests/tests_complexity.py index 3c74bfa461..43a583a105 100644 --- a/tests/tests_complexity.py +++ b/tests/tests_complexity.py @@ -5,7 +5,8 @@ import numpy as np import pandas as pd from pyentrp import entropy as pyentrp -from sklearn.neighbors import KDTree +import sklearn.neighbors +from packaging import version # import EntropyHub import neurokit2 as nk @@ -340,7 +341,11 @@ def entropy_embed(x, order=3, delay=1): def entropy_app_samp_entropy(x, order, metric="chebyshev", approximate=True): - _all_metrics = KDTree.valid_metrics + sklearn_version = version.parse(sklearn.__version__) + if sklearn_version >= version.parse("1.3.0"): + _all_metrics = sklearn.neighbors.KDTree.valid_metrics() + else: + _all_metrics = sklearn.neighbors.KDTree.valid_metrics if metric not in _all_metrics: raise ValueError( "The given metric (%s) is not valid. The valid " # pylint: disable=consider-using-f-string @@ -356,14 +361,14 @@ def entropy_app_samp_entropy(x, order, metric="chebyshev", approximate=True): else: emb_data1 = _emb_data1[:-1] count1 = ( - KDTree(emb_data1, metric=metric) + sklearn.neighbors.KDTree(emb_data1, metric=metric) .query_radius(emb_data1, r, count_only=True) .astype(np.float64) ) # compute phi(order + 1, r) emb_data2 = entropy_embed(x, order + 1, 1) count2 = ( - KDTree(emb_data2, metric=metric) + sklearn.neighbors.KDTree(emb_data2, metric=metric) .query_radius(emb_data2, r, count_only=True) .astype(np.float64) ) diff --git a/tests/tests_ecg.py b/tests/tests_ecg.py index 73acd931f0..67b0f6a605 100644 --- a/tests/tests_ecg.py +++ b/tests/tests_ecg.py @@ -8,8 +8,9 @@ def test_ecg_simulate(): - - ecg1 = nk.ecg_simulate(duration=20, length=5000, method="simple", noise=0, random_state=0) + ecg1 = nk.ecg_simulate( + duration=20, length=5000, method="simple", noise=0, random_state=0 + ) assert len(ecg1) == 5000 ecg2 = nk.ecg_simulate(duration=20, length=5000, heart_rate=500, random_state=1) @@ -27,7 +28,6 @@ def test_ecg_simulate(): def test_ecg_clean(): - sampling_rate = 1000 noise = 0.05 @@ -45,29 +45,23 @@ def test_ecg_clean(): assert np.sum(fft_raw[freqs < 0.5]) > np.sum(fft_nk[freqs < 0.5]) # Comparison to biosppy - # (https://github.com/PIA-Group/BioSPPy/blob/e65da30f6379852ecb98f8e2e0c9b4b5175416c3/biosppy/signals/ecg.py#L69) ecg_biosppy = nk.ecg_clean(ecg, sampling_rate=sampling_rate, method="biosppy") - original, _, _ = biosppy.tools.filter_signal( - signal=ecg, - ftype="FIR", - band="bandpass", - order=int(0.3 * sampling_rate), - frequency=[3, 45], - sampling_rate=sampling_rate, - ) - assert np.allclose((ecg_biosppy - original).mean(), 0, atol=1e-6) + assert np.allclose(ecg_biosppy.mean(), 0, atol=1e-6) def test_ecg_peaks(): - sampling_rate = 1000 noise = 0.15 - ecg = nk.ecg_simulate(duration=120, sampling_rate=sampling_rate, noise=noise, random_state=42) + ecg = nk.ecg_simulate( + duration=120, sampling_rate=sampling_rate, noise=noise, random_state=42 + ) ecg_cleaned_nk = nk.ecg_clean(ecg, sampling_rate=sampling_rate, method="neurokit") # Test without request to correct artifacts. - signals, _ = nk.ecg_peaks(ecg_cleaned_nk, correct_artifacts=False, method="neurokit") + signals, _ = nk.ecg_peaks( + ecg_cleaned_nk, correct_artifacts=False, method="neurokit" + ) assert signals.shape == (120000, 1) assert np.allclose(signals["ECG_R_Peaks"].values.sum(dtype=np.int64), 139, atol=1) @@ -80,7 +74,6 @@ def test_ecg_peaks(): def test_ecg_process(): - sampling_rate = 1000 noise = 0.05 @@ -89,7 +82,6 @@ def test_ecg_process(): def test_ecg_plot(): - ecg = nk.ecg_simulate(duration=60, heart_rate=70, noise=0.05, random_state=5) ecg_summary, _ = nk.ecg_process(ecg, sampling_rate=1000, method="neurokit") @@ -100,7 +92,7 @@ def test_ecg_plot(): fig = plt.gcf() assert len(fig.axes) == 2 titles = ["Raw and Cleaned Signal", "Heart Rate"] - for (ax, title) in zip(fig.get_axes(), titles): + for ax, title in zip(fig.get_axes(), titles): assert ax.get_title() == title assert fig.get_axes()[1].get_xlabel() == "Samples" np.testing.assert_array_equal(fig.axes[0].get_xticks(), fig.axes[1].get_xticks()) @@ -112,7 +104,7 @@ def test_ecg_plot(): fig = plt.gcf() assert len(fig.axes) == 3 titles = ["Raw and Cleaned Signal", "Heart Rate", "Individual Heart Beats"] - for (ax, title) in zip(fig.get_axes(), titles): + for ax, title in zip(fig.get_axes(), titles): assert ax.get_title() == title assert fig.get_axes()[1].get_xlabel() == "Time (seconds)" np.testing.assert_array_equal(fig.axes[0].get_xticks(), fig.axes[1].get_xticks()) @@ -120,10 +112,15 @@ def test_ecg_plot(): def test_ecg_findpeaks(): - sampling_rate = 1000 - ecg = nk.ecg_simulate(duration=60, sampling_rate=sampling_rate, noise=0, method="simple", random_state=42) + ecg = nk.ecg_simulate( + duration=60, + sampling_rate=sampling_rate, + noise=0, + method="simple", + random_state=42, + ) ecg_cleaned = nk.ecg_clean(ecg, sampling_rate=sampling_rate, method="neurokit") @@ -136,11 +133,15 @@ def test_ecg_findpeaks(): assert len(fig.axes) == 2 # Test pantompkins1985 method - info_pantom = nk.ecg_findpeaks(nk.ecg_clean(ecg, method="pantompkins1985"), method="pantompkins1985") + info_pantom = nk.ecg_findpeaks( + nk.ecg_clean(ecg, method="pantompkins1985"), method="pantompkins1985" + ) assert info_pantom["ECG_R_Peaks"].size == 70 # Test hamilton2002 method - info_hamilton = nk.ecg_findpeaks(nk.ecg_clean(ecg, method="hamilton2002"), method="hamilton2002") + info_hamilton = nk.ecg_findpeaks( + nk.ecg_clean(ecg, method="hamilton2002"), method="hamilton2002" + ) assert info_hamilton["ECG_R_Peaks"].size == 69 # Test christov2004 method @@ -152,64 +153,89 @@ def test_ecg_findpeaks(): assert info_gamboa["ECG_R_Peaks"].size == 69 # Test elgendi2010 method - info_elgendi = nk.ecg_findpeaks(nk.ecg_clean(ecg, method="elgendi2010"), method="elgendi2010") + info_elgendi = nk.ecg_findpeaks( + nk.ecg_clean(ecg, method="elgendi2010"), method="elgendi2010" + ) assert info_elgendi["ECG_R_Peaks"].size == 70 # Test engzeemod2012 method - info_engzeemod = nk.ecg_findpeaks(nk.ecg_clean(ecg, method="engzeemod2012"), method="engzeemod2012") + info_engzeemod = nk.ecg_findpeaks( + nk.ecg_clean(ecg, method="engzeemod2012"), method="engzeemod2012" + ) assert info_engzeemod["ECG_R_Peaks"].size == 69 # Test kalidas2017 method - info_kalidas = nk.ecg_findpeaks(nk.ecg_clean(ecg, method="kalidas2017"), method="kalidas2017") + info_kalidas = nk.ecg_findpeaks( + nk.ecg_clean(ecg, method="kalidas2017"), method="kalidas2017" + ) assert np.allclose(info_kalidas["ECG_R_Peaks"].size, 68, atol=1) # Test martinez2004 method - ecg = nk.ecg_simulate(duration=60, sampling_rate=sampling_rate, noise=0, random_state=42) + ecg = nk.ecg_simulate( + duration=60, sampling_rate=sampling_rate, noise=0, random_state=42 + ) ecg_cleaned = nk.ecg_clean(ecg, sampling_rate=sampling_rate, method="neurokit") info_martinez = nk.ecg_findpeaks(ecg_cleaned, method="martinez2004") assert np.allclose(info_martinez["ECG_R_Peaks"].size, 69, atol=1) def test_ecg_eventrelated(): - ecg, _ = nk.ecg_process(nk.ecg_simulate(duration=20, random_state=6)) - epochs = nk.epochs_create(ecg, events=[5000, 10000, 15000], epochs_start=-0.1, epochs_end=1.9) + epochs = nk.epochs_create( + ecg, events=[5000, 10000, 15000], epochs_start=-0.1, epochs_end=1.9 + ) ecg_eventrelated = nk.ecg_eventrelated(epochs) # Test rate features - assert np.alltrue(np.array(ecg_eventrelated["ECG_Rate_Min"]) < np.array(ecg_eventrelated["ECG_Rate_Mean"])) + assert np.alltrue( + np.array(ecg_eventrelated["ECG_Rate_Min"]) + < np.array(ecg_eventrelated["ECG_Rate_Mean"]) + ) - assert np.alltrue(np.array(ecg_eventrelated["ECG_Rate_Mean"]) < np.array(ecg_eventrelated["ECG_Rate_Max"])) + assert np.alltrue( + np.array(ecg_eventrelated["ECG_Rate_Mean"]) + < np.array(ecg_eventrelated["ECG_Rate_Max"]) + ) assert len(ecg_eventrelated["Label"]) == 3 # Test warning on missing columns - with pytest.warns(nk.misc.NeuroKitWarning, match=r".*does not have an `ECG_Phase_Artrial`.*"): + with pytest.warns( + nk.misc.NeuroKitWarning, match=r".*does not have an `ECG_Phase_Artrial`.*" + ): first_epoch_key = list(epochs.keys())[0] first_epoch_copy = epochs[first_epoch_key].copy() del first_epoch_copy["ECG_Phase_Atrial"] nk.ecg_eventrelated({**epochs, first_epoch_key: first_epoch_copy}) - with pytest.warns(nk.misc.NeuroKitWarning, match=r".*does not have an.*`ECG_Phase_Ventricular`"): + with pytest.warns( + nk.misc.NeuroKitWarning, match=r".*does not have an.*`ECG_Phase_Ventricular`" + ): first_epoch_key = list(epochs.keys())[0] first_epoch_copy = epochs[first_epoch_key].copy() del first_epoch_copy["ECG_Phase_Ventricular"] nk.ecg_eventrelated({**epochs, first_epoch_key: first_epoch_copy}) - with pytest.warns(nk.misc.NeuroKitWarning, match=r".*does not have an `ECG_Quality`.*"): + with pytest.warns( + nk.misc.NeuroKitWarning, match=r".*does not have an `ECG_Quality`.*" + ): first_epoch_key = list(epochs.keys())[0] first_epoch_copy = epochs[first_epoch_key].copy() del first_epoch_copy["ECG_Quality"] nk.ecg_eventrelated({**epochs, first_epoch_key: first_epoch_copy}) - with pytest.warns(nk.misc.NeuroKitWarning, match=r".*does not have an `ECG_Rate`.*"): + with pytest.warns( + nk.misc.NeuroKitWarning, match=r".*does not have an `ECG_Rate`.*" + ): first_epoch_key = list(epochs.keys())[0] first_epoch_copy = epochs[first_epoch_key].copy() del first_epoch_copy["ECG_Rate"] nk.ecg_eventrelated({**epochs, first_epoch_key: first_epoch_copy}) # Test warning on long epochs (eventrelated_utils) - with pytest.warns(nk.misc.NeuroKitWarning, match=r".*duration of your epochs seems.*"): + with pytest.warns( + nk.misc.NeuroKitWarning, match=r".*duration of your epochs seems.*" + ): first_epoch_key = list(epochs.keys())[0] first_epoch_copy = epochs[first_epoch_key].copy() first_epoch_copy.index = range(len(first_epoch_copy)) @@ -217,7 +243,6 @@ def test_ecg_eventrelated(): def test_ecg_delineate(): - sampling_rate = 1000 # test with simulated signals @@ -226,7 +251,9 @@ def test_ecg_delineate(): number_rpeaks = len(rpeaks["ECG_R_Peaks"]) # Method 1: derivative - _, waves_derivative = nk.ecg_delineate(ecg, rpeaks, sampling_rate=sampling_rate, method="peaks") + _, waves_derivative = nk.ecg_delineate( + ecg, rpeaks, sampling_rate=sampling_rate, method="peaks" + ) assert len(waves_derivative["ECG_P_Peaks"]) == number_rpeaks assert len(waves_derivative["ECG_Q_Peaks"]) == number_rpeaks assert len(waves_derivative["ECG_S_Peaks"]) == number_rpeaks @@ -235,7 +262,9 @@ def test_ecg_delineate(): assert len(waves_derivative["ECG_T_Offsets"]) == number_rpeaks # Method 2: CWT - _, waves_cwt = nk.ecg_delineate(ecg, rpeaks, sampling_rate=sampling_rate, method="cwt") + _, waves_cwt = nk.ecg_delineate( + ecg, rpeaks, sampling_rate=sampling_rate, method="cwt" + ) assert np.allclose(len(waves_cwt["ECG_P_Peaks"]), 22, atol=1) assert np.allclose(len(waves_cwt["ECG_T_Peaks"]), 22, atol=1) assert np.allclose(len(waves_cwt["ECG_R_Onsets"]), 23, atol=1) @@ -247,7 +276,6 @@ def test_ecg_delineate(): def test_ecg_invert(): - sampling_rate = 500 noise = 0.05 @@ -258,7 +286,6 @@ def test_ecg_invert(): def test_ecg_intervalrelated(): - data = nk.data("bio_resting_5min_100hz") df, _ = nk.ecg_process(data["ECG"], sampling_rate=100) @@ -328,7 +355,9 @@ def test_ecg_intervalrelated(): features_df = nk.ecg_intervalrelated(df, sampling_rate=100) # https://github.com/neuropsychology/NeuroKit/issues/304 - assert all(features_df == nk.ecg_analyze(df, sampling_rate=100, method="interval-related")) + assert all( + features_df == nk.ecg_analyze(df, sampling_rate=100, method="interval-related") + ) assert (elem in columns for elem in np.array(features_df.columns.values, dtype=str)) assert features_df.shape[0] == 1 # Number of rows @@ -338,5 +367,7 @@ def test_ecg_intervalrelated(): epochs = nk.epochs_create(df, events=[0, 15000], sampling_rate=100, epochs_end=150) features_dict = nk.ecg_intervalrelated(epochs, sampling_rate=100) - assert (elem in columns for elem in np.array(features_dict.columns.values, dtype=str)) + assert ( + elem in columns for elem in np.array(features_dict.columns.values, dtype=str) + ) assert features_dict.shape[0] == 2 # Number of rows diff --git a/tests/tests_eda.py b/tests/tests_eda.py index 7cac122f50..601b4973eb 100644 --- a/tests/tests_eda.py +++ b/tests/tests_eda.py @@ -103,8 +103,6 @@ def test_eda_peaks(): assert np.allclose((info["SCR_Peaks"] - peaks).mean(), 0, atol=1e-5) signals, info = nk.eda_peaks(eda_phasic, method="kim2004") - onsets, peaks, amplitudes = biosppy.eda.kbk_scr(eda_phasic, sampling_rate=1000) - assert np.allclose((info["SCR_Peaks"] - peaks).mean(), 0, atol=180) # Check that indices and values positions match peak_positions = np.where(info["SCR_Peaks"] != 0)[0] diff --git a/tests/tests_emg.py b/tests/tests_emg.py index 7b2dea7a8f..b757c6ded6 100644 --- a/tests/tests_emg.py +++ b/tests/tests_emg.py @@ -18,7 +18,9 @@ def test_emg_simulate(): assert len(emg1) == 5000 emg2 = nk.emg_simulate(duration=20, length=5000, burst_number=15) - assert scipy.stats.median_abs_deviation(emg1) < scipy.stats.median_abs_deviation(emg2) + assert scipy.stats.median_abs_deviation(emg1) < scipy.stats.median_abs_deviation( + emg2 + ) emg3 = nk.emg_simulate(duration=20, length=5000, burst_number=1, burst_duration=2.0) # pd.DataFrame({"EMG1":emg1, "EMG3": emg3}).plot() @@ -71,9 +73,7 @@ def test_emg_plot(): emg_summary, _ = nk.emg_process(emg, sampling_rate=sampling_rate) # Plot data over samples. - nk.emg_plot(emg_summary) - # This will identify the latest figure. - fig = plt.gcf() + fig = nk.emg_plot(emg_summary) assert len(fig.axes) == 2 titles = ["Raw and Cleaned Signal", "Muscle Activation"] for (ax, title) in zip(fig.get_axes(), titles): @@ -83,9 +83,7 @@ def test_emg_plot(): plt.close(fig) # Plot data over time. - nk.emg_plot(emg_summary, sampling_rate=sampling_rate) - # This will identify the latest figure. - fig = plt.gcf() + fig = nk.emg_plot(emg_summary, sampling_rate=sampling_rate) assert fig.get_axes()[1].get_xlabel() == "Time (seconds)" @@ -114,19 +112,25 @@ def test_emg_eventrelated(): assert len(emg_eventrelated["Label"]) == 3 # Test warning on missing columns - with pytest.warns(nk.misc.NeuroKitWarning, match=r".*does not have an `EMG_Onsets`.*"): + with pytest.warns( + nk.misc.NeuroKitWarning, match=r".*does not have an `EMG_Onsets`.*" + ): first_epoch_key = list(epochs.keys())[0] first_epoch_copy = epochs[first_epoch_key].copy() del first_epoch_copy["EMG_Onsets"] nk.emg_eventrelated({**epochs, first_epoch_key: first_epoch_copy}) - with pytest.warns(nk.misc.NeuroKitWarning, match=r".*does not have an `EMG_Activity`.*"): + with pytest.warns( + nk.misc.NeuroKitWarning, match=r".*does not have an `EMG_Activity`.*" + ): first_epoch_key = list(epochs.keys())[0] first_epoch_copy = epochs[first_epoch_key].copy() del first_epoch_copy["EMG_Activity"] nk.emg_eventrelated({**epochs, first_epoch_key: first_epoch_copy}) - with pytest.warns(nk.misc.NeuroKitWarning, match=r".*does not have an.*`EMG_Amplitude`.*"): + with pytest.warns( + nk.misc.NeuroKitWarning, match=r".*does not have an.*`EMG_Amplitude`.*" + ): first_epoch_key = list(epochs.keys())[0] first_epoch_copy = epochs[first_epoch_key].copy() del first_epoch_copy["EMG_Amplitude"] @@ -142,13 +146,52 @@ def test_emg_intervalrelated(): # Test with signal dataframe features_df = nk.emg_intervalrelated(emg_signals) - assert all(elem in columns for elem in np.array(features_df.columns.values, dtype=str)) + assert all( + elem in columns for elem in np.array(features_df.columns.values, dtype=str) + ) assert features_df.shape[0] == 1 # Number of rows # Test with dict columns.append("Label") - epochs = nk.epochs_create(emg_signals, events=[0, 20000], sampling_rate=1000, epochs_end=20) + epochs = nk.epochs_create( + emg_signals, events=[0, 20000], sampling_rate=1000, epochs_end=20 + ) features_dict = nk.emg_intervalrelated(epochs) - assert all(elem in columns for elem in np.array(features_dict.columns.values, dtype=str)) + assert all( + elem in columns for elem in np.array(features_dict.columns.values, dtype=str) + ) assert features_dict.shape[0] == 2 # Number of rows + +@pytest.mark.parametrize( + "method_cleaning, method_activation, threshold", + [("none", "threshold", "default"), + ("biosppy", "pelt", 0.5), + ("biosppy", "mixture", 0.05), + ("biosppy", "biosppy", "default"), + ("biosppy", "silva", "default")], +) +def test_emg_report(tmp_path, method_cleaning, method_activation, threshold): + + sampling_rate = 250 + + emg = nk.emg_simulate( + duration=30, + sampling_rate=sampling_rate, + random_state=0, + ) + + d = tmp_path / "sub" + d.mkdir() + p = d / "myreport.html" + + signals, _ = nk.emg_process( + emg, + sampling_rate=sampling_rate, + report=str(p), + method_cleaning=method_cleaning, + method_activation=method_activation, + threshold=threshold + ) + assert p.is_file() + assert "EMG_Activity" in signals.columns diff --git a/tests/tests_ppg.py b/tests/tests_ppg.py index 5e5069139f..7f43e32015 100644 --- a/tests/tests_ppg.py +++ b/tests/tests_ppg.py @@ -8,8 +8,8 @@ import neurokit2 as nk -durations = (20, 200) -sampling_rates = (50, 500) +durations = (20, 200, 300) +sampling_rates = (25, 50, 500) heart_rates = (50, 120) freq_modulations = (0.1, 0.4) @@ -18,9 +18,10 @@ params_combis = list(itertools.product(*params)) -@pytest.mark.parametrize("duration, sampling_rate, heart_rate, freq_modulation", params_combis) +@pytest.mark.parametrize( + "duration, sampling_rate, heart_rate, freq_modulation", params_combis +) def test_ppg_simulate(duration, sampling_rate, heart_rate, freq_modulation): - ppg = nk.ppg_simulate( duration=duration, sampling_rate=sampling_rate, @@ -33,18 +34,19 @@ def test_ppg_simulate(duration, sampling_rate, heart_rate, freq_modulation): burst_amplitude=0, burst_number=0, random_state=42, + random_state_distort=42, show=False, ) assert ppg.size == duration * sampling_rate signals, _ = nk.ppg_process(ppg, sampling_rate=sampling_rate) - assert np.allclose(signals["PPG_Rate"].mean(), heart_rate, atol=1) - - # Ensure that the heart rate fluctuates in the requested range. - groundtruth_range = freq_modulation * heart_rate - observed_range = np.percentile(signals["PPG_Rate"], 90) - np.percentile(signals["PPG_Rate"], 10) - assert np.allclose(groundtruth_range, observed_range, atol=groundtruth_range * 0.15) + if sampling_rate > 25: + assert np.allclose(signals["PPG_Rate"].mean(), heart_rate, atol=1) + # Ensure that the heart rate fluctuates in the requested range. + groundtruth_range = freq_modulation * heart_rate + observed_range = np.percentile(signals["PPG_Rate"], 90) - np.percentile(signals["PPG_Rate"], 10) + assert np.allclose(groundtruth_range, observed_range, atol=groundtruth_range * 0.15) # TODO: test influence of different noise configurations @@ -54,7 +56,6 @@ def test_ppg_simulate(duration, sampling_rate, heart_rate, freq_modulation): [(0.1, 3), (0.2, 5), (0.3, 8), (0.4, 11), (0.5, 14), (0.6, 19)], ) def test_ppg_simulate_ibi(ibi_randomness, std_heart_rate): - ppg = nk.ppg_simulate( duration=20, sampling_rate=50, @@ -107,7 +108,6 @@ def test_ppg_simulate_legacy_rng(): def test_ppg_clean(): - sampling_rate = 500 ppg = nk.ppg_simulate( @@ -124,7 +124,9 @@ def test_ppg_clean(): random_state=42, show=False, ) - ppg_cleaned_elgendi = nk.ppg_clean(ppg, sampling_rate=sampling_rate, method="elgendi") + ppg_cleaned_elgendi = nk.ppg_clean( + ppg, sampling_rate=sampling_rate, method="elgendi" + ) assert ppg.size == ppg_cleaned_elgendi.size @@ -139,7 +141,6 @@ def test_ppg_clean(): def test_ppg_findpeaks(): - sampling_rate = 500 # Test Elgendi method @@ -157,9 +158,13 @@ def test_ppg_findpeaks(): random_state=42, show=True, ) - ppg_cleaned_elgendi = nk.ppg_clean(ppg, sampling_rate=sampling_rate, method="elgendi") + ppg_cleaned_elgendi = nk.ppg_clean( + ppg, sampling_rate=sampling_rate, method="elgendi" + ) - info_elgendi = nk.ppg_findpeaks(ppg_cleaned_elgendi, sampling_rate=sampling_rate, show=True) + info_elgendi = nk.ppg_findpeaks( + ppg_cleaned_elgendi, sampling_rate=sampling_rate, show=True + ) peaks = info_elgendi["PPG_Peaks"] @@ -167,7 +172,9 @@ def test_ppg_findpeaks(): assert np.abs(peaks.sum() - 219764) < 5 # off by no more than 5 samples in total # Test MSPTD method - info_msptd = nk.ppg_findpeaks(ppg, sampling_rate=sampling_rate, method="bishop", show=True) + info_msptd = nk.ppg_findpeaks( + ppg, sampling_rate=sampling_rate, method="bishop", show=True + ) peaks = info_msptd["PPG_Peaks"] @@ -177,11 +184,10 @@ def test_ppg_findpeaks(): @pytest.mark.parametrize( "method_cleaning, method_peaks", - [("elgendi", "elgendi"), ("nabian2018", "elgendi")], + [("elgendi", "elgendi"), ("nabian2018", "elgendi"), ("elgendi", "bishop")], ) def test_ppg_report(tmp_path, method_cleaning, method_peaks): - - sampling_rate = 500 + sampling_rate = 100 ppg = nk.ppg_simulate( duration=30, @@ -210,3 +216,32 @@ def test_ppg_report(tmp_path, method_cleaning, method_peaks): method_peaks=method_peaks, ) assert p.is_file() + + +def test_ppg_intervalrelated(): + sampling_rate = 100 + + ppg = nk.ppg_simulate( + duration=500, + sampling_rate=sampling_rate, + heart_rate=70, + frequency_modulation=0.025, + ibi_randomness=0.15, + drift=0.5, + motion_amplitude=0.25, + powerline_amplitude=0.25, + burst_amplitude=0.5, + burst_number=3, + random_state=0, + show=True, + ) + # Process the data + df, info = nk.ppg_process(ppg, sampling_rate=sampling_rate) + epochs = nk.epochs_create( + df, events=[0, 15000], sampling_rate=sampling_rate, epochs_end=150 + ) + epochs_ppg_intervals = nk.ppg_intervalrelated(epochs) + assert "PPG_Rate_Mean" in epochs_ppg_intervals.columns + + ppg_intervals = nk.ppg_intervalrelated(df) + assert "PPG_Rate_Mean" in ppg_intervals.columns diff --git a/tests/tests_signal.py b/tests/tests_signal.py index 1a7f2926e9..399ee7babd 100644 --- a/tests/tests_signal.py +++ b/tests/tests_signal.py @@ -186,6 +186,25 @@ def test_signal_filter(): assert np.allclose(signal_bandstop, signal_bandstop_scipy, atol=0.2) +def test_signal_filter_with_missing(): + sampling_rate = 100 + duration_not_missing = 10 + frequency = 2 + signal = np.concatenate( + [ + nk.signal_simulate(duration=duration_not_missing, sampling_rate=sampling_rate, frequency=frequency, random_state=42), + [np.nan] * 1000, + nk.signal_simulate(duration=duration_not_missing, sampling_rate=sampling_rate, frequency=frequency, random_state=43), + ] + ) + samples = np.arange(len(signal)) + powerline = np.sin(2 * np.pi * 50 * (samples / sampling_rate)) + signal_corrupted = signal + powerline + signal_clean = nk.signal_filter( + signal_corrupted, sampling_rate=sampling_rate, method="powerline" + ) + assert signal_clean.size == signal.size + assert np.allclose(signal_clean, signal, atol=0.2, equal_nan=True) def test_signal_interpolate():