ENH: Use data-based padding instead of "odd" padding when filtering in raw.plot (#13183)

Co-authored-by: autofix-ci[bot] <114827586+autofix-ci[bot]@users.noreply.github.com>

Author: larsoner

azure-pipelines.yml

@@ -113,7 +113,7 @@ stages:
           - bash: |
               set -e
               python -m pip install --progress-bar off --upgrade pip
-              python -m pip install --progress-bar off "mne-qt-browser[opengl] @ git+https://github.com/mne-tools/mne-qt-browser.git@main" pyvista scikit-learn python-picard qtpy nibabel sphinx-gallery "PySide6!=6.8.0,!=6.8.0.1,!=6.8.1.1" pandas neo pymatreader antio defusedxml
+              python -m pip install --progress-bar off "mne-qt-browser[opengl] @ git+https://github.com/mne-tools/mne-qt-browser.git" pyvista scikit-learn python-picard qtpy nibabel sphinx-gallery "PySide6!=6.8.0,!=6.8.0.1,!=6.8.1.1" pandas neo pymatreader antio defusedxml
               python -m pip uninstall -yq mne
               python -m pip install --progress-bar off --upgrade -e .[test]
             displayName: 'Install dependencies with pip'

doc/changes/devel/13183.bugfix.rst

@@ -0,0 +1 @@
+Fixed bug with filter padding type in :func:`mne.io.Raw.plot` and related functions to reduce edge ringing during data display, by `Eric Larson`_.

mne/cuda.py

@@ -374,10 +374,9 @@ def _smart_pad(x, n_pad, pad="reflect_limited"):
     elif (n_pad < 0).any():
         raise RuntimeError("n_pad must be non-negative")
     if pad == "reflect_limited":
-        # need to pad with zeros if len(x) <= npad
         l_z_pad = np.zeros(max(n_pad[0] - len(x) + 1, 0), dtype=x.dtype)
         r_z_pad = np.zeros(max(n_pad[1] - len(x) + 1, 0), dtype=x.dtype)
-        return np.concatenate(
+        out = np.concatenate(
             [
                 l_z_pad,
                 2 * x[0] - x[n_pad[0] : 0 : -1],
@@ -387,4 +386,8 @@ def _smart_pad(x, n_pad, pad="reflect_limited"):
             ]
         )
     else:
-        return np.pad(x, (tuple(n_pad),), pad)
+        kwargs = dict()
+        if pad == "reflect":
+            kwargs["reflect_type"] = "odd"
+        out = np.pad(x, (tuple(n_pad),), pad, **kwargs)
+    return out

mne/filter.py

@@ -545,16 +545,21 @@ def _iir_filter(x, iir_params, picks, n_jobs, copy, phase="zero"):
         padlen = min(iir_params["padlen"], x.shape[-1] - 1)
         if "sos" in iir_params:
             fun = partial(
-                signal.sosfiltfilt, sos=iir_params["sos"], padlen=padlen, axis=-1
+                _iir_pad_apply_unpad,
+                func=signal.sosfiltfilt,
+                sos=iir_params["sos"],
+                padlen=padlen,
+                padtype="reflect_limited",
             )
             _check_coefficients(iir_params["sos"])
         else:
             fun = partial(
-                signal.filtfilt,
+                _iir_pad_apply_unpad,
+                func=signal.filtfilt,
                 b=iir_params["b"],
                 a=iir_params["a"],
                 padlen=padlen,
-                axis=-1,
+                padtype="reflect_limited",
             )
             _check_coefficients((iir_params["b"], iir_params["a"]))
     else:
@@ -2937,3 +2942,15 @@ def _filt_update_info(info, update_info, l_freq, h_freq):
         ):
             with info._unlock():
                 info["highpass"] = float(l_freq)
+
+
+def _iir_pad_apply_unpad(x, *, func, padlen, padtype, **kwargs):
+    x_out = np.reshape(x, (-1, x.shape[-1])).copy()
+    for this_x in x_out:
+        x_ext = this_x
+        if padlen:
+            x_ext = _smart_pad(x_ext, (padlen, padlen), padtype)
+        x_ext = func(x=x_ext, axis=-1, padlen=0, **kwargs)
+        this_x[:] = x_ext[padlen : len(x_ext) - padlen]
+    x_out.shape = x.shape
+    return x_out

mne/tests/test_filter.py

@@ -97,82 +97,80 @@ def test_estimate_ringing():
         assert estimate_ringing_samples(butter(4, 0.00001)) == 100000
 
 
-def test_1d_filter():
+@pytest.mark.parametrize("n_signal", (1, 2, 3, 5, 10, 20, 40))
+@pytest.mark.parametrize("n_filter", (1, 2, 3, 5, 10, 11, 20, 21, 40, 41, 100, 101))
+@pytest.mark.parametrize("filter_type", ("identity", "random"))
+def test_1d_filter(n_signal, n_filter, filter_type):
     """Test our private overlap-add filtering function."""
     # make some random signals and filters
     rng = np.random.RandomState(0)
-    for n_signal in (1, 2, 3, 5, 10, 20, 40):
-        x = rng.randn(n_signal)
-        for n_filter in (1, 2, 3, 5, 10, 11, 20, 21, 40, 41, 100, 101):
-            for filter_type in ("identity", "random"):
-                if filter_type == "random":
-                    h = rng.randn(n_filter)
-                else:  # filter_type == 'identity'
-                    h = np.concatenate([[1.0], np.zeros(n_filter - 1)])
-                # ensure we pad the signal the same way for both filters
-                n_pad = n_filter - 1
-                x_pad = _smart_pad(x, (n_pad, n_pad))
-                for phase in ("zero", "linear", "zero-double"):
-                    # compute our expected result the slow way
-                    if phase == "zero":
-                        # only allow zero-phase for odd-length filters
-                        if n_filter % 2 == 0:
-                            pytest.raises(
-                                RuntimeError,
-                                _overlap_add_filter,
-                                x[np.newaxis],
-                                h,
-                                phase=phase,
-                            )
-                            continue
-                        shift = (len(h) - 1) // 2
-                        x_expected = np.convolve(x_pad, h)
-                        x_expected = x_expected[shift : len(x_expected) - shift]
-                    elif phase == "zero-double":
-                        shift = len(h) - 1
-                        x_expected = np.convolve(x_pad, h)
-                        x_expected = np.convolve(x_expected[::-1], h)[::-1]
-                        x_expected = x_expected[shift : len(x_expected) - shift]
-                        shift = 0
-                    else:
-                        shift = 0
-                        x_expected = np.convolve(x_pad, h)
-                        x_expected = x_expected[: len(x_expected) - len(h) + 1]
-                    # remove padding
-                    if n_pad > 0:
-                        x_expected = x_expected[n_pad : len(x_expected) - n_pad]
-                    assert len(x_expected) == len(x)
-                    # make sure we actually set things up reasonably
-                    if filter_type == "identity":
-                        out = x_pad.copy()
-                        out = out[shift + n_pad :]
-                        out = out[: len(x)]
-                        out = np.concatenate((out, np.zeros(max(len(x) - len(out), 0))))
-                        assert len(out) == len(x)
-                        assert_allclose(out, x_expected)
-                    assert len(x_expected) == len(x)
-
-                    # compute our version
-                    for n_fft in (None, 32, 128, 129, 1023, 1024, 1025, 2048):
-                        # need to use .copy() b/c signal gets modified inplace
-                        x_copy = x[np.newaxis, :].copy()
-                        min_fft = 2 * n_filter - 1
-                        if phase == "zero-double":
-                            min_fft = 2 * min_fft - 1
-                        if n_fft is not None and n_fft < min_fft:
-                            pytest.raises(
-                                ValueError,
-                                _overlap_add_filter,
-                                x_copy,
-                                h,
-                                n_fft,
-                                phase=phase,
-                            )
-                        else:
-                            x_filtered = _overlap_add_filter(
-                                x_copy, h, n_fft, phase=phase
-                            )[0]
-                            assert_allclose(x_filtered, x_expected, atol=1e-13)
+    x = rng.randn(n_signal)
+    if filter_type == "random":
+        h = rng.randn(n_filter)
+    else:  # filter_type == 'identity'
+        h = np.concatenate([[1.0], np.zeros(n_filter - 1)])
+    # ensure we pad the signal the same way for both filters
+    n_pad = n_filter
+    x_pad = _smart_pad(x, (n_pad, n_pad))
+    for phase in ("zero", "linear", "zero-double"):
+        # compute our expected result the slow way
+        if phase == "zero":
+            # only allow zero-phase for odd-length filters
+            if n_filter % 2 == 0:
+                pytest.raises(
+                    RuntimeError,
+                    _overlap_add_filter,
+                    x[np.newaxis],
+                    h,
+                    phase=phase,
+                )
+                continue
+            shift = (len(h) - 1) // 2
+            x_expected = np.convolve(x_pad, h)
+            x_expected = x_expected[shift : len(x_expected) - shift]
+        elif phase == "zero-double":
+            shift = len(h) - 1
+            x_expected = np.convolve(x_pad, h)
+            x_expected = np.convolve(x_expected[::-1], h)[::-1]
+            x_expected = x_expected[shift : len(x_expected) - shift]
+            shift = 0
+        else:
+            shift = 0
+            x_expected = np.convolve(x_pad, h)
+            x_expected = x_expected[: len(x_expected) - len(h) + 1]
+        # remove padding
+        if n_pad > 0:
+            x_expected = x_expected[n_pad : len(x_expected) - n_pad]
+        assert len(x_expected) == len(x)
+        # make sure we actually set things up reasonably
+        if filter_type == "identity":
+            out = x_pad.copy()
+            out = out[shift + n_pad :]
+            out = out[: len(x)]
+            out = np.concatenate((out, np.zeros(max(len(x) - len(out), 0))))
+            assert len(out) == len(x)
+            assert_allclose(out, x_expected)
+        assert len(x_expected) == len(x)
+
+        # compute our version
+        for n_fft in (None, 32, 128, 129, 1023, 1024, 1025, 2048):
+            # need to use .copy() b/c signal gets modified inplace
+            x_copy = x[np.newaxis, :].copy()
+            min_fft = 2 * n_filter - 1
+            if phase == "zero-double":
+                min_fft = 2 * min_fft - 1
+            if n_fft is not None and n_fft < min_fft:
+                pytest.raises(
+                    ValueError,
+                    _overlap_add_filter,
+                    x_copy,
+                    h,
+                    n_fft,
+                    phase=phase,
+                )
+            else:
+                x_filtered = _overlap_add_filter(x_copy, h, n_fft, phase=phase)[0]
+                assert_allclose(x_filtered, x_expected, atol=1e-13)
 
 
 def test_iir_stability():
@@ -1107,3 +1105,32 @@ def test_filter_minimum_phase_bug():
     dB_min_half = 20 * np.log10(np.abs(H_min_half[mask]))
     assert_array_less(dB_min_half, -20)
     assert not (dB_min_half < -30).all()
+
+
+@pytest.mark.parametrize("dc", (0, 100))
+@pytest.mark.parametrize("sfreq", (1000.0, 999.0))
+def test_smart_pad(dc, sfreq):
+    """Test that smart pad does what it should."""
+    f = 1.0
+    t = (np.arange(0, int(round(sfreq)) + 1)) / sfreq
+    x = np.sin(2 * np.pi * f * t) + dc
+    x_want = np.r_[x, x, x]
+    padlen = (len(x),) * 2
+    x_pad = np.pad(x, padlen, mode="reflect", reflect_type="odd")
+    assert_allclose(x_pad, x_want, err_msg="np.pad", atol=0.1)  # slight DC shift
+    this_x_pad = _smart_pad(x, padlen, "reflect")
+    # slight DC shift from out "x_want" (doubled sample at each end)
+    assert_allclose(this_x_pad, x_want, atol=0.1, err_msg="_smart_pad reflect x_want")
+    assert_allclose(
+        this_x_pad, x_pad, atol=1e-6, err_msg="_smart_pad reflect vs np.pad"
+    )
+    this_x_pad = _smart_pad(x, padlen, "reflect_limited")
+    x_want = np.r_[0, x[:-1], x, x[1:], 0]
+    assert_allclose(
+        this_x_pad, x_want, atol=1e-7, err_msg="_smart_pad reflect_limited x_want"
+    )
+    # reflect_limited uses one fewer sample so ends up a little bit different
+    # with even more padding
+    x_want = np.r_[np.zeros_like(x), x_want, np.zeros_like(x)]
+    x_pad = _smart_pad(x, (len(x) * 2,) * 2, "reflect_limited")
+    assert_allclose(x_pad, x_want, atol=0.1, err_msg="reflect_limited with zeros")

mne/viz/_figure.py

@@ -5,6 +5,7 @@
 # Copyright the MNE-Python contributors.
 
 import importlib
+import inspect
 from abc import ABC, abstractmethod
 from collections import OrderedDict
 from contextlib import contextmanager
@@ -307,25 +308,29 @@ def _make_butterfly_selections_dict(self):
     def _get_start_stop(self):
         # update time
         start_sec = self.mne.t_start - self.mne.first_time
-        stop_sec = start_sec + self.mne.duration
         if self.mne.is_epochs:
             start, stop = np.round(
-                np.array([start_sec, stop_sec]) * self.mne.info["sfreq"]
+                np.array([start_sec, start_sec + self.mne.duration])
+                * self.mne.info["sfreq"]
             ).astype(int)
         else:
+            # ensure our end time includes the last sample
+            disp_duration = (
+                np.ceil(self.mne.duration * self.mne.info["sfreq"])
+                / self.mne.info["sfreq"]
+            )
+            stop_sec = start_sec + disp_duration
             start, stop = self.mne.inst.time_as_index((start_sec, stop_sec))
 
         return start, stop
 
     def _load_data(self, start=None, stop=None):
         """Retrieve the bit of data we need for plotting."""
         if "raw" in (self.mne.instance_type, self.mne.ica_type):
-            # Add additional sample to cover the case sfreq!=1000
-            # when the shown time-range wouldn't correspond to duration anymore
             if stop is None:
                 return self.mne.inst[:, start:]
             else:
-                return self.mne.inst[:, start : stop + 2]
+                return self.mne.inst[:, start:stop]
         else:
             # subtract one sample from tstart before searchsorted, to make sure
             # we land on the left side of the boundary time (avoid precision
@@ -362,9 +367,11 @@ def _apply_filter(self, data, start, stop, picks):
                 )
             data[_picks, _start:_stop] = this_data
 
-    def _process_data(self, data, start, stop, picks, thread=None):
+    def _process_data(self, data, start, stop, picks, thread=None, *, time_slice=None):
         """Update self.mne.data after user interaction."""
         # apply projectors
+        if time_slice is None:
+            time_slice = slice(None)
         if self.mne.projector is not None:
             # thread is the loading-thread only available in Qt-backend
             if thread:
@@ -376,12 +383,13 @@ def _process_data(self, data, start, stop, picks, thread=None):
         if self.mne.remove_dc:
             if thread:
                 thread.processText.emit("Removing DC...")
-            data -= np.nanmean(data, axis=1, keepdims=True)
+            data -= np.nanmean(data[..., time_slice], axis=1, keepdims=True)
         # apply filter
         if self.mne.filter_coefs is not None:
             if thread:
                 thread.processText.emit("Apply Filter...")
             self._apply_filter(data, start, stop, picks)
+        data = data[..., time_slice]
         # scale the data for display in a 1-vertical-axis-unit slot
         if thread:
             thread.processText.emit("Scale Data...")
@@ -400,12 +408,41 @@ def _process_data(self, data, start, stop, picks, thread=None):
 
         return data
 
+    @property
+    def _has_time_slice(self):
+        # check that mne-qt-browser is new enough to support time_slice
+        specs = inspect.getfullargspec(self._process_data)
+        return "time_slice" in specs.kwonlyargs or specs.varkw
+
     def _update_data(self):
         start, stop = self._get_start_stop()
-        # get the data
-        data, times = self._load_data(start, stop)
+        # get the data, with padding if necessary
+        kwargs = dict()
+        padlen = None
+        if isinstance(self.mne.filter_coefs, dict) and self._has_time_slice:  # IIR
+            padlen = self.mne.filter_coefs["padlen"]
+            use_start = max(0, start - padlen)
+            use_stop = min(self.mne.n_times, stop + padlen)
+            # now during filt step, only pad as much as needed
+            self.mne.filter_coefs["padlen"] = max(
+                padlen - (use_stop - stop), padlen - (start - use_start)
+            )
+            time_slice = slice(start - use_start, start - use_start + (stop - start))
+            kwargs["time_slice"] = time_slice
+        else:
+            use_start, use_stop = start, stop
+            time_slice = slice(None)
+
+        data, times = self._load_data(use_start, use_stop)
+        assert data.ndim >= 2 and data.shape[-1] == (use_stop - use_start)
         # process the data
-        data = self._process_data(data, start, stop, self.mne.picks)
+        data = self._process_data(
+            data, use_start, use_stop, picks=self.mne.picks, **kwargs
+        )
+        if padlen is not None:
+            self.mne.filter_coefs["padlen"] = padlen
+        times = times[time_slice]
+        assert data.ndim >= 2 and data.shape[-1] == (stop - start)
         # set the data as attributes
         self.mne.data = data
         self.mne.times = times

tools/circleci_dependencies.sh

@@ -7,6 +7,7 @@ python -m pip install --upgrade --progress-bar off \
     "git+https://github.com/pyvista/pyvista.git" \
     "git+https://github.com/sphinx-gallery/sphinx-gallery.git" \
     "git+https://github.com/mne-tools/mne-bids.git" \
+    "git+https://github.com/mne-tools/mne-qt-browser.git" \
     \
     alphaCSC autoreject bycycle conpy emd fooof meggie \
     mne-ari mne-bids-pipeline mne-faster mne-features \

tools/install_pre_requirements.sh

@@ -44,9 +44,7 @@ echo "OpenMEEG"
 python -m pip install $STD_ARGS --only-binary ":all:" --extra-index-url "https://test.pypi.org/simple" "openmeeg>=2.6.0.dev4"
 
 echo "nilearn"
-# TODO: Revert once settled:
-# https://github.com/scikit-learn/scikit-learn/pull/30268#issuecomment-2479701651
-python -m pip install $STD_ARGS "git+https://github.com/larsoner/nilearn@sklearn"
+python -m pip install $STD_ARGS "git+https://github.com/nilearn/nilearn"
 
 echo "VTK"
 python -m pip install $STD_ARGS --only-binary ":all:" --extra-index-url "https://wheels.vtk.org" vtk