Streaming Async Despreader¶

import numpy as np

from doppler.detection import (
    det_cn0_to_snr,
    det_n_noncoh,
    det_threshold_noncoherent,
    det_verify_count,
)
from doppler.impairment import DopplerChannel
from doppler.source import AWGN, LO
from doppler.track import Dll
from doppler.wfm import Gold, Synth, wfm_awgn_amplitude

# ── Setup: this scenario's numbers ─────────────────────────────────────────
SF = 1023  # Gold code period, chips (CCSDS)
CHIP_RATE = 3.069e6  # chips per second
SPC = 2  # samples per chip
FS = CHIP_RATE * SPC  # sample rate, 6.138 MSa/s
SYM_RATE = 2700.0  # BPSK symbols per second, asynchronous to the code
TE = SF * SPC  # one code period, samples: one push
EPOCH_S = TE / FS  # seconds of signal per push
CARRIER_HZ = 2.5e9  # the RF carrier the motion shifts
CARRIER_OFFSET_HZ = 1500.0  # this link's carrier offset
LINK_PPM = 1e6 * CARRIER_OFFSET_HZ / CARRIER_HZ  # the same motion: 0.6 ppm
INFLATED_PPM = 5.0  # inflated, to show what the aid is for
ESN0_DB = 10.0  # Es/N0 at the receiver input
CN0_DBHZ = ESN0_DB + 10.0 * np.log10(SYM_RATE)  # C/N0, dB-Hz
DURATION_S = 0.6  # seconds of signal per run

GOLD = np.asarray(Gold().generate(SF), np.uint8)

#: The code NCO's rate step, chips per chip: one phase increment of its
#: 32-bit accumulator, which turns once per code period.
NCO_STEP = TE / 2**32


def receive(doppler_ppm, seed=2):
    """The link at a Doppler: the emitter through a channel, then noise.

    The emitter is :class:`~doppler.wfm.Synth` in continuous DSSS mode on
    nominal frequency, with its own PRBS payload. The
    :class:`~doppler.impairment.DopplerChannel` applies the motion: it
    dilates every clock in the stream by ``1 + doppler_ppm * 1e-6`` and
    shifts the carrier by ``CARRIER_HZ * doppler_ppm * 1e-6``, from the one
    parameter. The noise is :class:`~doppler.source.AWGN` at the Es/N0.

    Returns
    -------
    x : NDArray[np.complex64]
        The received stream.
    channel : DopplerChannel
        The channel. Its ``offset_hz`` stands in for the carrier estimate an
        upstream carrier loop hands over, and its ``delay_samples`` for the
        code phase acquisition hands over.
    """
    n = int(DURATION_S / EPOCH_S) * TE
    tx = np.asarray(
        Synth(
            type="dsss",
            fs=FS,
            sps=SPC,
            data_code=GOLD.tobytes(),
            symbol_rate=SYM_RATE,
        ).steps(n),
        np.complex64,
    )
    channel = DopplerChannel(
        fs=FS, carrier_hz=CARRIER_HZ, doppler_ppm=doppler_ppm
    )
    # The channel takes at most 65536 samples a call.
    rx = np.concatenate(
        [channel.execute(tx[i : i + 65536]) for i in range(0, n, 65536)]
    )
    amp = float(
        wfm_awgn_amplitude(ESN0_DB - 10.0 * np.log10(FS / SYM_RATE), 1.0)
    )
    noise = AWGN(seed=seed, amplitude=amp).generate(len(rx))
    return (rx + noise).astype(np.complex64), channel

BN = 0.002  # narrow, for a quiet code phase: what the aid lets the loop keep
SEGMENTS = (
    4  # despread samples per code period, for an asynchronous data clock
)
PFA = 1e-3  # the lock detector's false-alarm rate per decision
PD = 0.99  # its detection rate per decision, at the link's C/N0

# One row per push, read from the loop's own properties.
ROW = np.dtype(
    [
        ("locked", bool),
        ("lock_stat", float),
        ("filter_out", float),
        ("samples", int),
    ]
)


def despreader(init_chip):
    """A Dll for this link, its lock detector sized for the link.

    Given the symbol period, the detector's looks are coherent over the
    transition-free part of a symbol rather than one despread sample.
    ``det_n_noncoh`` gives the looks a decision needs for ``PD`` at the C/N0,
    and ``det_verify_count`` the consecutive misses that hold false drops
    under 1e-6 per decision. ``AsyncDsssReceiver`` sizes its own the same
    way, because at its defaults, one look per despread sample, the flag
    drops on a loop that holds the code: measured on the receiver, a few
    times a second at Es/N0 10.7 dB and 96% of the time at 5.7 dB
    (``docs/design/async-dsss-receiver-measurements.md`` section 12.3).

    Returns the Dll and its lock threshold.
    """
    dll = Dll(GOLD, sps=SPC, init_chip=init_chip, bn=BN, segments=SEGMENTS)
    dll.set_symbol_period(SEGMENTS * CHIP_RATE / (SF * SYM_RATE))
    look = dll.symbol_window * (TE // SEGMENTS)  # samples per look
    n_looks = det_n_noncoh(det_cn0_to_snr(CN0_DBHZ, FS), look, PD, PFA, 4000)
    dll.configure_lock(PFA, n_looks)
    dll.set_lock_verify(
        det_verify_count(PFA, PFA * 1e-3), det_verify_count(1.0 - PD, 1e-6)
    )
    return dll, det_threshold_noncoherent(PFA, n_looks)


def track(x, channel, aided):
    """Wipe the carrier, then despread with a Dll, one code period per push.

    One carrier estimate does two jobs. It wipes the carrier, and divided by
    the carrier frequency it is the code-rate aid. The loop starts at the code
    phase acquisition reports, to one sample. The code phase runs opposite to
    the delay, so a code that arrives ``L`` chips late starts the loop at
    ``SF - L``.

    Returns one :data:`ROW` per push, the lock threshold, and the despread
    output: ``SEGMENTS`` samples per code period, the data with the code
    removed.
    """
    carrier_estimate_hz = channel.offset_hz
    wiped = x * LO(-carrier_estimate_hz / FS).steps(len(x))

    lag_chips = np.floor(channel.delay_samples) / SPC
    dll, threshold = despreader(SF - lag_chips)
    if aided:
        dll.set_rate_aid(carrier_estimate_hz / CARRIER_HZ)

    rows = np.zeros(len(wiped) // TE, ROW)
    output = []
    for k in range(len(rows)):
        despread = dll.steps(wiped[k * TE : (k + 1) * TE])
        output.append(despread)
        rows[k] = (
            dll.locked,
            dll.lock_stat,
            dll.code_rate - 1.0,  # the loop filter output, chips per chip
            len(despread),
        )
    return rows, threshold, np.concatenate(output)

def check(runs, outputs):
    """Assert what the runs show, and return what this run measured.

    ``runs`` maps ``(doppler_ppm, aided)`` to a run's rows, and ``outputs``
    to its despread output.

    Three bounds come from the object and the link, not the draw. The
    despread output is ``SEGMENTS`` samples per code period. A locked loop's
    output is the data along one axis, at whatever carrier phase the data
    sits: its power along that axis stands above its power across it by up
    to the link budget's ``1 + 2 * SNR`` per output sample, and reads below
    it because a sample straddling a data flip loses amplitude. The check
    takes half the budget as the line, below which the output is
    noise. And the NCO truncates its
    command, so a locked loop's settled filter output reads between the rate
    the loop carries and one NCO step above it. The check widens that by a
    quarter step either side, an allowance for the noise left in a mean over
    the settled half that is not part of the bound. The rest describes this
    run: first lock within a generous 50 ms, which loops held lock, and when
    the unaided loop lost it.
    """
    ms = 1e3 * EPOCH_S
    got = {}
    first = {}
    for (ppm, aided), rows in runs.items():
        assert abs(rows["samples"].sum() - SEGMENTS * len(rows)) <= SEGMENTS, (
            "the despread output is not SEGMENTS samples per code period"
        )
        assert rows["locked"].any(), f"{ppm} ppm aided={aided}: never locked"
        first[ppm, aided] = int(np.argmax(rows["locked"]))
        assert first[ppm, aided] * ms < 50.0, (
            f"{ppm} ppm aided={aided}: first lock at "
            f"{first[ppm, aided] * ms:.0f} ms"
        )

    # This run: aided, the loop holds lock from its first lock at both
    # Dopplers; unaided it holds this link's, and loses the inflated one.
    for key in ((LINK_PPM, True), (INFLATED_PPM, True), (LINK_PPM, False)):
        assert runs[key]["locked"][first[key] :].all(), (
            f"this run: {key[0]} ppm aided={key[1]} lost lock"
        )
    unaided = runs[INFLATED_PPM, False]["locked"]
    assert not unaided[len(unaided) // 2 :].any(), (
        "this run: the unaided loop held the inflated Doppler"
    )
    got["inflated_unaided_lost_ms"] = (
        len(unaided) - int(np.argmax(unaided[::-1]))
    ) * ms

    # The output: the data along one axis while the loop holds the code,
    # noise once it has lost it. The two eigenvalues of the output's I/Q
    # second moments are its power along its strongest axis and across it,
    # whatever carrier phase the data sits at. Per output sample the SNR is
    # Es/N0 over the samples a symbol spans.
    per_symbol = SEGMENTS * CHIP_RATE / (SF * SYM_RATE)
    budget = 1.0 + 2.0 * 10.0 ** (ESN0_DB / 10.0) / per_symbol
    for (ppm, aided), out in outputs.items():
        settled = out[len(out) // 2 :]
        iq = np.stack([settled.real, settled.imag]).astype(np.float64)
        across, along = np.linalg.eigvalsh(iq @ iq.T / iq.shape[1])
        ratio = float(along / across)
        holds = (ppm, aided) != (INFLATED_PPM, False)
        assert (ratio > budget / 2) == holds, (
            f"{ppm} ppm aided={aided}: output power along its axis against "
            f"across it {ratio:.2f}, against the budget's {budget:.2f}"
        )
        got[f"{ppm}_{aided}_axis"] = ratio
    got["axis_budget"] = budget

    # Physical: a locked loop's settled filter output is the rate it carries,
    # read high by up to one NCO step. Unaided, the loop carries the
    # dilation. Aided, the aid carries it and the loop carries nothing.
    for ppm, aided in (
        (LINK_PPM, False),
        (LINK_PPM, True),
        (INFLATED_PPM, True),
    ):
        rows = runs[ppm, aided]
        carried = 0.0 if aided else ppm * 1e-6
        mean = float(rows["filter_out"][len(rows) // 2 :].mean())
        steps = (mean - carried) / NCO_STEP
        assert -0.25 <= steps <= 1.25, (
            f"{ppm} ppm aided={aided}: filter output {1e6 * mean:+.2f} ppm, "
            f"{steps:+.2f} NCO steps from the {1e6 * carried:.2f} ppm the "
            "loop carries"
        )
        got[f"{ppm}_{aided}_filter_ppm"] = 1e6 * mean
    return got