Keep invalid samples invalid in smooth_frames (#608)

When a signal has more than one sample per frame, `smooth_frames`
averages each frame in digital units. It doesn't check for the invalid
sample value, so an invalid sample is averaged with the others and gives
a large value that looks valid:

```python
import numpy as np, wfdb

# fmt 16, 2 samples per frame, gain 200. -32768 is the invalid sample value.
np.array([0, 1, -1, -2, -32768, 100], dtype="<i2").tofile("r.dat")
open("r.hea", "w").write("r 1 100 3\nr.dat 16x2 200/mV 16 0 0 0 0 sig\n")

wfdb.rdrecord("r", smooth_frames=False).e_p_signal[0]
# [ 0.     0.005 -0.005 -0.01     nan  0.5  ]
wfdb.rdrecord("r").p_signal[:, 0]
# [ 0.    -0.005 -81.67]   (last frame should be NaN)
```

The WFDB library returns an invalid value for any frame that contains an
invalid sample (`meansamp()` in `lib/signal.c`), so this does the same:
after averaging a chunk of frames in digital units, frames with an
invalid sample are set to the invalid value for the channel's format.
That then becomes NaN when converted to physical units. Physical
smoothing (`sigtype="physical"`) already gave NaN, since NaN propagates
through the sum.

This doesn't change the rounding of valid frames, which is #335.

Tests: a fmt 16 record with two samples per frame and one invalid sample
reads back as `[3, -32768, 7]` (digital) and `[0.015, nan, 0.035]`
(physical). This fails on `main`. `pytest -m "not network"` passes (97
tests), and the docs build passes.

Author: tompollard

docs/changes.rst

@@ -17,6 +17,9 @@ Unreleased
 **Add timeouts to requests for remote files**
   Requests for files on PhysioNet (and other web servers) now time out if the server does not respond (10 seconds to connect, and 60 seconds to wait for data), instead of waiting indefinitely. A ``requests.exceptions.Timeout`` (``ConnectTimeout`` or ``ReadTimeout``) is raised in this case.
 
+**Fix smoothing of frames with invalid samples**
+  When a signal with more than one sample per frame is read with ``smooth_frames=True`` (the default), a frame that contains an invalid sample now gives an invalid value (NaN in physical units), as in the WFDB library. The invalid sample value used to be averaged with the other samples, which gave a large but valid-looking value.
+
 **Build with uv_build**
   The package is now built with ``uv_build`` instead of ``hatchling``. The version number is specified in ``pyproject.toml`` only, and ``wfdb.__version__`` is read from the installed package metadata (it is ``"unknown"`` if the package is run from a source tree without being installed).
 

src/wfdb/io/_signal.py

@@ -1060,15 +1060,20 @@ def smooth_frames(self, sigtype="physical"):
             if spf[ch] == 1:
                 signal[:, ch] = expanded_signal[ch]
             else:
+                # As in the WFDB library, a frame with an invalid sample gives
+                # an invalid value. Physical samples are already NaN.
+                invalid = None
+                if sigtype == "digital" and self.fmt and self.fmt[ch]:
+                    invalid = _digi_nan(self.fmt[ch])
                 frames = expanded_signal[ch].reshape(-1, spf[ch])
                 for chunk_start in range(0, sig_len, CHUNK_SIZE):
                     chunk_end = chunk_start + CHUNK_SIZE
-                    signal_sum = np.sum(
-                        frames[chunk_start:chunk_end, :],
-                        axis=1,
-                        dtype=intermediate_dtype,
-                    )
+                    chunk = frames[chunk_start:chunk_end, :]
+                    signal_sum = np.sum(chunk, axis=1, dtype=intermediate_dtype)
                     signal[chunk_start:chunk_end, ch] = signal_sum / spf[ch]
+                    if invalid is not None:
+                        has_invalid = np.any(chunk == invalid, axis=1)
+                        signal[chunk_start:chunk_end, ch][has_invalid] = invalid
 
         return signal
 

tests/test_record.py

@@ -857,6 +857,22 @@ def test_write_smoothed(self, temp_path):
         )
         np.testing.assert_array_equal(record.d_signal, record2.d_signal)
 
+    def test_smooth_frames_invalid_sample(self, temp_path):
+        """A frame with an invalid sample is invalid after smoothing."""
+        # Two samples per frame. -32768 is the invalid sample value for fmt 16.
+        raw = np.array([2, 4, -32768, 100, 6, 8], dtype="<i2")
+        raw.tofile(os.path.join(temp_path, "invalid.dat"))
+        with open(os.path.join(temp_path, "invalid.hea"), "w") as f:
+            f.write("invalid 1 100 3\n")
+            f.write("invalid.dat 16x2 200/mV 16 0 0 0 0 sig\n")
+        record_name = os.path.join(temp_path, "invalid")
+
+        record = wfdb.rdrecord(record_name, physical=False)
+        np.testing.assert_array_equal(record.d_signal[:, 0], [3, -32768, 7])
+
+        record = wfdb.rdrecord(record_name)
+        np.testing.assert_allclose(record.p_signal[:, 0], [0.015, np.nan, 0.035])
+
     def test_to_dataframe(self):
         record = wfdb.rdrecord("sample-data/test01_00s")
         df = record.to_dataframe()