Use the last RR interval in compute_hr (#601)

`compute_hr` loops over `range(len(qrs_inds) - 2)` and then copies the
rate at the last beat forward (`heart_rate[qrs_inds[-1]:] =
heart_rate[qrs_inds[-1]]`). So the interval between the last two beats
is never used:

- `compute_hr(14, [0, 2, 6, 12], fs=60)` gives 900 bpm after sample 12
(from the 2-sample interval before it) instead of 600 bpm (from the last
interval, 6 samples).
- `compute_hr(10, [2, 6], fs=1)` returns only NaN, although the
`len(qrs_inds) < 2` check implies two beats are enough.

The loop now runs over every RR interval, and the last one fills up to
the end of the signal. Everything before the last beat is unchanged. For
record 100 with the reference annotations, the only differences are the
8 samples after the last beat.

Two tests are added in `TestComputeHr`, and both fail on `main`. `ruff
check` is clean and `pytest -m "not network"` passes.

Related to #185, but that report was mostly caused by non-beat
annotations (see the comments there), so this doesn't close it.

Author: tompollard

docs/changes.rst

@@ -29,6 +29,9 @@ Unreleased
 **Fix parameter handling in csv_to_wfdb**
   ``csv_to_wfdb`` now builds the same ``Record`` whether it writes files or returns the record (``record_only=True``), so all header parameters (such as ``dat_file_name``, ``counter_freq``, ``base_counter``, ``skew``, ``byte_offset``, ``adc_res`` and ``adc_zero``) are written instead of being ignored. Single values of any numeric type (including NumPy scalars and floats) are used for all channels, and an explicit ``0`` is no longer replaced by the default. The initial value and checksum are computed from the digital samples, and the default ADC resolution now follows ``fmt`` in both modes. The ``samps_per_frame``, ``init_value`` and ``checksum`` arguments are deprecated and ignored. The parameter descriptions in the docstring have been corrected.
 
+**Use the last RR interval in compute_hr**
+  ``compute_hr`` now uses the interval between the last two beats for the heart rate after the last beat, instead of repeating the previous rate. With exactly two beats it now returns their rate after the second beat instead of only NaN.
+
 **Use Ruff for formatting and linting**
   The ``dev`` extra now includes ``ruff`` instead of ``black`` and ``pylint``, and the code is formatted with a line length of 88 (the default).
 

src/wfdb/processing/hr.py

@@ -9,6 +9,10 @@
 def compute_hr(sig_len, qrs_inds, fs):
     """Compute instantaneous heart rate from peak indices.
 
+    The rate from each RR interval is assigned to the samples after its second
+    beat, up to and including the next beat. After the last beat, the rate
+    from the last RR interval is used up to the end of the signal.
+
     Parameters
     ----------
     sig_len : int
@@ -31,15 +35,15 @@ def compute_hr(sig_len, qrs_inds, fs):
     if len(qrs_inds) < 2:
         return heart_rate
 
-    for i in range(0, len(qrs_inds) - 2):
+    # The rate from each RR interval holds from just after its second beat up
+    # to the next beat, or to the end of the signal after the last one.
+    for i in range(0, len(qrs_inds) - 1):
         a = qrs_inds[i]
         b = qrs_inds[i + 1]
-        c = qrs_inds[i + 2]
+        end = qrs_inds[i + 2] + 1 if i + 2 < len(qrs_inds) else sig_len
         rr = (b - a) * (1.0 / fs) * 1000
         hr = 60000.0 / rr
-        heart_rate[b + 1 : c + 1] = hr
-
-    heart_rate[qrs_inds[-1] :] = heart_rate[qrs_inds[-1]]
+        heart_rate[b + 1 : end] = hr
 
     return heart_rate
 

tests/test_processing.py

@@ -174,6 +174,37 @@ def test_correct_peaks(self):
         assert np.array_equal(yz, expected_peaks)
 
 
+class TestComputeHr:
+    """Tests for compute_hr."""
+
+    def test_uses_last_rr_interval(self):
+        """The rate after the last beat comes from the last RR interval."""
+        hr = processing.compute_hr(sig_len=14, qrs_inds=[0, 2, 6, 12], fs=60)
+
+        expected = np.full(14, np.nan)
+        expected[3:7] = 60 * 60 / 2
+        expected[7:13] = 60 * 60 / 4
+        expected[13:] = 60 * 60 / 6
+        np.testing.assert_allclose(hr, expected)
+
+    def test_two_beats(self):
+        """Two beats give one RR interval, so the rate is known after them."""
+        hr = processing.compute_hr(sig_len=10, qrs_inds=[2, 6], fs=1)
+
+        expected = np.full(10, np.nan)
+        expected[7:] = 60 / 4
+        np.testing.assert_allclose(hr, expected)
+
+    def test_last_beat_past_end(self):
+        """A last beat beyond the signal doesn't raise an IndexError."""
+        hr = processing.compute_hr(sig_len=8, qrs_inds=[0, 2, 6, 12], fs=60)
+
+        expected = np.full(8, np.nan)
+        expected[3:7] = 60 * 60 / 2
+        expected[7:] = 60 * 60 / 4
+        np.testing.assert_allclose(hr, expected)
+
+
 class TestQrs:
     """Testing QRS detectors"""