Claude Fable5

pd_pfd_emulation.py

"""
Emulation: average output of a phase detector (XOR) vs a phase-frequency
detector (tri-state PFD) when the two inputs differ in frequency.

Setup (open loop, no feedback):
    input A = reference square wave at f_a = 1
    input B = "VCO" square wave at f_b = f_a - df

Claims demonstrated:
  1. XOR PD: time-average -> 0 for every df != 0, identical for +df and -df
     -> the average carries no frequency information.
  2. Tri-state PFD: steady-state average = sign(df) * (1 - f_min / (2 f_max)),
     so |avg| >= 1/2 and sign(avg) == sign(df) always.
"""
import numpy as np
import matplotlib
# matplotlib.use("Agg")
import matplotlib.pyplot as plt

F_REF = 1.0
OUT = "./"

# ----------------------- tri-state PFD, event driven -----------------------
def pfd_run(fa, fb, t_end, pha=0.31831, phb=0.123456):
    """Ideal tri-state PFD driven by rising edges.
    state: +1 = UP, 0 = idle, -1 = DN.
    A edge: state = min(state+1, +1);  B edge: state = max(state-1, -1).
    Returns piecewise-constant trajectory (seg_t, seg_s): state on
    [seg_t[i], seg_t[i+1]) is seg_s[i]."""
    ea = np.arange(pha / fa, t_end, 1.0 / fa)
    eb = np.arange(phb / fb, t_end, 1.0 / fb)
    times = np.concatenate([ea, eb])
    kind = np.concatenate([np.ones(ea.size, np.int8), -np.ones(eb.size, np.int8)])
    o = np.argsort(times, kind="stable")
    times, kind = times[o], kind[o]
    seg_t = np.empty(times.size + 2)
    seg_s = np.empty(times.size + 2, np.int8)
    seg_t[0], seg_s[0] = 0.0, 0
    s = 0
    for i in range(times.size):
        s = min(s + 1, 1) if kind[i] > 0 else max(s - 1, -1)
        seg_t[i + 1], seg_s[i + 1] = times[i], s
    seg_t[-1], seg_s[-1] = t_end, s
    return seg_t, seg_s

def piecewise_avg(seg_t, seg_s, t0, t1):
    """Time-average of the piecewise-constant state over [t0, t1]."""
    t = np.clip(seg_t, t0, t1)
    return float(np.sum(seg_s[:-1] * np.diff(t)) / (t1 - t0))

def pfd_steady_avg(fa, fb, t_end, n_phase=8):
    """Steady-state average, ensemble-averaged over the initial B-edge
    alignment (theory's <gap> = T_A/2 assumes equidistributed alignments;
    rational f_a/f_b makes a single run alignment-dependent)."""
    vals = []
    for phb in np.linspace(0.04, 0.96, n_phase):
        st, ss = pfd_run(fa, fb, t_end, phb=phb)
        vals.append(piecewise_avg(st, ss, t_end / 2, t_end))  # discard transient
    return float(np.mean(vals))

def pfd_theory(fa, fb):
    if fa == fb:
        return 0.0
    hi, lo = max(fa, fb), min(fa, fb)
    return float(np.sign(fa - fb) * (1.0 - lo / (2.0 * hi)))

def running_avg_curve(seg_t, seg_s):
    cum = np.cumsum(seg_s[:-1] * np.diff(seg_t))
    return seg_t[1:], cum / seg_t[1:]

# ----------------------------- XOR PD, sampled -----------------------------
def xor_out(fa, fb, t, pha=0.7):
    """Bipolar XOR of two 50%-duty square waves: +1 when levels differ."""
    return -np.sign(np.sin(2 * np.pi * fa * t + pha)) * np.sign(np.sin(2 * np.pi * fb * t))

def xor_avg(fa, fb, df):
    n_beats = int(np.clip(4000 * abs(df), 10, 400))   # integer # of beat periods
    t_end = n_beats / abs(df)
    dt = 1.0 / (64 * max(fa, fb))
    t = np.arange(0.0, t_end, dt)
    return float(xor_out(fa, fb, t).mean())

# ------------------------------- experiments -------------------------------
if __name__ == "__main__":
    # ---- sweep df and record averages ----
    pos = np.concatenate([np.arange(0.005, 0.05, 0.005), np.arange(0.05, 0.301, 0.025)])
    deltas = np.concatenate([-pos[::-1], pos])
    xor_pts, pfd_pts, thr = [], [], []
    for df in deltas:
        fa, fb = F_REF, F_REF - df
        T = max(40.0 / abs(df), 1000.0)
        xor_pts.append(xor_avg(fa, fb, df))
        pfd_pts.append(pfd_steady_avg(fa, fb, T))
        thr.append(pfd_theory(fa, fb))

    print(f"{'df/f_ref':>9} | {'XOR PD avg':>11} | {'PFD avg (emul)':>14} | {'PFD theory':>10}")
    print("-" * 56)
    for df in [0.02, 0.05, 0.20, -0.02, -0.05, -0.20]:
        fa, fb = F_REF, F_REF - df
        T = max(40.0 / abs(df), 1000.0)
        print(f"{df:>+9.3f} | {xor_avg(fa, fb, df):>+11.5f} | "
              f"{pfd_steady_avg(fa, fb, T):>+14.4f} | {pfd_theory(fa, fb):>+10.4f}")

    # ---- figure 1: average output vs frequency offset ----
    fig, ax = plt.subplots(figsize=(7.2, 4.6))
    ax.axhline(0, color="0.65", lw=0.8)
    ax.axvline(0, color="0.85", lw=0.8)
    ax.plot(deltas, thr, "k--", lw=1.1,
            label=r"PFD theory: $\pm\,(1 - f_{min}/2f_{max})$")
    ax.plot(deltas, pfd_pts, "o", ms=5.5, color="#1D9E75",
            label="PFD, emulated (steady state)")
    ax.plot(deltas, xor_pts, "s", ms=4.5, color="#7F77DD",
            label="XOR PD, emulated")
    ax.annotate("sign(avg) = sign($\\Delta f$),\n|avg| $\\geq$ 1/2",
                xy=(0.12, 0.62), fontsize=9, color="#0F6E56")
    ax.annotate("avg $\\approx$ 0 for every $\\Delta f$\n(same for $\\pm\\Delta f$)",
                xy=(-0.28, 0.10), fontsize=9, color="#534AB7")
    ax.set_xlabel(r"frequency offset  $\Delta f / f_{ref}$   ($f_B = f_{ref} - \Delta f$)")
    ax.set_ylabel("time-averaged output  (UP $-$ DN, normalized)")
    ax.set_title("Average detector output under a pure frequency offset (open loop)")
    ax.set_ylim(-1.08, 1.08)
    ax.legend(loc="lower right", fontsize=9)
    fig.tight_layout()
    fig.savefig(f"{OUT}/fig1_avg_output_vs_frequency_offset.png", dpi=160)

    # ---- figure 2: time domain at df = +5% ----
    df = 0.05
    fa, fb = F_REF, F_REF - df
    st, ss = pfd_run(fa, fb, 1200.0)
    dn_times = st[1:][ss[1:] == -1]
    n_dn = dn_times.size
    last_dn = dn_times.max() if n_dn else float("nan")
    print(f"\ndf = +0.05:  PFD DN pulses emitted = {n_dn}, "
          f"last DN pulse at t = {last_dn:.2f}  (beat period 1/df = {1/df:.0f})")
    print(f"             PFD steady-state avg = {piecewise_avg(st, ss, 600, 1200):+.4f}"
          f"  vs theory {pfd_theory(fa, fb):+.4f}")

    fig, axs = plt.subplots(2, 1, figsize=(7.2, 5.6),
                            gridspec_kw={"height_ratios": [1, 1.5]})
    axs[0].step(st, ss, where="post", color="#0F6E56", lw=1.1)
    axs[0].set_xlim(0, 60)
    axs[0].set_ylim(-1.4, 1.4)
    axs[0].set_yticks([-1, 0, 1], labels=["DN", "idle", "UP"])
    axs[0].axvline(last_dn, color="#E24B4A", lw=0.9, ls=":")
    axs[0].annotate("last DN pulse ever", xy=(last_dn + 1, -1.1),
                    fontsize=9, color="#A32D2D")
    axs[0].set_title(r"PFD state vs time,  $\Delta f = +5\%$:  DN dies out within ~one beat period")
    axs[0].set_xlabel("time  (periods of $f_{ref}$)")

    tr, vr = running_avg_curve(st, ss)
    dt = 1.0 / (64 * fa)
    tt = np.arange(dt, 600.0, dt)
    vx = np.cumsum(xor_out(fa, fb, tt)) * dt / tt
    axs[1].axhline(0, color="0.65", lw=0.8)
    axs[1].axhline(pfd_theory(fa, fb), color="k", ls="--", lw=0.9)
    axs[1].annotate(rf"theory: $1 - f_B/2f_A = {pfd_theory(fa, fb):.3f}$",
                    xy=(330, pfd_theory(fa, fb) + 0.04), fontsize=9)
    axs[1].plot(tr, vr, color="#1D9E75", lw=1.4, label="PFD running average")
    axs[1].plot(tt, vx, color="#7F77DD", lw=1.4, label="XOR PD running average")
    axs[1].set_xlim(0, 600)
    axs[1].set_ylim(-0.6, 0.8)
    axs[1].set_xlabel("averaging window length  (periods of $f_{ref}$)")
    axs[1].set_ylabel("running average")
    axs[1].legend(loc="center right", fontsize=9)
    fig.tight_layout()
    fig.savefig(f"{OUT}/fig2_time_domain_df_plus5pct.png", dpi=160)
    print("\nfigures written to", OUT)
    plt.show()