Incorrect HR interpolation using systole.utils.heart_rate() when specifying sfreq different from 1000 Hz #75
Description
Lines 203 to 207 in e515392
Dear Nico (and Systole maintainers),
While using the function systole.utils.heart_rate() to interpolate the continuous heart rate time series from RR intervals, I encountered an issue whenever the argument sfreq is explicitly set to any value other than the default 1000 Hz. Although specifying sfreq is not essential for interpolating HR from RR intervals, using a different value without knowing can lead to incorrectly interpolated HR values.
This issue arises from the code lines reported above, mainly the calculation in the function:
heartrate = (rr / sfreq) * 1000
By this point in the function, rr should already be converted into RR intervals in ms, regardless of the input type (peaks, peaks_idx, rr_s or rr_ms). Dividing by sfreq incorrectly scales the HR values, causing them to be proportional to sfreq itself.
Reproducing the issue:
You can replicate this issue using the provided RR intervals in ms from import_rr().
from systole import import_rr
rr_ms = import_rr()
rr_ms = rr_ms.iloc[:,1]
rr_interpol_bpm, rr_interpol_time = heart_rate(x=rr_ms, unit='bpm', sfreq=500, kind='cubic', input_type='rr_ms')
Setting the argument sfreq=500 leads to incorrect HR (in bpm) values that are exactly half of the expected values. For example, the first RR interval is 992 ms, which should correspond to subsequently interpolated HR near to 60.48 bpm. Instead, the function returns rr_interpol_bpm to be:
array([30.24193548, 30.23581587, 30.22972094, 30.22365068, 30.21760504,
30.21158398, 30.20558748, 30.19961549, 30.19366797, 30.1877449 ,
30.18184623, 30.17597194, 30.17012197, 30.1642963 , 30.15849489,
30.15271771, 30.14696472, 30.14123587, 30.13553115, 30.1298505 ])
Of course, when sfreq=1000 (or when it is left to its default), the issue does not occur because sfreq cancels out with the multiplication by 1000, effectively returning heartrate = rr.
Proposed solution:
Removing the sfreq scaling when defining heartrate should potentially solve the issue:
# R-R intervals (in miliseconds)
heartrate = rr
if unit == "bpm":
# Beats per minutes
heartrate = 60000 / rr
Looking forward to your thoughts on this.
Thanks for maintaining this great package!