2aa8739f25
An SSB capture was tuned to the centroid of the detected energy, which is what every other mode wants and the one thing SSB cannot use. Its demodulator is a filter that opens at the suppressed carrier, so centring on the middle of the voice filtered off its lower half and shifted the rest down by the error. Measured against known transmitters that error was +2445 Hz on 20 m and -2486 Hz on 80 m, against 35 Hz for NFM and 9 Hz for AM, which do not care either way. On real speech the difference is a clean transcript versus nothing recognisable at all. Speech puts most of its power just above the carrier, so the occupied band leans towards it. ssb_alignment() reads that lean: it locates the carrier to within about 100 Hz and names the sideband at the same time, without recourse to any convention. The offset is applied inside the demodulator at the IF rate, where it costs a fraction of what shifting the full-rate stream would, and the reported frequency becomes the carrier -- the frequency an operator would dial in. The sideband was also decided by "LSB below 10 MHz, USB above", which overrode a band plan that already knew better and demodulated 60 m as LSB. The measurement decides now, the band plan when the signal has no lean to read, and the convention only when neither has anything to say. The simulator was transmitting SSB unfiltered, several times wider than anything on the air, because it applied no transmit audio filter -- and that filter is what makes a signal single-sideband. Its absence hid the whole problem from the tests. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
1071 lines
45 KiB
Python
Executable file
1071 lines
45 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"]
|
|
|
|
|
|
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
|
|
|
|
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))])
|
|
|
|
|
|
def _cyclic_line(feature: np.ndarray, fs: float,
|
|
lo_hz: float = 40.0, 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
|
|
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)
|
|
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 ---------------------------------------------------
|
|
# For FSK the transition magnitude of the discriminator carries the line;
|
|
# for linear modulations the squared envelope does.
|
|
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))))
|
|
cand.append((_cyclic_line(env ** 2, sample_rate), env ** 2))
|
|
(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 ""
|
|
|
|
|
|
# ---------------------------------------------------------------------------
|
|
# Entry point
|
|
# ---------------------------------------------------------------------------
|
|
|
|
def classify(x: np.ndarray, sample_rate: float, freq_hz: float = 0.0,
|
|
snr_db: float = 0.0, analyse_audio: bool = True) -> 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.
|
|
"""
|
|
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")
|
|
|
|
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}")
|
|
|
|
# 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,
|
|
)
|
|
|
|
|
|
# --------------------------------------------------------------------------
|
|
# 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))
|