diff --git a/sfm/shape_metrics.py b/sfm/shape_metrics.py new file mode 100644 index 0000000..05c9dd5 --- /dev/null +++ b/sfm/shape_metrics.py @@ -0,0 +1,58 @@ +"""Waveform-shape metrics for false-trigger detection. + +A false trigger is an isolated impulse (quiet → spike → quiet); a real event +rings for many cycles. Two numbers separate them: crest factor (how far the +peak stands above the typical sample) and how many samples sit near the peak. +""" +from __future__ import annotations +import numpy as np + +_GEO_CHANNELS = ("Tran", "Vert", "Long") +NEAR_PEAK_FRACTION = 0.5 # a sample "near the peak" is >= this * peak amplitude + + +def channel_shape(x) -> dict | None: + x = np.asarray(x, dtype=float) + if x.size < 2: + return None + peak = float(np.max(np.abs(x))) + if peak <= 0: + return None + rms = float(np.sqrt(np.mean(x ** 2))) + if rms <= 0: + return None + near = int(np.sum(np.abs(x) >= NEAR_PEAK_FRACTION * peak)) + return {"crest_factor": peak / rms, "near_peak_count": near, + "sample_count": int(x.size)} + + +def shape_from_samples(chans: dict) -> dict | None: + best_axis, best_peak, best_x = None, -1.0, None + for ax in _GEO_CHANNELS: + x = chans.get(ax) + if x is None: + continue + x = np.asarray(x, dtype=float) + if x.size < 2: + continue + p = float(np.max(np.abs(x))) + if p > best_peak: + best_axis, best_peak, best_x = ax, p, x + if best_axis is None: + return None + s = channel_shape(best_x) + if s is None: + return None + s["axis"] = best_axis + return s + + +def shape_from_h5(path) -> dict | None: + import h5py + try: + with h5py.File(path, "r") as f: + chans = {ax: f[f"samples/{ax}"][:] for ax in _GEO_CHANNELS + if f"samples/{ax}" in f} + except Exception: + return None + return shape_from_samples(chans) diff --git a/tests/test_shape_metrics.py b/tests/test_shape_metrics.py new file mode 100644 index 0000000..ae23475 --- /dev/null +++ b/tests/test_shape_metrics.py @@ -0,0 +1,31 @@ +import numpy as np +from sfm.shape_metrics import channel_shape, shape_from_samples + +def test_needle_spike_high_crest_few_near_peak(): + x = np.zeros(1024); x[500] = 1.0 # one isolated spike + s = channel_shape(x) + assert s["sample_count"] == 1024 + assert s["crest_factor"] > 15 # peak towers over rms + assert s["near_peak_count"] <= 3 # almost nothing near the peak + +def test_ringing_low_crest_many_near_peak(): + t = np.arange(1024) + x = np.sin(2*np.pi*t/32) * np.exp(-t/4000) # decaying oscillation + s = channel_shape(x) + assert s["crest_factor"] < 6 + assert s["near_peak_count"] > 30 # many samples near the peak + +def test_channel_shape_none_for_unusable(): + assert channel_shape(np.zeros(1024)) is None # flat / all-zero + assert channel_shape(np.array([1.0])) is None # too short + +def test_shape_from_samples_picks_max_peak_axis(): + chans = {"Tran": np.zeros(1024), "Vert": np.zeros(1024), "Long": np.zeros(1024)} + chans["Long"][10] = 0.5 + chans["Vert"] = np.sin(np.arange(1024)/5) * 0.01 + s = shape_from_samples(chans) + assert s["axis"] == "Long" # Long has the biggest peak + assert s["near_peak_count"] <= 3 + +def test_shape_from_samples_none_when_no_geo(): + assert shape_from_samples({"MicL": np.ones(1024)}) is None