istrain-public/scripts/eot/eot_scan.py

#!/usr/bin/env python3
# EOT scan/decode (numpy). FSK audio (rtl_fm .s16 @ 48k, or a wav) -> vectorized tone-energy discriminator
# -> 1200-baud clock recovery -> PyEOT frame-sync + parser + BCH check. Pure analysis — does not touch a radio.
#
# TWO passes (2026-07-05, the O-11 crack):
#   1. FAST/EXACT — the proven original: whole file, fixed tone pairs, exact sync. Catches on-frequency
#      EOTs (proven against PyEOT's demo3eot.wav, 3 packets, and weeks of live passages).
#   2. ADAPTIVE/FUZZY escalation — for the off-frequency emitters the fast path is deaf to. The 07-05
#      ground-truthed dwelling train (unit 76029, CUT at the crossing) transmitted ~600-800 Hz HIGH of
#      our bin: its audio tones sat ~(2000,2400), nowhere near the fixed pairs, in short ~100 ms bursts
#      drowned by clip noise. Fix: find the burst by FM-quieting (carrier = LOW audio rms), estimate its
#      tone hump by spectral centroid, build tone pairs around it, and allow <=2 bit errors in the sync
#      word (frames still gate on the PyEOT/BCH validity check, so a fuzzy sync never invents a packet).
#      Validated on 07-04/07-05 dwell clips: 0/357 -> 26/357 morning decodes; both dwell flavors read
#      (psi ~62 held = solid stopped train vs psi 0 held = cut/air-dumped). See OBSERVATIONS O-11.
#
# Packets carry dq: "exact" | "fuzzy" (decode-quality provenance — a lone fuzzy decode can be a BCH
# collision; corroboration [>=2 records or a head-addr match] is enforced at the card layer, not here).
# Consensus: the returned list is majority-vote ordered, so callers taking pk[0] get the agreed packet.
#   CLI:  eot_scan.py <file.wav|file.s16> ...
#   API:  scan_file(path) -> [ {unit,pressure,motion,marker,turbine,batt,charge,arm,dq}, ... ]  (empty = not EOT)
import sys, os, wave
from collections import Counter
import numpy as np
sys.path.insert(0, os.path.join(os.path.dirname(os.path.abspath(__file__)), "PyEOT"))
from eot_decoder import EOT_decode

SYNC = '10101011100010010'; BAUD = 1200
_INV = str.maketrans('01', '10')
_SYNC_A = np.array([int(c) for c in SYNC], np.int8)
TONE_PAIRS = [(1200, 1800), (1300, 1700), (1100, 1900)]
FUZZY_MAX_ERR = 2          # sync-word bit errors tolerated on the escalation path (BCH still gates)

def load(path):
    if path.lower().endswith('.wav'):
        w = wave.open(path, 'rb'); sw = w.getsampwidth(); fr = w.getframerate(); raw = w.readframes(w.getnframes()); w.close()
        if sw == 1: return (np.frombuffer(raw, np.uint8).astype(np.float64) - 128) / 128.0, fr
        return np.frombuffer(raw, np.int16).astype(np.float64) / 32768.0, fr
    return np.frombuffer(open(path, 'rb').read(), np.int16).astype(np.float64) / 32768.0, 48000

def discriminator(x, fr, flo, fhi, W):
    n = np.arange(len(x))
    def energy(f):
        c = np.concatenate(([0j], np.cumsum(x * np.exp(-2j * np.pi * f * n / fr))))
        return np.abs(c[W:] - c[:-W]) ** 2
    return energy(fhi) - energy(flo)

def bitstring(d, fr, thresh=None):
    t = np.median(d) if thresh is None else thresh
    sl = (d > t).astype(np.int8)
    spb = fr / BAUD; out = []; phase = spb / 2.0; prev = int(sl[0])
    for i in range(1, len(sl)):
        b = int(sl[i])
        if b != prev: phase = spb / 2.0; prev = b
        else:
            phase -= 1.0
            if phase <= 0: out.append(b); phase += spb
    return ''.join('1' if b else '0' for b in out)

def _decode_at(bb, i, seen, out, dq):
    for off in range(0, 13):
        seg = bb[i + off:i + off + 80]
        if len(seg) < 74: break
        try:
            e = EOT_decode(seg)
            if e.valid:
                key = (e.unit_addr, e.pressure, int(e.motion))
                if key not in seen:
                    seen.add(key)
                    out.append(dict(unit=e.unit_addr, pressure=e.pressure, motion=int(e.motion),
                                    marker=int(e.mkr_light), turbine=int(e.turbine),
                                    batt=e.batt_cond_text, charge=e.batt_charge, arm=e.arm_status, dq=dq))
        except Exception: pass

def _packets(bs, fuzzy=False, seen=None, out=None):
    seen = set() if seen is None else seen
    out = [] if out is None else out
    dq = "fuzzy" if fuzzy else "exact"
    for bb in (bs, bs.translate(_INV)):
        if fuzzy and len(bb) >= len(SYNC):
            a = np.array([int(c) for c in bb], np.int8)
            win = np.lib.stride_tricks.sliding_window_view(a, len(SYNC))
            for i in np.where((win != _SYNC_A).sum(axis=1) <= FUZZY_MAX_ERR)[0]:
                _decode_at(bb, int(i), seen, out, dq)
        else:
            i = bb.find(SYNC)
            while i != -1:
                _decode_at(bb, i, seen, out, dq)
                i = bb.find(SYNC, i + 1)
    return out

def _quiet_segs(x, fr, thresh):
    """Carrier segments found by FM-quieting: demodulated NOISE is loud, a keyed carrier is quiet.
    Returns [(a,b)] sample ranges (>=60 ms quiet runs, padded), so the decoder works the burst alone
    instead of drowning a ~100 ms frame in a 0.7 s clip's noise (the clock recovery killer)."""
    win = int(fr * 0.010)
    n = len(x) // win
    if n < 8: return []
    r = np.array([x[i*win:(i+1)*win].std() for i in range(n)])
    med = np.median(r)
    if med <= 0: return []
    quiet = r < med * thresh
    segs = []; cur = 0; start = 0
    for i, v in enumerate(quiet):
        if v:
            if cur == 0: start = i
            cur += 1
        else:
            if cur >= 6: segs.append((start, start + cur))
            cur = 0
    if cur >= 6: segs.append((start, start + cur))
    return [(max(0, (a - 2) * win), (b + 2) * win) for a, b in segs]

def _hump_center(seg, fr):
    """Spectral centroid of the burst's tone hump (600-3500 Hz) — where this emitter's FSK audio
    actually sits. An off-frequency carrier shifts the whole hump (O-11: +600-800 Hz)."""
    w = np.hanning(len(seg))
    F = np.abs(np.fft.rfft(seg * w)) ** 2
    f = np.fft.rfftfreq(len(seg), 1 / fr)
    m = (f > 600) & (f < 3500)
    P = F[m]
    return float((f[m] * P).sum() / P.sum()) if P.sum() > 0 else None

def _consensus(pk):
    """Majority-vote order so pk[0] is the agreed packet (fuzzy sync can add a stray BCH collision)."""
    if len(pk) <= 1: return pk
    votes = Counter((d["unit"], d["pressure"], d["motion"]) for d in pk)
    return sorted(pk, key=lambda d: -votes[(d["unit"], d["pressure"], d["motion"])])

def scan_file(path):
    x, fr = load(path); W = max(8, int(round(fr / BAUD)))
    # pass 1 — FAST/EXACT (the proven original, unchanged behavior)
    for flo, fhi in TONE_PAIRS:
        bs = bitstring(discriminator(x, fr, flo, fhi, W), fr)
        if SYNC in bs or SYNC in bs.translate(_INV):
            pk = _packets(bs)
            if pk: return _consensus(pk)
    # pass 2 — ADAPTIVE/FUZZY escalation (off-frequency / short-burst emitters)
    for thresh in (0.65, 0.8):
        for a, b in _quiet_segs(x, fr, thresh):
            seg = x[a:b]
            if len(seg) < int(fr * 0.05):        # <50 ms can't hold a frame
                continue
            c = _hump_center(seg, fr)
            cands = []
            if c:
                for cc in (c, c - 150, c + 150):
                    for half in (200, 250, 300, 350):
                        if cc - half > 300: cands.append((cc - half, cc + half))
            cands += TONE_PAIRS + [(2000, 2400)]
            for flo, fhi in cands:
                d = discriminator(seg, fr, flo, fhi, W)
                for thr in (None, d.mean()):     # median + mean slicers (skewed short bursts)
                    pk = _packets(bitstring(d, fr, thr), fuzzy=True)
                    if pk: return _consensus(pk)
    return []

def run(p):
    label = os.path.basename(p)
    try: pk = scan_file(p)
    except Exception as e: print(f"[{label}] err: {e}"); return 0
    for d in pk:
        print(f"  ✅ [{label}] unit={d['unit']} pressure={d['pressure']}psig motion={d['motion']} "
              f"marker={d['marker']} turbine={d['turbine']} batt={d['batt']}/{d['charge']} arm={d['arm']} ({d['dq']})")
    if not pk: print(f"[{label}] no valid EOT")
    return len(pk)

if __name__ == '__main__':
    t = sum(run(p) for p in sys.argv[1:]); print(f"\n=== {t} valid EOT packet(s) ===")