istrain
559caead36
Public community release. The complete working system — DSP + decoders (scripts/), web dashboard + API (dashboard/), container stack (docker/), config templates (config/) — plus APOCALYPSE-EDITION.md, the full build guide for a human or a coding agent. GPLv3 (PyEOT dependency). Scrubbed of all secrets, internal IPs, hostnames, and location detail. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
315 lines
15 KiB
Python
315 lines
15 KiB
Python
#!/usr/bin/env python3
|
||
# iq_channelize.py — split ONE RTL-SDR IQ dwell into several sub-channels by frequency offset, detect
|
||
# presence per sub-channel, and (on a burst) save a decode-ready clip. This is the IQ swap that lets dongle
|
||
# 1001 see mid-train DPU distinctly: DPU rides ±12.5 kHz off the HOT/EOT centers (452.925/.950, 457.925/.950),
|
||
# which FM-demod audio blurs into one unrecoverable channel. Raw IQ keeps the whole slice; we digitally
|
||
# down-convert each sub-channel and judge it on its own.
|
||
#
|
||
# 1001 stays a 2-position hop (452.9375 / 457.9375 — 5 MHz apart, wider than the tuner can grab at once).
|
||
# This script only changes what each dwell EXTRACTS:
|
||
# 452 dwell -> BOT (center) + MID (-12.5k, +12.5k)
|
||
# 457 dwell -> EOT (center) + MID (-12.5k, +12.5k)
|
||
# So MID is fed from BOTH dwells (a train uses one DPU channel, picked by its lead loco's road number).
|
||
#
|
||
# PRESENCE = in-channel carrier power above an adaptive noise floor (NOT FM-audio envelope) — the correct,
|
||
# squelch-polarity-proof presence metric, and what makes clean channel isolation testable. We STILL emit a
|
||
# discriminated s16 audio clip per burst so eot_scan.py / listening keep working unchanged.
|
||
#
|
||
# Reads CU8 IQ on stdin: rtl_sdr -d 1001 -f 452.9375M -s 240000 -g 40 - | iq_channelize.py --dwell 452
|
||
# RUNTIME helper — built by the maintainer, you run it. numpy only.
|
||
# No radio needed to prove the DSP: python3 iq_channelize.py --selftest
|
||
#
|
||
# Detection mirrors vumon.py's adaptive floor+margin (promote-candidate: share one detector module later).
|
||
import os, sys, time, json, math
|
||
import numpy as np
|
||
|
||
FS_IN = int(os.environ.get("IQ_FS", "240000")) # rtl_sdr sample rate (Hz)
|
||
DECIM = int(os.environ.get("IQ_DECIM", "5")) # 240k -> 48k channel rate
|
||
LP_CUT = float(os.environ.get("IQ_LPCUT", "6000")) # per-channel low-pass cutoff (Hz); rejects the ±12.5k neighbor
|
||
NTAPS = int(os.environ.get("IQ_TAPS", "151")) # FIR length; ~53 dB stopband -> clean sub-channel isolation
|
||
MARGIN = float(os.environ.get("IQ_MARGIN", "6")) # fire this many dB above the rolling floor
|
||
# Dial 2 (2026-07-10): only SAVE mid .s16 audio for bursts this many dB above the floor. The booster
|
||
# floor is characterized furniture; its audio adds nothing. 0 = save all (the pre-dial behavior).
|
||
MID_CLIP_MIN_OVER = float(os.environ.get("MID_CLIP_MIN_OVER", "6"))
|
||
FLOOR_A = 0.02 # noise-floor EMA rate (updated only while quiet)
|
||
CHUNK = int(os.environ.get("IQ_CHUNK", "24000")) # IQ samples per read (~0.1s @ 240k); multiple of DECIM
|
||
FS_AUD = FS_IN // DECIM
|
||
HOME = os.path.expanduser("~")
|
||
IST = os.path.join(HOME, "istrain")
|
||
LEVEL = os.path.join(IST, "level.json") # live meter the dashboard polls (vumon-compatible shape)
|
||
|
||
_last_level = 0.0
|
||
def write_level(now, db, active, label, peak, floor):
|
||
"""Throttled (~5 Hz) atomic write of the instantaneous level — same shape vumon writes, so the
|
||
dashboard meter works unchanged. In IQ mode we report the HOTTEST sub-channel (label = bot/eot/mid)."""
|
||
global _last_level
|
||
if now - _last_level < 0.2:
|
||
return
|
||
_last_level = now
|
||
try:
|
||
with open(LEVEL + ".tmp", "w") as f:
|
||
json.dump({"ts": round(now, 2), "db": round(db, 1), "active": bool(active),
|
||
"label": label, "peak": (round(peak, 1) if (active and peak is not None) else None),
|
||
"floor": (round(floor, 1) if floor is not None else None)}, f)
|
||
os.replace(LEVEL + ".tmp", LEVEL)
|
||
except OSError:
|
||
pass
|
||
|
||
|
||
def design_lpf(cut, fs, ntaps):
|
||
"""Windowed-sinc low-pass FIR (Hamming). Linear phase, ~53 dB stopband at 151 taps."""
|
||
n = np.arange(ntaps) - (ntaps - 1) / 2.0
|
||
fc = cut / fs # normalized cutoff (cycles/sample)
|
||
h = 2 * fc * np.sinc(2 * fc * n)
|
||
h *= np.hamming(ntaps)
|
||
return (h / h.sum()).astype(np.float64)
|
||
|
||
|
||
class SubChannel:
|
||
"""One frequency-offset sub-channel: NCO mix -> LPF -> decimate -> power-presence + discriminated clip."""
|
||
def __init__(self, offset_hz, label, jsonl, clips_dir, sub=None):
|
||
self.offset = float(offset_hz)
|
||
self.label = label
|
||
self.sub = sub # sub-channel freq tag (e.g. "457925") — de-collides + labels mid clips
|
||
self.jsonl = jsonl
|
||
self.clips = clips_dir
|
||
if clips_dir:
|
||
os.makedirs(clips_dir, exist_ok=True)
|
||
self.h = design_lpf(LP_CUT, FS_IN, NTAPS)
|
||
# NCO: a fixed per-chunk exponential * a carried start-phase scalar (unit magnitude, renormalized)
|
||
self._base = np.exp(-2j * math.pi * self.offset * np.arange(CHUNK) / FS_IN)
|
||
self._step = complex(np.exp(-2j * math.pi * self.offset * CHUNK / FS_IN))
|
||
self._ph = 1 + 0j
|
||
self._hist = np.zeros(NTAPS - 1, dtype=np.complex128) # filter state (overlap)
|
||
self._lastc = 0 + 0j # FM discriminator state (prev decimated sample)
|
||
# detector state
|
||
self.floor = None
|
||
self.active = False
|
||
self.since = 0.0
|
||
self.peak = -999.0
|
||
self.last_db = -999.0
|
||
self.last_on = False
|
||
self._clipbuf = [] # discriminated audio (int16) for the live burst
|
||
self._preroll = [] # ~0.5s pre-roll of audio chunks
|
||
|
||
def _mix_filter_decimate(self, iq):
|
||
"""iq: complex chunk @ FS_IN -> decimated complex @ FS_AUD for this sub-channel."""
|
||
mixed = iq * self._base * self._ph
|
||
self._ph *= self._step
|
||
if abs(self._ph) > 1e-6:
|
||
self._ph /= abs(self._ph) # kill slow magnitude drift
|
||
x = np.concatenate((self._hist, mixed))
|
||
y = np.convolve(x, self.h, "valid") # len == len(mixed)
|
||
self._hist = mixed[-(NTAPS - 1):]
|
||
return y[::DECIM]
|
||
|
||
def _discriminate(self, yc):
|
||
"""Polar FM discriminator (same core as rtl_fm -M fm) -> int16 audio for clips."""
|
||
prev = np.empty_like(yc)
|
||
prev[0] = self._lastc
|
||
prev[1:] = yc[:-1]
|
||
self._lastc = yc[-1]
|
||
ang = np.angle(yc * np.conj(prev)) # [-pi, pi]
|
||
return np.clip(ang / math.pi * 32767, -32768, 32767).astype(np.int16)
|
||
|
||
def process(self, iq, now):
|
||
yc = self._mix_filter_decimate(iq)
|
||
pw = float(np.mean(np.abs(yc) ** 2)) + 1e-12
|
||
db = 10.0 * math.log10(pw)
|
||
audio = self._discriminate(yc)
|
||
|
||
if self.floor is None:
|
||
self.floor = db
|
||
on = db > self.floor + MARGIN
|
||
if not on:
|
||
self.floor += FLOOR_A * (db - self.floor) # track floor ONLY while quiet
|
||
self.last_db = db; self.last_on = on # for the live meter (hottest-channel pick)
|
||
|
||
if on and not self.active:
|
||
self.active = True; self.since = now; self.peak = db
|
||
if self.clips:
|
||
self._clipbuf = list(self._preroll); self._clipbuf.append(audio)
|
||
elif on:
|
||
self.peak = max(self.peak, db)
|
||
if self.clips:
|
||
self._clipbuf.append(audio)
|
||
elif self.active:
|
||
self.active = False
|
||
self._emit(now)
|
||
elif self.clips:
|
||
self._preroll.append(audio)
|
||
if len(self._preroll) > 5: # ~0.5s pre-roll
|
||
self._preroll.pop(0)
|
||
return db, on
|
||
|
||
def _emit(self, now):
|
||
# height above the tracked noise/booster floor = how "train-shaped" this burst is (0 ≈ floor).
|
||
over = round(self.peak - self.floor, 1) if self.floor is not None else None
|
||
rec = {"ts": round(self.since, 2), "dur": round(now - self.since, 2),
|
||
"peak": round(self.peak, 1), "over": over}
|
||
if self.sub: # which ±12.5k sub-channel (457925 = the booster)
|
||
rec["sub"] = self.sub
|
||
if self.jsonl: # the durable ROW always writes — the record is cheap
|
||
try:
|
||
with open(self.jsonl, "a") as f:
|
||
f.write(json.dumps(rec) + "\n")
|
||
except OSError:
|
||
pass
|
||
# Dial 2: don't SAVE floor-level mid AUDIO (the booster). Only gates the mid label's .s16 write;
|
||
# eot/bot always save. Reversible via MID_CLIP_MIN_OVER=0. The row above still logged the burst.
|
||
if self.label == "mid" and over is not None and over < MID_CLIP_MIN_OVER:
|
||
self._clipbuf = []
|
||
return
|
||
if self.clips and self._clipbuf:
|
||
try:
|
||
# sub tag in the filename: de-collides the two mid sub-channels (they could share a
|
||
# start-second and overwrite each other) AND sorts clips by frequency. "mid457925_…".
|
||
prefix = self.label + (self.sub or "")
|
||
fn = os.path.join(self.clips, "%s_%s.s16" % (prefix,
|
||
time.strftime("%Y%m%d_%H%M%S", time.localtime(self.since))))
|
||
np.concatenate(self._clipbuf).tofile(fn)
|
||
except OSError:
|
||
pass
|
||
self._clipbuf = []
|
||
return rec
|
||
|
||
|
||
def dwell_channels(dwell):
|
||
"""Map a dwell (452|457) to its sub-channels. MID writes midtrain.jsonl from either dwell.
|
||
Each mid sub-channel is tagged with its absolute frequency in kHz (dwell center ±12.5), so the
|
||
fixed 457.925 booster carrier separates in the log + filenames from any real DPU frequency."""
|
||
if dwell == "452":
|
||
center = ("bot", os.path.join(IST, "bot.jsonl"), os.path.join(IST, "bot_clips"))
|
||
cen_khz = 452937.5 # AAR 452.9375 MHz (head-end / BOT)
|
||
elif dwell == "457":
|
||
center = ("eot", os.path.join(IST, "eot.jsonl"), os.path.join(IST, "eot_clips"))
|
||
cen_khz = 457937.5 # AAR 457.9375 MHz (EOT); −12.5 = 457.925 booster
|
||
else:
|
||
sys.exit("iq_channelize: --dwell must be 452 or 457")
|
||
midlog = os.path.join(IST, "midtrain.jsonl")
|
||
midclip = os.path.join(IST, "mid_clips")
|
||
sub = lambda off: str(int(cen_khz + off / 1000.0)) # e.g. 457937.5 + (-12.5) = 457925
|
||
return [
|
||
SubChannel(0, center[0], center[1], center[2]),
|
||
SubChannel(-12500, "mid", midlog, midclip, sub=sub(-12500)),
|
||
SubChannel(+12500, "mid", midlog, midclip, sub=sub(+12500)),
|
||
]
|
||
|
||
|
||
def run_stream(channels):
|
||
"""Read CU8 IQ from stdin, drive each sub-channel. CU8 = interleaved uint8 I,Q centered at 127.5."""
|
||
nbytes = CHUNK * 2
|
||
fd = sys.stdin.buffer
|
||
sys.stderr.write("[iq_channelize] %d ch, fs=%d decim=%d -> %d Hz, lpf=%.0f, %d taps | LEVEL=%s\n"
|
||
% (len(channels), FS_IN, DECIM, FS_AUD, LP_CUT, NTAPS, LEVEL)); sys.stderr.flush()
|
||
while True:
|
||
raw = fd.read(nbytes)
|
||
if not raw or len(raw) < nbytes:
|
||
break
|
||
u = np.frombuffer(raw, dtype=np.uint8).astype(np.float64)
|
||
iq = (u[0::2] - 127.5) + 1j * (u[1::2] - 127.5)
|
||
iq /= 127.5
|
||
now = time.time()
|
||
for ch in channels:
|
||
ch.process(iq, now)
|
||
hot = max(channels, key=lambda c: c.last_db) # meter shows the hottest sub-channel
|
||
write_level(now, hot.last_db, any(c.active for c in channels),
|
||
hot.label, (hot.peak if hot.active else None), hot.floor)
|
||
|
||
|
||
# ---- synthetic self-test: prove routing + isolation without a radio ----------------------------------
|
||
def _fsk_burst(t, f_off, on_lo, on_hi, amp=0.6):
|
||
"""A crude FFSK-ish burst at carrier offset f_off, active only in [on_lo, on_hi]."""
|
||
base = 2 * math.pi * f_off * t
|
||
sub = 1200.0 * np.sign(np.sin(2 * math.pi * 90 * t)) # toy ±1200 Hz audio FSK
|
||
phase = base + np.cumsum(2 * math.pi * sub / FS_IN)
|
||
mask = ((t >= on_lo) & (t < on_hi)).astype(np.float64)
|
||
return amp * mask * np.exp(1j * phase)
|
||
|
||
|
||
def selftest():
|
||
dur = 1.2
|
||
n = int(FS_IN * dur)
|
||
t = np.arange(n) / FS_IN
|
||
rng = np.random.default_rng(7)
|
||
noise = 0.05 * (rng.standard_normal(n) + 1j * rng.standard_normal(n))
|
||
sig = noise.copy()
|
||
sig += _fsk_burst(t, 0.0, 0.10, 0.25) # CENTER (bot/eot) -> 0.10..0.25
|
||
sig += _fsk_burst(t, +12500.0, 0.45, 0.60) # MID (+12.5k) -> 0.45..0.60
|
||
sig += _fsk_burst(t, -12500.0, 0.80, 0.95) # MID (-12.5k) -> 0.80..0.95
|
||
|
||
# build channels with NO file output (capture detections in memory)
|
||
chans = [SubChannel(0, "center", None, None),
|
||
SubChannel(-12500, "mid_lo", None, None),
|
||
SubChannel(+12500, "mid_hi", None, None)]
|
||
hits = {c.label: [] for c in chans}
|
||
orig = {c.label: c._emit for c in chans}
|
||
|
||
def make_cap(c):
|
||
def cap(now):
|
||
hits[c.label].append((round(c.since, 2), round(now - c.since, 2), round(c.peak, 1)))
|
||
return orig[c.label](now)
|
||
return cap
|
||
for c in chans:
|
||
c._emit = make_cap(c)
|
||
|
||
# convert absolute sample stream to CU8 and feed chunk by chunk (using a fake monotonic clock)
|
||
u = np.empty(2 * n, dtype=np.float64)
|
||
u[0::2] = np.real(sig) * 127.5 + 127.5
|
||
u[1::2] = np.imag(sig) * 127.5 + 127.5
|
||
cu8 = np.clip(u, 0, 255).astype(np.uint8).tobytes()
|
||
nbytes = CHUNK * 2
|
||
pos = 0; k = 0
|
||
while pos + nbytes <= len(cu8):
|
||
raw = cu8[pos:pos + nbytes]; pos += nbytes
|
||
arr = np.frombuffer(raw, dtype=np.uint8).astype(np.float64)
|
||
iq = ((arr[0::2] - 127.5) + 1j * (arr[1::2] - 127.5)) / 127.5
|
||
now = k * CHUNK / FS_IN # synthetic clock (seconds)
|
||
for c in chans:
|
||
c.process(iq, now)
|
||
k += 1
|
||
for c in chans: # flush any still-active burst
|
||
if c.active:
|
||
c._emit((k * CHUNK) / FS_IN)
|
||
|
||
def near(hitlist, lo, hi):
|
||
return any(lo - 0.15 <= h[0] <= hi + 0.05 for h in hitlist)
|
||
|
||
ok = True
|
||
def check(name, cond):
|
||
nonlocal ok
|
||
print((" PASS " if cond else " FAIL ") + name)
|
||
ok = ok and cond
|
||
|
||
print("self-test — synthetic bursts: center@0.10-0.25, mid+12.5k@0.45-0.60, mid-12.5k@0.80-0.95")
|
||
print("detections:")
|
||
for c in chans:
|
||
print(" %-7s %s" % (c.label, hits[c.label]))
|
||
check("center detects its burst", near(hits["center"], 0.10, 0.25))
|
||
check("mid_hi (+12.5k) detects its burst", near(hits["mid_hi"], 0.45, 0.60))
|
||
check("mid_lo (-12.5k) detects its burst", near(hits["mid_lo"], 0.80, 0.95))
|
||
check("center ISOLATED from the two mid bursts", not near(hits["center"], 0.45, 0.60) and not near(hits["center"], 0.80, 0.95))
|
||
check("mid_hi ISOLATED from center + mid_lo", not near(hits["mid_hi"], 0.10, 0.25) and not near(hits["mid_hi"], 0.80, 0.95))
|
||
check("mid_lo ISOLATED from center + mid_hi", not near(hits["mid_lo"], 0.10, 0.25) and not near(hits["mid_lo"], 0.45, 0.60))
|
||
print("RESULT:", "PASS — routing + isolation hold" if ok else "FAIL")
|
||
return 0 if ok else 1
|
||
|
||
|
||
def main():
|
||
a = sys.argv[1:]
|
||
if "--selftest" in a:
|
||
sys.exit(selftest())
|
||
dwell = None
|
||
while a:
|
||
t = a.pop(0)
|
||
if t == "--dwell": dwell = a.pop(0)
|
||
else: sys.exit("iq_channelize: unknown arg %r" % t)
|
||
if not dwell:
|
||
sys.exit("usage: rtl_sdr ... - | iq_channelize.py --dwell {452|457} (or --selftest)")
|
||
os.makedirs(IST, exist_ok=True)
|
||
run_stream(dwell_channels(dwell))
|
||
|
||
|
||
if __name__ == "__main__":
|
||
main()
|