Merge pull request #680 from danibene/feature/hrv_utils_interpolation…

…_rate

[Feature] Option to change interpolation rate of interbeat intervals for HRV frequency

neurokit2/hrv/hrv_frequency.py

@@ -24,6 +24,7 @@ def hrv_frequency(
     silent=True,
     normalize=True,
     order_criteria=None,
+    interpolation_rate=100,
     **kwargs
 ):
     """**Computes frequency-domain indices of Heart Rate Variability (HRV)**
@@ -61,8 +62,7 @@ def hrv_frequency(
         Can be a list of indices or the output(s) of other functions such as :func:`.ecg_peaks`,
         :func:`.ppg_peaks`, :func:`.ecg_process` or :func:`.bio_process`.
     sampling_rate : int, optional
-        Sampling rate (Hz) of the continuous cardiac signal in which the peaks occur. Should be at
-        least twice as high as the highest frequency in vhf. By default 1000.
+        Sampling rate (Hz) of the continuous cardiac signal in which the peaks occur.
     ulf : tuple, optional
         Upper and lower limit of the ultra-low frequency band. By default (0, 0.0033).
     vlf : tuple, optional
@@ -86,6 +86,9 @@ def hrv_frequency(
     order_criteria : str
         The criteria to automatically select order in parametric PSD (only used for autoregressive
         (AR) methods such as ``"burg"``). Defaults to ``None``.
+    interpolation_rate : int, optional
+        Sampling rate (Hz) of the interpolated interbeat intervals. Should be at least twice as 
+        high as the highest frequency in vhf. By default 100. To replicate Kubios defaults, set to 4. 
     **kwargs
         Additional other arguments.
 
@@ -160,9 +163,9 @@ def hrv_frequency(
     if isinstance(peaks, tuple):  # Detect actual sampling rate
         peaks, sampling_rate = peaks[0], peaks[1]
 
-    # Compute R-R intervals (also referred to as NN) in milliseconds (interpolated at 1000 Hz by default)
+    # Compute R-R intervals (also referred to as NN) in milliseconds (interpolated at 4 Hz by default)
     rri, sampling_rate = _hrv_get_rri(
-        peaks, sampling_rate=sampling_rate, interpolate=True, **kwargs
+        peaks, sampling_rate=sampling_rate, interpolate=True, interpolation_rate=interpolation_rate, **kwargs
     )
 
     frequency_band = [ulf, vlf, lf, hf, vhf]

neurokit2/hrv/hrv_utils.py

@@ -1,32 +1,49 @@
 # -*- coding: utf-8 -*-
+from warnings import warn
+
 import numpy as np
 import pandas as pd
 
 from ..signal import signal_interpolate
+from ..misc import NeuroKitWarning
 
 
-def _hrv_get_rri(peaks=None, sampling_rate=1000, interpolate=False, **kwargs):
+def _hrv_get_rri(peaks=None, sampling_rate=1000, interpolate=False, interpolation_rate=100, **kwargs):
 
     rri = np.diff(peaks) / sampling_rate * 1000
 
     if interpolate is False:
-        sampling_rate = None
+        interpolation_rate = None
 
     else:
-
-        # Sanitize minimum sampling rate for interpolation to 10 Hz
-        sampling_rate = max(sampling_rate, 10)
-
-        # Compute length of interpolated heart period signal at requested sampling rate.
-        desired_length = int(np.rint(peaks[-1]))
+        # Rate should be at least 1 Hz (due to Nyquist & frequencies we are interested in)
+        # We considered an interpolation rate 4 Hz by default to match Kubios
+        # but in case of some applications with high heart rates we decided to make it 100 Hz
+        # See https://github.com/neuropsychology/NeuroKit/pull/680 for more information 
+        # and if you have any thoughts to contribute, please let us know!
+        if interpolation_rate < 1:
+                warn(
+                    "The interpolation rate of the R-R intervals is too low for "
+                    " computing the frequency-domain features."
+                    " Consider increasing the interpolation rate to at least 1 Hz.",
+                    category=NeuroKitWarning,
+                )
+
+        # Compute x-values of interpolated heart period signal at requested sampling rate.
+        x_new = np.arange(
+            start=peaks[1] / sampling_rate,
+            stop=peaks[-1] / sampling_rate + 1 / interpolation_rate,
+            step=1 / interpolation_rate,
+        )
+        rri_time = peaks[1:] / sampling_rate # Skip first peak since it has no corresponding element in heart_period
 
         rri = signal_interpolate(
-            peaks[1:],  # Skip first peak since it has no corresponding element in heart_period
+            rri_time,
             rri,
-            x_new=np.arange(desired_length),
+            x_new=x_new,
             **kwargs
         )
-    return rri, sampling_rate
+    return rri, interpolation_rate
 
 
 def _hrv_sanitize_input(peaks=None):