diff --git a/src/ezmsg/sigproc/gaussiansmoothing.py b/src/ezmsg/sigproc/gaussiansmoothing.py index 6a697962..f0f01507 100644 --- a/src/ezmsg/sigproc/gaussiansmoothing.py +++ b/src/ezmsg/sigproc/gaussiansmoothing.py @@ -14,10 +14,13 @@ class GaussianSmoothingSettings(FilterBaseSettings): - sigma: float | None = 1.0 + sigma: float | None = 0.01 """ sigma : float - Standard deviation of the Gaussian kernel. + Standard deviation of the Gaussian kernel, in **seconds**. Converted + to samples using the sampling rate of the first message. + The -3 dB corner frequency is sqrt(ln 2) / (2 * pi * sigma); the + default of 0.01 s is equivalent to a ~13.2 Hz low-pass. """ width: int | None = 4 @@ -38,6 +41,8 @@ def gaussian_smoothing_filter_design( width: int = 4, kernel_size: int | None = None, ) -> 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.""" # Parameter checks if sigma <= 0: raise ValueError(f"sigma must be positive. Received: {sigma}") @@ -61,6 +66,12 @@ def gaussian_smoothing_filter_design( "The kernel may be truncated." ) + if kernel_size == 1: + warnings.warn( + f"kernel_size=1 (sigma={sigma} samples, width={width}) yields an " + "identity (single-tap) kernel: no smoothing will be performed." + ) + from scipy.signal.windows import gaussian b = gaussian(kernel_size, std=sigma) @@ -74,10 +85,9 @@ class GaussianSmoothingFilterTransformer(FilterByDesignTransformer[GaussianSmoot def get_design_function( self, ) -> Callable[[float], BACoeffs]: - # Create a wrapper function that ignores fs parameter since gaussian smoothing doesn't need it def design_wrapper(fs: float) -> BACoeffs: return gaussian_smoothing_filter_design( - sigma=self.settings.sigma, + sigma=self.settings.sigma * fs, # settings.sigma is in seconds width=self.settings.width, kernel_size=self.settings.kernel_size, ) diff --git a/tests/unit/test_gaussian_smoothing_filter.py b/tests/unit/test_gaussian_smoothing_filter.py index 1d701008..5ab2c2d7 100644 --- a/tests/unit/test_gaussian_smoothing_filter.py +++ b/tests/unit/test_gaussian_smoothing_filter.py @@ -36,7 +36,7 @@ def test_gaussian_smoothing_filter_function(axis, sigma, width, kernel_size): def test_gaussian_smoothing_settings_defaults(): """Test the GaussianSmoothingSettings class with default values.""" settings = GaussianSmoothingSettings() - assert settings.sigma == 1.0 + assert settings.sigma == 0.01 # seconds; ~13.2 Hz low-pass (-3 dB) assert settings.width == 4 assert settings.kernel_size is None @@ -136,9 +136,9 @@ def test_gaussian_smoothing_filter_process(data_shape): msg = AxisArray(data=data, dims=dims, axes=axes, key="test_gaussian_smoothing") - # Instantiate transformer (not unit) + # Instantiate transformer (not unit). sigma is in seconds: 0.02 s @ 100 Hz = 2 samples. transformer = GaussianSmoothingFilterTransformer( - settings=GaussianSmoothingSettings(axis="time", sigma=2.0, width=4) + settings=GaussianSmoothingSettings(axis="time", sigma=0.02, width=4) ) # Process message using __call__ method @@ -152,8 +152,9 @@ def test_gaussian_smoothing_filter_process(data_shape): def test_gaussian_smoothing_edge_cases(): """Test edge cases for gaussian smoothing filter.""" - # Test with very small sigma - coefs_small = gaussian_smoothing_filter_design(sigma=0.01) + # Sub-sample sigma degenerates to an identity kernel and warns. + with pytest.warns(UserWarning, match="identity"): + coefs_small = gaussian_smoothing_filter_design(sigma=0.01) b_small, a_small = coefs_small assert len(b_small) > 0 assert np.isclose(np.sum(b_small), 1.0) @@ -206,10 +207,11 @@ def _calc_smoothing_effect(msg): original_variance = _calc_smoothing_effect(msg_in) - # Initialize filter with small sigma (minimal smoothing) + # Initialize filter with small sigma (minimal smoothing). + # sigma is in seconds: 0.0025 s @ 200 Hz = 0.5 samples. proc = GaussianSmoothingFilterTransformer( axis="time", - sigma=0.5, + sigma=0.0025, width=4, coef_type="ba", ) @@ -221,8 +223,8 @@ def _calc_smoothing_effect(msg): # Small sigma should have minimal effect assert np.allclose(variance1, original_variance, rtol=0.1) - # Update settings - change to larger sigma (more smoothing) - proc.update_settings(sigma=3.0) + # Update settings - change to larger sigma (more smoothing): 3 samples @ 200 Hz + proc.update_settings(sigma=0.015) # Process the same message with new settings result2 = proc(msg_in) @@ -234,7 +236,7 @@ def _calc_smoothing_effect(msg): # Test update_settings with complete new settings object new_settings = GaussianSmoothingSettings( axis="time", - sigma=5.0, # Even larger sigma + sigma=0.025, # Even larger sigma: 5 samples @ 200 Hz width=6, kernel_size=None, coef_type="ba", @@ -248,6 +250,40 @@ def _calc_smoothing_effect(msg): assert np.all(variance3 < variance2) +def test_gaussian_sigma_is_in_seconds(): + """The same settings must yield the same temporal smoothing at any fs: + kernel length in samples scales with the sampling rate.""" + + def _kernel_len(fs: float) -> int: + proc = GaussianSmoothingFilterTransformer(GaussianSmoothingSettings(sigma=0.02, axis="time")) + msg = AxisArray( + data=np.zeros((int(fs), 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_sigma_seconds", + ) + _ = proc(msg) + b, _a = proc.state.filter.settings.coefs + return len(b) + + len_100 = _kernel_len(100.0) # sigma = 2 samples + len_1000 = _kernel_len(1000.0) # sigma = 20 samples + assert len_100 == int(2 * 4 * 0.02 * 100.0 + 1) + assert len_1000 == int(2 * 4 * 0.02 * 1000.0 + 1) + + +def test_gaussian_identity_kernel_warns_and_passes_through(): + """An explicit single-tap kernel warns and leaves the data unchanged.""" + proc = GaussianSmoothingFilterTransformer(GaussianSmoothingSettings(sigma=0.02, kernel_size=1, axis="time")) + msg = make_msg() + with pytest.warns(UserWarning): + result = proc(msg) + assert np.allclose(result.data, msg.data) + + def test_gaussian_empty_after_init(): from ezmsg.sigproc.gaussiansmoothing import GaussianSmoothingFilterTransformer, GaussianSmoothingSettings