DsssBurstReceiver — the Composed Burst Chain¶

# The frame, as a DESCRIPTION: fields in wire order, and the stage that
# fills the CRC over the span it covers. Written once as JSON -- the file
# `wfmgen --frame` reads and the value a scene's "frame" key holds -- and
# once as the Python object, by the same names. Property 0 below is what
# proves the two are one description: they render byte-identical bursts.
FRAME = {
    "fields": [
        {"name": "sync", "spec": SYNC_SPEC},  # found, not decoded
        {"name": "payload", "spec": PAYLOAD_SPEC},
        {"name": "crc", "bits": 16, "derived_by": 1},  # stage 0 fills it
    ],
    "stages": [{"kind": "crc16", "first_field": 1, "n_fields": 2}],
}

desc = FrameDesc()
desc.add_field("sync", field_bits(SYNC_SPEC))
desc.add_field("payload", field_bits(PAYLOAD_SPEC))
desc.add_derived("crc", 16)
desc.add_stage_over(STAGE_CRC16, "payload", "crc")  # CRC-16 over both
desc.build()

FRAME_SYMS = desc.nbits  # the frame's length is the description's to say

# One dsss segment, played N_BURSTS times: the engine tiles the preamble,
# spreads the described frame by the data code, sizes the on-time itself,
# and draws fresh noise per burst. `gap_noise` stays at its default, so the
# floor runs through the gaps, as it does in a real capture.
#
# Python: the description is the FrameDesc object itself.
python_face = Composer(
    [
        Segment(
            type="dsss",
            fs=FS,
            snr=ESN0_DB,
            snr_mode="esno",  # Es/N0 of the DATA_SF-chip data symbol
            seed=SEED,
            sps=SPC,  # samples per CHIP
            acq_code=acq_code,
            acq_reps=REPS,
            data_code=data_code,
            frame=desc,
            repeats=N_BURSTS,
            off_samples=GAP,
        )
    ]
).compose()

# A scene: the same segment as JSON, the description under "frame", every
# code as its Field text. A scene has no `acq_reps`: the repetition rides
# in the preamble's Field, as `*REPS`.
scene = {
    "version": 1,
    "segments": [
        {
            "type": "dsss",
            "fs": FS,
            "snr": ESN0_DB,
            "snr_mode": "esno",
            "seed": SEED,
            "sps": SPC,
            "acq_code": f"{ACQ_SPEC}*{REPS}",
            "data_code": DATA_SPEC,
            "frame": FRAME,
            "repeats": N_BURSTS,
            "off_samples": GAP,
        }
    ],
}
scene_face = Composer.from_json(json.dumps(scene)).compose()

# The command line: the description is a file, and `--frame` reads it.
argv = [
    "wfmgen", "--type", "dsss", "--fs", f"{FS:.0f}",
    "--snr", f"{ESN0_DB:g}", "--snr-mode", "esno",
    "--seed", f"{SEED}", "--sps", f"{SPC}",
    "--acq-code", f"{ACQ_SPEC}*{REPS}", "--data-code", DATA_SPEC,
    "--frame", "spread.json",
    "--repeats", f"{N_BURSTS}", "--off", f"{GAP}",
    "-o", "train.cf32",
]  # fmt: skip
with tempfile.TemporaryDirectory() as tmp:
    Path(tmp, "spread.json").write_text(json.dumps(FRAME, indent=2))
    # `cli._runnable()` is the one locator for the bundled C binary.
    subprocess.run([cli._runnable(), *argv[1:]], cwd=tmp, check=True)
    cli_face = np.fromfile(Path(tmp, "train.cf32"), dtype=np.complex64)

assert np.array_equal(python_face, scene_face), "Python vs scene differ"
assert np.array_equal(cli_face, python_face), "CLI vs Python differ"
capture = python_face

def decode(block_size):
    """Push the capture in blocks of `block_size`; return valid starts+bits.

    A frame is kept only when its event says `frame_valid`: the receiver
    stops at decisions, so a false alarm in the noise (at the configured
    pfa) still returns a frame -- and that flag is how a caller tells it
    from a burst.
    """
    rx = receiver()
    starts, bits = [], []
    for off in range(0, capture.size, block_size):
        out = np.asarray(rx.push(capture[off : off + block_size]))
        for k, e in enumerate(rx.events()):
            if e["frame_valid"]:
                bits.append(out[k * FRAME_SYMS : (k + 1) * FRAME_SYMS])
                starts.append(int(e["preamble_start"]))
    return (
        starts,
        (np.concatenate(bits) if bits else np.empty(0, np.uint8)),
        rx,
    )

cut = truth[0] + BURST_LEN // 2
rx = receiver()
first = np.asarray(rx.push(capture[:cut]))
held = rx.pending
second = np.asarray(rx.push(capture[cut:]))
print(
    f"  fed {cut} samples (mid-burst): {first.size // FRAME_SYMS} frame(s), "
    f"pending={held}"
)
print(
    f"  fed the rest:                  {second.size // FRAME_SYMS} frame(s), "
    f"pending={rx.pending}"
)
assert first.size == 0, "a half-arrived burst must not be emitted"
assert held == 1, "pending must report the burst being held"
assert rx.pending == 0, "pending must clear once the burst is emitted"
assert np.array_equal(second[:FRAME_SYMS], desc.bits(1)), (
    "the split burst is not exact"
)
print("  -> held, then returned whole. Read pending before you stop feeding.")

# The dead air between bursts is READ from the receiver: `min_gap`, edge to
# edge, is the gap it guarantees. The demo places its bursts AT it, with no
# margin, because the guarantee is the thing being shown -- the assertions
# below are what would catch it failing. (`refine_span` is a start-to-start
# reach, not a gap; spacing by it is how this demo used to be wrong.)
probe = receiver()
MIN_GAP, RETAIN_SPAN = probe.min_gap, probe.retain_span
GAP = MIN_GAP
SPACING = BURST_LEN + GAP

# The scalar properties describe the LAST burst only — one push() can
# complete several and one set of scalars cannot speak for all of them.
# `events()` hands back the same fields per burst, and between them they
# are the whole diagnostic surface: where the burst was, what the
# search grid thought its Doppler was and how wide that grid is, what C/N0
# the hit implies, what the estimator refined the residual frequency and
# chirp rate to and how confident it was, how far the winning preamble beat
# its runner-up, and whether the CRC checked out.
_, _, rx_all = results[capture.size]
events = np.asarray(rx_all.events(rx_all.events_max_out()))

# What the scene says each read-back must be. Derived from the geometry
# above, not from a previous run of this script: `acq` transforms sf*spc
# verbatim, so the bin width is fs over that, and a segment generated at
# Es/N0 with a DATA_SF-chip symbol carries C/N0 = Es/N0 + 10log10(Rs).
BIN_HZ = FS / (ACQ_SF * SPC)
CN0_DBHZ_TRUE = ESN0_DB + 10.0 * np.log10(CHIP_RATE / DATA_SF)

print("\nevery read-back, per burst (events(), one row per decoded burst):")
print(
    f"  {'#':>2} {'start':>7} {'dopp_hz':>8} {'res_hz':>8} {'cn0_dBHz':>9} "
    f"{'freq_hz':>8} {'rate_hz':>8} {'conf_dB':>8}"
)
for i, e in enumerate(events):
    print(
        f"  {i:>2} {int(e['preamble_start']):>7} "
        f"{e['doppler_hz_est']:>8.1f} {e['doppler_res_hz']:>8.1f} "
        f"{e['cn0_dbhz_est']:>9.2f} {e['est_freq_hz']:>8.2f} "
        f"{e['est_rate_hz']:>8.2f} {e['demod_cn0_dbhz']:>8.2f}"
    )
print(
    f"  scene: bin {BIN_HZ:.1f} Hz, C/N0 {CN0_DBHZ_TRUE:.2f} dB-Hz "
    f"(Es/N0 {ESN0_DB:.0f} dB at {CHIP_RATE / DATA_SF / 1e3:.1f} ksym/s), "
    "true offset 0 Hz, no chirp"
)

assert events.size == N_BURSTS, "one event per decoded burst, always"
for i, e in enumerate(events):
    where = f"burst {i}"
    # The record must stand on its own: its start is the burst's start.
    assert int(e["preamble_start"]) == truth[i], (
        f"{where}: event says {int(e['preamble_start'])}, burst starts at "
        f"{truth[i]}"
    )
    # The search grid: its width is derived, and its estimate must sit in
    # the bin that contains the truth (0 Hz).
    assert e["doppler_res_hz"] == BIN_HZ, (
        f"{where}: bin width {e['doppler_res_hz']} != fs/(sf*spc) {BIN_HZ}"
    )
    assert abs(e["doppler_hz_est"]) <= BIN_HZ / 2, (
        f"{where}: coarse Doppler {e['doppler_hz_est']:.1f} Hz is outside "
        "the bin containing the true 0 Hz"
    )
    # ...and the whole reason the chain does not stop at acquisition: the
    # refined residual is finer than the grid by orders of magnitude.
    assert abs(e["est_freq_hz"]) < BIN_HZ / 100, (
        f"{where}: refined residual {e['est_freq_hz']:.2f} Hz is no better "
        f"than a hundredth of the {BIN_HZ:.0f} Hz search bin — refine and "
        "demod are not adding anything over the grid"
    )
    # A zero here is a CONFIGURATION fact, not a measurement: max_rate=0
    # switches the chirp axis off. Pass a non-zero max_rate to ask for one.
    assert e["est_rate_hz"] == 0.0, (
        f"{where}: chirp rate {e['est_rate_hz']} Hz/s from a receiver built "
        "with max_rate=0, which does not search that axis"
    )
    # The C/N0 estimate is a LOWER BOUND, so it may not run hot; the scene
    # says what it is bounding.
    assert e["cn0_dbhz_est"] <= CN0_DBHZ_TRUE + 1.5, (
        f"{where}: C/N0 estimate {e['cn0_dbhz_est']:.2f} dB-Hz runs hot "
        f"against the scene's {CN0_DBHZ_TRUE:.2f} dB-Hz — it is documented "
        "as a lower bound, and per-burst estimator spread here is a few "
        "tenths of a dB"
    )
    assert e["cn0_dbhz_est"] >= CN0_DBHZ_TRUE - 3.0, (
        f"{where}: C/N0 estimate {e['cn0_dbhz_est']:.2f} dB-Hz is more than "
        f"3 dB under the scene's {CN0_DBHZ_TRUE:.2f} dB-Hz"
    )
    # `demod_cn0_dbhz` IS the same quantity as the scene's C/N0, measured
    # on the decoded symbols instead of bounded from the acquisition hit, so
    # it is checked against the truth from BOTH sides. Its predecessor
    # `est_snr_db` could only be checked against a bare threshold: it was the
    # preamble estimator's peak-to-mean, carrying the coherent processing
    # gain, and comparing it with a link budget was a category error
    # (doppler#1304).
    assert abs(e["demod_cn0_dbhz"] - CN0_DBHZ_TRUE) < 3.0, (
        f"{where}: demod C/N0 {e['demod_cn0_dbhz']:.2f} dB-Hz is more than "
        f"3 dB off the scene's {CN0_DBHZ_TRUE:.2f} dB-Hz"
    )
    assert abs(e["demod_timing_chips"]) < 0.5, (
        f"{where}: start error {e['demod_timing_chips']:+.3f} chip sits at "
        "the half-chip search bound, so the estimate is a floor"
    )

# The scalars are not a second source: they ARE the last row.
last = events[-1]
for name in events.dtype.names:
    assert getattr(rx_all, name) == last[name], (
        f"scalar {name} disagrees with the last event row — the scalars "
        "describe the most recent burst and nothing else"
    )
print(
    f"  -> all {N_BURSTS} rows check out against the scene; the scalar "
    "read-backs equal row -1, which is all they ever claim to be"
)

# The receiver stopped at decisions: `push()` handed back FRAME BITS and no
# opinion about them (doppler#1022). The verdict and the payload need the
# frame's description -- and it is the SAME object the transmitter was
# given. `check()` reverses every checking stage and reports; `deframe()`
# hands back the corrected frame, where the payload is a slice found by name.
off = desc.field_off(desc.field_index("payload"))
deframed_ok, payloads = [], []
for k in range(N_BURSTS):
    frame = ref_bits[k * FRAME_SYMS : (k + 1) * FRAME_SYMS]
    verdict = desc.check(frame)
    deframed_ok.append(verdict.passed == verdict.checked == 1)
    got = np.asarray(desc.deframe(frame))
    payloads.append(got[off : off + PAYLOAD])

sent = field_bits(PAYLOAD_SPEC)
exact = sum(int(np.array_equal(p, sent)) for p in payloads)
print("\nchecked against the description (the receiver has no opinion):")
print(
    f"  {sum(deframed_ok)}/{N_BURSTS} frames check out, "
    f"{exact}/{N_BURSTS} payloads bit-exact"
)
assert all(deframed_ok), "every frame's CRC must check out"
for k, got in enumerate(payloads):
    assert np.array_equal(got, sent), f"burst {k} payload is not exact"

# ...and a check that cannot fail is not a check: one flipped payload bit
# turns the verdict over.
bad = ref_bits[:FRAME_SYMS].copy()
bad[off] ^= 1
assert desc.check(bad).passed == 0, "a corrupted frame must fail its CRC"
print("  -> one description: transmitted from, decided on, checked against.")

bits = rx.push(x)                    # hard, per burst, frame rows
llr = rx.llrs(rx.llrs_max_out(1))    # soft, same order, same rows

BASE_EBN0_DB = 6.0
FRAME = FrameDesc()  # [sync | data]: no CRC, so nothing follows the payload
FRAME.add_field("sync", SYNC)
FRAME.add_data("payload", len(PAYLOAD))
rx = receiver()  # ONE receiver for the whole sweep: reset() between trials
# One dsss source is the whole scene, so the Composer takes it directly.
# snr_mode="ebno" makes `snr` the Eb/N0 of one data bit (a dsss source's
# symbol is the despread data bit), and Plan's at() takes the SAME number:
# no hand conversion to the over-fs SNR anywhere below.
scene = Composer(
    type="dsss",
    fs=FS,
    snr=BASE_EBN0_DB,
    snr_mode="ebno",
    seed=11,
    sps=SPC,  # samples per CHIP
    acq_code=ACQ,
    acq_reps=REPS,
    data_code=DATA,
    frame=FRAME,
    data=PAYLOAD,  # PN data, and nothing after it
    delay_samples=LEAD,
    # the receiver holds a burst until it has seen refine_span past its
    # end, so that much trailing noise is what lets it emit
    off_samples=rx.refine_span,
    gap_noise="auto",
)
plan = prepare(scene)  # the burst's signal is rendered ONCE, here

# The anchor contract: at the scene's own Eb/N0 and the anchor seed, a Plan
# draw IS a full compose of the same scene, byte for byte.
anchor = plan.at(BASE_EBN0_DB, plan.anchor_seed)
assert anchor.tobytes() == scene.compose().tobytes(), (
    "plan.at(base_ebn0, anchor_seed) is not byte-identical to compose()"
)

# Building a receiver solves its detection design (the Pd/Pfa threshold)
# every time; reset() keeps the design and clears the stream. A reused
# receiver must decode a trial exactly as a fresh one does -- pinned here,
# on a draw pushed AFTER the reused one has already decoded another.
rx.push(plan.at(BASE_EBN0_DB, 1))
rx.reset()
x = plan.at(BASE_EBN0_DB, 2)
fresh = receiver()
assert np.array_equal(np.asarray(rx.push(x)), np.asarray(fresh.push(x)))
assert rx.events().tobytes() == fresh.events().tobytes(), (
    "a reset() receiver's events differ from a fresh receiver's"
)
rx.reset()

def measure(ebn0_db: float, seed0: int, meter: BerMeter) -> dict:
    """One Eb/N0 point: trials until TARGET_ERRORS (and MIN_TRIALS bursts).

    Each trial is a fresh noise draw of the same cached burst, pushed whole
    through the one receiver, reset between trials (output-identical to a
    fresh one, asserted above). The window scored is the one the receiver
    decoded at the burst's true start; the payload is the frame after the
    sync word, compared bit for bit with the Field's own bits.
    """
    ideal = ber_theory_ber(2, 10.0 ** (ebn0_db / 10.0))
    max_trials = max(
        MIN_TRIALS, int(BUDGET * TARGET_ERRORS / (ideal * PAYLOAD.size))
    )
    errors = bits = misses = trials = 0
    draws = set()  # every trial's noise must be its own
    while (errors < TARGET_ERRORS or trials < MIN_TRIALS) and (
        trials < max_trials
    ):
        x = plan.at(ebn0_db, seed0 + trials)
        trials += 1
        draws.add(hash(x.tobytes()))
        rx.reset()
        out = np.asarray(rx.push(x))
        hit = [
            k
            for k, e in enumerate(rx.events())
            if int(e["preamble_start"]) == LEAD
        ]
        if not hit:
            misses += 1
            continue
        frame = out[hit[0] * FRAME_SYMS :][:FRAME_SYMS]
        got = frame[len(SYNC) :]
        errors += int(np.count_nonzero(got != PAYLOAD))
        bits += PAYLOAD.size
    return {
        "ebn0": ebn0_db,
        "trials": trials,
        "misses": misses,
        "distinct": len(draws),
        "ci": meter.interval(errors, bits),
        "quoted": errors >= TARGET_ERRORS and trials >= MIN_TRIALS,
    }

def sync_phase_limited_ber(ebn0_db: float) -> float:
    """BPSK BER for THIS receiver's carrier-phase estimate, not a fit.

    BurstDemod is feedforward: the residual carrier phase comes from the
    complex peak of the sync-word correlation (burst_demod_core.h, step 4),
    a sum over len(SYNC) = 13 despread symbols each at Es/N0 = Eb/N0. That
    estimate's error is Gaussian with variance 1 / (2 * 13 * Eb/N0), and a
    phase error phi scales the decision amplitude by cos(phi), so the
    expected BER is  E_phi[ Q(sqrt(2 Eb/N0) cos phi) ].  Averaged by
    Gauss-Hermite quadrature over doppler's own Q, it is ideal BPSK plus
    the loss this estimator must cost: ~0.2 dB at 0 dB, falling as 1/Eb/N0.
    """
    g = 10.0 ** (ebn0_db / 10.0)
    sigma = np.sqrt(1.0 / (2.0 * len(SYNC) * g))
    x, w = np.polynomial.hermite_e.hermegauss(40)
    ber = [ber_theory_ber(2, g * np.cos(sigma * xi) ** 2) for xi in x]
    return float(np.dot(w, ber) / w.sum())

    quoted = 0
    for p in points:
        ci, eb = p["ci"], p["ebn0"]
        ideal = ber_theory_ber(2, 10.0 ** (eb / 10.0))
        limited = sync_phase_limited_ber(eb)
        # Stated only when the point stopped on errors: a point that ran
        # out of trials has not resolved its rate and is not quoted.
        loss = eb - ber_esn0_db_for_ser(2, ci.p_hat) if p["quoted"] else None
        print(
            f"{eb:6.1f} {p['trials']:6d} {p['misses']:4d} {ci.errors:6d} "
            f"{ci.p_hat:9.3e} [{ci.lo:9.3e}, {ci.hi:9.3e}] {ideal:9.3e} "
            f"{limited:9.3e} "
            + (f"{loss:+6.2f}dB" if loss is not None else "  (n/q)")
        )
        # Independent Monte Carlo: every trial drew its own noise.
        assert p["distinct"] == p["trials"], (
            f"{eb} dB: {p['trials']} trials but {p['distinct']} distinct "
            "draws -- the Plan's seed is not varying the noise"
        )
        # Acquisition met its design: the miss rate's interval admits 1-pd.
        miss = meter.interval(p["misses"], p["trials"])
        assert miss.lo <= 1.0 - PD, (
            f"{eb} dB: {p['misses']}/{p['trials']} bursts missed, more than "
            f"a pd={PD} acquisition allows"
        )
        quoted += bool(p["quoted"])
        # Both limits hold for EVERY point, quoted or not: an unresolved
        # point's interval is wide, not wrong. The receiver cannot beat
        # coherent BPSK ...
        assert ci.hi >= ideal, (
            f"{eb} dB: BER interval [{ci.lo:.3e}, {ci.hi:.3e}] lies below "
            f"ideal BPSK {ideal:.3e} -- the Eb/N0 axis is wrong"
        )
        # ... and must be no worse than its own sync-phase estimate costs.
        assert ci.lo <= limited, (
            f"{eb} dB: BER interval [{ci.lo:.3e}, {ci.hi:.3e}] lies above "
            f"the sync-phase-limited {limited:.3e} -- a loss the receiver's "
            "phase estimate does not account for"
        )
    assert quoted == len(points), (
        f"only {quoted}/{len(points)} points reached {TARGET_ERRORS} errors"
    )