Add a causal half-Gaussian option to the Gaussian smoothing filter

src/ezmsg/sigproc/gaussiansmoothing.py

@@ -33,37 +33,76 @@ class GaussianSmoothingSettings(FilterBaseSettings):
     """
     kernel_size : int | None
         Length of the kernel in samples. If provided, overrides automatic calculation.
+        In causal mode this is the number of *causal* taps, i.e. the kernel spans
+        ``kernel_size`` samples into the past rather than ``kernel_size // 2``.
+    """
+
+    causal: bool = False
+    """
+    causal : bool
+        If False (default), the kernel is a symmetric Gaussian of
+        ``2 * width * sigma + 1`` taps. Filtering is applied causally (``lfilter``),
+        so the acausal half of the kernel manifests purely as group delay of
+        ``(kernel_size - 1) / 2 == width * sigma`` samples -- 4 * sigma at the
+        default ``width=4``.
+
+        If True, the kernel is the causal half of that Gaussian (peak at lag 0,
+        tail extending only into the past), renormalized to unit sum. Its group
+        delay is the centroid of a half-Gaussian, ``sigma * sqrt(2 / pi)``
+        (~0.8 * sigma), i.e. roughly a factor of 5 less lag than the symmetric
+        kernel at the same sigma.
+
+        The two modes are **not** interchangeable at equal sigma: halving the
+        kernel also halves the effective averaging window, so the causal kernel
+        smooths less and its stopband rolls off less steeply (-12 dB/octave
+        versus the symmetric kernel's much sharper Gaussian rolloff) for a given
+        sigma. Compare them at matched white-noise variance reduction
+        (``sum(b ** 2)``) rather than at matched sigma; on that footing the
+        causal kernel reaches the same noise gain at roughly a third of the lag.
+        For example, at 100 Hz a symmetric sigma of 20 ms gives a noise gain of
+        0.141 for 80 ms of delay, while a causal sigma of 38 ms gives the same
+        0.141 for 27 ms.
     """
 
 
 def gaussian_smoothing_filter_design(
     sigma: float = 1.0,
     width: int = 4,
     kernel_size: int | None = None,
+    causal: bool = False,
 ) -> BACoeffs | None:
     """Design a normalized Gaussian FIR kernel. ``sigma`` is in **samples**;
-    callers with a time-domain sigma must scale by the sampling rate first."""
+    callers with a time-domain sigma must scale by the sampling rate first.
+
+    If ``causal`` is True, only the causal half of the Gaussian is kept -- the
+    peak sits at lag 0 and the tail extends into the past -- and ``kernel_size``
+    counts causal taps. See :class:`GaussianSmoothingSettings` for the group
+    delay of each mode.
+    """
     # Parameter checks
     if sigma <= 0:
         raise ValueError(f"sigma must be positive. Received: {sigma}")
 
     if width <= 0:
         raise ValueError(f"width must be positive. Received: {width}")
 
+    # A symmetric kernel spans ``width`` sigmas either side of the peak; a causal
+    # kernel spans them on the past side only.
+    expected_kernel_size = int(width * sigma + 1) if causal else int(2 * width * sigma + 1)
+
     if kernel_size is not None:
         if kernel_size < 1:
             raise ValueError(f"kernel_size must be >= 1. Received: {kernel_size}")
     else:
-        kernel_size = int(2 * width * sigma + 1)
+        kernel_size = expected_kernel_size
 
     # Warn if kernel_size is smaller than recommended but don't fail
-    expected_kernel_size = int(2 * width * sigma + 1)
     if kernel_size < expected_kernel_size:
         ## TODO: Either add a warning or determine appropriate kernel size and raise an error
         warnings.warn(
             f"Provided kernel_size {kernel_size} is smaller than recommended "
-            f"size {expected_kernel_size} for sigma={sigma} and width={width}. "
-            "The kernel may be truncated."
+            f"size {expected_kernel_size} for sigma={sigma}, width={width} and "
+            f"causal={causal}. The kernel may be truncated."
         )
 
     if kernel_size == 1:
@@ -74,8 +113,15 @@ def gaussian_smoothing_filter_design(
 
     from scipy.signal.windows import gaussian
 
-    b = gaussian(kernel_size, std=sigma)
-    b /= np.sum(b)  # Ensure normalization
+    if causal:
+        # Take the peak and everything to its right from a symmetric kernel of
+        # 2 * kernel_size - 1 taps, then reverse-normalize: lfilter convolves
+        # b[0] with the newest sample, so index 0 is the peak (lag 0) and
+        # increasing index reaches further into the past.
+        b = gaussian(2 * kernel_size - 1, std=sigma)[kernel_size - 1 :]
+    else:
+        b = gaussian(kernel_size, std=sigma)
+    b = b / np.sum(b)  # Ensure normalization
     a = np.array([1.0])
 
     return b, a
@@ -98,6 +144,7 @@ def design_wrapper(fs: float) -> BACoeffs | None:
                 sigma=self.settings.sigma * fs,  # settings.sigma is in seconds
                 width=self.settings.width,
                 kernel_size=self.settings.kernel_size,
+                causal=self.settings.causal,
             )
 
         return design_wrapper

tests/unit/test_gaussian_smoothing_filter.py

@@ -42,6 +42,7 @@ def test_gaussian_smoothing_settings_defaults():
     assert settings.sigma == 0.01  # seconds; ~13.2 Hz low-pass (-3 dB)
     assert settings.width == 4
     assert settings.kernel_size is None
+    assert settings.causal is False  # default preserves the symmetric kernel
 
 
 def test_gaussian_smoothing_settings_custom():
@@ -278,6 +279,144 @@ def _kernel_len(fs: float) -> int:
     assert len_1000 == int(2 * 4 * 0.02 * 1000.0 + 1)
 
 
+def _dc_group_delay(b: np.ndarray) -> float:
+    """Group delay at DC, in samples, of the FIR kernel ``b``."""
+    from scipy.signal import group_delay
+
+    _w, gd = group_delay((b, np.array([1.0])), w=[0.0])
+    return float(gd[0])
+
+
+@pytest.mark.parametrize("causal", [False, True])
+@pytest.mark.parametrize("sigma", [1.0, 2.0, 5.0, 20.0])
+def test_gaussian_kernel_unit_sum(causal, sigma):
+    """Both modes yield a unit-sum kernel, so DC gain is exactly 1."""
+    b, a = gaussian_smoothing_filter_design(sigma=sigma, causal=causal)
+    assert np.isclose(np.sum(b), 1.0)
+    assert np.all(b > 0)
+    assert len(a) == 1 and a[0] == 1.0
+
+
+@pytest.mark.parametrize("sigma", [1.0, 2.0, 5.0, 20.0])
+def test_gaussian_causal_kernel_shape(sigma):
+    """The causal kernel peaks at lag 0 and decays monotonically into the past."""
+    b, _a = gaussian_smoothing_filter_design(sigma=sigma, width=4, causal=True)
+
+    assert len(b) == int(4 * sigma + 1)  # width * sigma taps, one side only
+    assert np.argmax(b) == 0  # peak at lag 0: no acausal half
+    assert np.all(np.diff(b) < 0)  # strictly decreasing tail
+
+    # It is exactly the causal half of the symmetric kernel of the same sigma,
+    # up to renormalization.
+    b_sym, _ = gaussian_smoothing_filter_design(sigma=sigma, width=4)
+    half = b_sym[len(b_sym) // 2 :]
+    assert np.allclose(b, half / np.sum(half))
+
+
+@pytest.mark.parametrize("sigma", [2.0, 5.0, 10.0, 20.0])
+def test_gaussian_causal_group_delay_matches_half_gaussian_centroid(sigma):
+    """Measured DC group delay of the causal kernel is the half-Gaussian
+    centroid, sigma * sqrt(2 / pi). Discrete truncation biases it low by a
+    fixed ~1/pi tap, so allow one tap of slack."""
+    b, _a = gaussian_smoothing_filter_design(sigma=sigma, causal=True)
+    assert _dc_group_delay(b) == pytest.approx(sigma * np.sqrt(2 / np.pi), abs=1.0)
+
+
+@pytest.mark.parametrize("sigma", [2.0, 5.0, 10.0])
+def test_gaussian_symmetric_group_delay_is_half_the_kernel(sigma):
+    """The symmetric kernel's delay is pure group delay of (n_taps - 1) / 2,
+    which is the lag the causal option exists to avoid."""
+    b, _a = gaussian_smoothing_filter_design(sigma=sigma, width=4)
+    assert _dc_group_delay(b) == pytest.approx((len(b) - 1) / 2, abs=1e-6)
+
+
+@pytest.mark.parametrize(
+    "sym_sigma,causal_sigma",
+    [
+        (2.0, 3.8),  # noise gain ~0.141: 8 samples of lag vs ~2.7
+        (5.0, 9.8),  # noise gain ~0.056: 20 samples of lag vs ~7.5
+    ],
+)
+def test_gaussian_causal_wins_on_lag_at_matched_noise_gain(sym_sigma, causal_sigma):
+    """Matched on white-noise variance reduction -- not on nominal sigma -- the
+    causal kernel achieves the same noise gain at substantially lower delay."""
+    b_sym, _ = gaussian_smoothing_filter_design(sigma=sym_sigma, width=4)
+    b_causal, _ = gaussian_smoothing_filter_design(sigma=causal_sigma, width=4, causal=True)
+
+    gain_sym = np.sum(b_sym**2)
+    gain_causal = np.sum(b_causal**2)
+    assert gain_causal == pytest.approx(gain_sym, rel=0.02)  # matched on noise gain
+
+    delay_sym = _dc_group_delay(b_sym)
+    delay_causal = _dc_group_delay(b_causal)
+    assert delay_causal < delay_sym / 2  # and the causal one is much cheaper in lag
+
+
+def test_gaussian_causal_kernel_size_counts_causal_taps():
+    """kernel_size overrides the automatic length in causal mode too, and is
+    interpreted as the number of causal taps."""
+    b, _a = gaussian_smoothing_filter_design(sigma=2.0, width=4, kernel_size=9, causal=True)
+    assert len(b) == 9
+    assert np.isclose(np.sum(b), 1.0)
+    assert np.argmax(b) == 0
+
+
+def test_gaussian_causal_small_kernel_size_warns():
+    """The too-small-kernel warning applies to the causal branch, against the
+    one-sided recommended length (width * sigma + 1)."""
+    with pytest.warns(UserWarning, match="smaller than recommended"):
+        b, _a = gaussian_smoothing_filter_design(sigma=5.0, width=4, kernel_size=5, causal=True)
+    assert len(b) == 5
+
+    # ...and the one-sided length itself does not trip it, though the same
+    # length would be undersized for a symmetric kernel.
+    with warnings.catch_warnings():
+        warnings.simplefilter("error")
+        gaussian_smoothing_filter_design(sigma=5.0, width=4, kernel_size=21, causal=True)
+    with pytest.warns(UserWarning, match="smaller than recommended"):
+        gaussian_smoothing_filter_design(sigma=5.0, width=4, kernel_size=21, causal=False)
+
+
+def test_gaussian_causal_identity_kernel_warns():
+    """A single causal tap is still an identity kernel."""
+    with pytest.warns(UserWarning, match="identity"):
+        gaussian_smoothing_filter_design(sigma=2.0, kernel_size=1, causal=True)
+
+
+def test_gaussian_causal_transformer_end_to_end():
+    """The causal setting reaches the designed kernel through design_wrapper,
+    with sigma still interpreted in seconds and scaled by fs."""
+    fs = 100.0
+    sigma_s = 0.05
+
+    def _kernel(causal: bool) -> np.ndarray:
+        proc = GaussianSmoothingFilterTransformer(
+            GaussianSmoothingSettings(axis="time", sigma=sigma_s, width=4, causal=causal)
+        )
+        msg = AxisArray(
+            data=np.random.randn(200, 2),
+            dims=["time", "ch"],
+            axes={
+                "time": AxisArray.TimeAxis(fs=fs, offset=0),
+                "ch": AxisArray.CoordinateAxis(data=np.arange(2).astype(str), dims=["ch"]),
+            },
+            key="test_gaussian_causal",
+        )
+        out = proc(msg)
+        assert np.isfinite(out.data).all()
+        b, _a = proc.state.filter.settings.coefs
+        return b
+
+    b_causal = _kernel(True)
+    b_sym = _kernel(False)
+
+    # sigma = 5 samples at 100 Hz
+    assert len(b_causal) == int(4 * sigma_s * fs + 1)
+    assert len(b_sym) == int(2 * 4 * sigma_s * fs + 1)
+    assert np.argmax(b_causal) == 0
+    assert _dc_group_delay(b_causal) < _dc_group_delay(b_sym)
+
+
 def test_gaussian_identity_kernel_design_warns():
     """The standalone design function still warns for a single-tap kernel."""
     with pytest.warns(UserWarning, match="identity"):