Arbitrary-Rate Symbol Recovery¶

import numpy as np

from doppler.ber import ber_evm_db, ber_settle_syms
from doppler.track import RateSync
from doppler.wfm import rrc_h

SPS = 17.33389  # non-integer samples per symbol -- the whole point
BETA = 0.35
SPAN = 8
BN = 0.005  # loop bandwidth, symbol-rate normalised
CLOCK_PPM = 200.0  # the true rate is this far off the nominal
ES_N0_DB = 15.0  # symbol energy to noise density -- NOT a per-sample SNR

# Long enough that a SETTLED window exists. `ber_settle_syms` is the
# library's own answer to where steady state may start -- 2*(5/Bn) for the
# one running loop, which is 2000 symbols at Bn = 0.005. A record pinned to
# a round number instead (this demo used 1500) is shorter than its own
# settling budget, so "the last quarter" was measuring the tail of the
# acquisition transient and calling it steady state.
NSYM = 3 * int(ber_settle_syms(BN, 0.0))


def make_signal(sps, tau=0.37, seed=7, es_n0_db=ES_N0_DB):
    """RRC-shaped BPSK at a fractional sps, offset by tau, in AWGN.

    Noise is set from **Es/N0**, the symbol energy over the noise density,
    because "SNR" alone is ambiguous here: at 17.33 samples per symbol a
    per-sample SNR is ~12 dB below the Es/N0 the receiver actually sees, so
    the same number means two very different links.

    The noise is **complex**, N0 total (N0/2 per dimension), because a real
    receiver's baseband is complex and its Q channel carries noise even when
    the modulation is real. Injecting real-only noise instead is the same
    thing as discarding Q, and it flatters the reported EVM by 3 dB against
    the convention everyone else quotes -- EVM is measured on the I/Q plane
    unless it says otherwise.

    Symbols are presented at the CONTRACTED unit amplitude, and that is the
    whole of the level story. RateSync carries no AGC by design -- a
    composing receiver already levels in its own front-end cascade -- so the
    caller owns the input level, and the level to hit is not a tuned number:
    the TED normalises by its own construct-time slope, computed for the
    reference the matched bank already defines
    (``RateConverter_agc_ref_db()`` = ``10*log10(bank_e0/bank_sps)``, ~0 dB
    because the bank normalises by its own pulse energy).

    This demo used to scale the shaped stream to ``0.25`` of its PEAK for
    "CIC headroom", and it is worth naming what that cost. An RRC stream
    peaks at ~1.582x its symbol amplitude, so symbols arrived at ~0.158
    against a contract written in unit amplitude -- and a Gardner detector's
    slope goes as A^2, so the loop ran ~40x under its stated bandwidth and
    failed its own lock assertion with nothing pointing at the level. The
    peak backoff was also a hand-rolled AGC: a peak detector with no loop
    and no time constant.

    That ``A^2`` is **Gardner's exponent, not the object's**: the two
    detectors do not share an amplitude law. Gardner's raw error is a
    product of two signal samples and carries ``A^2``; DTTL's multiplies
    one signal sample by a difference of hard-decision SIGNS, which is
    amplitude-free, so it carries ``A^1``. Both are measured -- fitted
    2.000 and 1.000 -- in section 2.6b of the validation report. The level
    CONTRACT is the same either way, but the penalty for missing it is
    squared for the default detector and linear for the other, so the same
    0.158 mistake above would have cost ~6x rather than ~40x under DTTL.

    The CIC's input bound is **2.0**, not 1.0 -- ``CIC_PAPR_HEADROOM``
    reserves exactly the 6 dB a unit-amplitude RRC's 1.582 peak needs, so
    the contracted level fits with margin. What does not fit is the NOISE:
    15 dB Es/N0 at 17.33 samples per symbol is only ~2.6 dB per-sample SNR,
    so the composite crest runs past the bound on noise alone. See the note
    beside the assertions below, and gh-668.
    """
    rng = np.random.default_rng(seed)
    syms = np.where(rng.integers(0, 2, NSYM) > 0, 1.0, -1.0)
    n = int(NSYM * sps) + 64
    idx = np.arange(n, dtype=float)
    x = np.zeros(n)
    for k, a in enumerate(syms):
        t = (idx - (k + SPAN) * sps) / sps - tau
        near = np.abs(t) <= SPAN
        x[near] += a * rrc_h(t[near], BETA)
    # Average symbol energy from the STEADY-STATE span only. Including the
    # ramp-up and the tail padding underestimates the power, which would
    # quietly set a higher Es/N0 than advertised -- and the giveaway is an
    # EVM that beats the matched-filter bound, which nothing can do.
    core = x[int(SPAN * sps) : -int(SPAN * sps)]
    es = np.mean(core**2) * sps
    n0 = es / 10 ** (es_n0_db / 10)
    noise = (rng.standard_normal(n) + 1j * rng.standard_normal(n)) * np.sqrt(
        n0 / 2
    )
    return (x + noise).astype(np.complex64), syms


def evm_floor_db(es_n0_db=ES_N0_DB):
    """Matched-filter EVM bound on the I/Q plane.

    At the matched-filter output the error vector is the complex noise, of
    total variance N0, against a reference of energy Es -- so
    EVM^2 = N0/Es and the bound in dB is simply -(Es/N0). (The familiar
    factor of two belongs to an I-only measurement, which discards the Q
    channel; EVM is a plane quantity unless stated otherwise.)
    """
    return -es_n0_db


# The receiver: one object, one call. `sps` is the NOMINAL rate -- the loop
# finds the true one, which is 200 ppm away.
true_sps = SPS * (1.0 + CLOCK_PPM * 1e-6)
rx, tx_syms = make_signal(true_sps)

sync = RateSync(sps=SPS, pulse="rrc", beta=BETA, span=SPAN, m=2, bn=BN)
symbols = np.asarray(sync.steps(rx))

from doppler.resample import MatchedRateConverter  # noqa: E402


# Why a high input rate is nearly free: the cascade's HB/CIC stages do the
# bulk decimation at no multiplies, so the matched filter is sized by the
# POST-decimation rate. Applied at the input rate it would grow with sps.
def bank_cost(sps_values, span=SPAN, m=2):
    """(taps/arm on the cascade, taps/arm at the input rate) per sps."""
    cascade, at_input = [], []
    for s in sps_values:
        rc = MatchedRateConverter(
            rate=m / s,
            compensate=1,
            pulse="rrc",
            span=span,
            pulse_sps=float(m),
        )
        cascade.append(rc.bank_shape[1])
        at_input.append(int(np.ceil((2 * span + 1.0 / m) * s)) + 1)
    return np.array(cascade), np.array(at_input)

def compare_detectors(rx_signal=None):
    """Recover the same stream with each timing-error detector.

    `ted` is the one knob on this object whose two settings are both
    correct, so it is a CHOICE rather than a default -- and a demo that
    only ever builds the default leaves a caller with no way to make it.

    Returns ``{name: (evm_db, lock_stat, locked)}``.
    """
    rx_signal = rx if rx_signal is None else rx_signal
    out = {}
    for name in ("gardner", "dttl"):
        obj = RateSync(
            sps=SPS, pulse="rrc", beta=BETA, span=SPAN, m=2, bn=BN, ted=name
        )
        y = np.asarray(obj.steps(rx_signal))
        out[name] = (_evm_db(y), float(obj.lock_stat), bool(obj.locked))
    return out


# How to choose, given both lock and (at this Es/N0) both land on the
# matched-filter floor:
#
#   gardner  blind. Works for any constellation, and its `bn` means the
#            same loop bandwidth at every roll-off -- measured flat across
#            beta 0.1 to 0.9. The default, and the safe one.
#   dttl     decision-directed. Lower self-noise near lock: on a NOISELESS
#            stream it rests ~6x closer to the eye centre, because Gardner
#            pays a self-noise cost on every non-transition symbol and DTTL
#            gates those out. Costs: BPSK/QPSK only (it assumes rectangular
#            I/Q decision boundaries), and its effective loop bandwidth
#            currently varies ~8x across the roll-off range, so `bn` does
#            NOT mean one bandwidth for it -- doppler#669.
#
# The EVM the two reach HERE is the same, and that is the point worth
# taking away: at a realistic Es/N0 the noise dominates and the detector
# choice does not move the number. DTTL's advantage is a self-noise
# advantage, so it shows where self-noise is what is left.