diff --git a/src/ezmsg/sigproc/gaussiansmoothing.py b/src/ezmsg/sigproc/gaussiansmoothing.py index 25a0b76..8735e4d 100644 --- a/src/ezmsg/sigproc/gaussiansmoothing.py +++ b/src/ezmsg/sigproc/gaussiansmoothing.py @@ -33,6 +33,35 @@ 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. """ @@ -40,9 +69,16 @@ 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}") @@ -50,20 +86,23 @@ def gaussian_smoothing_filter_design( 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 diff --git a/tests/unit/test_gaussian_smoothing_filter.py b/tests/unit/test_gaussian_smoothing_filter.py index 29491f9..7217c58 100644 --- a/tests/unit/test_gaussian_smoothing_filter.py +++ b/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"):