Merge branch 'dev' into master

ThomasGebhardt · Feb 25, 2023 · 2c84e22 · 2c84e22
2 parents 575b148 + 7a916bd
commit 2c84e22
Show file tree

Hide file tree

Showing 26 changed files with 760 additions and 426 deletions.
diff --git a/AUTHORS.rst b/AUTHORS.rst
@@ -30,6 +30,7 @@ Core contributors
 Contributors
 -------------
 
+* `Gansheng Tan <https://github.com/GanshengT>`_ *(Washington University, USA)*
 * `Hung Pham <https://github.com/hungpham2511>`_ *(Eureka Robotics, Singapore)*
 * `Christopher Schölzel <https://github.com/CSchoel>`_ *(THM University of Applied Sciences, Germany)*
 * `Duy Le <https://github.com/duylp>`_ *(Hubble, Singapore)*
@@ -53,7 +54,7 @@ Contributors
 * `Marek Sokol <https://github.com/sokolmarek>`_ *(Faculty of Biomedical Engineering of the CTU in Prague, Czech Republic)*
 
 
-Thanks also to `Gansheng Tan <https://github.com/GanshengT>`_, `Chuan-Peng Hu <https://github.com/hcp4715>`_, `@ucohen <https://github.com/ucohen>`_, `Anthony Gatti <https://github.com/gattia>`_, `Julien Lamour <https://github.com/lamourj>`_, `@renatosc <https://github.com/renatosc>`_, `Nicolas Beaudoin-Gagnon <https://github.com/Fegalf>`_ and `@rubinovitz <https://github.com/rubinovitz>`_ for their contribution in `NeuroKit 1 <https://github.com/neuropsychology/NeuroKit.py>`_.
+Thanks also to `Chuan-Peng Hu <https://github.com/hcp4715>`_, `@ucohen <https://github.com/ucohen>`_, `Anthony Gatti <https://github.com/gattia>`_, `Julien Lamour <https://github.com/lamourj>`_, `@renatosc <https://github.com/renatosc>`_, `Nicolas Beaudoin-Gagnon <https://github.com/Fegalf>`_ and `@rubinovitz <https://github.com/rubinovitz>`_ for their contribution in `NeuroKit 1 <https://github.com/neuropsychology/NeuroKit.py>`_.
 
 
 

diff --git a/NEWS.rst b/NEWS.rst
@@ -1,7 +1,16 @@
 News
 =====
 
+0.2.4
+-------------------
+Fixes
++++++++++++++
 
+* `eda_sympathetic()` has been reviewed: low-pass filter and resampling have been added to be in
+  line with the original paper
+* `eda_findpeaks()` using methods proposed in nabian2018 is reviewed and improved. Differentiation
+  has been added before smoothing. Skin conductance response criteria have been revised based on
+  the original paper.
 
 
 

diff --git a/docs/functions/eda.rst b/docs/functions/eda.rst
@@ -31,14 +31,6 @@ Preprocessing
 """""""""""""""""""""
 .. autofunction:: neurokit2.eda.eda_phasic
 
-*eda_autocor()*
-"""""""""""""""
-.. autofunction:: neurokit2.eda.eda_autocor
-
-*eda_changepoints()*
-"""""""""""""""""""""
-.. autofunction:: neurokit2.eda.eda_changepoints
-
 *eda_peaks()*
 """""""""""""""""""""
 .. autofunction:: neurokit2.eda.eda_peaks
@@ -47,9 +39,6 @@ Preprocessing
 """""""""""""""""""""
 .. autofunction:: neurokit2.eda.eda_fixpeaks
 
-*eda_sympathetic()*
-"""""""""""""""""""""
-.. autofunction:: neurokit2.eda.eda_sympathetic
 
 
 Analysis
@@ -62,6 +51,18 @@ Analysis
 """""""""""""""""""""""""
 .. autofunction:: neurokit2.eda.eda_intervalrelated
 
+*eda_sympathetic()*
+"""""""""""""""""""""
+.. autofunction:: neurokit2.eda.eda_sympathetic
+
+*eda_autocor()*
+"""""""""""""""
+.. autofunction:: neurokit2.eda.eda_autocor
+
+*eda_changepoints()*
+"""""""""""""""""""""
+.. autofunction:: neurokit2.eda.eda_changepoints
+
 
 
 Miscellaneous

diff --git a/docs/functions/hrv.rst b/docs/functions/hrv.rst
@@ -46,5 +46,5 @@ Intervals
 
 .. automodule:: neurokit2.hrv
    :members:
-   :exclude-members: hrv, hrv_time, hrv_frequency, hrv_nonlinear, hrv_rqa, hrv_rsa
+   :exclude-members: hrv, hrv_time, hrv_frequency, hrv_nonlinear, hrv_rqa, hrv_rsa, intervals_process, intervals_to_peaks
 
diff --git a/neurokit2/__init__.py b/neurokit2/__init__.py
@@ -32,7 +32,7 @@
 from .video import *
 
 # Info
-__version__ = "0.2.3"
+__version__ = "0.2.4"
 
 
 # Maintainer info

diff --git a/neurokit2/ecg/ecg_delineate.py b/neurokit2/ecg/ecg_delineate.py
@@ -164,6 +164,11 @@ def ecg_delineate(
             "NeuroKit error: ecg_delineate(): 'method' should be one of 'peak'," "'cwt' or 'dwt'."
         )
 
+    # Ensure that all indices are not larger than ECG signal indices
+    for _, value in waves.items():
+        if value[-1] >= len(ecg_cleaned):
+            value[-1] = np.nan
+
     # Remove NaN in Peaks, Onsets, and Offsets
     waves_noNA = waves.copy()
     for feature in waves_noNA.keys():
@@ -947,11 +952,6 @@ def _ecg_delineator_peak(ecg, rpeaks=None, sampling_rate=1000):
         "ECG_T_Offsets": T_offsets,
     }
 
-    # Ensure that all indices are not larger than ECG signal indices
-    for _, value in info.items():
-        if value[-1] >= len(ecg):
-            value[-1] = np.nan
-
     # Return info dictionary
     return info
 

diff --git a/neurokit2/ecg/ecg_segment.py b/neurokit2/ecg/ecg_segment.py
@@ -59,6 +59,11 @@ def ecg_segment(ecg_cleaned, rpeaks=None, sampling_rate=1000, show=False):
         epochs_end=epochs_end,
     )
 
+    # pad last heartbeat with nan so that segments are equal length
+    last_heartbeat_key = str(np.max(np.array(list(heartbeats.keys()), dtype=int)))
+    after_last_index = heartbeats[last_heartbeat_key]["Index"] < len(ecg_cleaned)
+    heartbeats[last_heartbeat_key].loc[after_last_index, "Signal"] = np.nan
+
     if show:
         heartbeats_plot = epochs_to_df(heartbeats)
         heartbeats_pivoted = heartbeats_plot.pivot(index="Time", columns="Label", values="Signal")

diff --git a/neurokit2/eda/eda_autocor.py b/neurokit2/eda/eda_autocor.py
@@ -7,7 +7,7 @@
 def eda_autocor(eda_cleaned, sampling_rate=1000, lag=4):
     """**EDA Autocorrelation**
 
-    Compute autocorrelation measure of raw EDA signal i.e., the correlation between the time
+    Compute the autocorrelation measure of raw EDA signal i.e., the correlation between the time
     series data and a specified time-lagged version of itself.
 
     Parameters
@@ -44,9 +44,10 @@ def eda_autocor(eda_cleaned, sampling_rate=1000, lag=4):
 
     References
     -----------
-    - Halem, S., van Roekel, E., Kroencke, L., Kuper, N., & Denissen, J. (2020). Moments That Matter?
-      On the Complexity of Using Triggers Based on Skin Conductance to Sample Arousing Events Within
-      an Experience Sampling Framework. European Journal of Personality.
+    * van Halem, S., Van Roekel, E., Kroencke, L., Kuper, N., & Denissen, J. (2020). Moments that
+      matter? On the complexity of using triggers based on skin conductance to sample arousing
+      events within an experience sampling framework. European Journal of Personality, 34(5),
+      794-807.
 
     """
     # Sanity checks

diff --git a/neurokit2/eda/eda_changepoints.py b/neurokit2/eda/eda_changepoints.py
@@ -53,9 +53,10 @@ def eda_changepoints(eda_cleaned, penalty=10000, show=False):
 
     References
     -----------
-    * Halem, S., van Roekel, E., Kroencke, L., Kuper, N., & Denissen, J. (2020). Moments That
-      Matter? On the Complexity of Using Triggers Based on Skin Conductance to Sample Arousing
-      Events Within an Experience Sampling Framework. European Journal of Personality.
+    * van Halem, S., Van Roekel, E., Kroencke, L., Kuper, N., & Denissen, J. (2020). Moments that
+      matter? On the complexity of using triggers based on skin conductance to sample arousing
+      events within an experience sampling framework. European Journal of Personality, 34(5),
+      794-807.
 
     """
     # Sanity checks

diff --git a/neurokit2/eda/eda_findpeaks.py b/neurokit2/eda/eda_findpeaks.py
@@ -114,7 +114,6 @@ def eda_findpeaks(eda_phasic, sampling_rate=1000, method="neurokit", amplitude_m
 
 
 def _eda_findpeaks_neurokit(eda_phasic, amplitude_min=0.1):
-
     peaks = signal_findpeaks(eda_phasic, relative_height_min=amplitude_min, relative_max=True)
 
     info = {
@@ -129,27 +128,9 @@ def _eda_findpeaks_neurokit(eda_phasic, amplitude_min=0.1):
 def _eda_findpeaks_vanhalem2020(eda_phasic, sampling_rate=1000):
     """Follows approach of van Halem et al. (2020).
 
-    A peak is considered when there is a consistent increase of 0.5 seconds following a consistent decrease
-    of 0.5 seconds.
-
-    Parameters
-    ----------
-    eda_phasic : array
-        Input filterd EDA signal.
-    sampling_rate : int
-        Sampling frequency (Hz). Defaults to 1000Hz.
+    A peak is considered when there is a consistent increase of 0.5 seconds following a consistent
+    decrease of 0.5 seconds.
 
-    Returns
-    -------
-    onsets : array
-        Indices of the SCR onsets.
-    peaks : array
-        Indices of the SRC peaks.
-    amplitudes : array
-        SCR pulse amplitudes.
-
-    References
-    ----------
     * van Halem, S., Van Roekel, E., Kroencke, L., Kuper, N., & Denissen, J. (2020).
       Moments That Matter? On the Complexity of Using Triggers Based on Skin Conductance to Sample
       Arousing Events Within an Experience Sampling Framework. European Journal of Personality.
@@ -170,7 +151,8 @@ def _eda_findpeaks_vanhalem2020(eda_phasic, sampling_rate=1000):
     threshold = 0.5 * sampling_rate
 
     # Define each peak as a consistent increase of 0.5s
-    peaks = peaks[info["Width"] > threshold]
+    increase = info["Peaks"] - info["Onsets"]
+    peaks = peaks[increase > threshold]
     idx = np.where(peaks[:, None] == info["Peaks"][None, :])[1]
 
     # Check if each peak is followed by consistent decrease of 0.5s
@@ -184,32 +166,16 @@ def _eda_findpeaks_vanhalem2020(eda_phasic, sampling_rate=1000):
     info = {
         "SCR_Onsets": info["Onsets"][idx],
         "SCR_Peaks": info["Peaks"][idx],
-        "SCR_Height": info["Height"][idx],
+        "SCR_Height": eda_phasic[info["Peaks"][idx]],
     }
 
     return info
 
 
 def _eda_findpeaks_gamboa2008(eda_phasic):
-    """Basic method to extract Skin Conductivity Responses (SCR) from an EDA signal following the approach in the thesis
-    by Gamboa (2008).
+    """Basic method to extract Skin Conductivity Responses (SCR) from an EDA signal following the
+    approach in the thesis by Gamboa (2008).
 
-    Parameters
-    ----------
-    eda_phasic : array
-        Input filterd EDA signal.
-
-    Returns
-    -------
-    onsets : array
-        Indices of the SCR onsets.
-    peaks : array
-        Indices of the SRC peaks.
-    amplitudes : array
-        SCR pulse amplitudes.
-
-    References
-    ----------
     * Gamboa, H. (2008). Multi-modal behavioral biometrics based on hci and electrophysiology.
       PhD Thesis Universidade.
 
@@ -255,26 +221,6 @@ def _eda_findpeaks_kim2004(eda_phasic, sampling_rate=1000, amplitude_min=0.1):
     """KBK method to extract Skin Conductivity Responses (SCR) from an EDA signal following the approach by Kim et
     al.(2004).
 
-    Parameters
-    ----------
-    eda_phasic : array
-        Input filterd EDA signal.
-    sampling_rate : int
-        Sampling frequency (Hz). Defaults to 1000Hz.
-    amplitude_min : float
-        Minimum treshold by which to exclude SCRs. Defaults to 0.1.
-
-    Returns
-    -------
-    onsets : array
-        Indices of the SCR onsets.
-    peaks : array
-        Indices of the SRC peaks.
-    amplitudes : array
-        SCR pulse amplitudes.
-
-    References
-    ----------
     * Kim, K. H., Bang, S. W., & Kim, S. R. (2004). Emotion recognition system using short-term
       monitoring of physiological signals. Medical and biological engineering and computing, 42(3),
       419-427.
@@ -322,40 +268,36 @@ def _eda_findpeaks_kim2004(eda_phasic, sampling_rate=1000, amplitude_min=0.1):
 
 
 def _eda_findpeaks_nabian2018(eda_phasic):
-    """Basic method to extract Skin Conductivity Responses (SCR) from an EDA signal following the approach by Nabian et
-    al. (2018).
+    """Basic method to extract Skin Conductivity Responses (SCR) from an EDA signal following the
+    approach by Nabian et al. (2018). The amplitude of the SCR is obtained by finding the maximum
+    value between these two zero-crossings, and calculating the difference between the initial zero
+    crossing and the maximum value. Detected SCRs with amplitudes smaller than 10 percent of the
+    maximum SCR amplitudes that are already detected on the differentiated signal will be
+    eliminated. It is crucial that artifacts are removed before finding peaks.
 
-    Parameters
-    ----------
-    eda_phasic : array
-        Input filterd EDA signal.
-
-    Returns
-    -------
-    onsets : array
-        Indices of the SCR onsets.
-    peaks : array
-        Indices of the SRC peaks.
-    amplitudes : array
-        SCR pulse amplitudes.
-
-    References
-    ----------
-    - Nabian, M., Yin, Y., Wormwood, J., Quigley, K. S., Barrett, L. F., & Ostadabbas, S. (2018). An
-    Open-Source Feature Extraction Tool for the Analysis of Peripheral Physiological Data. IEEE
-    journal of translational engineering in health and medicine, 6, 2800711.
-    https://doi.org/10.1109/JTEHM.2018.2878000
+    * Nabian, M., Yin, Y., Wormwood, J., Quigley, K. S., Barrett, L. F., & Ostadabbas, S. (2018). An
+      Open-Source Feature Extraction Tool for the Analysis of Peripheral Physiological Data. IEEE
+      journal of translational engineering in health and medicine, 6, 2800711.
+      https://doi.org/10.1109/JTEHM.2018.2878000
 
     """
 
+    # differentiation
+    eda_phasic_diff = np.diff(eda_phasic)
+
     # smooth
-    eda_phasic = signal_smooth(eda_phasic, kernel="bartlett", size=20)
+    eda_phasic_smoothed = signal_smooth(eda_phasic_diff, kernel="bartlett", size=20)
 
     # zero crossings
-    pos_crossings = signal_zerocrossings(eda_phasic, direction="positive")
-    neg_crossings = signal_zerocrossings(eda_phasic, direction="negative")
+    pos_crossings = signal_zerocrossings(eda_phasic_smoothed, direction="positive")
+    neg_crossings = signal_zerocrossings(eda_phasic_smoothed, direction="negative")
 
+    # if negative crossing happens before the positive crossing
+    # delete first negative crossing because we want to identify peaks
+    if neg_crossings[0] < pos_crossings[0]:
+        neg_crossings = neg_crossings[1:]
     # Sanitize consecutive crossings
+
     if len(pos_crossings) > len(neg_crossings):
         pos_crossings = pos_crossings[0 : len(neg_crossings)]
     elif len(pos_crossings) < len(neg_crossings):
@@ -366,15 +308,33 @@ def _eda_findpeaks_nabian2018(eda_phasic):
     amps_list = []
     for i, j in zip(pos_crossings, neg_crossings):
         window = eda_phasic[i:j]
-        amp = np.max(window)
+        # The amplitude of the SCR is obtained by finding the maximum value
+        # between these two zero-crossings and calculating the difference
+        # between the initial zero crossing and the maximum value.
+        # amplitude defined in neurokit2
+        amp = np.nanmax(window)
 
         # Detected SCRs with amplitudes less than 10% of max SCR amplitude will be eliminated
-        diff = amp - eda_phasic[i]
-        if not diff < (0.1 * amp):
+        # we append the first SCR
+        if len(amps_list) == 0:
+            # be careful, if two peaks have the same amplitude, np.where will return a list
             peaks = np.where(eda_phasic == amp)[0]
-            peaks_list.append(peaks)
+            # make sure that the peak is within the window
+            peaks = [peak for peak in [peaks] if peak > i and peak < j]
+            peaks_list.append(peaks[0])
             onsets_list.append(i)
             amps_list.append(amp)
+        else:
+            # we have a list of peaks
+            # amplitude defined in the paper
+            diff = amp - eda_phasic[i]
+            if not diff < (0.1 * max(amps_list)):
+                peaks = np.where(eda_phasic == amp)[0]
+                # make sure that the peak is within the window
+                peaks = [peak for peak in [peaks] if peak > i and peak < j]
+                peaks_list.append(peaks[0])
+                onsets_list.append(i)
+                amps_list.append(amp)
 
     # output
     info = {