#!/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()