4a272eb1d5
A trunked system keeps one frequency transmitting a data stream around the clock so its radios know where each conversation has been put. There is no speech on it and it never stops, which makes it the strongest and most useless signal in the band: the scanner parked on 856.561 MHz for the full record limit, saved four minutes of buzzing, and found it again on the next sweep. Five signatures, matched against a constant-envelope stream that never pauses: 3600 baud two-level (Motorola SMARTNET/SmartZone), 9600 (EDACS), 1200 (MPT-1327), 4800 four-level (P25 or DMR Tier III), 2400 (NXDN). The first two are believed at once -- nothing else sends at those rates without pausing. The rest share their shape with a digital voice call on the same system, so they wait for the carrier to run unbroken past --control-seconds, longer than a conversation goes without a breath. Being in a trunked allocation raises confidence but is never required; trunking is licensed on business pairs all over the spectrum. One is named on screen, abandoned within a second or so, and its capture deleted. --keep-control records them for a decoder; --lockout-control writes them into the lock-out list. Three things had to be fixed to get there. The simulator's "pseudo-random" symbols were a counter: multiplying the symbol index by an odd constant and taking it modulo the level count returns the low bits, so two-level FSK came out 0,1,0,1. Every FSK test in the suite was measuring a tone. Its FSK is now shaped the way GFSK and C4FM shape a stream, too, square-edged keying being a signal no licensed transmitter would radiate. The symbol-rate estimator locked onto harmonics -- 3600 baud read as 18000 -- because a transition impulse train is a comb of equal lines; it now walks down to the fundamental. The squared envelope is no longer a candidate: it is not a transition signal, and its DC lobe made every random OOK signal measure ninety baud. The search starts at 200 Hz rather than 40, below which it was reading drift, which is how a bare carrier was awarded a symbol rate. And a clean two-level signal counted zero discriminator levels, because its modes land in the first and last histogram bin, where find_peaks cannot see them. Separately: locking out a frequency wrote to the settings file even under --no-config, which has no settings file by definition. It now writes only where it read from, and --simulate never writes at all -- an invented frequency would sit in a real config for ever, skipping whatever genuine signal happened to land near it. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_016PsWPTweCT6pwxKngvVxcg
1232 lines
52 KiB
Python
Executable file
1232 lines
52 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
|
|
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))
|