diff --git a/docs/changes.rst b/docs/changes.rst
index 8c63c074..445341cf 100644
--- a/docs/changes.rst
+++ b/docs/changes.rst
@@ -41,6 +41,9 @@ Unreleased
 **Fix resample_sig for signals starting with NaN**
   ``resample_sig`` (and the functions that use it, such as ``resample_singlechan`` and ``resample_multichan``) returned only NaN values if the signal started with NaN. Leading NaN values now take the first valid value before resampling, as trailing ones already took the last.
 
+**Keep missing values missing in resample_sig**
+  ``resample_sig`` (and ``resample_singlechan`` and ``resample_multichan``) now return NaN for output samples with a NaN input value within half an output sample spacing (or within one input sample, when upsampling). NaN values are still filled before resampling, so the rest of the signal is unaffected, but the filled values are no longer returned as if they were real data.
+
 **Reject aux notes that can't be written in wrann**
   ``wrann`` now raises a ``ValueError`` if an ``aux_note`` is longer than 255 characters or contains characters outside Latin-1, instead of writing a corrupt annotation file. A longer note used to be read back as many extra annotations, and other characters were silently replaced.
 
diff --git a/src/wfdb/processing/basic.py b/src/wfdb/processing/basic.py
index 67ad8fdf..5f8eb42d 100644
--- a/src/wfdb/processing/basic.py
+++ b/src/wfdb/processing/basic.py
@@ -42,7 +42,9 @@ def resample_sig(x, fs, fs_target):
     Returns
     -------
     resampled_x : ndarray
-        Array of the resampled signal values.
+        Array of the resampled signal values. Samples with a NaN value in
+        `x` within half an output sample spacing (or within one input
+        sample, when upsampling) are NaN.
     resampled_t : ndarray
         Array of the resampled signal locations.
 
@@ -55,12 +57,23 @@ def resample_sig(x, fs, fs_target):
     new_length = int(x.shape[0] * fs_target / fs)
     # Fill NaN values before resampling. Leading and trailing NaN values take
     # the nearest valid value, since interpolate() leaves leading ones as NaN.
-    if np.isnan(x).any():
+    nan_mask = np.isnan(x)
+    if nan_mask.any():
         x = pd.Series(x.reshape((-1,))).interpolate(limit_direction="both").values
     resampled_x, resampled_t = signal.resample(x, num=new_length, t=t)
     assert resampled_x.shape == resampled_t.shape and resampled_x.shape[0] == new_length
     assert np.all(np.diff(resampled_t) > 0)
 
+    # The filled values aren't real data, so mark output samples with a NaN
+    # input sample within half an output spacing (and at least one input
+    # sample, so neighbours count when upsampling) as NaN again.
+    if nan_mask.any():
+        half = max(fs / fs_target / 2, 1)
+        n_nan = np.concatenate([[0], np.cumsum(nan_mask.reshape((-1,)))])
+        lo = np.clip(np.floor(resampled_t - half).astype(int) + 1, 0, len(t))
+        hi = np.clip(np.ceil(resampled_t + half).astype(int), 0, len(t))
+        resampled_x[n_nan[hi] - n_nan[lo] > 0] = np.nan
+
     return resampled_x, resampled_t
 
 
diff --git a/tests/test_processing.py b/tests/test_processing.py
index 60bb7bd0..50990271 100644
--- a/tests/test_processing.py
+++ b/tests/test_processing.py
@@ -38,26 +38,51 @@ def test_resample_multi(self):
         assert new_sig.shape[1] == sig.shape[1]
 
     def test_resample_leading_nan(self):
-        """Leading NaN values take the first valid value before resampling."""
+        """Leading NaN values are filled before resampling, then marked as NaN."""
         x = np.array([np.nan, np.nan, 1.0, 2.0, 3.0, 4.0])
         filled = np.array([1.0, 1.0, 1.0, 2.0, 3.0, 4.0])
 
         new_x, _ = processing.resample_sig(x, fs=2, fs_target=4)
         expected, _ = processing.resample_sig(filled, fs=2, fs_target=4)
 
-        assert np.isfinite(new_x).all()
+        # Output positions 0 to 1.5 are next to a NaN input sample
+        expected[:4] = np.nan
         np.testing.assert_allclose(new_x, expected)
 
     def test_resample_trailing_and_inner_nan(self):
-        """Inner NaN values are interpolated, and trailing ones take the last."""
+        """Inner and trailing NaN values are filled, then marked as NaN."""
         x = np.array([1.0, np.nan, 3.0, 4.0, np.nan, np.nan])
         filled = np.array([1.0, 2.0, 3.0, 4.0, 4.0, 4.0])
 
         new_x, _ = processing.resample_sig(x, fs=2, fs_target=4)
         expected, _ = processing.resample_sig(filled, fs=2, fs_target=4)
 
+        # Positions 0.5 to 1.5 and 3.5 onwards are next to a NaN input sample
+        expected[[1, 2, 3]] = np.nan
+        expected[7:] = np.nan
         np.testing.assert_allclose(new_x, expected)
 
+    def test_resample_marks_nan(self):
+        """Output samples at or next to NaN input samples are NaN."""
+        x = np.arange(12.0)
+        x[6] = np.nan
+
+        new_x, new_t = processing.resample_sig(x, fs=4, fs_target=8)
+
+        assert np.isnan(new_x).tolist() == ((new_t > 5) & (new_t < 7)).tolist()
+
+    @pytest.mark.parametrize(
+        "gap, expected", [([500], [50]), (list(range(501, 509)), [50, 51])]
+    )
+    def test_resample_marks_nan_downsampling(self, gap, expected):
+        """Gaps between two output samples are still marked when downsampling."""
+        x = np.sin(np.arange(1000) / 50)
+        x[gap] = np.nan
+
+        new_x, _ = processing.resample_sig(x, fs=1000, fs_target=100)
+
+        assert np.flatnonzero(np.isnan(new_x)).tolist() == expected
+
     def test_normalize_bound(self):
         sig, _ = wfdb.rdsamp("sample-data/100")
         lb = -5