Skip to content

Arbitrary-Rate Symbol Recovery

Arbitrary-rate matched filtering and timing recovery

track.RateSync recovers symbols from a stream whose sample clock has no integer relationship to the symbol clock — here 17.33389 samples per symbol, and 200 ppm off even that.

It gets there by inverting the usual arrangement. Where SymbolSync runs a matched FIR and then a separate Farrow interpolator steered by an integer NCO, RateSync owns a MatchedRateConverter whose terminal stage carries the pulse. The cascade's last dot product is the matched filter, and the polyphase arm that dot product selects is the fractional timing delay. One filter, no Farrow, no separate matched-filtering pass — and because that stage's accumulator is a double, sps is a double.

What you're seeing

Top — Recovered symbols. The real part per symbol index: a short pull-in while the loop acquires, then clean rails — at a fractional sps, with the clock 200 ppm off nominal.

The noise is quoted as Es/N0, not as an "SNR": at 17.33 samples per symbol those differ by more than 12 dB, so the same number would describe two very different links. The noise is complex, N0 total — a real receiver's baseband is complex and its Q channel carries noise even when the modulation is real.

At the matched-filter output the error vector is that complex noise, of total variance N0, against a reference of energy Es. So on the I/Q plane

\[ \\mathrm{EVM}^2 = \\frac{N_0}{E_s} \\qquad\\Longleftrightarrow\\qquad \\mathrm{EVM},[\\mathrm{dB}] = -\\bigl(E_s/N_0\\bigr)[\\mathrm{dB}] \]

(The familiar factor of two belongs to an I-only measurement, which discards the Q channel. EVM is a plane quantity unless it says otherwise.)

Measured, ten seeds per point:

Es/N0 bound measured offset
10 dB −10.0 dB −10.02 dB −0.02 dB
15 dB −15.0 dB −15.01 dB −0.01 dB
20 dB −20.0 dB −19.97 dB +0.03 dB

So the fused matched filter is on the bound, across a 10 dB span — fusing it into the resampler's polyphase bank costs nothing in detection performance. It is the same filter, evaluated at a steered instant.

(One 375-symbol window carries ~0.3 dB of estimator noise, so a single seed lands either side. The script asserts ±1.5 dB at three operating points rather than pretending one measurement is tighter than it is.)

Middle — Tracked clock. RateSync.rate pulling off the nominal it was constructed with (dotted) onto the stream's true rate (dashed). The loop is tracking an asynchronous, non-integer clock; rate reports the truth, which is what a caller disciplining a clock reads. It is taken from the loop integrator, not the instantaneous control — the integrator is the rate memory, and reading the noisy total instead biases the estimate.

Bottom — Matched-filter cost. Taps per arm against input samples per symbol. Applied at the input rate (red) a root-raised-cosine matched filter grows linearly with sps — 4225 taps per arm at 256 samples/symbol, some 35 MB of bank. Riding the cascade's terminal stage (blue) it is flat: the HB/CIC stages ahead of it have already done the bulk decimation, at no multiplies, so the filter is sized by the post-decimation rate. The single 34 → 40 step is the CIC droop compensator folding into the bank once the planner picks a CIC; a halfband cascade has no droop to correct.

How it works

The receiver is one object and one call. sps is the nominal rate — the loop finds the true one:

import numpy as np

from doppler.track import RateSync

SPS = 17.33389  # non-integer samples per symbol -- the whole point
BETA = 0.35
SPAN = 8
NSYM = 1500
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


def rrc_h(t, beta=BETA):
    """Analytic RRC pulse at arbitrary (non-grid) times, in symbol periods."""
    t = np.asarray(t, float)
    out = np.empty_like(t)
    pt = np.pi * t
    at_zero = np.abs(t) < 1e-9
    at_sing = (
        np.abs(np.abs(t) - 1.0 / (4 * beta)) < 1e-9
        if beta > 0
        else np.zeros(t.shape, bool)
    )
    gen = ~(at_zero | at_sing)
    out[at_zero] = 1.0 - beta + 4.0 * beta / np.pi
    if at_sing.any():
        a = np.pi / (4 * beta)
        out[at_sing] = (beta / np.sqrt(2)) * (
            (1 + 2 / np.pi) * np.sin(a) + (1 - 2 / np.pi) * np.cos(a)
        )
    tg, ptg = t[gen], pt[gen]
    num = np.sin(ptg * (1 - beta)) + 4 * beta * tg * np.cos(ptg * (1 + beta))
    out[gen] = num / (ptg * (1 - (4 * beta * tg) ** 2))
    return out


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.

    Amplitude is kept well inside +-1.0: the planned cascade contains a CIC,
    which bounds its input and clips silently past that (RateSync.clipped is
    the only signal, and this demo asserts it stays clear).
    """
    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])
    x *= 0.25 / np.max(np.abs(x))  # headroom for the CIC's +-1.0 bound
    # 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=0.005)
symbols = np.asarray(sync.steps(rx))

RateSync asks its MatchedRateConverter for rate = m/sps and lets the planner decide the shape. At sps = 17.33389 with m = 2 that is CIC(8) followed by a Resampler(0.923, rrc): the CIC throws away the bulk of the rate for free, and the fractional remainder lands on the stage that carries the pulse. Ask for sps = 64 instead and the plan becomes CIC(32) + Resampler(1.0, rrc) — the same bank, a 32× cheaper front end.

That terminal stage always exists when a pulse is selected, even when the rate divides exactly. It is not there to correct the rate; it is the matched filter and the timing element, and a cascade that ended in a bare CIC(32) would have nothing steerable at the end.

Every m-th terminal output is the on-time strobe and the output m/2 back is the Gardner transition gate, so the oversampled stream and the symbol stream come out of the same dot products. A half-symbol error in that assignment is an equilibrium of the detector — but an unstable one: each parity's S-curve has one zero at the eye centre with negative slope and one at the T/2 point with positive slope, so the loop runs away from the wrong one unaided. No eye-sign detector, no counter flip, and no second bank.

The cost claim in the third panel is measured, not asserted:

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)

Reading the object honestly

Two habits worth forming, both of which cost real debugging time to learn:

Judge lock by lock_stat / locked, not by an error-vector magnitude. A single cycle slip during acquisition drags a windowed EVM by 20 dB while the eye is wide open at +0.75. The eye statistic is the honest indicator; EVM is only meaningful once the loop has settled.

Check clipped at least once against real input. The cascade inherits its CIC's ±1.0 input bound, and an overdriven front end costs about 25 dB of EVM with a perfectly healthy lock — nothing in the timing metrics reveals it. That is exactly why the flag exists.

One tuning note: use m >= 4 with pulse="iandd". The rectangle is one symbol wide, so at m = 2 its matched filter is a two-tap sum and the eye barely opens (lock_stat −0.34 at m=2 against +0.95 at m=4 on the same NRZ stream). The RRC spans many symbols and is unaffected.

SymbolSync is unchanged and remains the right answer when the matched filter is not one this family builds, or when the front end is already at a small integer sps and a Farrow interpolator is the cheaper shape.

Reproduce

uv run python src/doppler/examples/ratesync_demo.py ratesync_demo.png

Source: src/doppler/examples/ratesync_demo.py