[Bug] RSP Rate Calculation gives strange result - Potentially buggy `_rsp_findpeaks_outliers`

[DerAndereJohannes]

Describe the bug
Using the nk.rsp_rate function on both the trough and xcorr methods gives strange results, especially at the edges of the signal. This leads to the potential bug on _rsp_findpeaks_outliers as this seems to remove the last peak (or trough) regardless as it bases computing vertical distance on a difference function which naturally removes a value.

To Reproduce
Please see below a code snippet and resulting image of how the results of the functions look:

import sys
import neurokit2 as nk
import numpy as np
import matplotlib.pyplot as plt

sampling_rate = 1000
rsp_rate = 10
rsp = nk.rsp_simulate(
    duration=60, sampling_rate=sampling_rate, respiratory_rate=rsp_rate,
    random_state=42, noise=0
)

rsp = nk.rsp_clean(rsp, sampling_rate=sampling_rate)

_, info = nk.rsp_peaks(rsp, sampling_rate=1000)

print(info["RSP_Troughs"])

# Testing alt. test method for troughs (only check the rate between actual data points)
# start, end = info['RSP_Troughs'][0], info['RSP_Troughs'][-1]
# rate_trough = nk.rsp_rate(rsp, troughs=info['RSP_Troughs'], sampling_rate=sampling_rate, method='trough')
# print(rate_trough.shape, rsp.size)
# print(rate_trough[start:end].mean(), rsp_rate, rate_trough[start:end].mean() - rsp_rate)
# print(np.abs(rate_trough[start:end].mean()-rsp_rate))

# Test Troughs
rate_trough = nk.rsp_rate(rsp, troughs=info['RSP_Troughs'], sampling_rate=sampling_rate, method='trough')
print(rate_trough.shape, rsp.size)
print(rate_trough.mean(), rsp_rate, rate_trough.mean() - rsp_rate)
print(np.abs(rate_trough.mean()-rsp_rate))

# Test xcorr
rate_xcorr = nk.rsp_rate(rsp, troughs=info['RSP_Troughs'], sampling_rate=sampling_rate, method='xcorr')
print(rate_xcorr.shape, rsp.size)
print(rate_xcorr.mean(), rsp_rate, rate_xcorr.mean() - rsp_rate)
print(np.abs(rate_xcorr.mean()-rsp_rate))

fig, axs = plt.subplots(2)
fig.suptitle("Comparison between rsp_rate methods")
axs[0].set_title("RSP Signal")
axs[0].plot(rsp)
axs[0].scatter(x=info['RSP_Troughs'], y = [rsp[i] for i in info["RSP_Troughs"]], color='red')
axs[0].scatter(x=info['RSP_Peaks'], y = [rsp[i] for i in info["RSP_Peaks"]], color='green')
axs[0].legend(["rsp signal", "troughs", "peals"], loc="best")
axs[1].set_title("Comparison of Rate Computation methods")
axs[1].plot(rate_xcorr)
axs[1].plot(rate_trough)
axs[1].legend(["xcorr", "trough"], loc="best")
plt.show()

Expected behaviour
The result of RSP Rate should be closer to the actual rsp rate value and the edge breaths should also be accounted for

System Specifications

- OS: Linux (ELF 64bit) 
- Python: 3.13.2 
- NeuroKit2: 0.2.11 

- NumPy: 2.1.3 
- Pandas: 2.2.3 
- SciPy: 1.15.2 
- sklearn: 1.6.1 
- matplotlib: 3.10.1

[welcome]

Hi 👋 Thanks for reaching out and opening your first issue here! We'll try to come back to you as soon as possible. ❤️

[me-pic]

Hello ! @DerAndereJohannes I'm curious to know if you've been working on that issue since I'm also seeing the same behavior in my data (last peak/trough missing). I've also noticed that the first peak is also systematically missing in my data too. Would be happy to collaborate on that if you need help :)

[DerAndereJohannes]

Hi me-pic! This issue mainly stems from trying to get code checks to work and me digging a bit of a rabbit hole when we were working on the 2.11 release (see here: #1022 (comment)).

Since then I have not had the time to look at the problem in more detail, as I have been pretty knee deep with other priorities. I also fear for it to be not so trivial to solve as it also needs to remain consistent as detailed in the response to my comment above. It may also help to explain that missing the first peak seems to be intentional.

I would absolutely be open to collaborate and discuss this issue with you!

[me-pic]

@DerAndereJohannes Thank you for your reply ! I'll first try to familiarize myself a bit more with the rsp_findpeaks.py code :)

I have a question regarding the intentionality of the missing first peak. I'm just trying to understand better the motivation behind it, and also how it is implemented in the code.

For the code portion (just trying to consolidate my understanding here), it seems like _rsp_findpeaks_extrema(), which is used by _rsp_findpeaks_khodadad() and _rsp_findpeaks_biosppy only ignores the first peak if the amplitude of the signal starts above 0, since the function uses the zero crossing points. In the case where the signal crosses zero at the beginning, the first peak is identified but further removed by the _rsp_findpeaks_sanitize() function. While the scipy function applied on the same data might just ignored the first peak depending on the specified prominence value, since it is comparing the neighboring values.

Regarding the motivation behind that choice: I'm certainly missing some background for it , so I'll definitely like to understand it better 🤓 Is there a reason for not computing the breathing amplitude from the first peak to the last trough ? Edit: @DerAndereJohannes just saw your comment on the PR #1022, and I'm also wondering if the _rsp_findpeaks_sanitize() is even needed for rsp_rate() !

Thank you for your time 🚀

[DerAndereJohannes]

Hi, sorry for the late response. It has been some hectic weeks.

My thoughts on missing the first peak is to start at a point where the breath starts (at the trough) rather than relying on a breath that is half way through and we are unsure of when it started. This actually does, after some more thought, seem quite logical to me.

My guess that sanitize is used to make each method consistent and comparable to each other. This does probably change the result a bit compared to the pure implementation, however.

Both these reasons feel to me like a genuine good argument to keep sanitize in the rsp_rate function, however I do feel like it is for these reasons that the rsp_rate also starts creating crazy values outside of the breathing peak and trough ranges (as seen in the first image at the top). I would be interested as to your thoughts on this.

I feel like the first priority would be to look at _rsp_findpeaks_outliers actually as I remember that it had quite a bit of suspicious code in it.

[DerAndereJohannes]

Here is a proposal fix for _rsp_findpeaks_outliers, It seems to fix it for the test that I created in the original post at the top:

def _rsp_findpeaks_outliers(rsp_cleaned, extrema, amplitude_min=0.3):

    while True:
        vertical_diff = np.abs(np.diff(rsp_cleaned[extrema]))
        min_diff = np.median(vertical_diff) * amplitude_min

        if vertical_diff.min() < min_diff:
            extrema = np.delete(extrema, np.array([vertical_diff.argmin(), vertical_diff.argmin() + 1]))
        else:
            break

    amplitudes = rsp_cleaned[extrema]

    return extrema, amplitudes

The new function looks at the extremas in pairs rather than individually and simply removes any pair that is below. The iterative style also makes it less likely to filter out real extremas as we recompute the vertical_diff on every removal. It does not completely fix the rsp_rate results though as you can see below (but at least the peaks / troughs look good now):

Rate trough: 10.718
Rate xcorr: 11.88

[stale]

Is this still relevant? If so, what is blocking it? Is there anything you can do to help move it forward?

This issue has been automatically marked as stale because it has not had recent activity. It will be closed if no further activity occurs.