[Docs] window size calculation for locreg detrending

[vasileermicioi]

Description

This PR aims to fix docs for detrending.

Proposed Changes

Doc comments were updated.

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 targeted 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. (not applicable)
• I have added the newly added features to News.rst (not applicable)

[welcome]

Thanks for opening this pull request! We'll make sure it's perfect before merging 🤗
Make sure to read the contributing guide. Also, if you think that your contribution is worthy of it, you can consider adding yourself to the Contributors list (feel free to ask us if you have any doubts).

[codecov-commenter]

Codecov Report

Base: 53.35% // Head: 53.36% // Increases project coverage by +0.00% 🎉

Coverage data is based on head (92058c0) compared to base (a430771).
Patch has no changes to coverable lines.

Additional details and impacted files

@@           Coverage Diff           @@
##              dev     #734   +/-   ##
=======================================
  Coverage   53.35%   53.36%           
=======================================
  Files         283      283           
  Lines       12848    12848           
=======================================
+ Hits         6855     6856    +1     
+ Misses       5993     5992    -1     

Impacted Files	Coverage Δ
neurokit2/signal/signal_detrend.py	39.68% <ø> (ø)
neurokit2/rsp/rsp_simulate.py	94.64% <0.00%> (+0.89%)	⬆️

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.

[vasileermicioi]

The bug was in using a window of 0.5 * sample_rate, which corresponds to a frequency of 2Hz, not to a frequency of 0.5Hz as desired.

A few comments to clarify

1. Any average method is a lowpass filter including simple averages, moving averages, Savitzki-Golay, etc.

2. Formula is window=(sample_rate/frequency).
   E.g. if the sample rate is 1000Hz and we want to keep frequencies below 8Hz, then the window size will be 1000/8=125

3. Heart rate information is in the range of 0.5Hz-5Hz, and, in order not to distort the signal, only filters that remove frequencies below 0.5Hz are used (lowpass filters). If the heart rate is known, then the perfect detrend is below the heart rate frequency (heart rate frequency = heart rate / 60).

4. Any lowpass filter can be used as a detrending method, including Butterworth and Chebyshev, and after the bandpass, the signal is already detrended.
   For e.g. applying a bandpass filter of 0.5hz-15hz will detrend (frequencies below 0.5Hz are removed) and smooth (frequencies above 15Hz are removed) the signal

[vasileermicioi]

Here 2 examples of PPG signals with the heart rate of 60bpm and 120bpm respectively, detrended with the locreg method at 0.5Hz and 2Hz.
You can see that the 2Hz filter distorted the signal with bpm=60,
and the 0.5Hz filter on the signal with bpm=120 left didn't remove the trend completely.

the 60 bpm signal

the 120 bpm signal

[danibene]

Thank you for this @vasileermicioi ! I added some comments above.

• Heart rate information is in the range of 0.5Hz-5Hz, and, in order not to distort the signal, only filters that remove frequencies below 0.5Hz are used (lowpass filters). If the heart rate is known, then the perfect detrend is below the heart rate frequency (heart rate frequency = heart rate / 60).

Do you happen to have a reference for this? I agree that we would not filter out frequencies above 0.5Hz for interbeat interval detrending, but for preprocessing a raw heartbeat signal like ECG, I've seen higher cutoff frequencies e.g. Pan-Tompkins uses 5 Hz:

NeuroKit/neurokit2/ecg/ecg_clean.py

Lines 191 to 196 in e66a732

	def _ecg_clean_pantompkins(ecg_signal, sampling_rate=1000):
	"""Adapted from https://github.com/PIA-
	Group/BioSPPy/blob/e65da30f6379852ecb98f8e2e0c9b4b5175416c3/biosppy/signals/ecg.py#L69."""
	f1 = 5 / (0.5 * sampling_rate)
	f2 = 15 / (0.5 * sampling_rate)

No worries if not, but how the appropriate cutoff frequency is determined is something I've wondered myself!

[vasileermicioi]

heart rate information is in 0.5Hz-5Hz because the heart rate is in 30-300bpm (multiply frequency by 60 :) ),
but usual high-cut for PPG is higher than 5Hz like 8Hz, 10Hz or 15Hz in order to preserve some details like diastolic peaks,
(https://www.researchgate.net/publication/328256630_Sinus_or_not_a_new_beat_detection_algorithm_based_on_a_pulse_morphology_quality_index_to_extract_normal_sinus_rhythm_beats_from_wrist-worn_photoplethysmography_recordings)

for ECG it is 40Hz, see https://www.medteq.net/article/2017/4/1/ecg-filters

[danibene]

Thank you for the references @vasileermicioi !

heart rate information is in 0.5Hz-5Hz because the heart rate is in 30-300bpm (multiply frequency by 60 :) ), but usual high-cut for PPG is higher than 5Hz like 8Hz, 10Hz or 15Hz in order to preserve some details like diastolic peaks, (https://www.researchgate.net/publication/328256630_Sinus_or_not_a_new_beat_detection_algorithm_based_on_a_pulse_morphology_quality_index_to_extract_normal_sinus_rhythm_beats_from_wrist-worn_photoplethysmography_recordings)

for ECG it is 40Hz, see https://www.medteq.net/article/2017/4/1/ecg-filters

I meant that the low-cut of Pan-Tompkins, an algorithm used to process ECG signals, is 5 Hz (the high-cut is 15 Hz), so I was curious if you had any sources for digitally filtering out frequencies below 5 Hz being necessarily problematic?

[vasileermicioi]

I meant that the low-cut of Pan-Tompkins, an algorithm used to process ECG signals, is 5 Hz (the high-cut is 15 Hz), so I was curious if you had any sources for digitally filtering out frequencies below 5 Hz being necessarily problematic?

you are right about ECG having a different range for bandpass than PPG, my experience and example are with PPG, and my statement is true for PPG only,
https://www.elgendi.net/papers/QRS_Bands_final.pdf
this article compares ECG bandpass bands and suggests 8Hz-20Hz being optimal,
while for PPG, it is usually 0.5Hz-8Hz

other than that the PR is still correct as

• docs for the window calculation for the locreg method is computed incorrectly, it should be window = sample_rate / freq, and not window=sample_rate * freq
• for consistency as the other 2 methods (polynomial and standardized) are using 0.5Hz for lowcut, so it should locreg
• using lowcut frequencies above 0.5Hz for PPG will distort the signal (see my example)

PS:
At first, I thought the polynomial is wrong too, but using an additional order for every 2 seconds (in the code is length / 2), has the same effect as locreg with window=2*sample_rate (which is the same as window=sample_rate/0.5, just in case :) )

[danibene]

you are right about ECG having a different range for bandpass than PPG, my experience and example are with PPG, and my statement is true for PPG only, https://www.elgendi.net/papers/QRS_Bands_final.pdf this article compares ECG bandpass bands and suggests 8Hz-20Hz being optimal, while for PPG, it is usually 0.5Hz-8Hz

Ah interesting, thank you :)

other than that the PR is still correct as

docs for the window calculation for the locreg method is computed incorrectly, it should be window = sample_rate / freq, and not window=sample_rate * freq
Were you able to compare with the following implementation? https://github.com/sappelhoff/pyprep/blob/master/pyprep/removeTrend.py#L92 They use window = (1.5/freq)*sample_rate not window = sample_rate / freq

for consistency as the other 2 methods (polynomial and standardized) are using 0.5Hz for lowcut, so it should locreg

Agreed, just want to make sure we compute that correctly, and I think since sampling_rate is now an argument for signal_detrend() we should include sampling_rate as an input rather than computing the window and stepsize with sampling_rate.

[danibene]

@vasileermicioi update: I tried out various definitions of the window with a simulated signal: https://gist.github.com/danibene/e3ed1200a75ea72a3c52bfb0ac248a2c

Now I'm also wondering if the pyprep implementation that I linked earlier is incorrect. Any insights you have here would be appreciated!

In any case, I agree that the documentation should be updated so that it's clear that the lowcut is not equivalent to the window.

[vasileermicioi]

Now I'm also wondering if the pyprep implementation that I linked earlier is incorrect. Any insights you have here would be appreciated!

yes, it is wrong too, because taking a window with 50% larger will have the effect of reducing the frequency by 1/3, for example if we want to remove all frequencies below 0.5Hz, it will remove all below 0.33Hz instead, and if we want to remove below 2Hz, it will remove only below 1.34Hz

[vasileermicioi]

@danibene now the window parameter is not a window anymore, but is the detrend frequency, as (sample_rate / freq) / sample_rate is the same as freq, so perhaps the parameter name (window) should be changed to detrend_frequency

[danibene]

@danibene now the window parameter is not a window anymore, but is the detrend frequency, as (sample_rate / freq) / sample_rate is the same as freq, so perhaps the parameter name (window) should be changed to detrend_frequency

Wouldn't it be 1/freq? But we could change the implementation so that it's freq

[vasileermicioi]

@danibene now the window parameter is not a window anymore, but is the detrend frequency, as (sample_rate / freq) / sample_rate is the same as freq, so perhaps the parameter name (window) should be changed to detrend_frequency

Wouldn't it be 1/freq? But we could change the implementation so that it's freq

yes, you are right

[danibene]

@DominiqueMakowski any thoughts on whether we should deviate from the pyprep implementation?

[DominiqueMakowski]

pyprep is supposedly pretty robust, perhaps it's worth opening an issue on their repo to get their thoughts?

[danibene]

pyprep is supposedly pretty robust, perhaps it's worth opening an issue on their repo to get their thoughts?

So I think we can hold off on merging for a bit to see if they respond here: sappelhoff/pyprep#125

[vasileermicioi]

@danibene
I am not sure I understand completely the answer, but to me, it sounds like "you are right, but 0.67Hz is a good frequency for detrending"
VisLab/EEG-Clean-Tools#33

[danibene]

@vasileermicioi I'm away for the next week but I can take a look when I get back! Otherwise @DominiqueMakowski feel free to address this without my input

[DominiqueMakowski]

you are right, but 0.67Hz is a [reasonable] frequency for EEG detrending

That's my understanding as well..., but we can wait 'til Dani comes back :)

[danibene]

@vasileermicioi based on the responses we received: VisLab/EEG-Clean-Tools#33 (comment) it looks like you're correct, but they don't plan to change the implementation.

How about we keep the window parameter as is so that our implementation is still based off of prep's, but just update our documentation that window should be 1/cutoff? If you agree could you please make those changes? Then I think we can merge.

[DominiqueMakowski]

Yeah that seem like a good resolve, this way we don't change the default as compared to prep but provide more context and info for users to potentially fine tune it

[DominiqueMakowski]

Maybe we can make the docs change and close the PR (I'd like to merge the dev branch soon :)

[DominiqueMakowski]

Maybe we can make the docs change and close the PR

I was thinking "merge" the PR actually 😢 @vasileermicioi could you reopen one?