[Fix] EDA peak detection and sympathetic activity

[GanshengT]

Description

This PR aims at fix the issue mentioned at #778, #780 and #779

Proposed Changes

if the user chooses to 'posada' to calculate Sympathetic Nervous System Index, EDA signal will be resampled at 400Hz. A low-pass filter is also added according to 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.

In the eda.find_peaks sub-module, _eda_findpeaks_vanhalem2020 has been revised. The algorithm is changed based on the orginal paper.
We specified a peak when a consistent increase of 0.5 seconds would follow a consistent decrease of 0.5 seconds.
However, if eda signals ends before the peak recovers, the offset will not be detected.

Checklist

Here are some things to check before creating the PR. If you encounter any issues, do let us know :)

• [+ ] I have read the CONTRIBUTING file.
• [ +] My PR is targetted at the dev branch (and not towards the master branch).
• [ +] I ran the CODE CHECKS on the files I added or modified and fixed the errors.
• [+ ] I have added the newly added features to News.rst (if applicable)

[codecov-commenter]

Codecov Report

Base: 53.28% // Head: 53.77% // Increases project coverage by +0.49% 🎉

Coverage data is based on head (d1e6815) compared to base (49b7d22).
Patch coverage: 81.25% of modified lines in pull request are covered.

Additional details and impacted files

@@            Coverage Diff             @@
##              dev     #783      +/-   ##
==========================================
+ Coverage   53.28%   53.77%   +0.49%     
==========================================
  Files         293      293              
  Lines       13276    13295      +19     
==========================================
+ Hits         7074     7150      +76     
+ Misses       6202     6145      -57     

Impacted Files	Coverage Δ
neurokit2/ecg/ecg_segment.py	66.66% <ø> (ø)
neurokit2/eda/eda_peaks.py	89.79% <ø> (ø)
neurokit2/eda/eda_sympathetic.py	66.03% <75.00%> (+43.81%)	⬆️
neurokit2/eda/eda_intervalrelated.py	77.77% <75.86%> (-5.56%)	⬇️
neurokit2/eda/eda_findpeaks.py	90.16% <94.73%> (+33.40%)	⬆️
neurokit2/signal/signal_fixpeaks.py	73.76% <0.00%> (+0.76%)	⬆️
neurokit2/rsp/rsp_simulate.py	94.64% <0.00%> (+0.89%)	⬆️
neurokit2/eda/eda_eventrelated.py	83.67% <0.00%> (+2.04%)	⬆️
neurokit2/signal/signal_zerocrossings.py	100.00% <0.00%> (+22.22%)	⬆️
... and 2 more

Help us with your feedback. Take ten seconds to tell us how you rate us. Have a feature suggestion? Share it here.

☔ View full report at Codecov.
📢 Do you have feedback about the report comment? Let us know in this issue.

[DominiqueMakowski]

Thanks :) It seems like it throws an error when a short signal is passed. Should we try to adjust nperseg or the overlap so that it doesn't fail in that case? If not,skip when shorter than some minimal length?

import neurokit2 as nk

eda = nk.data('bio_resting_8min_100hz')['EDA']

indexes_posada = nk.eda_sympathetic(eda[0:2000], sampling_rate=100, method='posada', show=True)

ValueError: noverlap must be less than nperseg.

[GanshengT]

Thanks :) It seems like it throws an error when a short signal is passed. Should we try to adjust nperseg or the overlap so that it doesn't fail in that case? If not,skip when shorter than some minimal length?

import neurokit2 as nk

eda = nk.data('bio_resting_8min_100hz')['EDA']

indexes_posada = nk.eda_sympathetic(eda[0:2000], sampling_rate=100, method='posada', show=True)

ValueError: noverlap must be less than nperseg.

The algorithm is proposed for data that is more than $\frac{nperseg = 128}{sampling frequency = 2} = 64 seconds$. I suggest that we skip or maybe return a warning "Input signal should be more than 64 seconds when using posada method".

[DominiqueMakowski]

Seems like there is bug in the nabian method that makes the docstring examples fail:

File D:\a\NeuroKit\NeuroKit\neurokit2\eda\eda_findpeaks.py:373, in _eda_findpeaks_nabian2018(eda_phasic)
    368 window = eda_phasic[i:j]
    369 # The amplitude of the SCR is obtained by finding the maximum value
    370 # between these two zero-crossings and calculating the difference
    371 # between the initial zero crossing and the maximum value.
    372 # amplitude defined in neurokit2
--> 373 amp = np.max(window)
    375 # Detected SCRs with amplitudes less than 10% of max SCR amplitude will be eliminated
    376 # we append the first SCR
    377 if len(amps_list) == 0:
    378     # be careful, if two peaks have the same amplitude, np.where will return a list

File <__array_function__ internals>:200, in amax(*args, **kwargs)

File C:\hostedtoolcache\windows\Python\3.10.10\x64\lib\site-packages\numpy\core\fromnumeric.py:2820, in amax(a, axis, out, keepdims, initial, where)
   2703 @array_function_dispatch(_amax_dispatcher)
   2704 def amax(a, axis=None, out=None, keepdims=np._NoValue, initial=np._NoValue,
   2705          where=np._NoValue):
   2706     """
   2707     Return the maximum of an array or maximum along an axis.
   2708 
   (...)
   2818     5
   2819     """
-> 2820     return _wrapreduction(a, np.maximum, 'max', axis, None, out,
   2821                           keepdims=keepdims, initial=initial, where=where)

File C:\hostedtoolcache\windows\Python\3.10.10\x64\lib\site-packages\numpy\core\fromnumeric.py:86, in _wrapreduction(obj, ufunc, method, axis, dtype, out, **kwargs)
     83         else:
     84             return reduction(axis=axis, out=out, **passkwargs)
---> 86 return ufunc.reduce(obj, axis, dtype, out, **passkwargs)

ValueError: zero-size array to reduction operation maximum which has no identity

[GanshengT]

Could you tell me which example failed?
I tried to run the following lines, no error is reported.
eda_signal = nk.eda_simulate(duration=30, scr_number=5, drift=0.1, noise=0)
eda_cleaned = nk.eda_clean(eda_signal)
eda = nk.eda_phasic(eda_cleaned)
eda_phasic = eda["EDA_Phasic"].values
nabian2018 = nk.eda_findpeaks(eda_phasic, method="nabian2018")

It looks like assert all([i in features_df.columns.values for i in columns]) in test_eda_intervalrelated() is failed
features_df.columns.values: ['SCR_Peaks_N' 'SCR_Peaks_Amplitude_Mean']

columns: ['SCR_Peaks_N', 'SCR_Peaks_Amplitude_Mean', 'Label']

[DominiqueMakowski]

Changing max to nanmax seemed to have fixed the issue 🤷‍♂️

I've also taken the liberty of adding you as a contributor, do let me know if it's okay

[DominiqueMakowski]

Actually this still does fail:

      import neurokit2 as nk

      # Get phasic component
      eda_signal = nk.eda_simulate(duration=30, scr_number=5, drift=0.1, noise=0, sampling_rate=100)
      eda_cleaned = nk.eda_clean(eda_signal, sampling_rate=100)
      eda = nk.eda_phasic(eda_cleaned, sampling_rate=100)
      eda_phasic = eda["EDA_Phasic"].values
      _, nabian2018 = nk.eda_peaks(eda_phasic, sampling_rate=100, method="nabian2018")

Reason is that the windows are empty ([]), not sure how to best fix it

[GanshengT]

The problem happens because a negative crossing happens before the first positive crossing. In this case, delete the first negative crossing would work.

[DominiqueMakowski]

Could you add that or show me how to add?

EDIT: oops I missed that you made a commit already

[DominiqueMakowski]

All good it seems, let me know if it's good to merge :)

[GanshengT]

I agree to merge.

[DominiqueMakowski]

Thanks a lot for diving into these functions, much appreciated