Timing Synchronization & Timing Error Detectors (TED)

Est. read time: 8 minutes | Last updated: July 27, 2026 by John Gentile

Contents

Timing Error Detector (TED)
- Gardner TED
- Mueller & Muller (M&M) TED
Polyphase Matched Filter for Timing Synchronization
References

import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
from rfproto import filter, measurements, modulation, pi_filter, plot, sig_gen

Timing Error Detector (TED)

# simulate random binary input values
num_symbols  = 2400
sym_rate     = 1e6 # Baseband symbol rate
# Generate random QPSK symbols
rand_symbols = np.random.randint(0, 4, num_symbols)

L  = 2               # Upsample ratio (sim 2x Samples per Symbol)
fs = L * sym_rate    # Output sample rate (Hz)

rolloff          = 0.5 # Alpha of RRC
num_filt_symbols = 6    # Symbol length of RRC matched filter

qpsk_tx_filtered_2x = sig_gen.gen_mod_signal(
    "QPSK",
    rand_symbols,
    fs,
    sym_rate,
    "RRC",
    rolloff,
    num_filt_symbols,
)

# Add timing offset (0.3 samples delay) via fractional delay filter
timing_offset = 0.3
frac_N = 21 # number of taps in fractional filter (should be odd)
frac_idx = np.arange(-(frac_N-1)//2, frac_N//2+1)
frac_h = np.sinc(frac_idx - timing_offset) # calc filter taps
frac_h *= np.hamming(frac_N) # window filter to make sure it decays to 0 on both sides
frac_h /= np.sum(frac_h) # normalize to get unity gain filter weights
samples_rx = np.convolve(qpsk_tx_filtered_2x, frac_h)

# Add AWGN noise
#samples_rx += 0.05 * (np.random.randn(len(samples_rx)) + 1j * np.random.randn(len(samples_rx)))

plot.IQ(samples_rx, alpha=0.1)
plt.show()

png

plot.spec_an(samples_rx, fs=fs, fft_shift=True, show_SFDR=False, y_unit="dB", title=f"QPSK ({L}x SPS)")
plt.show()

png

plt.plot(samples_rx[20:40].real, '-o', linewidth=0.5)
plt.plot(qpsk_tx_filtered_2x[20:40].real, '-o', linewidth=0.5)
plt.show()

png

Gardner TED

From this DSP StackExchange Post

e(n) = \mathbb{R}[\overline{y_{n}}(y_{n+1} - y_{n-1})]

early = prompt = late = 0.0 * 1j*0.0
gardner_ted_out = []

# Perform TED error calculation -> real{conj(y[n]) * (y[n+1] - y[n-1])}
#  early (y[2n - 1]), prompt (y[2n]), late (y[2n + 1])
def gardner_ted(early, prompt, late):
    return (np.real(prompt)*(np.real(late)-np.real(early))) + (np.imag(prompt)*(np.imag(late)-np.imag(early)))

for sample in samples_rx[:50]:
    # Shift in samples at 2x SPS rate to TED early/prompt/late registers
    early = prompt # y[2n - 1]
    prompt = late  # y[2n]
    late = sample  # y[2n + 1]
    
    gardner_ted_out.append(gardner_ted(early, prompt, late))

plt.plot(gardner_ted_out, linewidth=0.5)
plt.show()

png

Mueller & Muller (M&M) TED

e(n) = \mathbb{R}[(y_{n} - y_{n-1}) * \overline{\text{sign}(y_{n-1})}] - \mathbb{R}[(\text{sign}(y_{n}) - \text{sign}(y_{n-1})) * \overline{y_{n-1}}]

rrc_coef = filter.RootRaisedCosine(L, 1, rolloff, 2 * num_filt_symbols * L + 1)
samples_rx_post_mf = signal.lfilter(rrc_coef, 1, samples_rx)
# downsample to 1x SPS
samples_rx_post_mf = samples_rx_post_mf[::2]

plot.IQ(samples_rx_post_mf, alpha=0.1)
plt.show()

png

# current (y[n]) and previous (y[n-1])
def mm_ted(curr, prev):
    mm_real = (np.real(curr) * np.sign(np.real(prev))) - (np.real(prev) * np.sign(np.real(curr)))
    mm_imag = (np.imag(curr) * np.sign(np.imag(prev))) - (np.imag(prev) * np.sign(np.imag(curr)))
    return mm_real + mm_imag

mm_out = []
for i in range(1, len(samples_rx_post_mf)):
    mm_out.append(mm_ted(samples_rx_post_mf[i], samples_rx_post_mf[i-1]))

plt.plot(mm_out, linewidth=0.5)
plt.show()

print(f"Mean timing error: {np.mean(mm_out):.4f}")
print(f"Error std dev: {np.std(mm_out):.4f}")

png

Mean timing error: 0.2344 Error std dev: 0.1271

Polyphase Matched Filter for Timing Synchronization

From the Multirate DSP page on Polyphase Filters, we know that each leg of a polyphase filter is an all-pass filter, with each leg providing a fractional amount of signal delay ( $1/M$ samples of delay, where $M$ is the number of polyphase branches), spread equally over the duration of a symbol period. We can exploit this feature by coopting the commutator- traditionally used to cycle through each polyphase leg to perform the interpolation or decimation operation- to act as the delay selector in a timing error correction block!

This also leads to implementation optimizations, as only one FIR filter structure need be instantiated, with the commutator indexing a Look-Up Table (LUT) of each polyphase leg’s filter coefficients at a given time.

In a simple simulated system above, we were able to generate a QPSK signal (no frequency or phase offset present) at 2x Samples per Symbol (SPS). Below, we will see the ideal sampling location (minimal deviation from ideal constellation points in I/Q plane) can be found after matched filtering with a polyphase filter implementation, selecting a timing offset based off commutator/branch, and final downsampling to 1x SPS.

For the Polyphase Interpolation Filter case, we start by generating a similar Root Raised Cosine (RRC) filter as before, however the prototype filter is scaled by the number of polyphase filter legs- as done in other polyphase prototype designs- to maintain the correct matched filter response.

# The number of taps of this filter is based on how long you expect the channel to be; that is, 
# how many symbols do you want to combine to get the current symbols energy back, usually 5 to 10+
taps = 2 * num_filt_symbols * L + 1
# With 32 filters, you get a good enough resolution in the phase to produce very small, almost 
# unnoticeable, ISI. Going to 64 filters can reduce this more, but after that there is very little 
# gain for the extra complexity. Total prototype filter taps = taps * num_filters, since we're 
# instantiating segments of these taps into the filterbanks in such a way that each bank now 
# represents the filter at different phases, equally spaced at 2pi/N, where N is the number of filters.
num_filters = 31
polyphase_rrc_coef = filter.RootRaisedCosine(num_filters * L, 1, rolloff, taps * num_filters)

plot.filter_response(polyphase_rrc_coef, title=f"{num_filters} x {taps}-tap RRC Polyphase Prototype Filter")
plot.plt.show()

png

h_poly_rrc = polyphase_rrc_coef.reshape(len(polyphase_rrc_coef)//num_filters, num_filters).T
print(np.shape(h_poly_rrc))

(31, 25)

We can then simulate sweeping of the commutator and filtering input samples across each static leg of the RRC polyphase filter to show the eye opening and closing based on timing correction:

from matplotlib.animation import FuncAnimation

fig, ax = plt.subplots()
line, = ax.plot([], [], ".", alpha=0.3)
plt.axvline(x=0, color="orange")
plt.axhline(y=0, color="orange")
plt.margins(x=0)
plt.grid(True, linestyle="--")
plt.minorticks_on()
plt.tick_params(labelsize=8)
plt.xlabel("In-Phase (I)", fontsize=12)
plt.ylabel("Quadrature (Q)", fontsize=12)
ax.set_yticklabels([])
ax.set_xticklabels([])
iq_vhlim = 0.04
ax.set_xlim([-iq_vhlim, iq_vhlim])
ax.set_ylim([-iq_vhlim, iq_vhlim])

def update_plot(frame):
    poly_match_out = signal.lfilter(h_poly_rrc[frame], 1, samples_rx)
    # NOTE: perform final decimation-by-2 here to bring output to 1x SPS
    line.set_xdata(np.real(poly_match_out[::2]))
    line.set_ydata(np.imag(poly_match_out[::2]))
    ax.set_title(f"Polyphase RRC Leg: {frame}")
    return line,

anim = FuncAnimation(fig, update_plot, frames=num_filters, interval=100, blit=True, repeat=True)
anim.save('iq_polyphase_timing.gif', writer='pillow')
plt.close()

TED S-Curves: Measuring Detector Gain

Sweeping the commutator through each static leg and averaging the TED output over many symbols yields the detector’s S-curve: mean timing error versus leg index. This measurement provides the two numbers a proper loop design needs: the stable lock point (the zero crossing with falling slope), and the detector gain $K_d$ , the S-curve slope at that crossing, in error units per leg, which calibrates the PI loop filter in the timing loops below. Since this bank only spans $\pm \frac{1}{2}$ of an input sample of delay, both S-curves below are nearly linear across all legs, making the slope estimate robust.

One trap when measuring a Gardner S-curve at 2x SPS: the Gardner S-curve is antiperiodic in half a symbol, $S(\tau + T/2) = -S(\tau)$ , so if the prompt sample steps through both 2x-SPS sample parities, the two contributions cancel to numerical zero (measured: the parities correlate at exactly $-1.0$ , leaving only $\sim 10^{-9}$ floating-point residue that plots as pure noise). The prompt must stay on one parity. The sweep below steps the prompt index by 2:

s_curve = []
for i in range(num_filters):
    poly_match_out = signal.lfilter(h_poly_rrc[i], 1, samples_rx)
    temp = 0.0
    # Gardner's S-curve is antiperiodic in T/2, so summing the TED with the
    # prompt on BOTH 2x-SPS parities cancels it to ~0. Keep the prompt on one
    # parity (the symbol strobes) by stepping j by 2:
    for j in range(2, len(poly_match_out) - 1, 2):
        temp += gardner_ted(poly_match_out[j-1], poly_match_out[j], poly_match_out[j+1])
    temp /= len(poly_match_out) // 2
    s_curve.append(temp)
plt.plot(s_curve)
plt.show()

png

s_curve = []
for i in range(num_filters):
    poly_match_out = signal.lfilter(h_poly_rrc[i], 1, samples_rx)
    poly_match_out = poly_match_out[::2]
    temp = 0.0
    for j in range(1, len(poly_match_out)):
        temp += mm_ted(poly_match_out[j], poly_match_out[j-1])
    temp /= len(poly_match_out)
    s_curve.append(temp)
plt.plot(s_curve)
plt.show()

png

Closing the Loop: Type-2 Tracking of the Commutator Index

The tempting way to wire this up — TED error into the PI loop filter, then comm_idx = round(pi_out) % num_filters — is a trap that costs orders of magnitude in performance. The PI loop filter derives its constants

K_p = \frac{2 \zeta \omega_n}{K_d}, \qquad K_i = \frac{\omega_n^2}{K_d}

assuming the element it drives contributes a pole at $z = 1$ (an NCO/accumulator — see the derivation note in pi_filter.py itself), forming the classic type-2 digital PLL with open-loop transfer function:

G_{OL}(z) = K_d \underbrace{\left( K_p + \frac{K_i}{1 - z^{-1}} \right)}_{\text{PI loop filter}} \underbrace{\left(\frac{1}{1 - z^{-1}}\right)}_{\text{NCO}}

Using the PI accumulator itself as the timing phase deletes the NCO’s pole, leaving a type-1 loop with two damaging consequences:

Convergence degrades from $\propto 1/\omega_n$ to $\propto 1/\omega_n^2$ . In the type-1 arrangement, a static timing offset must be held on the PI integrator, whose pull-in pole sits at $\approx 1 - \omega_n^2$ per symbol, versus the type-2 loop’s dominant poles at $\approx 1 - \zeta\omega_n$ . Measured with identical detector gain (DMF TED): halving the loop bandwidth quadruples the type-1 loop’s convergence time ( $511 \rightarrow 2166 \rightarrow$ off the end of a 6000-symbol record) but only doubles the type-2 loop’s ( $51 \rightarrow 110 \rightarrow 248$ symbols). This is exactly the “lower the bandwidth to quiet the leg index, and now it never locks” wall.
The proportional path dumps raw TED noise onto the commutator. The type-1 output is $\text{acc} + K_p e[n]$ , so every strobe the leg index gets kicked by $K_p e[n]$ before rounding. With the gains this loop needs ( $K_p \sim 10^3$ ), even small per-strobe TED self-noise becomes $\gtrsim 1$ leg of index jitter that no downstream averaging can remove. In the type-2 loop that same product is merely an NCO increment — integrated, hence low-pass filtered — so the index walks smoothly.

The fix is structural, not a tuning knob: leave the PI accumulator at zero and treat its output $v[n]$ as a timing frequency in legs/symbol (when tracking a sample-clock offset, $v$ settles at that offset and reads it out directly), then integrate it in a separate arm-phase accumulator that is the commutator:

\phi[n+1] = \phi[n] + v[n], \qquad k[n] = \mathrm{round}(\phi[n])

Three more details matter to get full performance from an implementation:

Calibrate $K_d$ from the measured S-curve (the slope at the stable crossing, in error units per leg, from the sweep above). A mis-set $K_d$ scales both loop constants: a value $81\times$ too small (as this page previously used) runs the loop $81\times$ hotter than the designed bandwidth and pushes the proportional feedthrough $K_p g$ toward the $-1$ sign-flip instability boundary.
A wrap of the arm phase must slip the symbol strobe by one input sample. The bank spans exactly one input sample of delay, so wrapping past leg $M-1$ or leg $0$ means total timing crossed a sample boundary; applying modulo alone snaps timing back by a full sample (half a symbol at 2x SPS) and the loop falls out of lock under any clock-frequency offset. Instead, stretch or shrink the strobe countdown by one input sample for that one strobe ( $L \pm 1$ ). Notably, the delay line and dot product need no special handling for ordinary leg changes — the line always shifts once per input sample, and each strobe’s MAC is a consistent snapshot at the currently-selected delay. Also note the matched filter output is only needed at the symbol strobes (1x SPS), so computing it every input sample doubles the MAC rate for nothing.
Hold the loop off until the delay line fills (and start the commutator from the seeded arm phase), so the first several strobes don’t integrate garbage error into the accumulators.

With these changes, the M&M loop locks in $\approx 27$ symbols at frac_loop_bw=0.02 — versus $\approx 1400$ symbols previously at a nominal frac_loop_bw=0.001 — using the exact same PiFilter class.

Measuring Constellation Quality with EVM

To judge each loop’s output we use Error Vector Magnitude (EVM): the RMS distance from each received symbol to its ideal constellation point, expressed as a percentage of the reference RMS. The wrinkle here is that the RRC $\cdot$ RRC polyphase matched-filter chain applies an unknown, signal-dependent gain (unit-sum leg coefficients, decimation, the fractional-delay filter), so the ideal QPSK points do not sit at a fixed, known magnitude — hard-coding $\pm 1 \pm 1j$ would fold that gain error into the EVM.

Instead we estimate the ideal amplitude from the data. For QPSK the ideal points are

A(\pm 1 \pm 1j)

, and since $$\mathbb{E}[\,

\,] = \mathbb{E}[\,

\,] = A

, the mean of the rectified I and Q recovers

directly. Each sample is then sliced to its nearest ideal point

A\,(\text{sign}\,I + j\,\text{sign}\,Q)

to form the per-symbol reference, and we hand both to `rfproto`'s RMS-normalized `measurements.EVM`. Estimating

A$$ this way makes the metric scale-invariant — whatever gain the chain applies, the reported EVM is unchanged — which is exactly what we want when comparing loops whose outputs live at different amplitudes:

def qpsk_evm(iq_samples: np.ndarray, transient: int = 0) -> float:
    """Error Vector Magnitude (%) of QPSK symbols vs a dynamically estimated
    ideal constellation.

    The RRC + polyphase matched-filter chain applies an unknown, signal-dependent
    gain, so the ideal QPSK points are not at a fixed magnitude. We recover the
    per-axis amplitude A from the data itself: the ideal points are A*(+/-1 +/- 1j),
    where mean(|I|) = mean(|Q|) = A. Each sample is then sliced to its nearest
    ideal point A*(sign(I) + j*sign(Q)) to form the reference, and rfproto's
    RMS-normalized EVM reports the error as a percentage. Scale-invariant by
    construction, since A tracks whatever gain the chain applies.

    Args:
        iq_samples: complex symbol-rate (1x SPS) samples to measure
        transient:  number of leading symbols to skip (e.g. loop acquisition)

    Returns:
        EVM in percent
    """
    iq = np.asarray(iq_samples)[transient:]
    # dynamic amplitude estimate: mean(|I|) == mean(|Q|) == A for ideal QPSK
    A = 0.5 * (np.mean(np.abs(iq.real)) + np.mean(np.abs(iq.imag)))
    ideal = A * (np.sign(iq.real) + 1j * np.sign(iq.imag))
    return 100.0 * measurements.EVM(iq, ideal)

# Calibrate the loop from the measured S-curve above: the PI filter's
# detector_gain must equal the actual TED slope (error units per polyphase arm)
# at the stable zero crossing for the designed bandwidth/damping to be realized.
s = np.asarray(s_curve)
k0 = np.where((s[:-1] > 0) & (s[1:] <= 0))[0][0]  # falling = stable crossing
mm_ted_gain = -np.polyfit(np.arange(k0 - 2, k0 + 3), s[k0 - 2 : k0 + 3], 1)[0]
print(f"Measured M&M TED gain: {mm_ted_gain:.3e} err/arm (crossing near arm {k0})")

in_delay = np.zeros(taps, dtype=complex)
prev = 0.0 + 0.0j
filt_out, ted_out, freq_out, arm_out = [], [], [], []

# PiFilter derives Kp/Ki assuming its output drives a plant with a pole at z=1
# (an NCO). Using the PI accumulator itself as the arm index makes a type-1 loop
# whose settle time grows as 1/BW^2 and feeds Kp*err jitter straight into the
# commutator. Instead, treat the PI output as a timing *frequency* (arms/symbol)
# and integrate it in a separate arm-phase accumulator (the NCO) -> type-2 loop.
loop_filt = pi_filter.PiFilter(frac_loop_bw=0.02, detector_gain=mm_ted_gain)
arm_phase = float(num_filters // 2)  # NCO state: fractional arm index
comm_idx = round(arm_phase)
strobe_count = L - 1                 # first symbol strobe lands on sample 0
strobe_period = L                    # L +/- 1 when the arm phase wraps

for idx, x in enumerate(samples_rx):
    # shift newest sample into the delay line (always at the input rate)
    in_delay[1:] = in_delay[:-1]
    in_delay[0] = x

    strobe_count += 1
    if strobe_count < strobe_period:
        continue
    # symbol strobe: the MF output is only needed at 1x SPS, halving the MACs
    strobe_count, strobe_period = 0, L
    mac = np.dot(in_delay, h_poly_rrc[comm_idx])
    filt_out.append(mac)
    if idx < taps:
        continue  # hold the loop off until the delay line has filled

    err = mm_ted(mac, prev)
    prev = mac
    ted_out.append(err)

    v = loop_filt.Step(err)  # timing frequency (arms per symbol)
    freq_out.append(v)
    arm_phase += v           # NCO integrates frequency -> timing phase

    # A wrap of the arm phase means the total timing crossed one full input
    # sample, so slip the next strobe by a sample to keep delay continuous
    if arm_phase > num_filters - 0.5:
        arm_phase -= num_filters
        strobe_period = L + 1  # stuff: strobe one sample later
    elif arm_phase < -0.5:
        arm_phase += num_filters
        strobe_period = L - 1  # skip: strobe one sample earlier
    comm_idx = round(arm_phase)
    arm_out.append(arm_phase)

plot.IQ(np.asarray(filt_out)[100:], alpha=0.4, title="Timing Correction Output")
plt.show()

plot.samples(ted_out, title="TED Output")
plt.show()

plot.samples(freq_out, title="PI Filter Output (timing frequency, arms/symbol)")
plt.show()

plot.samples(arm_out, title="Arm-Phase NCO (selected polyphase leg)")
plt.show()

# skip the acquisition transient, then measure output constellation quality
print(f"M&M timing loop output EVM: {qpsk_evm(filt_out, transient=100):.2f}%")

Measured M&M TED gain: 8.190e-04 err/arm (crossing near arm 24)

png

M&M timing loop output EVM: 0.89%

The residual $\pm 1$ leg chatter visible at lock is a quantization limit cycle, not loop noise: the true optimum sits between legs (here at $\phi \approx 24.3$ ), so a nearest-leg loop must forever hunt between legs 24 and 25 — the integrator walks $\phi$ across the rounding boundary, the selected leg flips, the error sign flips, and it walks back. When that matters, blend the two adjacent legs’ coefficients by the fractional arm phase, replacing the single mac dot product above:

kf = int(np.floor(arm_phase)) % num_filters
frac = arm_phase - np.floor(arm_phase)
h_blend = (1.0 - frac) * h_poly_rrc[kf] + frac * h_poly_rrc[(kf + 1) % num_filters]
mac = np.dot(in_delay, h_blend)

This makes the applied delay continuous rather than quantized to $1/M$ sample steps: measured on this signal it eliminates leg toggling entirely, settles $\phi$ at exactly the S-curve zero crossing (24.30), and drops the output EVM from 0.90% to 0.31% — essentially the open-loop floor of an ideally-placed fractional delay (0.28% here). It also answers “how many polyphase legs are enough?”: with blending, the leg count stops mattering once the S-curve is linear between adjacent legs. (One caveat: at a wrap boundary, leg $M-1$ blends with leg $0$ of the adjacent input sample, a transient $1/M$ -sample error for a single strobe — irrelevant once locked mid-bank.)

Harris Polyphase Derivative Matched Filter TED

The M&M loop above derives its timing error from the shape of a single matched-filter output stream. The classic fred harris approach is more direct: run a matched filter (MF) polyphase bank alongside a derivative matched filter (DMF) bank, with both banks driven by the same commutator / arm-select index. At the maximum-likelihood (peak) sampling instant the MF output is maximized, so its derivative — the DMF output — passes through zero. The DMF output, sign-corrected by the current symbol decision, is therefore a direct, near-linear estimate of the timing error:

e[n] = \text{sign}\big(\Re\{y_{MF}[n]\}\big)\,\Re\{y_{DMF}[n]\} + \text{sign}\big(\Im\{y_{MF}[n]\}\big)\,\Im\{y_{DMF}[n]\}

We build the DMF bank by taking the central-difference derivative of the same RRC polyphase prototype (np.gradient) and partitioning it into the same $M = 31$ arms — see the Derivative Matched Filter section of the Pulse-Shaping notebook for why the cheap finite-difference derivative is indistinguishable from an ideal $j\omega$ differentiator here (the prototype is heavily oversampled by its polyphase arms, so the signal band sits far below Nyquist). Both banks convolve the same input delay line, indexed by the same arm select, and the loop uses the same type-2 structure as the M&M loop above: the measured DMF S-curve slope calibrates $K_d$ , the PI output (timing frequency) is integrated by the arm-phase NCO, and arm wraps slip the symbol strobe — so unlike a mod-only implementation, this loop genuinely tracks a sample-clock-frequency offset, slipping the integer sample index at each wrap while the PI output reads out the clock offset in arms/symbol.

One tradeoff to respect: at the symbol strobes the composite MF cascade (RRC $\cdot$ RRC = raised cosine) is Nyquist, so the M&M TED’s decision-directed error self-cancels there — but the derivative of the composite pulse is not zero at the neighboring symbol instants, so every DMF strobe carries data-dependent self-noise (measured here at $\sim 70\times$ the M&M TED’s, input-referred to arms of timing). Self-noise cannot be removed by the loop filter, only averaged, so this loop acquires at a moderate bandwidth and then gear-shifts down ( $0.01 \rightarrow 0.003$ ) once locked — the standard acquire-wide/track-narrow modem practice. The DMF’s compensating virtue is an S-curve that stays linear across the entire bank, giving wide, hang-up-free pull-in; “acquire with the DMF, track with M&M” is a legitimate hybrid in demanding designs.

# --- Derivative matched-filter (DMF) polyphase bank ---
dproto_rrc = np.gradient(polyphase_rrc_coef)
h_poly_drrc = dproto_rrc.reshape(len(dproto_rrc) // num_filters, num_filters).T


def dmf_ted(y_mf, y_dmf):
    return np.sign(np.real(y_mf)) * np.real(y_dmf) + np.sign(np.imag(y_mf)) * np.imag(y_dmf)


# Measure the DMF S-curve the same way as the M&M one above and take its slope
# at the stable crossing as the loop's detector gain
dmf_s_curve = []
for i in range(num_filters):
    y = signal.lfilter(h_poly_rrc[i], 1, samples_rx)[::2]
    yd = signal.lfilter(h_poly_drrc[i], 1, samples_rx)[::2]
    dmf_s_curve.append(np.mean(dmf_ted(y, yd)))
sd = np.asarray(dmf_s_curve)
k0 = np.where((sd[:-1] > 0) & (sd[1:] <= 0))[0][0]
dmf_ted_gain = -np.polyfit(np.arange(k0 - 2, k0 + 3), sd[k0 - 2 : k0 + 3], 1)[0]
print(f"Measured DMF TED gain: {dmf_ted_gain:.3e} err/arm (crossing near arm {k0})")

in_delay = np.zeros(taps, dtype=complex)
mf_out, ted_out, freq_out, arm_out = [], [], [], []

# Same type-2 structure as the M&M loop. The DMF TED's per-strobe self-noise is
# ~70x the M&M's (input-referred), so acquire at a moderate bandwidth then
# gear-shift down once locked to keep the arm index quiet.
loop_filt = pi_filter.PiFilter(frac_loop_bw=0.01, detector_gain=dmf_ted_gain)
arm_phase = float(num_filters // 2)
comm_idx = round(arm_phase)
strobe_count = L - 1
strobe_period = L

for idx, x in enumerate(samples_rx):
    in_delay[1:] = in_delay[:-1]
    in_delay[0] = x

    strobe_count += 1
    if strobe_count < strobe_period:
        continue
    strobe_count, strobe_period = 0, L
    y_mf = np.dot(in_delay, h_poly_rrc[comm_idx])
    mf_out.append(y_mf)
    if idx < taps:
        continue

    err = dmf_ted(y_mf, np.dot(in_delay, h_poly_drrc[comm_idx]))
    ted_out.append(err)
    if len(ted_out) == 100:
        loop_filt.frac_BW = 0.003  # gear-shift: acquire wide, track narrow

    v = loop_filt.Step(err)
    freq_out.append(v)
    arm_phase += v
    if arm_phase > num_filters - 0.5:
        arm_phase -= num_filters
        strobe_period = L + 1
    elif arm_phase < -0.5:
        arm_phase += num_filters
        strobe_period = L - 1
    comm_idx = round(arm_phase)
    arm_out.append(arm_phase)

sym_out = np.asarray(mf_out)
plot.IQ(sym_out[250:], alpha=0.4, title="Harris PFB TED: Timing Correction Output")
plt.show()

plot.samples(ted_out, title="Derivative Matched Filter TED Output")
plt.show()

plot.samples(freq_out, title="PI Filter Output (timing frequency, arms/symbol)")
plt.show()

plot.samples(arm_out, title="Arm-Phase NCO (selected polyphase leg)")
plt.show()

# skip the (longer) DMF acquisition transient, then measure constellation quality
print(f"DMF timing loop output EVM: {qpsk_evm(sym_out, transient=250):.2f}%")