M
muq
93ae08bc91
The OSS threshold_qualifier was last synced from utils on Feb 23 2026, before the major robustness improvements landed in utils on Mar 19 2026. This commit brings it fully up to date. Changes ported from utils: - Multi-pass detection (Pass 1 strong burst, Pass 2 weak residual, Pass 3 sustained faint burst via macro-window averaging) - Noise floor estimation via percentile instead of simple max*threshold - Dynamic range ratio guard (early exit on low-contrast captures) - Improved _find_ranges, _expand_and_filter_ranges, _merge_ranges helpers - Spectral smoothing in _estimate_spectral_bounds for wideband bursts - Minimum duration filter expressed in absolute time (5ms) not sample count Also includes the Pass 2 hysteresis spillover fix: - Pass 2 expansion now runs against residual_power (masked) instead of smoothed_power, preventing it from walking into Pass 1 territory - Pass 2 mask now has a window_size guard band around Pass 1 ranges, matching the guard already used in Pass 3 Only change from utils: import swapped to ria_toolkit_oss.datatypes.
360 lines
15 KiB
Python
360 lines
15 KiB
Python
"""
|
|
Temporal signal detection and boundary refinement via Hysteresis Thresholding.
|
|
|
|
Provides methods to detect signal bursts in the time domain by triggering on
|
|
smoothed power peaks and expanding boundaries to capture the full energy envelope.
|
|
|
|
This module implements a **dual-threshold trigger** to solve the 'chatter'
|
|
problem in noisy environments, ensuring that signal annotations encapsulate
|
|
the entire rise and fall of a burst rather than just the peak.
|
|
|
|
**Key Design Decisions**:
|
|
|
|
1. **Hysteresis Logic (Dual-Threshold)**:
|
|
- **Trigger**: High threshold (`threshold * max_power`) ensures high confidence
|
|
in signal presence.
|
|
- **Boundary**: Low threshold (`0.5 * trigger`) allows the annotation to
|
|
"crawl" outward, capturing the lower-energy start and end of the burst
|
|
often missed by simple single-threshold detectors.
|
|
|
|
2. **Temporal Smoothing**: Uses a moving average window (`window_size`) prior
|
|
- to thresholding. This prevents high-frequency noise spikes from causing
|
|
fragmented annotations and provides a more stable estimate of the
|
|
signal's power envelope.
|
|
|
|
3. **Spectral Profiling**: Once a temporal segment is isolated, the module
|
|
- performs an automated FFT analysis. It identifies the **90% spectral
|
|
occupancy** to define the frequency boundaries (`f_min`, `f_max`),
|
|
allowing the detector to work on narrowband and wideband signals without
|
|
manual frequency tuning.
|
|
|
|
4. **Baseband/RF Mapping**: Automatically handles the conversion from
|
|
- relative FFT bin frequencies to absolute RF frequencies by referencing
|
|
`recording.metadata["center_frequency"]`.
|
|
|
|
5. **False Positive Mitigation**: Implements a hard minimum duration check
|
|
- (10ms) to ignore transient hardware spikes or noise floor fluctuations
|
|
that do not constitute a valid signal burst.
|
|
|
|
The module is designed to be the primary "first-pass" detector for pulsed
|
|
waveforms (like ADS-B, Lora, or bursty FSK) before passing them to
|
|
classification or demodulation stages.
|
|
"""
|
|
|
|
import json
|
|
from typing import Optional
|
|
|
|
import numpy as np
|
|
|
|
from ria_toolkit_oss.datatypes import Annotation, Recording
|
|
|
|
|
|
def _find_ranges(indices, max_gap):
|
|
"""
|
|
Groups individual indices into continuous temporal ranges.
|
|
|
|
Args:
|
|
indices: Array of indices where the signal exceeded a threshold.
|
|
max_gap: Maximum gap allowed between indices to consider them part
|
|
of the same range.
|
|
|
|
Returns:
|
|
A list of (start, stop) tuples representing detected signal segments.
|
|
"""
|
|
|
|
if len(indices) == 0:
|
|
return []
|
|
|
|
start = indices[0]
|
|
prev = indices[0]
|
|
ranges = []
|
|
|
|
for i in range(1, len(indices)):
|
|
if indices[i] - prev > max_gap:
|
|
ranges.append((start, prev))
|
|
start = indices[i]
|
|
prev = indices[i]
|
|
|
|
ranges.append((start, prev))
|
|
|
|
return ranges
|
|
|
|
|
|
def _expand_and_filter_ranges(
|
|
smoothed_power: np.ndarray,
|
|
initial_ranges: list[tuple[int, int]],
|
|
boundary_val: float,
|
|
min_duration_samples: int,
|
|
) -> list[tuple[int, int]]:
|
|
"""Apply hysteresis expansion and minimum-duration filtering."""
|
|
out: list[tuple[int, int]] = []
|
|
n = len(smoothed_power)
|
|
for start, stop in initial_ranges:
|
|
if (stop - start) < min_duration_samples:
|
|
continue
|
|
|
|
true_start = start
|
|
while true_start > 0 and smoothed_power[true_start] > boundary_val:
|
|
true_start -= 1
|
|
|
|
true_stop = stop
|
|
while true_stop < n - 1 and smoothed_power[true_stop] > boundary_val:
|
|
true_stop += 1
|
|
|
|
if (true_stop - true_start) >= min_duration_samples:
|
|
out.append((true_start, true_stop))
|
|
return out
|
|
|
|
|
|
def _merge_ranges(ranges: list[tuple[int, int]], max_gap: int) -> list[tuple[int, int]]:
|
|
"""Merge overlapping or near-adjacent ranges."""
|
|
if not ranges:
|
|
return []
|
|
ranges = sorted(ranges, key=lambda r: r[0])
|
|
merged = [ranges[0]]
|
|
for s, e in ranges[1:]:
|
|
last_s, last_e = merged[-1]
|
|
if s <= last_e + max_gap:
|
|
merged[-1] = (last_s, max(last_e, e))
|
|
else:
|
|
merged.append((s, e))
|
|
return merged
|
|
|
|
|
|
def _estimate_noise_floor(power: np.ndarray, quantile: float = 20.0) -> float:
|
|
"""Estimate baseline from the quieter portion of the envelope."""
|
|
return float(np.percentile(power, quantile))
|
|
|
|
|
|
def _estimate_group_gap(sample_rate: float) -> int:
|
|
"""Use a fixed temporal grouping gap instead of reusing the smoothing window."""
|
|
return max(1, int(0.001 * sample_rate))
|
|
|
|
|
|
def _estimate_spectral_bounds(signal_segment: np.ndarray, sample_rate: float) -> tuple[float, float]:
|
|
"""Estimate occupied bandwidth from a smoothed magnitude spectrum."""
|
|
if len(signal_segment) == 0:
|
|
return -sample_rate / 4, sample_rate / 4
|
|
|
|
window = np.hanning(len(signal_segment))
|
|
windowed = signal_segment * window
|
|
|
|
fft_data = np.abs(np.fft.fftshift(np.fft.fft(windowed)))
|
|
fft_freqs = np.fft.fftshift(np.fft.fftfreq(len(signal_segment), 1 / sample_rate))
|
|
|
|
# Smooth the spectrum so noise-like wideband bursts form a contiguous mask
|
|
# instead of thousands of tiny isolated runs.
|
|
spectral_smooth_bins = max(5, min(257, (len(signal_segment) // 512) | 1))
|
|
spectral_kernel = np.ones(spectral_smooth_bins, dtype=np.float64) / spectral_smooth_bins
|
|
smoothed_fft = np.convolve(fft_data, spectral_kernel, mode="same")
|
|
|
|
spectral_floor = float(np.percentile(smoothed_fft, 20))
|
|
spectral_peak = float(np.max(smoothed_fft))
|
|
spectral_ratio = spectral_peak / max(spectral_floor, 1e-12)
|
|
|
|
if spectral_ratio < 1.2:
|
|
return -sample_rate / 4, sample_rate / 4
|
|
|
|
spectral_thresh = spectral_floor + 0.1 * (spectral_peak - spectral_floor)
|
|
sig_indices = np.where(smoothed_fft > spectral_thresh)[0]
|
|
|
|
if len(sig_indices) == 0:
|
|
peak_idx = int(np.argmax(smoothed_fft))
|
|
bin_hz = sample_rate / len(signal_segment)
|
|
half_bins = max(1, int(np.ceil(10_000.0 / bin_hz)))
|
|
lo_idx = max(0, peak_idx - half_bins)
|
|
hi_idx = min(len(smoothed_fft) - 1, peak_idx + half_bins)
|
|
else:
|
|
runs = _find_ranges(sig_indices, max_gap=max(1, spectral_smooth_bins // 2))
|
|
peak_idx = int(np.argmax(smoothed_fft))
|
|
lo_idx, hi_idx = min(
|
|
runs,
|
|
key=lambda run: 0 if run[0] <= peak_idx <= run[1] else min(abs(run[0] - peak_idx), abs(run[1] - peak_idx)),
|
|
)
|
|
|
|
# Prevent extremely narrow tone boxes from collapsing to just a few bins.
|
|
min_total_bw_hz = 20_000.0
|
|
min_half_bins = max(1, int(np.ceil((min_total_bw_hz / 2) / (sample_rate / len(signal_segment)))))
|
|
center_idx = int(round((lo_idx + hi_idx) / 2))
|
|
lo_idx = max(0, min(lo_idx, center_idx - min_half_bins))
|
|
hi_idx = min(len(smoothed_fft) - 1, max(hi_idx, center_idx + min_half_bins))
|
|
|
|
return float(fft_freqs[lo_idx]), float(fft_freqs[hi_idx])
|
|
|
|
|
|
def threshold_qualifier(
|
|
recording: Recording,
|
|
threshold: float,
|
|
window_size: Optional[int] = None,
|
|
label: Optional[str] = None,
|
|
annotation_type: Optional[str] = "standalone",
|
|
channel: int = 0,
|
|
) -> Recording:
|
|
"""
|
|
Annotate a recording with bounding boxes for regions above a threshold.
|
|
Threshold is defined as a fraction of the maximum sample magnitude.
|
|
This algorithm searches for samples above the threshold and combines them into ranges if they
|
|
are within window_size of each other.
|
|
Detects and annotates signals using energy thresholding and spectral analysis.
|
|
|
|
The algorithm follows these steps:
|
|
1. Smooths power data using a moving average.
|
|
2. Identifies 'peak' regions exceeding a high trigger threshold.
|
|
3. Uses hysteresis to expand boundaries until power drops below a lower threshold.
|
|
4. Performs an FFT on each segment to determine frequency occupancy.
|
|
|
|
Args:
|
|
recording: The Recording object containing IQ or real signal data.
|
|
threshold: Sensitivity multiplier (0.0 to 1.0) applied to max power.
|
|
window_size: Size of the smoothing filter in samples. Defaults to 1ms worth of samples.
|
|
label: Custom string label for annotations.
|
|
annotation_type: Metadata string for the 'type' field in the annotation.
|
|
channel: Index of the channel to annotate. Defaults to 0.
|
|
|
|
Returns:
|
|
A new Recording object populated with detected Annotations.
|
|
"""
|
|
# Extract signal and metadata
|
|
sample_data = recording.data[channel]
|
|
sample_rate = recording.metadata["sample_rate"]
|
|
center_frequency = recording.metadata.get("center_frequency", 0)
|
|
|
|
if window_size is None:
|
|
window_size = max(64, int(sample_rate * 0.001))
|
|
|
|
# --- 1. SIGNAL CONDITIONING ---
|
|
# Convert to power (Magnitude squared)
|
|
power_data = np.abs(sample_data) ** 2
|
|
smoothing_window = np.ones(window_size) / window_size
|
|
smoothed_power = np.convolve(power_data, smoothing_window, mode="same")
|
|
group_gap_samples = _estimate_group_gap(sample_rate)
|
|
|
|
# Define thresholds using peak relative to baseline.
|
|
max_power = np.max(smoothed_power)
|
|
noise_floor = _estimate_noise_floor(smoothed_power)
|
|
dynamic_range_ratio = max_power / max(noise_floor, 1e-12)
|
|
|
|
# Soft early exit: keep a guard for low-contrast noise, but compute it from
|
|
# the quieter tail of the envelope so burst-heavy captures are not rejected.
|
|
if dynamic_range_ratio < 1.5:
|
|
return Recording(data=recording.data, metadata=recording.metadata, annotations=recording.annotations)
|
|
|
|
trigger_val = noise_floor + threshold * (max_power - noise_floor)
|
|
boundary_val = noise_floor + 0.5 * threshold * (max_power - noise_floor)
|
|
|
|
# --- 2. INITIAL DETECTION ---
|
|
# Enforce an explicit minimum duration in seconds; this is stable across
|
|
# varying capture lengths and avoids over-fitting to recording length.
|
|
min_duration_samples = max(1, int(0.005 * sample_rate))
|
|
annotations = []
|
|
|
|
# Pass 1: Detect stronger bursts.
|
|
indices = np.where(smoothed_power > trigger_val)[0]
|
|
pass1_initial = _find_ranges(indices=indices, max_gap=group_gap_samples)
|
|
pass1_ranges = _expand_and_filter_ranges(
|
|
smoothed_power=smoothed_power,
|
|
initial_ranges=pass1_initial,
|
|
boundary_val=boundary_val,
|
|
min_duration_samples=min_duration_samples,
|
|
)
|
|
|
|
# Pass 2: Recover weaker bursts on residual power not already covered.
|
|
# This improves recall in mixed-amplitude captures.
|
|
# Expand each Pass-1 range by the smoothing window on both sides so the
|
|
# smoothing skirts of a strong burst are not re-detected as a weak burst
|
|
# immediately adjacent to it (mirrors the guard used in Pass 3).
|
|
mask = np.ones_like(smoothed_power, dtype=np.float32)
|
|
pass2_mask_expand = window_size
|
|
for s, e in pass1_ranges:
|
|
mask[max(0, s - pass2_mask_expand) : min(len(mask), e + pass2_mask_expand)] = 0.0
|
|
residual_power = smoothed_power * mask
|
|
|
|
residual_max = float(np.max(residual_power))
|
|
residual_ratio = residual_max / max(noise_floor, 1e-12)
|
|
|
|
pass2_ranges: list[tuple[int, int]] = []
|
|
if residual_ratio >= 2.0:
|
|
weak_threshold = max(0.3, threshold * 0.7)
|
|
weak_trigger = noise_floor + weak_threshold * (residual_max - noise_floor)
|
|
weak_boundary = noise_floor + 0.5 * weak_threshold * (residual_max - noise_floor)
|
|
weak_indices = np.where(residual_power > weak_trigger)[0]
|
|
pass2_initial = _find_ranges(indices=weak_indices, max_gap=group_gap_samples)
|
|
pass2_ranges = _expand_and_filter_ranges(
|
|
smoothed_power=residual_power,
|
|
initial_ranges=pass2_initial,
|
|
boundary_val=weak_boundary,
|
|
min_duration_samples=min_duration_samples,
|
|
)
|
|
|
|
# Pass 3: Detect sustained faint bursts via macro-window averaging.
|
|
# Targets bursts whose peak power is near the trigger level but whose
|
|
# *average* power is consistently elevated above the noise floor — these
|
|
# are missed by peak-based detection because only a few short spikes exceed
|
|
# the trigger, all too brief to pass the minimum-duration filter.
|
|
#
|
|
# The mask is applied to power_data *before* convolving so that bright
|
|
# burst energy does not bleed through the long window into adjacent regions,
|
|
# which would inflate macro_residual_max and push the trigger above the
|
|
# faint burst's average power.
|
|
macro_window_size = max(window_size * 16, int(sample_rate * 0.02))
|
|
macro_kernel = np.ones(macro_window_size, dtype=np.float64) / macro_window_size
|
|
# Expand each annotated range by half the macro window on both sides so that
|
|
# the long convolution cannot "see" the leading/trailing edges of already-
|
|
# annotated bursts, which would produce spurious short fragments in Pass 3.
|
|
macro_expand = macro_window_size * 2
|
|
masked_power_for_macro = power_data.copy()
|
|
n = len(masked_power_for_macro)
|
|
for s, e in pass1_ranges + pass2_ranges:
|
|
masked_power_for_macro[max(0, s - macro_expand) : min(n, e + macro_expand)] = 0.0
|
|
macro_residual = np.convolve(masked_power_for_macro, macro_kernel, mode="same")
|
|
|
|
macro_residual_max = float(np.max(macro_residual))
|
|
|
|
pass3_ranges: list[tuple[int, int]] = []
|
|
if macro_residual_max / max(noise_floor, 1e-12) >= 1.3:
|
|
macro_trigger = noise_floor + threshold * (macro_residual_max - noise_floor)
|
|
macro_boundary = noise_floor + 0.5 * threshold * (macro_residual_max - noise_floor)
|
|
macro_indices = np.where(macro_residual > macro_trigger)[0]
|
|
macro_initial = _find_ranges(indices=macro_indices, max_gap=group_gap_samples)
|
|
pass3_ranges = _expand_and_filter_ranges(
|
|
smoothed_power=macro_residual,
|
|
initial_ranges=macro_initial,
|
|
boundary_val=macro_boundary,
|
|
min_duration_samples=min_duration_samples,
|
|
)
|
|
|
|
all_ranges = _merge_ranges(pass1_ranges + pass2_ranges + pass3_ranges, max_gap=group_gap_samples)
|
|
|
|
for true_start, true_stop in all_ranges:
|
|
|
|
# --- 4. SPECTRAL ANALYSIS (Frequency Detection) ---
|
|
signal_segment = sample_data[true_start:true_stop]
|
|
f_min, f_max = _estimate_spectral_bounds(signal_segment, sample_rate)
|
|
|
|
# --- 5. ANNOTATION GENERATION ---
|
|
ann_label = label if label is not None else f"{int(threshold*100)}%"
|
|
|
|
# Pack metadata for the UI/Downstream processing
|
|
comment_data = {
|
|
"type": annotation_type,
|
|
"generator": "threshold_qualifier",
|
|
"params": {
|
|
"threshold": threshold,
|
|
"window_size": window_size,
|
|
},
|
|
}
|
|
|
|
anno = Annotation(
|
|
sample_start=true_start,
|
|
sample_count=true_stop - true_start,
|
|
freq_lower_edge=center_frequency + f_min,
|
|
freq_upper_edge=center_frequency + f_max,
|
|
label=ann_label,
|
|
comment=json.dumps(comment_data),
|
|
detail={"generator": "hysteresis_qualifier"},
|
|
)
|
|
annotations.append(anno)
|
|
|
|
# Return a new Recording object including the new annotations
|
|
return Recording(data=recording.data, metadata=recording.metadata, annotations=recording.annotations + annotations)
|