homelab-radio/sdr/discover/discover.py

#!/usr/bin/env python3
"""Signal-source discovery for the SDR scanner (HackRF, hackrf_sweep).

Run by sdr-discover.service (root) every 30 min while sdr-scanner is stopped (~70 s):
  1. hackrf_sweep over voice bands (same gains as the scanner), ~SWEEP_S seconds;
  2. per 5 kHz bin: local noise floor (rolling median over +-250 kHz), excess per sweep;
  3. clusters of bins >= HIT_DB above floor -> sources (peak freq snapped to channel grid,
     bandwidth estimate, activity = fraction of sweeps with the signal, level variance);
  4. merged into /srv/sdr/discovery/sources.json (first/last seen, runs, hits ...).
Persistent carriers that never change level (HackRF clone spurs, beacons, data links) are
learned and flagged kind="carrier" so the UI can hide them.
"""
import json
import os
import subprocess
import sys
import time

import numpy as np

# `python3 -I` implies -P (no script dir on sys.path): put the package root (the parent of discover/) there ourselves.
sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.abspath(__file__))))

from radio.config import CHANNELS_FILE, DISC_DIR, SOURCES_FILE  # noqa: E402
from radio.util import read_json  # noqa: E402

RUNS = os.path.join(DISC_DIR, "runs.jsonl")
SWEEP_S = int(os.environ.get("SWEEP_S", "60"))
BIN_HZ = 5000
HIT_DB = 10.0          # excess (block-averaged) over local floor counted as "signal present"
EDGE_DB = 6.0          # bandwidth edge
GAIN = ["-l", "32", "-g", "24", "-a", "0"]
# hackrf_sweep wants whole MHz; max 10 ranges
RANGES = [(118, 138), (144, 147), (148, 175), (420, 450)]


MAX_SWEEPS = 400


def sweep():
    """Run hackrf_sweep for SWEEP_S s; stream-parse into freqs (Hz) and M (n_sweeps x n_bins, dB)."""
    args = ["hackrf_sweep", "-w", str(BIN_HZ)] + GAIN
    for a, b in RANGES:
        args += ["-f", f"{a}:{b}"]
    p = subprocess.Popen(args, stdout=subprocess.PIPE, stderr=subprocess.DEVNULL, text=True)
    t_end = time.time() + SWEEP_S
    sweeps, curr = [], {}
    layout = None          # lo -> (bin width, nbins), fixed after the first full sweep
    try:
        for line in p.stdout:
            if time.time() > t_end:
                break
            parts = line.split(",")
            if len(parts) < 8:
                continue
            try:
                lo = int(parts[2])
                bw = float(parts[4])
                vals = np.array(parts[6:], dtype=np.float32)
            except ValueError:
                continue
            if lo in curr:               # chunk repeats -> a new sweep begins
                if layout is None:
                    layout = {k: (bw, len(v)) for k, (bw, v) in curr.items()}
                sweeps.append({k: v for k, (_, v) in curr.items()})
                curr = {}
                if len(sweeps) >= MAX_SWEEPS:
                    break
            curr[lo] = (bw, vals)
    finally:
        p.terminate()
        try:
            p.wait(timeout=10)
        except subprocess.TimeoutExpired:
            p.kill()
            p.wait()
    if layout is None or len(sweeps) < 2:
        return None, None
    los = sorted(layout)
    freqs, offs, o = [], {}, 0
    for lo in los:
        bw, nb = layout[lo]
        offs[lo] = o
        freqs.extend(int(lo + bw * (i + 0.5)) for i in range(nb))
        o += nb
    M = np.full((len(sweeps), o), np.nan, dtype=np.float32)
    for si, sw in enumerate(sweeps):
        for lo, v in sw.items():
            if lo in offs and len(v) == layout[lo][1]:
                M[si, offs[lo]:offs[lo] + len(v)] = v
    order = np.argsort(freqs)
    return np.array(freqs, dtype=np.int64)[order], M[:, order]


def floor_of(freqs, med):
    """Rolling median of the per-bin median over +-250 kHz (computed per contiguous segment)."""
    fl = np.empty_like(med)
    half = int(250000 / BIN_HZ)
    step = 10
    centers = list(range(0, len(med), step))
    vals = [np.nanpercentile(med[max(0, c - half):min(len(med), c + half + 1)], 40) for c in centers]
    fl[:] = np.interp(np.arange(len(med)), centers, vals)
    return fl


def grid(f):
    """Snap to the most plausible channel raster for the band."""
    mhz = f / 1e6
    if 118 <= mhz < 137:
        step, off = 25000, 0
    elif 446.0 <= mhz < 446.2:
        step, off = 12500, 6250
    elif 433.05 <= mhz < 434.8:
        step, off = 25000, 0
    else:
        step, off = 12500, 0
    return int(round((f - off) / step) * step + off)


def load_channels():
    try:
        with open(CHANNELS_FILE, encoding="utf-8") as fh:
            return {int(k): v for k, v in json.load(fh).get("channels", {}).items()}
    except (OSError, ValueError):
        return {}


# Community-sourced candidates (2026-10-20 research) that are NOT scanned: discovery labels them
# when they show up, so the user can promote them with one click.
CANDIDATES = {
    122.475: ("Левашово ATIS (непрерывно)", "am"),
    127.400: ("Пулково ATIS RU (непрерывно)", "am"),
    127.300: ("Пулково ATIS EN (непрерывно)", "am"),
    128.000: ("кандидат: Пулково резерв 128.0", "am"),
    129.000: ("кандидат: Пулково резерв 129.0", "am"),
    131.800: ("кандидат: Пулково Транзит 131.8", "am"),
    129.800: ("кандидат: Подход (стар.) 129.8", "am"),
    123.850: ("кандидат: воен. канал 3 123.85", "am"),
    135.000: ("кандидат: воен. канал 7 135.0", "am"),
    131.500: ("кандидат: Ржевка верт. 131.5", "am"),
    123.700: ("кандидат: Гостилицы «Бутон» 123.7", "am"),
    125.500: ("кандидат: «Саксофон» 125.5", "am"),
    122.125: ("кандидат: Таллин 122.125", "am"),
    131.725: ("ACARS (данные)", "am"),
    131.825: ("ACARS (данные)", "am"),
    438.750: ("UB1AAM DMR (цифра)", "nfm"),
    438.875: ("C4FM узел (цифра)", "nfm"),
    438.900: ("C4FM узел (цифра)", "nfm"),
    145.7375: ("RD1AS DMR (цифра)", "nfm"),
    144.800: ("APRS 144.800 (данные)", "nfm"),
}


def guess(f, bw, chans):
    mhz = f / 1e6
    for k, v in chans.items():
        if abs(k - f) <= 3000:
            return v["label"], v.get("mod", "nfm"), True
    for k, (lab, mod) in CANDIDATES.items():
        if abs(k * 1e6 - f) <= 3000:
            return lab, mod, False
    if 118 <= mhz < 137:
        return "Авиа (неизв.)", "am", False
    if 156.0 <= mhz <= 157.45:
        for n in list(range(1, 29)) + list(range(60, 89)):
            ship = 156.000 + 0.05 * n if n <= 28 else 156.025 + 0.05 * (n - 60)
            if abs(ship - mhz) < 0.004:
                return f"Мор. {n}", "nfm", False
        return "Морской диапазон", "nfm", False
    if 160.6 <= mhz <= 162.05:
        for n in list(range(1, 29)) + list(range(60, 89)):
            ship = 156.000 + 0.05 * n if n <= 28 else 156.025 + 0.05 * (n - 60)
            if abs(ship + 4.6 - mhz) < 0.004:
                return f"Мор. {n} (берег)", "nfm", False
        return "Морской (берег)", "nfm", False
    if 144 <= mhz < 146:
        return "Любит. 2 м", "nfm", False
    if 151.7 <= mhz < 152.1:
        return "ЖД (поездная?)", "nfm", False
    if 446.0 <= mhz < 446.2:
        return f"PMR {int(round((mhz - 446.00625) / 0.0125)) + 1}", "nfm", False
    if 433.05 <= mhz < 434.8:
        n = int(round((mhz - 433.075) / 0.025)) + 1
        return (f"LPD {n}" if 1 <= n <= 69 else "ISM 433"), "nfm", False
    if 430 <= mhz < 440:
        return "Любит. 70 см", "nfm", False
    if 146 <= mhz < 175:
        return "Служебная УКВ", "nfm", False
    if 420 <= mhz < 450:
        return "Служебная UHF", "nfm", False
    return "Неизвестно", "nfm", False


BLOCK = 8          # sweeps averaged per block (~0.7 s): per-bin noise std 5.6 dB -> ~1.6 dB
MIN_BLOCKS = 2     # signal must be present in >= 2 blocks (kills single-sweep noise spikes)
MERGE_HZ = 15000   # clusters closer than this are one source


def detect(freqs, M):
    inband = np.zeros(len(freqs), bool)
    for a, b in RANGES:
        inband |= (freqs >= a * 1e6) & (freqs < b * 1e6)
    freqs, M = freqs[inband], M[:, inband]
    L = np.power(10.0, M / 10.0)
    nb = L.shape[0] // BLOCK
    T = 10 * np.log10(np.nanmean(L[:nb * BLOCK].reshape(nb, BLOCK, -1), axis=1))
    med = np.nanmedian(T, axis=0)
    fl = floor_of(freqs, med)
    X = T - fl
    hit = X >= HIT_DB
    ok = hit.sum(axis=0) >= MIN_BLOCKS
    mx = np.nanmax(X, axis=0)
    out = []
    idx = np.where(ok)[0]
    groups, g = [], []
    for i in idx:
        if g and freqs[i] - freqs[g[-1]] > MERGE_HZ:
            groups.append(g)
            g = []
        g.append(i)
    if g:
        groups.append(g)
    for g in groups:
        seg = np.array(g)
        pk = int(seg[np.nanargmax(mx[seg])])
        if len(seg) < 2 and mx[pk] < 13:
            continue          # lone weak bin: almost always noise
        # bandwidth from the block where the peak occurred (max-over-time is too noisy)
        row = X[int(np.nanargmax(X[:, pk]))]
        edge = max(EDGE_DB, row[pk] - 20)
        a = b = pk
        while a - 1 >= 0 and row[a - 1] >= edge and freqs[pk] - freqs[a - 1] < 200000:
            a -= 1
        while b + 1 < len(freqs) and row[b + 1] >= edge and freqs[b + 1] - freqs[pk] < 200000:
            b += 1
        col = hit[:, seg].any(axis=1)
        lvl = T[:, pk]
        on = lvl[col] if col.any() else lvl
        out.append({
            "f": grid(int(freqs[pk])),
            "peak_hz": int(freqs[pk]),
            "bw_khz": round((b - a + 1) * BIN_HZ / 1000, 1),
            "act": round(float(col.mean()), 3),
            "max_db": round(float(mx[pk]), 1),
            "lvl_dbfs": round(float(np.nanmax(lvl)), 1),
            "std_db": round(float(np.nanstd(on)), 2),
            "hits": int(col.sum()),
        })
    return out, float(np.nanmedian(fl)), nb


def main():
    os.makedirs(DISC_DIR, exist_ok=True)
    t0 = time.time()
    freqs, M = sweep()
    if freqs is None or M.shape[0] < 5:
        print(f"too few sweeps ({0 if M is None else M.shape[0]}), skipping", flush=True)
        return 1
    if os.environ.get("DUMP"):
        np.savez_compressed(os.environ["DUMP"], freqs=freqs, M=M)
        return 0
    found, floor_db, nblocks = detect(freqs, M)
    chans = load_channels()
    db = read_json(SOURCES_FILE, {"sources": {}, "runs": 0})
    db["runs"] = db.get("runs", 0) + 1
    now = int(time.time())
    src = db.setdefault("sources", {})
    seen = set()
    for s in found:
        k = str(s["f"])
        if k in seen:     # two clusters snapped to the same channel: keep the stronger
            continue
        seen.add(k)
        e = src.get(k) or {"f": s["f"], "first": now, "runs": 0, "hits": 0, "sweeps": 0,
                           "max_db": 0, "act_sum": 0.0, "std_sum": 0.0}
        e["last"] = now
        e["runs"] += 1
        e["hits"] += s["hits"]
        e["sweeps"] += nblocks
        e["max_db"] = max(e["max_db"], s["max_db"])
        e["act_sum"] += s["act"]
        e["std_sum"] += s["std_db"]
        e["bw_khz"] = s["bw_khz"]
        e["last_db"] = s["max_db"]
        e["lvl_dbfs"] = s["lvl_dbfs"]
        lab, mod, known = guess(s["f"], s["bw_khz"], chans)
        e["label"], e["mod"], e["in_scan"] = lab, mod, known
        if s["bw_khz"] > 40:
            e["mod"] = "wide"
        src[k] = e
    # classify: persistent carrier = present in ~every sweep of most runs it was seen in,
    # with a steady level -> spur / beacon / data link, never voice. Intermittent -> voice candidate.
    for k, e in src.items():
        r = max(1, e["runs"])
        act = e["act_sum"] / r
        stdv = e["std_sum"] / r
        e["act"] = round(act, 3)
        e["std"] = round(stdv, 2)
        e["seen_frac"] = round(e["runs"] / db["runs"], 3)
        if act >= 0.9 and stdv < 2.0:
            e["kind"] = "carrier"      # steady unmodulated-looking carrier: spur / beacon / data
        elif act >= 0.9:
            e["kind"] = "continuous"   # always on but level varies: digital / broadcast / repeater
        elif e.get("mod") == "wide":
            e["kind"] = "wide"
        else:
            e["kind"] = "active"       # intermittent: voice candidate
    # forget one-off weak blips older than 7 days
    cutoff = now - 7 * 86400
    for k in [k for k, e in src.items() if e["last"] < cutoff and e["runs"] <= 1]:
        del src[k]
    db["last_run"] = now
    db["last_run_s"] = round(time.time() - t0, 1)
    db["last_sweeps"] = int(M.shape[0])
    db["floor_dbfs"] = round(floor_db, 1)
    tmp = SOURCES_FILE + ".tmp"
    with open(tmp, "w", encoding="utf-8") as fh:
        json.dump(db, fh, ensure_ascii=False, separators=(",", ":"))
    os.chmod(tmp, 0o644)
    os.replace(tmp, SOURCES_FILE)
    with open(RUNS, "a", encoding="utf-8") as fh:
        fh.write(json.dumps({"t": now, "sweeps": int(M.shape[0]), "found": len(found),
                             "active": [s for s in found if s["act"] < 0.9]}, ensure_ascii=False) + "\n")
    print(f"sweeps={M.shape[0]} found={len(found)} intermittent={sum(1 for s in found if s['act'] < 0.9)} "
          f"total_known={len(src)}", flush=True)
    return 0


if __name__ == "__main__":
    sys.exit(main())