b4718aa425
Most of what identifies itself on the air identifies itself in Morse. A repeater, a beacon, an unattended transmitter: four to six characters, over in a second or two, and no speech anywhere in the capture. Every one of those was being thrown away, in three separate places. The callsign book and the map were built inside the transcription branch, on the reasoning that callsigns come out of transcripts. They also come out of Morse and out of APRS headers, neither of which involves a speech recogniser -- so a receiver with none installed found none of them, and a CW ident reached the sidecar and stopped there. Both are now built whenever classification is on, and all three sources go through one place. The CW decoder only ran where the classifier had already said cw, ook or carrier. A two-second ident is a fraction of a capture named after whatever filled the rest of it. Every capture is offered to it now, once it has finished; a decode does not relabel a capture that plainly holds speech. And the decoder's own gates were written for a paragraph. Three characters, eight elements, and any repeated character refused -- which read VVV, DE, AR and K correctly and then discarded them. Short is the normal case now, on a second bar: perfect timing, nothing undecoded, and the keyed tone at least 20 dB over its band. That last is not decoration. With four elements the dot length is fitted to those very elements, so noise lands on the grid as neatly as keying does; a third of a second of white noise decodes as a perfectly timed V. Over 200 noise blocks the loudest bin never rose 13 dB above the median while keying at 3 dB SNR sits above 40. One keyed element is still refused: a single pulse is an E or a T whether a person sent it or the squelch opened on a click. 288 non-Morse cases, no false positives. Feeding that text to a callsign lookup made truncation matter. A capture opens when the squelch does, halfway through an element as often as not, and half a character is not a smaller reading -- a K missing its first dash is an A. So the sliced character is dropped, and so is the rest of its word, because what is left can read as a whole one: K1AA caught halfway through is K1A, which is somebody else. Across 1805 truncated captures that is 107 invented callsigns down to none, with 550 correct ones still found. Phonetics, which is the other half of the ask. A recogniser has never heard of the alphabet -- it writes what the words sounded like: Whiskey-One-Alpha-Whiskey hyphenated WhiskeyOneAlphaWhiskey run together Whiskey1AlphaWhiskey and half in digits wiskey one alfa whisky spelled the way it sounded whiskey one alpha, uh, whiskey with the hesitation written down All read back to W1AW now. A word is only taken apart when it is phonetic all the way through, which is what keeps it off "kilometre" and "victorious". Two bugs found on the way, both of which invented a callsign: - Nothing is joined across a slash any more. The beacon W1AW/B came back as W1AWB, which belongs to nobody, and W1AW-4 came back as nothing at all. - Nor across a gap the sender chose. A transcript's spacing is the recogniser's guess and may be closed up; a word gap in Morse is seven dot units, so "KU0W K" is a station signing off, not a longer callsign. Also here: - saunterbrowse gives Morse a panel of its own, with the licence under it, searchable with / and readable with t. - --simulate no longer looks anything up or writes a map. The demo band is invented but W1AW is the ARRL's own station, and it would have been pinned to the same map a real scan writes. - classify._psk_order took the logarithm of zero on a silent block. - The demo band has a repeater ident in it, and its Morse no longer runs one repeat into the next. - conftest refuses a real licence lookup from any test. 1161 tests, up from 1014. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_016PsWPTweCT6pwxKngvVxcg
1238 lines
53 KiB
Python
Executable file
1238 lines
53 KiB
Python
Executable file
"""Modulation and signal-type classification.
|
|
|
|
The classifier works on a block of complex baseband that has already been
|
|
centred on the signal and decimated to a rate a few times its bandwidth. It
|
|
extracts a feature vector, scores it against a rule set for the common
|
|
modulation families, then refines the answer with a table of known systems
|
|
keyed on frequency, bandwidth and symbol rate.
|
|
|
|
Nothing here decodes traffic; it names what the signal *is*.
|
|
"""
|
|
|
|
from __future__ import annotations
|
|
|
|
import math
|
|
from dataclasses import dataclass, field
|
|
|
|
import numpy as np
|
|
from scipy import signal as sps
|
|
|
|
from .dsp import (db, instantaneous_frequency,
|
|
occupied_bandwidth, spectral_flatness, welch_psd)
|
|
|
|
__all__ = ["classify", "Classification", "SignalFeatures", "extract_features",
|
|
"CTCSS_TONES", "detect_ctcss", "SSBAlignment", "ssb_alignment",
|
|
"ControlChannel", "trunk_control", "CONTROL_SIGNATURES"]
|
|
|
|
|
|
def _pow2_floor(n: int, cap: int = 1 << 16) -> int:
|
|
"""Largest power of two that is <= n (and <= cap)."""
|
|
n = int(min(n, cap))
|
|
return 1 << int(math.floor(math.log2(max(2, n))))
|
|
|
|
|
|
# EIA/TIA-603 standard CTCSS tones, Hz.
|
|
CTCSS_TONES = (
|
|
67.0, 69.3, 71.9, 74.4, 77.0, 79.7, 82.5, 85.4, 88.5, 91.5,
|
|
94.8, 97.4, 100.0, 103.5, 107.2, 110.9, 114.8, 118.8, 123.0, 127.3,
|
|
131.8, 136.5, 141.3, 146.2, 151.4, 156.7, 159.8, 162.2, 165.5, 167.9,
|
|
171.3, 173.8, 177.3, 179.9, 183.5, 186.2, 189.9, 192.8, 196.6, 199.5,
|
|
203.5, 206.5, 210.7, 218.1, 225.7, 229.1, 233.6, 241.8, 250.3, 254.1,
|
|
)
|
|
|
|
# Symbol rates worth naming when the cyclostationary estimator lands near one.
|
|
_KNOWN_BAUD = (
|
|
(300, "300 baud"), (512, "POCSAG-512"), (1200, "1200 baud"),
|
|
(1600, "FLEX-1600"), (2400, "2400 baud"), (3200, "FLEX-3200"),
|
|
(4800, "4800 baud"), (6400, "FLEX-6400"), (9600, "9600 baud"),
|
|
(18000, "TETRA"), (19200, "19.2 kbaud"), (36000, "36 kbaud"),
|
|
)
|
|
|
|
|
|
@dataclass
|
|
class SignalFeatures:
|
|
"""Everything the rules and the report are allowed to look at."""
|
|
|
|
sample_rate: float
|
|
n_samples: int
|
|
duration: float
|
|
|
|
# spectrum
|
|
bandwidth: float = 0.0 # span holding 99% of the power
|
|
bw_noise: float = 0.0 # span standing above the noise floor
|
|
am_depth: float = 0.0 # envelope modulation in the audio band
|
|
bw3: float = 0.0
|
|
bw20: float = 0.0
|
|
centre_offset: float = 0.0
|
|
flatness: float = 0.0
|
|
papr_spectral: float = 0.0 # dB, peak bin over median bin
|
|
carrier_ratio: float = 0.0 # fraction of power in the peak bin
|
|
symmetry: float = 0.0 # -1 all lower sideband, +1 all upper
|
|
snr_db: float = 0.0
|
|
|
|
# envelope
|
|
env_cv: float = 0.0 # std/mean of |x|
|
|
env_kurtosis: float = 0.0
|
|
ook_contrast_db: float = 0.0
|
|
ook_duty: float = 0.0
|
|
keying_regularity: float = 0.0 # do on/off runs fit a symbol grid?
|
|
is_bursty: bool = False
|
|
burst_rate_hz: float = 0.0
|
|
duty_cycle: float = 1.0
|
|
|
|
# frequency / phase
|
|
fdev_rms: float = 0.0
|
|
fdev_peak: float = 0.0
|
|
fdev_ratio: float = 0.0 # fdev_rms / bandwidth: separates AM from SSB/FM
|
|
ifreq_kurtosis: float = 0.0 # peaky (PSK) vs multimodal (FSK)
|
|
level_dwell: float = 0.0 # fraction of time the tone sits still
|
|
freq_modes: int = 0
|
|
mode_spacing: float = 0.0
|
|
psk_order: int = 0
|
|
psk_strength: float = 0.0
|
|
|
|
# timing
|
|
baud: float = 0.0
|
|
baud_strength: float = 0.0
|
|
baud_stability: float = 0.0 # does the symbol rate hold across the capture?
|
|
|
|
# audio-domain
|
|
ctcss_hz: float = 0.0
|
|
has_subaudible_data: bool = False
|
|
stereo_pilot: bool = False
|
|
afsk_1200: bool = False
|
|
audio_peak_hz: float = 0.0
|
|
|
|
extras: dict = field(default_factory=dict)
|
|
|
|
|
|
@dataclass
|
|
class Classification:
|
|
label: str
|
|
family: str
|
|
confidence: float
|
|
reasons: list[str] = field(default_factory=list)
|
|
alternatives: list[tuple[str, float]] = field(default_factory=list)
|
|
suggested_mode: str = "nfm"
|
|
features: SignalFeatures | None = None
|
|
control: "ControlChannel | None" = None
|
|
|
|
def summary(self) -> str:
|
|
pct = int(round(self.confidence * 100))
|
|
return f"{self.label} ({pct}%)"
|
|
|
|
|
|
# ---------------------------------------------------------------------------
|
|
# Feature extraction
|
|
# ---------------------------------------------------------------------------
|
|
|
|
def _otsu(values: np.ndarray, bins: int = 128) -> float:
|
|
"""Otsu threshold -- splits an on/off envelope into its two populations."""
|
|
hist, edges = np.histogram(values, bins=bins)
|
|
hist = hist.astype(np.float64)
|
|
total = hist.sum()
|
|
if total == 0:
|
|
return float(np.median(values))
|
|
centres = 0.5 * (edges[1:] + edges[:-1])
|
|
w0 = np.cumsum(hist)
|
|
w1 = total - w0
|
|
mu0 = np.cumsum(hist * centres) / np.maximum(w0, 1e-12)
|
|
grand = (hist * centres).sum()
|
|
mu1 = (grand - np.cumsum(hist * centres)) / np.maximum(w1, 1e-12)
|
|
var_between = w0 * w1 * (mu0 - mu1) ** 2
|
|
var_between[~np.isfinite(var_between)] = 0.0
|
|
return float(centres[int(np.argmax(var_between))])
|
|
|
|
|
|
# Nothing this program identifies keys slower than a few hundred baud -- the
|
|
# slowest named system is 512 baud POCSAG -- and below a couple of hundred
|
|
# hertz the transition signals are dominated by slow drift in gain and
|
|
# frequency, whose spectrum is a lobe running down to DC rather than a line.
|
|
# Searching from 40 Hz meant a bare carrier was routinely awarded a "symbol
|
|
# rate" of fifty-something baud, taken from the bottom edge of that lobe.
|
|
SYMBOL_RATE_FLOOR_HZ = 200.0
|
|
|
|
|
|
def _cyclic_line(feature: np.ndarray, fs: float,
|
|
lo_hz: float = SYMBOL_RATE_FLOOR_HZ,
|
|
hi_hz: float | None = None):
|
|
"""Find the strongest periodic line in a nonnegative feature signal.
|
|
|
|
Symbol transitions sit on a symbol-rate grid, so the transition-magnitude
|
|
signal carries a spectral line at the baud rate. Returns
|
|
``(frequency_hz, prominence_db)``.
|
|
"""
|
|
n = feature.size
|
|
if n < 256:
|
|
return 0.0, 0.0
|
|
hi_hz = hi_hz or fs / 2.5
|
|
x = feature.astype(np.float64)
|
|
x = x - x.mean()
|
|
if not np.any(x):
|
|
return 0.0, 0.0
|
|
nfft = _pow2_floor(n)
|
|
x = x[:nfft] * np.hanning(nfft)
|
|
spec = np.abs(np.fft.rfft(x, nfft))
|
|
freqs = np.fft.rfftfreq(nfft, 1.0 / fs)
|
|
|
|
band = (freqs >= lo_hz) & (freqs <= hi_hz)
|
|
if not np.any(band):
|
|
return 0.0, 0.0
|
|
sub = spec[band]
|
|
subf = freqs[band]
|
|
k = int(np.argmax(sub))
|
|
peak = sub[k]
|
|
med = np.median(sub) + 1e-12
|
|
if peak <= 0:
|
|
return 0.0, 0.0
|
|
|
|
# Symbol transitions are impulses on the symbol grid, and an impulse train
|
|
# has a *comb* of lines: the symbol rate and every multiple of it, all of
|
|
# comparable height. Taking the tallest therefore returns a harmonic as
|
|
# often as the fundamental -- 3600 baud came back as 18000. So walk back
|
|
# down the comb and report the lowest sub-multiple that still carries a
|
|
# line of its own.
|
|
df = float(freqs[1] - freqs[0])
|
|
for m in range(8, 1, -1):
|
|
want = subf[k] / m
|
|
if want < lo_hz:
|
|
continue
|
|
j = int(round((want - subf[0]) / df))
|
|
w = max(1, int(round(0.02 * want / df))) # +/- 2%, for a little drift
|
|
lo_i, hi_i = max(0, j - w), min(sub.size, j + w + 1)
|
|
if hi_i <= lo_i:
|
|
continue
|
|
window = sub[lo_i:hi_i]
|
|
height = float(window.max())
|
|
if height > 0.30 * peak and height > 6.0 * med:
|
|
i = lo_i + int(np.argmax(window))
|
|
return float(subf[i]), float(20.0 * math.log10(height / med))
|
|
return float(subf[k]), float(20.0 * math.log10(peak / med))
|
|
|
|
|
|
def _count_modes(values: np.ndarray, weights: np.ndarray | None = None,
|
|
bins: int = 96):
|
|
"""Histogram-based mode counting for FSK level detection."""
|
|
if values.size < 64:
|
|
return 0, 0.0, np.zeros(0)
|
|
lo, hi = np.percentile(values, [1.0, 99.0])
|
|
if hi <= lo:
|
|
return 0, 0.0, np.zeros(0)
|
|
# Leave room at both ends. find_peaks cannot return the first or last
|
|
# bin, and a clean two-level FSK signal puts its two modes exactly there:
|
|
# the levels *are* the 1st and 99th percentiles. Without the margin such
|
|
# a signal counts zero modes, and the cleaner it is the worse it gets.
|
|
margin = 0.10 * (hi - lo)
|
|
lo, hi = lo - margin, hi + margin
|
|
hist, edges = np.histogram(values, bins=bins, range=(lo, hi), weights=weights)
|
|
hist = hist.astype(np.float64)
|
|
if hist.sum() == 0:
|
|
return 0, 0.0, np.zeros(0)
|
|
# Light smoothing so shot noise does not create spurious modes.
|
|
kern = np.array([1.0, 3.0, 6.0, 8.0, 6.0, 3.0, 1.0])
|
|
kern /= kern.sum()
|
|
sm = np.convolve(hist, kern, mode="same")
|
|
centres = 0.5 * (edges[1:] + edges[:-1])
|
|
|
|
peaks, props = sps.find_peaks(sm, height=0.22 * sm.max(),
|
|
distance=max(3, bins // 16),
|
|
prominence=0.15 * sm.max())
|
|
if peaks.size == 0:
|
|
return 0, 0.0, np.zeros(0)
|
|
order = np.argsort(props["peak_heights"])[::-1][:8]
|
|
sel = np.sort(peaks[order])
|
|
locs = centres[sel]
|
|
if locs.size < 2:
|
|
return int(locs.size), 0.0, locs
|
|
|
|
diffs = np.diff(locs)
|
|
spacing = float(np.median(diffs))
|
|
# FSK levels are evenly spaced. If the peaks fit a uniform grid, report
|
|
# the grid size instead of the raw peak count -- one spurious shoulder
|
|
# should not turn 4-FSK into 5-FSK.
|
|
if spacing > 0:
|
|
grid = (locs - locs[0]) / spacing
|
|
if np.max(np.abs(grid - np.round(grid))) < 0.28:
|
|
n_levels = int(round(grid[-1])) + 1
|
|
if 2 <= n_levels <= 8:
|
|
return n_levels, spacing, locs
|
|
return int(locs.size), spacing, locs
|
|
|
|
|
|
def _psk_order(x: np.ndarray, fs: float):
|
|
"""Detect M-PSK by looking for the spectral line produced by x**M."""
|
|
if x.size < 1024:
|
|
return 0, 0.0
|
|
xn = x / (np.abs(x) + 1e-9) # constant-modulus, phase only
|
|
best = (0, 0.0)
|
|
for m in (2, 4, 8):
|
|
y = xn ** m
|
|
nfft = _pow2_floor(y.size, 1 << 15)
|
|
spec = np.abs(np.fft.fftshift(np.fft.fft(y[:nfft] * np.hanning(nfft), nfft)))
|
|
peak = spec.max()
|
|
med = np.median(spec) + 1e-12
|
|
if peak <= 0.0:
|
|
# Nothing in the window at all. A real receiver always has
|
|
# noise, so this only happens on a synthetic silence -- but the
|
|
# logarithm of zero is not a number and the caller was left with
|
|
# a math domain error where it expected a measurement.
|
|
continue
|
|
strength = 20.0 * math.log10(peak / med)
|
|
if strength > best[1]:
|
|
best = (m, strength)
|
|
return best
|
|
|
|
|
|
def detect_ctcss(audio: np.ndarray, fs: float):
|
|
"""Return ``(tone_hz, is_dcs_like)`` from a demodulated FM audio block."""
|
|
if audio.size < int(fs * 0.25):
|
|
return 0.0, False
|
|
n = _pow2_floor(audio.size)
|
|
x = audio[:n].astype(np.float64)
|
|
x = x - x.mean()
|
|
spec = np.abs(np.fft.rfft(x * np.hanning(n), n))
|
|
freqs = np.fft.rfftfreq(n, 1.0 / fs)
|
|
|
|
sub = (freqs >= 60.0) & (freqs <= 260.0)
|
|
voice = (freqs >= 300.0) & (freqs <= 3000.0)
|
|
if not np.any(sub):
|
|
return 0.0, False
|
|
sub_spec = spec[sub]
|
|
sub_f = freqs[sub]
|
|
k = int(np.argmax(sub_spec))
|
|
peak_f = float(sub_f[k])
|
|
peak_v = float(sub_spec[k])
|
|
floor = float(np.median(spec[voice])) + 1e-12 if np.any(voice) else 1e-12
|
|
|
|
if peak_v / floor < 6.0:
|
|
return 0.0, False
|
|
# A CTCSS tone is a single sharp line; DCS is a 134.4 bps square wave and
|
|
# spreads its energy across the whole subaudible region.
|
|
band_energy = float(np.sum(sub_spec ** 2))
|
|
tone_energy = float(np.sum(sub_spec[max(0, k - 2):k + 3] ** 2))
|
|
if tone_energy / max(band_energy, 1e-12) < 0.35:
|
|
return 0.0, True
|
|
|
|
nearest = min(CTCSS_TONES, key=lambda t: abs(t - peak_f))
|
|
if abs(nearest - peak_f) <= max(1.5, 0.02 * nearest):
|
|
return float(nearest), False
|
|
return 0.0, False
|
|
|
|
|
|
def _center_and_filter(x: np.ndarray, sample_rate: float, offset_hz: float,
|
|
bw_hz: float, wide_bw_hz: float = 0.0) -> np.ndarray:
|
|
"""Shift the signal to DC and low-pass it to its own occupied bandwidth.
|
|
|
|
Time-domain features (envelope, discriminator, phase) are meaningless when
|
|
they are dominated by noise from the rest of the IF, so every measurement
|
|
after the spectrum step runs on this filtered copy.
|
|
"""
|
|
# Be generous: a carrier-dominated signal (AM) has a small 99%-power
|
|
# bandwidth but sidebands well outside it, so take the wider of the two
|
|
# measures and leave headroom on top.
|
|
keep = max(bw_hz, wide_bw_hz) * 1.8
|
|
keep = max(keep, sample_rate / 200.0)
|
|
if keep >= sample_rate * 0.9 and abs(offset_hz) < sample_rate / 50.0:
|
|
return x # already occupies most of the band; filtering buys nothing
|
|
|
|
if abs(offset_hz) > sample_rate / 1000.0:
|
|
n = np.arange(x.size, dtype=np.float64)
|
|
x = (x * np.exp(-2j * math.pi * offset_hz * n / sample_rate)).astype(np.complex64)
|
|
|
|
# firwin's cutoff is in units of Nyquist; we want +-keep/2 around DC.
|
|
norm = (keep / 2.0) / (sample_rate / 2.0)
|
|
if norm >= 0.95:
|
|
return x
|
|
ntaps = 127
|
|
taps = sps.firwin(ntaps, norm).astype(np.float64)
|
|
y = sps.lfilter(taps, [1.0], x).astype(np.complex64)
|
|
# Drop the filter's start-up ramp: it looks exactly like a signal fading
|
|
# in, which would otherwise register as on-off keying.
|
|
return y[ntaps:] if y.size > 4 * ntaps else y
|
|
|
|
|
|
def _robust_floor(psd_db: np.ndarray) -> tuple[float, float]:
|
|
"""Noise floor and spread of one spectrum, by sigma clipping.
|
|
|
|
Repeatedly drops the bins that stand out until only the noise population
|
|
is left. A sliding percentile is the right tool while sweeping, where
|
|
signals are narrow slivers of a wide span, but not here: by this point
|
|
the capture is centred on one signal that may fill most of the analysis
|
|
band, and a windowed percentile would sit on the signal itself.
|
|
"""
|
|
v = np.asarray(psd_db, dtype=np.float64)
|
|
if v.size < 8:
|
|
return float(np.median(v)), 1.0
|
|
mask = np.ones(v.size, dtype=bool)
|
|
floor = float(np.median(v))
|
|
sigma = 1.0
|
|
for _ in range(6):
|
|
sel = v[mask]
|
|
if sel.size < max(8, v.size // 10):
|
|
break
|
|
floor = float(np.median(sel))
|
|
sigma = float(1.4826 * np.median(np.abs(sel - floor))) or 1.0
|
|
new_mask = v < floor + max(3.0 * sigma, 3.0)
|
|
if new_mask.sum() < max(8, v.size // 10):
|
|
break
|
|
if np.array_equal(new_mask, mask):
|
|
break
|
|
mask = new_mask
|
|
return floor, max(sigma, 0.3)
|
|
|
|
|
|
def _occupied_span(psd_db: np.ndarray, bin_hz: float,
|
|
margin_db: float = 0.0) -> tuple[float, float]:
|
|
"""Occupied bandwidth as an analyst reads it off a spectrum display.
|
|
|
|
Returns ``(bandwidth_hz, centre_offset_hz)`` for the contiguous run of
|
|
bins standing above the noise floor around the strongest peak.
|
|
|
|
A 99%-of-power measure cannot be used here: AM puts almost all of its
|
|
power in the carrier, so 99% of the power lives in a single bin and the
|
|
channel would be reported as tens of hertz wide.
|
|
"""
|
|
n = psd_db.size
|
|
if n < 8:
|
|
return bin_hz, 0.0
|
|
floor, sigma = _robust_floor(psd_db)
|
|
threshold = floor + (margin_db or max(6.0, 4.0 * sigma))
|
|
above = psd_db > threshold
|
|
if not np.any(above):
|
|
return bin_hz, 0.0
|
|
|
|
idx = np.flatnonzero(above)
|
|
# 2nd-to-98th percentile of the *positions* that stand above the noise.
|
|
# Counting positions rather than weighting them by power matters: a
|
|
# carrier holds so much more power than its sidebands that a weighted
|
|
# measure collapses onto the carrier bin and reports a few hertz.
|
|
lo = int(np.percentile(idx, 2))
|
|
hi = int(np.percentile(idx, 98))
|
|
if hi < lo:
|
|
lo, hi = hi, lo
|
|
w = np.maximum(psd_db[idx] - floor, 1e-9)
|
|
centroid = float(np.dot(idx.astype(np.float64), w) / w.sum())
|
|
offset = (centroid - (n - 1) / 2.0) * bin_hz
|
|
return max(bin_hz, float((hi - lo + 1) * bin_hz)), offset
|
|
|
|
|
|
def _run_spans(on: np.ndarray):
|
|
"""Yield ``(state, start, length)`` for each constant run of a mask."""
|
|
if on.size == 0:
|
|
return
|
|
change = np.flatnonzero(np.diff(on.astype(np.int8)))
|
|
bounds = np.concatenate(([0], change + 1, [on.size]))
|
|
for i in range(bounds.size - 1):
|
|
yield bool(on[bounds[i]]), int(bounds[i]), int(bounds[i + 1] - bounds[i])
|
|
|
|
|
|
def _longest_active_run(on: np.ndarray) -> tuple[int, int]:
|
|
"""Start and length of the longest contiguous key-down / active stretch."""
|
|
best = (0, 0)
|
|
for state, start, length in _run_spans(on):
|
|
if state and length > best[1]:
|
|
best = (start, length)
|
|
return best
|
|
|
|
|
|
def extract_features(x: np.ndarray, sample_rate: float,
|
|
snr_db: float = 0.0) -> SignalFeatures:
|
|
"""Compute the full feature vector for one captured block."""
|
|
x = np.asarray(x, dtype=np.complex64)
|
|
n = x.size
|
|
f = SignalFeatures(sample_rate=float(sample_rate), n_samples=n,
|
|
duration=n / float(sample_rate), snr_db=float(snr_db))
|
|
if n < 512:
|
|
return f
|
|
|
|
# ---- keying structure, measured across the whole capture ----------
|
|
env_full = np.abs(x).astype(np.float64)
|
|
smooth_n = max(4, int(sample_rate / 4000.0))
|
|
env_s = np.convolve(env_full, np.ones(smooth_n) / smooth_n, mode="same")
|
|
thr = _otsu(env_s)
|
|
on = env_s > thr
|
|
f.ook_duty = float(on.mean())
|
|
f.duty_cycle = f.ook_duty
|
|
if np.any(on) and np.any(~on):
|
|
hi = float(np.mean(env_s[on]))
|
|
lo = float(np.mean(env_s[~on])) + 1e-12
|
|
f.ook_contrast_db = float(20.0 * math.log10(hi / lo))
|
|
|
|
# Do the on and off runs land on a common grid? Keyed data and Morse
|
|
# both quantise to a symbol or dot length; a signal drifting across the
|
|
# squelch threshold produces runs of every length, which is what tells
|
|
# real keying apart from a fading carrier.
|
|
spans = [ln for _, _, ln in _run_spans(on)]
|
|
if len(spans) >= 6:
|
|
arr = np.array(spans[1:-1] if len(spans) > 8 else spans, dtype=np.float64)
|
|
unit = float(np.percentile(arr, 20))
|
|
# The unit has to be resolvable. When the envelope merely jitters
|
|
# across the threshold the runs are one or two samples long, and every
|
|
# length is then trivially an "integer multiple" of one sample -- a
|
|
# test that noise passes perfectly.
|
|
if unit >= 8.0:
|
|
ratios = arr / unit
|
|
usable = ratios[ratios <= 12.0]
|
|
if usable.size >= 4:
|
|
f.keying_regularity = float(
|
|
np.mean(np.abs(usable - np.round(usable)) < 0.2))
|
|
|
|
# Burst / TDMA structure: how often does the envelope gate on and off?
|
|
if 0.02 < f.ook_duty < 0.98:
|
|
edges = np.diff(on.astype(np.int8))
|
|
rises = np.where(edges > 0)[0]
|
|
if rises.size >= 3:
|
|
periods = np.diff(rises) / sample_rate
|
|
med = float(np.median(periods))
|
|
if med > 0:
|
|
f.burst_rate_hz = 1.0 / med
|
|
f.is_bursty = bool(np.std(periods) / med < 0.5)
|
|
|
|
# ---- pick the stretch to characterise ------------------------------
|
|
# Speech has pauses, and on SSB the carrier disappears with them. Judging
|
|
# modulation across the silence would describe the silence: an AM voice
|
|
# channel reads as a bare carrier, SSB voice reads as on-off keying. So
|
|
# measure the longest continuously-active stretch instead -- provided it
|
|
# is long enough to be a transmission rather than a data symbol.
|
|
seg = x
|
|
if 0.05 < f.ook_duty < 0.92:
|
|
s_start, s_len = _longest_active_run(on)
|
|
if s_len >= max(int(0.3 * sample_rate), 4096) and s_len < int(0.92 * n):
|
|
seg = x[s_start:s_start + s_len]
|
|
f.extras["analysed_seconds"] = round(s_len / sample_rate, 3)
|
|
f.extras["analysed_fraction"] = round(s_len / n, 3)
|
|
|
|
# ---- spectrum ------------------------------------------------------
|
|
nfft = min(4096, 1 << int(math.floor(math.log2(seg.size))))
|
|
freqs, psd = welch_psd(seg, nfft)
|
|
bin_hz = sample_rate / nfft
|
|
psd_db = db(psd)
|
|
|
|
# Primary bandwidth is the span holding 99% of the power. The span
|
|
# standing above the noise is kept alongside it because the two disagree
|
|
# in a useful way: AM puts nearly all its power in the carrier, so a large
|
|
# gap between them is itself evidence of a carrier-dominated signal.
|
|
f.bandwidth, f.centre_offset = occupied_bandwidth(psd, bin_hz, 0.99)
|
|
f.bw_noise, _ = _occupied_span(psd_db, bin_hz)
|
|
total = psd.sum()
|
|
peak_lin = psd.max()
|
|
f.carrier_ratio = float(peak_lin / total) if total > 0 else 0.0
|
|
f.papr_spectral = float(psd_db.max() - np.median(psd_db))
|
|
f.flatness = spectral_flatness(psd)
|
|
|
|
peak_db = psd_db.max()
|
|
for lvl, attr in ((3.0, "bw3"), (20.0, "bw20")):
|
|
above = psd_db >= (peak_db - lvl)
|
|
if np.any(above):
|
|
idx = np.where(above)[0]
|
|
setattr(f, attr, float((idx[-1] - idx[0] + 1) * bin_hz))
|
|
|
|
# Sideband asymmetry measured about the strongest bin and only across the
|
|
# occupied band -- comparing the two halves of the whole IF just measures
|
|
# where the noise sits.
|
|
pk = int(np.argmax(psd))
|
|
span = max(2, int(f.bandwidth / bin_hz))
|
|
lo_i, hi_i = max(0, pk - span), min(nfft, pk + span + 1)
|
|
lower = float(psd[lo_i:pk].sum())
|
|
upper = float(psd[pk + 1:hi_i].sum())
|
|
if lower + upper > 0:
|
|
f.symmetry = (upper - lower) / (upper + lower)
|
|
|
|
# Every time-domain measurement below runs on the signal alone.
|
|
# Size the analysis filter from the widest honest estimate. Using the
|
|
# 99%-power figure alone would band-limit an AM channel to its carrier and
|
|
# every later measurement would describe a dead carrier.
|
|
keep_bw = max(f.bandwidth, f.bw_noise, f.bw20)
|
|
xf = _center_and_filter(seg, sample_rate, f.centre_offset, keep_bw, f.bw20)
|
|
f.extras["_filtered"] = xf
|
|
seg = xf
|
|
|
|
# ---- envelope -----------------------------------------------------
|
|
env = np.abs(seg).astype(np.float64)
|
|
mean_env = float(env.mean())
|
|
if mean_env > 0:
|
|
f.env_cv = float(env.std() / mean_env)
|
|
centred = env - mean_env
|
|
var = float(centred.var())
|
|
if var > 0:
|
|
f.env_kurtosis = float(np.mean(centred ** 4) / (var ** 2) - 3.0)
|
|
|
|
# Does the envelope carry audio? A bare carrier's envelope is flat, an
|
|
# FM carrier's is flat by construction, and an AM or SSB voice channel's
|
|
# envelope *is* the speech. This is what tells a modulated AM channel
|
|
# apart from a dead carrier, which the 99%-power bandwidth cannot do
|
|
# because the carrier holds nearly all the power either way.
|
|
if mean_env > 0 and env.size > 1024:
|
|
ac = env - mean_env
|
|
dec = max(1, int(sample_rate / 16_000))
|
|
a = ac[::dec]
|
|
fs_a = sample_rate / dec
|
|
n_a = 1 << int(math.floor(math.log2(max(256, min(a.size, 1 << 15)))))
|
|
if a.size >= n_a:
|
|
spec = np.abs(np.fft.rfft(a[:n_a] * np.hanning(n_a), n_a)) ** 2
|
|
fr = np.fft.rfftfreq(n_a, 1.0 / fs_a)
|
|
band = (fr >= 100.0) & (fr <= min(4000.0, fs_a / 2.2))
|
|
if np.any(band):
|
|
f.am_depth = float(np.sqrt(np.sum(spec[band])) /
|
|
(mean_env * n_a / 2.0))
|
|
|
|
# ---- instantaneous frequency --------------------------------------
|
|
ifreq = instantaneous_frequency(seg, sample_rate)
|
|
if ifreq.size:
|
|
w = env[1:]
|
|
strong = w > (0.5 * mean_env) if mean_env > 0 else np.ones_like(w, bool)
|
|
sel = ifreq[strong] if strong.sum() > 64 else ifreq
|
|
f.fdev_rms = float(np.std(sel))
|
|
f.fdev_peak = float(np.percentile(np.abs(sel - np.mean(sel)), 99.0))
|
|
f.fdev_ratio = f.fdev_rms / max(f.bandwidth, 1.0)
|
|
centred = sel - np.mean(sel)
|
|
var = float(centred.var())
|
|
if var > 0:
|
|
# Negative => the discriminator sits on discrete levels (FSK);
|
|
# strongly positive => it is flat with impulses (PSK phase jumps).
|
|
f.ifreq_kurtosis = float(np.mean(centred ** 4) / (var ** 2) - 3.0)
|
|
|
|
# How much of the time does the tone hold still? FSK parks on a level
|
|
# for a whole symbol and jumps between them, so its slope is zero
|
|
# almost everywhere. Tone-modulated or voice FM sweeps continuously
|
|
# and is never flat for long. Without this, the two peaks that
|
|
# sinusoidal FM puts at the ends of its swing (a CTCSS tone during a
|
|
# speech pause, say) read as a pair of FSK levels.
|
|
slope = np.diff(sel)
|
|
if slope.size > 64:
|
|
scale = float(np.std(slope))
|
|
if scale > 0:
|
|
f.level_dwell = float(np.mean(np.abs(slope) < 0.25 * scale))
|
|
modes, spacing, _ = _count_modes(sel, weights=None)
|
|
f.freq_modes = modes
|
|
f.mode_spacing = spacing
|
|
|
|
# ---- phase --------------------------------------------------------
|
|
order, strength = _psk_order(seg, sample_rate)
|
|
f.psk_order, f.psk_strength = order, strength
|
|
|
|
# ---- symbol rate ---------------------------------------------------
|
|
# Transition magnitude: for FSK the discriminator's, for anything that
|
|
# keys its amplitude the envelope's. Both are impulse trains sitting on
|
|
# the symbol grid, which is what puts a line at the symbol rate.
|
|
#
|
|
# The squared envelope used to be offered as a third candidate. It is not
|
|
# a transition signal: for data-carrying NRZ its spectrum is a sinc lobe
|
|
# peaking at DC, so the lowest bin in the search band beat the band median
|
|
# by 40-odd dB and won every time. Every random OOK signal measured about
|
|
# 90 baud because of it.
|
|
cand = []
|
|
if ifreq.size > 512:
|
|
cand.append((_cyclic_line(np.abs(np.diff(ifreq)), sample_rate),
|
|
np.abs(np.diff(ifreq))))
|
|
cand.append((_cyclic_line(np.abs(np.diff(env)), sample_rate),
|
|
np.abs(np.diff(env))))
|
|
(baud, strength), winner = max(cand, key=lambda c: c[0][1])
|
|
if strength > 8.0:
|
|
f.baud, f.baud_strength = baud, strength
|
|
# A real symbol rate is a property of the transmission and holds for
|
|
# its whole length. The estimator always returns its best peak, so
|
|
# on noise it returns a different answer for each half of the same
|
|
# capture -- which is exactly how to tell the two apart.
|
|
half = winner.size // 2
|
|
if half > 1024:
|
|
b1, _ = _cyclic_line(winner[:half], sample_rate)
|
|
b2, _ = _cyclic_line(winner[half:], sample_rate)
|
|
if b1 > 0 and b2 > 0:
|
|
f.baud_stability = float(
|
|
1.0 - abs(b1 - b2) / max(b1, b2))
|
|
|
|
return f
|
|
|
|
|
|
def _analyse_fm_audio(x: np.ndarray, sample_rate: float, f: SignalFeatures) -> None:
|
|
"""Fill in the audio-domain features that need an FM demodulation."""
|
|
ifreq = instantaneous_frequency(x, sample_rate)
|
|
if ifreq.size < 1024:
|
|
return
|
|
|
|
# 19 kHz stereo pilot -> broadcast FM.
|
|
if sample_rate > 60_000:
|
|
n = _pow2_floor(ifreq.size)
|
|
spec = np.abs(np.fft.rfft((ifreq[:n] - ifreq[:n].mean()) * np.hanning(n), n))
|
|
freqs = np.fft.rfftfreq(n, 1.0 / sample_rate)
|
|
near = (freqs > 18_800) & (freqs < 19_200)
|
|
ref = (freqs > 22_000) & (freqs < 30_000)
|
|
if np.any(near) and np.any(ref):
|
|
f.stereo_pilot = bool(spec[near].max() > 8.0 * (np.median(spec[ref]) + 1e-12))
|
|
|
|
# Decimate the discriminator output to a voice rate for tone work.
|
|
dec = max(1, int(sample_rate // 16_000))
|
|
audio = sps.decimate(ifreq, dec, ftype="fir", zero_phase=False) if dec > 1 else ifreq
|
|
fs_a = sample_rate / dec
|
|
|
|
tone, dcs = detect_ctcss(audio, fs_a)
|
|
f.ctcss_hz = tone
|
|
f.has_subaudible_data = dcs
|
|
|
|
n = _pow2_floor(audio.size, 1 << 15)
|
|
if n >= 1024:
|
|
a = audio[:n].astype(np.float64)
|
|
a -= a.mean()
|
|
spec = np.abs(np.fft.rfft(a * np.hanning(n), n))
|
|
freqs = np.fft.rfftfreq(n, 1.0 / fs_a)
|
|
band = (freqs > 250.0) & (freqs < 3500.0)
|
|
if np.any(band):
|
|
sb, sf = spec[band], freqs[band]
|
|
f.audio_peak_hz = float(sf[int(np.argmax(sb))])
|
|
# AFSK1200 (APRS, Bell 202) sits on 1200 Hz and 2200 Hz marks.
|
|
def energy(target, width=90.0):
|
|
m = (sf > target - width) & (sf < target + width)
|
|
return float(sb[m].max()) if np.any(m) else 0.0
|
|
base = float(np.median(sb)) + 1e-12
|
|
e12, e22 = energy(1200.0), energy(2200.0)
|
|
f.afsk_1200 = bool(e12 > 5 * base and e22 > 5 * base)
|
|
|
|
|
|
# ---------------------------------------------------------------------------
|
|
# Rule engine
|
|
# ---------------------------------------------------------------------------
|
|
|
|
def _score_rules(f: SignalFeatures) -> list[tuple[str, str, float, str, str]]:
|
|
"""Return ``(label, family, score, reason, suggested_mode)`` candidates."""
|
|
out = []
|
|
|
|
def add(label, family, score, reason, mode="nfm"):
|
|
if score > 0:
|
|
out.append((label, family, float(score), reason, mode))
|
|
|
|
const_env = f.env_cv < 0.35
|
|
narrow = f.bandwidth < 30_000
|
|
very_narrow = f.bandwidth < 1_500
|
|
# A dead carrier occupies essentially no bandwidth, whatever the noise does
|
|
# to its measured phase jitter.
|
|
tone_bw = max(60.0, f.sample_rate / 400.0)
|
|
|
|
# -- unmodulated carrier / tone -------------------------------------
|
|
if f.bandwidth < tone_bw and f.env_cv < 0.20 \
|
|
and (f.ook_contrast_db < 6 or f.ook_duty > 0.93):
|
|
add("Unmodulated carrier", "carrier",
|
|
0.72 + min(0.2, f.papr_spectral / 200.0),
|
|
f"single steady line only {f.bandwidth:.0f} Hz wide, no keying or "
|
|
f"modulation sidebands", "cw")
|
|
|
|
# -- CW / on-off keying ---------------------------------------------
|
|
if f.ook_contrast_db > 10 and very_narrow and 0.05 < f.ook_duty < 0.85:
|
|
base = 0.55 + min(0.3, f.ook_contrast_db / 60.0)
|
|
if f.baud and f.baud < 60:
|
|
base += 0.1
|
|
add("CW / Morse (on-off keyed carrier)", "cw", base,
|
|
f"keyed carrier, {f.ook_contrast_db:.0f} dB on/off contrast, "
|
|
f"{f.bandwidth:.0f} Hz wide", "cw")
|
|
elif f.ook_contrast_db > 12 and f.bandwidth < 60_000 and 0.02 < f.ook_duty < 0.9:
|
|
add("OOK / ASK data burst", "ook",
|
|
0.5 + min(0.25, f.ook_contrast_db / 80.0),
|
|
f"on-off keying, {f.ook_contrast_db:.0f} dB contrast, "
|
|
f"{f.baud:.0f} baud" if f.baud else "on-off keying", "raw")
|
|
|
|
# -- FSK families ----------------------------------------------------
|
|
# Discrete levels, and the tone actually rests on them.
|
|
level_like = f.ifreq_kurtosis < 1.5 and f.level_dwell > 0.32
|
|
if const_env and f.freq_modes >= 2 and f.fdev_rms > 200 and level_like:
|
|
levels = f.freq_modes
|
|
if levels in (2, 3):
|
|
name, conf = "2-FSK (binary FSK)", 0.6
|
|
elif levels == 4:
|
|
name, conf = "4-FSK / C4FM", 0.68
|
|
elif levels in (5, 6, 7, 8):
|
|
name, conf = f"{levels}-level FSK", 0.5
|
|
else:
|
|
name, conf = "Multi-level FSK", 0.4
|
|
if f.baud_strength > 12:
|
|
conf += 0.12
|
|
reason = (f"constant envelope, {levels} discriminator levels "
|
|
f"{f.mode_spacing:.0f} Hz apart")
|
|
if f.baud:
|
|
reason += f", ~{f.baud:.0f} baud"
|
|
add(name, "fsk", conf, reason, "nfm")
|
|
|
|
# -- analogue FM -----------------------------------------------------
|
|
# A quiet FM channel carrying only a CTCSS tone has very little deviation,
|
|
# so this floor has to sit low.
|
|
if const_env and f.fdev_rms > 120 and not level_like:
|
|
if f.bandwidth > 100_000:
|
|
conf = 0.72 + (0.15 if f.stereo_pilot else 0.0)
|
|
label = "Wideband FM (broadcast)"
|
|
if f.stereo_pilot:
|
|
label = "Wideband FM broadcast (stereo, 19 kHz pilot)"
|
|
add(label, "wfm", conf,
|
|
f"{f.bandwidth/1e3:.0f} kHz wide, {f.fdev_rms/1e3:.1f} kHz rms deviation",
|
|
"wfm")
|
|
elif narrow:
|
|
conf = 0.6
|
|
if f.ctcss_hz:
|
|
conf += 0.2
|
|
if 300 < f.audio_peak_hz < 3200:
|
|
conf += 0.08
|
|
label = "Narrowband FM voice"
|
|
if f.ctcss_hz:
|
|
label += f" (CTCSS {f.ctcss_hz:.1f} Hz)"
|
|
elif f.has_subaudible_data:
|
|
label += " (DCS subaudible data)"
|
|
add(label, "nfm", conf,
|
|
f"{f.bandwidth/1e3:.1f} kHz wide, {f.fdev_rms/1e3:.1f} kHz rms deviation",
|
|
"nfm")
|
|
|
|
# -- AM ---------------------------------------------------------------
|
|
# AM keeps its carrier, so the phase hardly moves and the sidebands are
|
|
# mirror images. That is exactly what separates it from SSB.
|
|
# Real AM often runs at modest modulation depth, so the envelope only has
|
|
# to vary more than receiver noise alone would make it vary.
|
|
snr_lin = 10.0 ** (max(f.snr_db, 0.0) / 10.0)
|
|
env_noise = 1.0 / math.sqrt(max(2.0, snr_lin))
|
|
am_modulated = f.env_cv > max(0.045, 2.0 * env_noise)
|
|
if (am_modulated and narrow and f.bandwidth > 4 * tone_bw
|
|
and f.carrier_ratio > 0.08
|
|
and f.fdev_ratio < 0.18 and abs(f.symmetry) < 0.60):
|
|
conf = 0.55 + min(0.2, 2.0 * f.env_cv) - 0.2 * abs(f.symmetry)
|
|
add("AM (amplitude modulation)", "am", conf,
|
|
f"symmetric sidebands around a surviving carrier, "
|
|
f"{f.env_cv:.2f} envelope variation, {f.bandwidth/1e3:.1f} kHz wide", "am")
|
|
|
|
# -- SSB ---------------------------------------------------------------
|
|
if (300 < f.bandwidth < 6_500 and f.env_cv > 0.35 and f.fdev_ratio > 0.15
|
|
and f.ook_contrast_db < 14):
|
|
# Which sideband cannot be recovered once we have centred on the
|
|
# signal, so fall back on the HF convention (LSB below 10 MHz).
|
|
conf = 0.52 + min(0.2, f.fdev_ratio) + min(0.15, abs(f.symmetry) * 0.2)
|
|
add("SSB voice (suppressed carrier)", "ssb", conf,
|
|
f"no carrier line, phase swings across the full {f.bandwidth:.0f} Hz "
|
|
f"of audio bandwidth", "usb")
|
|
|
|
# -- PSK ---------------------------------------------------------------
|
|
impulsive_phase = f.ifreq_kurtosis > 3.0
|
|
fsk_like = f.freq_modes >= 2 and f.mode_spacing > 0 and level_like
|
|
# Raising to the Mth power must *create* the spectral line. An
|
|
# unmodulated carrier -- including the silent gaps between phrases on an
|
|
# FM channel -- already has a line at every power, and would otherwise
|
|
# look like textbook PSK.
|
|
psk_line_created = f.psk_strength > f.papr_spectral + 8.0
|
|
if const_env and f.psk_strength > 20 and f.psk_order and impulsive_phase \
|
|
and psk_line_created and f.carrier_ratio < 0.10 \
|
|
and not fsk_like and f.bandwidth > tone_bw:
|
|
m = f.psk_order
|
|
name = {2: "BPSK", 4: "QPSK / pi-4 DQPSK", 8: "8-PSK"}.get(m, f"{m}-PSK")
|
|
conf = 0.5 + min(0.25, (f.psk_strength - 20) / 60.0)
|
|
reason = f"x^{m} produces a spectral line ({f.psk_strength:.0f} dB)"
|
|
if f.baud:
|
|
reason += f", ~{f.baud:.0f} baud"
|
|
add(name, "psk", conf, reason, "raw")
|
|
|
|
# -- wideband digital / noise-like ------------------------------------
|
|
if f.flatness > 0.55 and f.bandwidth > 200_000 and f.carrier_ratio < 0.05:
|
|
add("Wideband digital carrier (OFDM/CDMA-like)", "digital",
|
|
0.5 + min(0.25, f.flatness - 0.55),
|
|
f"flat, noise-like spectrum {f.bandwidth/1e6:.2f} MHz wide", "raw")
|
|
|
|
# -- pulsed / radar-like ----------------------------------------------
|
|
if f.ook_duty < 0.05 and f.ook_contrast_db > 15 and f.bandwidth > 200_000:
|
|
add("Pulsed transmission (radar / DME / Mode S-like)", "pulse",
|
|
0.5 + min(0.2, f.ook_contrast_db / 100.0),
|
|
f"{f.ook_duty*100:.1f}% duty cycle, wide pulses", "raw")
|
|
|
|
if not out:
|
|
add("Unidentified signal", "unknown", 0.2,
|
|
f"{f.bandwidth/1e3:.1f} kHz wide, SNR {f.snr_db:.0f} dB", "nfm")
|
|
return out
|
|
|
|
|
|
# ---------------------------------------------------------------------------
|
|
# Known-system refinement
|
|
# ---------------------------------------------------------------------------
|
|
|
|
def _baud_near(f: SignalFeatures, target: float, tol: float = 0.12) -> bool:
|
|
return bool(f.baud) and abs(f.baud - target) <= tol * target
|
|
|
|
|
|
def _identify_system(freq_hz: float, f: SignalFeatures,
|
|
best_family: str) -> tuple[str, float, str] | None:
|
|
"""Name a specific system when frequency + shape + baud all agree."""
|
|
mhz = freq_hz / 1e6
|
|
|
|
def within(lo, hi):
|
|
return lo <= mhz <= hi
|
|
|
|
# Aviation surveillance
|
|
if within(1089, 1091) and f.ook_duty < 0.2:
|
|
return ("ADS-B / Mode S (1090 MHz extended squitter)", 0.9,
|
|
"1090 MHz, pulse-position keyed bursts")
|
|
if within(977, 979):
|
|
return ("UAT ADS-B / FIS-B (978 MHz)", 0.85, "978 MHz UAT channel")
|
|
if within(960, 1215) and f.ook_duty < 0.1:
|
|
return ("DME / TACAN pulse pairs", 0.7, "pulsed navigation band")
|
|
|
|
# Broadcast
|
|
if within(87.9, 108.1) and best_family in ("wfm", "fsk", "nfm"):
|
|
label = "FM broadcast station"
|
|
if f.stereo_pilot:
|
|
label += " (stereo)"
|
|
return (label, 0.9, "FM broadcast band, wideband FM")
|
|
if within(162.39, 162.56) and best_family == "nfm":
|
|
return ("NOAA Weather Radio", 0.9, "NWR channel, narrowband FM")
|
|
if within(136.9, 138.1) and f.bandwidth > 25_000:
|
|
return ("Weather satellite downlink (NOAA APT / Meteor)", 0.7,
|
|
"137 MHz satellite band")
|
|
|
|
# Marine / maritime
|
|
if within(161.96, 162.04) and _baud_near(f, 9600, 0.2):
|
|
return ("AIS ship transponder (9600 GMSK)", 0.88,
|
|
"AIS channel A/B, 9600 baud GMSK")
|
|
if within(156.0, 162.1) and best_family == "nfm":
|
|
ch = _marine_channel(freq_hz)
|
|
return (f"Marine VHF voice{ch}", 0.75, "marine VHF band, narrowband FM")
|
|
|
|
# Aviation voice / data
|
|
if within(118.0, 137.0):
|
|
if within(129.0, 137.0) and _baud_near(f, 2400, 0.2) and f.ook_duty < 0.6:
|
|
return ("ACARS datalink (2400 baud MSK)", 0.85,
|
|
"ACARS band, 2400 baud bursts")
|
|
if best_family == "am":
|
|
return ("VHF airband voice (AM)", 0.85, "118-137 MHz airband, AM")
|
|
if within(225.0, 400.0) and best_family == "am":
|
|
return ("Military UHF air voice (AM)", 0.7, "225-400 MHz UHF air band")
|
|
|
|
# Amateur
|
|
if within(144.38, 144.40) and (f.afsk_1200 or _baud_near(f, 1200, 0.2)):
|
|
return ("APRS packet (AFSK 1200 baud)", 0.88, "144.390 MHz APRS channel")
|
|
ham_hf = any(lo <= mhz <= hi for lo, hi in (
|
|
(1.8, 2.0), (3.5, 4.0), (5.33, 5.41), (7.0, 7.3), (10.1, 10.15),
|
|
(14.0, 14.35), (18.068, 18.168), (21.0, 21.45), (24.89, 24.99),
|
|
(28.0, 29.7)))
|
|
if ham_hf and best_family == "ssb":
|
|
return ("Amateur HF SSB voice", 0.75, "inside a US amateur HF phone band")
|
|
if ham_hf and best_family == "cw":
|
|
return ("Amateur HF CW (Morse)", 0.8, "inside a US amateur HF CW segment")
|
|
if within(26.965, 27.405) and best_family in ("am", "ssb"):
|
|
return ("CB radio (Citizens Band)", 0.75, "11 m CB channel")
|
|
if within(144.0, 148.0) and best_family == "nfm":
|
|
return ("2 m amateur FM", 0.7, "2 m band, narrowband FM")
|
|
if within(420.0, 450.0) and best_family == "nfm":
|
|
return ("70 cm amateur FM", 0.7, "70 cm band, narrowband FM")
|
|
|
|
# Land mobile digital voice
|
|
if best_family == "fsk" and f.freq_modes == 4:
|
|
if _baud_near(f, 4800) and 8_000 < f.bandwidth < 16_000:
|
|
if f.is_bursty and 25.0 < f.burst_rate_hz < 45.0:
|
|
return ("DMR digital voice (TDMA, 4800 baud C4FM)", 0.8,
|
|
"12.5 kHz 4-FSK with ~30 ms TDMA bursts")
|
|
return ("P25 Phase 1 C4FM digital voice (4800 baud)", 0.75,
|
|
"12.5 kHz 4-level FSK at 4800 baud")
|
|
if _baud_near(f, 2400) and f.bandwidth < 8_000:
|
|
return ("NXDN digital voice (2400 baud, 6.25 kHz)", 0.72,
|
|
"6.25 kHz 4-FSK at 2400 baud")
|
|
if best_family == "fsk" and f.freq_modes <= 3:
|
|
if _baud_near(f, 4800) and f.bandwidth < 8_000:
|
|
return ("D-STAR digital voice (4800 baud GMSK)", 0.65,
|
|
"6.25 kHz GMSK at 4800 baud")
|
|
for baud, name in ((512, "POCSAG 512"), (1200, "POCSAG 1200"),
|
|
(2400, "POCSAG 2400")):
|
|
if _baud_near(f, baud) and (within(929, 932) or within(150, 160)):
|
|
return (f"{name} pager traffic", 0.75,
|
|
f"paging band, {baud} baud 2-FSK")
|
|
for baud, name in ((1600, "FLEX 1600"), (3200, "FLEX 3200"),
|
|
(6400, "FLEX 6400")):
|
|
if _baud_near(f, baud) and within(929, 932):
|
|
return (f"{name} pager traffic", 0.72,
|
|
f"900 MHz paging, {baud} baud FLEX")
|
|
|
|
if best_family == "psk" and f.psk_order == 4 and _baud_near(f, 18000, 0.15):
|
|
return ("TETRA (pi/4-DQPSK, 18 kbaud)", 0.7, "25 kHz pi/4-DQPSK")
|
|
|
|
# ISM / short range devices
|
|
if within(314.5, 315.5) or within(433.0, 434.9) or within(389.5, 390.5):
|
|
if best_family in ("ook", "fsk"):
|
|
kind = "OOK" if best_family == "ook" else "FSK"
|
|
return (f"ISM short-range device ({kind}: TPMS / remote / sensor)",
|
|
0.7, f"ISM band burst, {kind}")
|
|
if within(902, 928) and f.bandwidth > 100_000:
|
|
return ("902-928 MHz ISM (FHSS / LoRa / smart meter)", 0.6,
|
|
"wideband ISM emission")
|
|
|
|
# Time signals
|
|
if abs(mhz - 2.5) < 0.005 or abs(mhz - 5.0) < 0.005 or \
|
|
abs(mhz - 10.0) < 0.005 or abs(mhz - 15.0) < 0.005 or \
|
|
abs(mhz - 20.0) < 0.005:
|
|
return ("WWV/WWVH standard time and frequency station", 0.8,
|
|
"exact WWV carrier frequency")
|
|
|
|
# Cellular
|
|
if (within(824, 894) or within(1850, 1990) or within(614, 698)) \
|
|
and f.bandwidth > 800_000:
|
|
return ("Cellular downlink (LTE/5G-NR)", 0.65,
|
|
"wide flat carrier in a cellular allocation")
|
|
return None
|
|
|
|
|
|
def _marine_channel(freq_hz: float) -> str:
|
|
"""Best-effort marine VHF channel label."""
|
|
known = {156_800_000: " (ch 16 distress)", 156_650_000: " (ch 13 bridge)",
|
|
157_100_000: " (ch 22A USCG)", 156_450_000: " (ch 9)",
|
|
156_600_000: " (ch 12)", 156_700_000: " (ch 14)"}
|
|
for hz, name in known.items():
|
|
if abs(freq_hz - hz) < 6_000:
|
|
return name
|
|
return ""
|
|
|
|
|
|
# ---------------------------------------------------------------------------
|
|
# Trunked-radio control channels
|
|
# ---------------------------------------------------------------------------
|
|
|
|
# A trunked system keeps one channel permanently transmitting a data stream
|
|
# that tells the radios in the fleet which channel each conversation has been
|
|
# assigned. There is no speech on it and it never stops, so a scanner that
|
|
# treats it as a signal parks on it for as long as the record limit allows and
|
|
# then finds it again on the next sweep. It is the single most common way for
|
|
# an unattended scan to fill a disk with nothing.
|
|
#
|
|
# Each entry is (baud, levels, name, confidence, ambiguous). ``ambiguous``
|
|
# marks a shape that digital *voice* also uses, which must therefore be held
|
|
# to a much longer run of unbroken carrier before it is called a control
|
|
# channel -- a P25 talkgroup and a P25 control channel look alike for the
|
|
# first few seconds, and only the control channel is still there a minute
|
|
# later.
|
|
CONTROL_SIGNATURES: tuple[tuple[float, int, str, float, bool], ...] = (
|
|
(3600.0, 2, "Motorola SMARTNET / SmartZone (Type I/II)", 0.90, False),
|
|
(9600.0, 2, "EDACS / ProVoice", 0.80, False),
|
|
(1200.0, 2, "MPT-1327", 0.62, True),
|
|
(4800.0, 4, "P25 or DMR Tier III", 0.75, True),
|
|
(2400.0, 4, "NXDN / NEXEDGE", 0.65, True),
|
|
)
|
|
|
|
# Allocations where trunked systems live in the US. Being inside one is
|
|
# corroboration, never a requirement: trunking turns up on licensed business
|
|
# pairs all over the spectrum, and the shape of the signal is the real
|
|
# evidence.
|
|
_TRUNKED_BANDS = (
|
|
(136.0, 174.0), (380.0, 400.0), (406.0, 420.0), (450.0, 470.0),
|
|
(470.0, 512.0), (758.0, 775.0), (788.0, 805.0), (806.0, 824.0),
|
|
(851.0, 869.0), (896.0, 902.0), (935.0, 941.0),
|
|
)
|
|
|
|
|
|
@dataclass
|
|
class ControlChannel:
|
|
"""A trunking control channel, and how sure we are of it."""
|
|
|
|
system: str
|
|
confidence: float
|
|
reason: str
|
|
baud: float = 0.0
|
|
ambiguous: bool = False
|
|
|
|
def describe(self) -> str:
|
|
return f"{self.system} trunking control channel"
|
|
|
|
|
|
def _continuous_data(f: SignalFeatures) -> bool:
|
|
"""True when the signal is an unbroken, constant-envelope data stream.
|
|
|
|
Every control channel looks like this. Bursty data -- a pager page, an
|
|
ACARS message, a packet frame -- does not, and neither does voice, which
|
|
swings the envelope around as the speaker pauses.
|
|
"""
|
|
return (f.env_cv < 0.30 and f.ook_contrast_db < 12.0
|
|
and f.duty_cycle > 0.75 and f.freq_modes >= 2
|
|
and f.mode_spacing > 500.0 and f.level_dwell > 0.30)
|
|
|
|
|
|
def trunk_control(f: SignalFeatures, freq_hz: float = 0.0,
|
|
continuous_for: float = 0.0,
|
|
min_seconds: float = 20.0) -> ControlChannel | None:
|
|
"""Identify a trunking control channel from its shape and symbol rate.
|
|
|
|
``continuous_for`` is how many seconds of unbroken carrier this capture
|
|
has seen. Signatures shared with digital voice are only believed once
|
|
that passes ``min_seconds``; the unambiguous symbol rates -- nothing but a
|
|
control channel sends 3600 baud two-level FSK without pause -- are
|
|
believed straight away.
|
|
"""
|
|
if f.baud <= 0 or f.baud_stability < 0.60 or f.baud_strength < 12.0:
|
|
return None
|
|
if not _continuous_data(f):
|
|
return None
|
|
|
|
levels = 4 if f.freq_modes >= 4 else 2
|
|
for baud, want_levels, name, conf, ambiguous in CONTROL_SIGNATURES:
|
|
if want_levels != levels or not _baud_near(f, baud, 0.06):
|
|
continue
|
|
if ambiguous and continuous_for < min_seconds:
|
|
return None
|
|
reason = (f"unbroken {baud:.0f} baud {want_levels}-level data with no "
|
|
f"speech and no gaps -- a trunking control channel")
|
|
if ambiguous:
|
|
reason += f", still transmitting after {continuous_for:.0f} s"
|
|
mhz = freq_hz / 1e6
|
|
if any(lo <= mhz <= hi for lo, hi in _TRUNKED_BANDS):
|
|
conf = min(0.95, conf + 0.05)
|
|
reason += " in a trunked allocation"
|
|
return ControlChannel(system=name, confidence=conf, reason=reason,
|
|
baud=f.baud, ambiguous=ambiguous)
|
|
return None
|
|
|
|
|
|
# ---------------------------------------------------------------------------
|
|
# Entry point
|
|
# ---------------------------------------------------------------------------
|
|
|
|
def classify(x: np.ndarray, sample_rate: float, freq_hz: float = 0.0,
|
|
snr_db: float = 0.0, analyse_audio: bool = True,
|
|
continuous_for: float = 0.0,
|
|
control_seconds: float = 20.0) -> Classification:
|
|
"""Identify what kind of signal ``x`` is.
|
|
|
|
``x`` should be complex baseband centred on the signal. ``freq_hz`` is the
|
|
real-world centre frequency and is used only for the known-system lookup.
|
|
``continuous_for`` is how long the carrier has been up without a break,
|
|
which is what separates a trunking control channel from the digital voice
|
|
call it otherwise resembles.
|
|
"""
|
|
f = extract_features(x, sample_rate, snr_db=snr_db)
|
|
filtered = f.extras.pop("_filtered", x)
|
|
if analyse_audio:
|
|
try:
|
|
_analyse_fm_audio(filtered, sample_rate, f)
|
|
except Exception:
|
|
pass # audio features are a bonus, never a hard failure
|
|
|
|
cands = _score_rules(f)
|
|
cands.sort(key=lambda c: c[2], reverse=True)
|
|
label, family, score, reason, mode = cands[0]
|
|
reasons = [reason]
|
|
|
|
if family == "ssb":
|
|
# Region 2 convention: LSB on 160/80/40 m, USB everywhere else.
|
|
lower = freq_hz > 0 and freq_hz < 10_000_000
|
|
mode = "lsb" if lower else "usb"
|
|
label = f"SSB voice ({mode.upper()})"
|
|
reasons.append(f"{mode.upper()} assumed from the band convention")
|
|
|
|
control = trunk_control(f, freq_hz, continuous_for=continuous_for,
|
|
min_seconds=control_seconds)
|
|
|
|
system = _identify_system(freq_hz, f, family)
|
|
if system:
|
|
sys_label, sys_conf, sys_reason = system
|
|
if sys_conf >= score:
|
|
label = sys_label
|
|
score = min(0.97, 0.5 * sys_conf + 0.5 * score + 0.15)
|
|
reasons.insert(0, sys_reason)
|
|
else:
|
|
reasons.append(f"also consistent with {sys_label}")
|
|
|
|
# A control channel outranks the generic name for its modulation: "2-FSK"
|
|
# is true but useless, and the whole point of spotting one is to say so.
|
|
if control is not None:
|
|
label = control.describe()
|
|
score = max(score, control.confidence)
|
|
reasons.insert(0, control.reason)
|
|
|
|
# Low SNR means low trust, whatever the rules said.
|
|
if f.snr_db < 8:
|
|
score *= 0.65
|
|
reasons.append(f"low SNR ({f.snr_db:.0f} dB) -- treat with caution")
|
|
elif f.snr_db < 15:
|
|
score *= 0.85
|
|
|
|
alts = [(c[0], round(min(0.99, c[2]), 2)) for c in cands[1:4]]
|
|
return Classification(
|
|
label=label, family=family, confidence=round(min(0.99, score), 3),
|
|
reasons=reasons, alternatives=alts, suggested_mode=mode, features=f,
|
|
control=control,
|
|
)
|
|
|
|
|
|
# --------------------------------------------------------------------------
|
|
# SSB alignment
|
|
# --------------------------------------------------------------------------
|
|
|
|
# How far below the first audio energy the suppressed carrier sits. Transmit
|
|
# filters start passing somewhere around 200-300 Hz, and our own SSB filter
|
|
# opens at 200 Hz, so guessing a little low costs nothing while guessing high
|
|
# clips the bottom of the voice.
|
|
CARRIER_GUARD_HZ = 250.0
|
|
|
|
|
|
@dataclass
|
|
class SSBAlignment:
|
|
"""Where an SSB signal's suppressed carrier is, and which sideband it is."""
|
|
|
|
lower_hz: float # low edge of the occupied band, about the centre
|
|
upper_hz: float # high edge
|
|
sideband: str # "usb" or "lsb", from which way the band leans
|
|
confidence: float # 0 = perfectly symmetric, 1 = all to one side
|
|
width_hz: float
|
|
|
|
def carrier_for(self, sideband: str) -> float:
|
|
"""Carrier offset in Hz, given a sideband -- measured or assumed."""
|
|
if sideband == "lsb":
|
|
return self.upper_hz + CARRIER_GUARD_HZ
|
|
return self.lower_hz - CARRIER_GUARD_HZ
|
|
|
|
|
|
def ssb_alignment(x: np.ndarray, sample_rate: float,
|
|
hint_bw: float = 0.0) -> SSBAlignment | None:
|
|
"""Find the suppressed carrier of an SSB signal centred near DC.
|
|
|
|
Every other mode tolerates a kilohertz of tuning error -- an FM
|
|
discriminator and an AM envelope detector do not care where in their
|
|
passband the signal sits. SSB does not: its demodulator is a filter that
|
|
opens at the carrier, so tuning to the middle of the voice throws away
|
|
everything below that point and shifts what is left down by the error.
|
|
|
|
Speech puts most of its power in the first few hundred hertz above the
|
|
carrier, so an SSB signal's occupied band is lopsided, and which way it
|
|
leans identifies the sideband without recourse to any convention: energy
|
|
bunched at the low edge is upper sideband, at the high edge lower.
|
|
Returns None when the band is too odd to read, leaving the caller to fall
|
|
back on what the band plan says.
|
|
"""
|
|
nfft = min(32768, _pow2_floor(int(sample_rate / 100.0)))
|
|
if x.size < nfft or nfft < 256:
|
|
return None
|
|
_, psd = welch_psd(x, nfft)
|
|
psd_db = db(psd)
|
|
freqs = np.fft.fftshift(np.fft.fftfreq(nfft, 1.0 / sample_rate))
|
|
win = np.abs(freqs) <= max(hint_bw * 2.5, 8_000.0)
|
|
if not np.any(win):
|
|
return None
|
|
floor = float(np.percentile(psd_db[win], 25.0))
|
|
peak = float(psd_db[win].max())
|
|
hot = win & (psd_db > max(floor + 8.0, peak - 20.0))
|
|
idx = np.flatnonzero(hot)
|
|
if idx.size < 8:
|
|
return None
|
|
|
|
f_hot = freqs[idx]
|
|
# Percentiles rather than the extremes: one stray bin in the skirt would
|
|
# otherwise put the carrier a kilohertz out, and the carrier is the whole
|
|
# point of the measurement.
|
|
lo = float(np.percentile(f_hot, 2.0))
|
|
hi = float(np.percentile(f_hot, 98.0))
|
|
width = hi - lo
|
|
if not (400.0 <= width <= 8_000.0):
|
|
return None
|
|
|
|
# Where the power actually sits inside that band, by weight rather than by
|
|
# count, so a wide quiet skirt cannot outvote the loud part of the voice.
|
|
w = np.maximum(psd[idx], 0.0)
|
|
total = float(w.sum())
|
|
if total <= 0:
|
|
return None
|
|
order = np.argsort(f_hot)
|
|
f_sorted, w_sorted = f_hot[order], w[order]
|
|
cumulative = np.cumsum(w_sorted) / total
|
|
median = float(f_sorted[int(np.searchsorted(cumulative, 0.5))])
|
|
median = min(max(median, lo), hi)
|
|
|
|
to_lo, to_hi = median - lo, hi - median
|
|
upper = to_lo < to_hi
|
|
confidence = abs(to_hi - to_lo) / width if width > 0 else 0.0
|
|
return SSBAlignment(lower_hz=lo, upper_hz=hi,
|
|
sideband="usb" if upper else "lsb",
|
|
confidence=float(confidence), width_hz=float(width))
|