Phase-locked loop (digital)

Deep dive

How it works

In: complex or real sample stream x[n], loop bandwidth Bn*T, damping · Out: locked phase and frequency estimate, derotated signal · State: NCO phase and the loop filter integrator · Cost: O(1) per sample: one phase detector, two multiply-adds, one sin/cos

 x[n] ──►(×)──► y[n] ──► phase detector ──► e[n] ──► loop filter ──► NCO ──┐
          ▲                                          alpha*e + integral    │
          └───────────────── e^{-j phase[n]} ◄───────────────────────────┘

1. Derotate: y = x * exp(-j*phase).
2. Detect: carrier present, e = atan2(Im y, Re y). Suppressed BPSK carrier (Costas): e = sign(Re y) * Im y. QPSK: Re*sign(Im) - Im*sign(Re) (not in the code).
3. Filter: freq += beta*e (integrator), then phase += freq + alpha*e.
4. Gains from the design bandwidth: theta = BnT / (zeta + 1/(4*zeta)), d = 1 + 2*zeta*theta + theta^2, alpha = 4*zeta*theta/d, beta = 4*theta^2/d. This assumes detector gain 1 and NCO gain 1 (rad per unit); rescale alpha, beta by 1/(Kd*K0) otherwise (the Costas error above has gain about the signal power, so normalize amplitude or divide).
5. Locked when |e| stays small; freq is then the offset estimate in rad/sample.

Code

use std::f64::consts::PI;
pub type C = (f64, f64); // (re, im)

/// Second-order digital PLL: PI loop filter + NCO (phase accumulator).
pub struct Pll { phase: f64, freq: f64, alpha: f64, beta: f64 }

impl Pll {
    /// `bn_t`: noise bandwidth x sample period (Bn/fs), e.g. 0.01. `zeta`: damping, 0.707.
    pub fn new(bn_t: f64, zeta: f64) -> Self {
        let theta = bn_t / (zeta + 1.0 / (4.0 * zeta));
        let d = 1.0 + 2.0 * zeta * theta + theta * theta;
        // Unit detector gain and unit NCO gain assumed (atan2 detector).
        Pll { phase: 0.0, freq: 0.0, alpha: 4.0 * zeta * theta / d, beta: 4.0 * theta * theta / d }
    }

    fn derotate(&self, x: C) -> C {
        let (s, c) = (-self.phase).sin_cos();            // NCO output e^{-j phase}
        (x.0 * c - x.1 * s, x.0 * s + x.1 * c)
    }

    /// One sample with a carrier present; returns (derotated sample, phase error).
    pub fn step(&mut self, x: C) -> (C, f64) {
        let y = self.derotate(x);
        let err = y.1.atan2(y.0);
        self.update(err);
        (y, err)
    }

    /// Costas variant for BPSK (suppressed carrier): 180-degree ambiguity.
    pub fn step_costas(&mut self, x: C) -> (C, f64) {
        let y = self.derotate(x);
        let err = y.0.signum() * y.1;
        self.update(err);
        (y, err)
    }

    fn update(&mut self, err: f64) {
        self.freq += self.beta * err;                    // integrator: tracks frequency offset
        self.phase += self.freq + self.alpha * err;      // proportional: pulls phase
        if self.phase >= PI { self.phase -= 2.0 * PI } else if self.phase < -PI { self.phase += 2.0 * PI }
    }
    pub fn freq(&self) -> f64 { self.freq }   // rad/sample
    pub fn phase(&self) -> f64 { self.phase } // phase used for the NEXT sample
}

Run with asserts: a tone at 0.05 rad/sample and 1 rad initial phase locked (error below 0.01 rad after 206 samples at BnT=0.01) and the frequency estimate matched to 1e-6; a frequency ramp of 1e-6 rad/sample^2 gave a constant error of 0.0028 rad, equal to r/beta; the Costas loop on random BPSK settled with either all or no decisions inverted.

Complexity and constants

• Per sample: one complex multiply, one detector, 2 multiply-adds, one sin_cos. State is 2 floats; no memory traffic.
• Measured here (scalar f64, one Apple-silicon core, pre-generated input): about 31 ns/sample with atan2 plus sin_cos, about 20 ns with the Costas sign*Q detector. The trigonometry dominates, not the loop arithmetic.
• Noise: for input amplitude A and total complex noise power sigma2, phase variance is about sigma2/(2A^2) * 2*BnT. Measured at 20 dB SNR: 5.4e-5, 1.05e-4, 2.08e-4 rad^2 at BnT = 0.005, 0.01, 0.02, against theory 5e-5, 1e-4, 2e-4.
• Acquisition, measured at BnT=0.01, initial phase 1 rad: locked in about 300 samples for offsets up to 0.02 rad/sample, 1,732 samples at 0.3, and not at 0.5 or above within 200k samples. At BnT=0.001: 2,220 samples at 0.01, 39k at 0.05, 162k at 0.1, no lock at 0.2.

Efficiency

• Bound: compute-bound on transcendental functions, not memory. State fits in registers, so the cost is arithmetic per sample.
• Cheaper detector: replacing atan2 with Im(y)/|y| or just Im(y) (error is proportional to the sine of the phase error, near-linear for small errors) saves about a third of the per-sample time in the measurement above (31 vs 20 ns, Costas detector); gain then depends on amplitude, so normalize or run an AGC ahead of it.
• Cheaper NCO: a rotation by e^{j*freq} per sample plus occasional renormalization, or a sine lookup table with a u32 phase accumulator, replaces sin_cos, which is most of the roughly 20 ns per sample measured above.
• Narrow loop = less noise, more work to acquire: jitter variance scales with BnT; pick the narrowest bandwidth that still tracks your worst-case frequency drift, then use a wider loop only during acquisition.
• Decimate first: run the loop after filtering and decimation to the lowest sufficient rate; work scales with samples per second.
• Cheaper alternative: for a single burst with a known preamble, a feed-forward estimate (angle of a correlation sum) needs no loop and no acquisition time.

Implementation notes

• Keep the phase wrapped in [-pi, pi) (floats) or use a wrapping u32 accumulator where 2^32 = one turn: wrap is free, and a frequency word of f/fs * 2^32 gives exact NCO frequency resolution fs/2^32.
• Order matters: update freq first, then phase (as above), or the loop gets an extra sample of delay and the damping changes.
• Clamp freq to the plausible offset range so a noise burst cannot send the NCO to the wrong lobe.
• Lock detector: low-pass |e| or the ratio Re(y)^2 / (Re(y)^2 + Im(y)^2) and compare with a threshold; widen the loop while unlocked, narrow it once locked.
• Damping 0.707 is the conventional compromise between settling time and overshoot; larger values overshoot less and settle slower.
• Keep BnT small (about 0.05 or less): the gain formulas approximate a continuous-time loop, so verify by simulation at larger values.

Pitfalls

• Wrong gain scaling: detector gain depends on amplitude (Costas) or modulation; alpha, beta are then off by that factor.
• Pulling in a large offset: the atan2 detector is a sawtooth, so large offsets cause cycle slips and slow or failed acquisition (see the measured numbers above); use a coarse estimate first.
• Cycle slips and ambiguity: Costas/PLL on suppressed-carrier signals lock at multiples of 180 or 90 degrees; resolve with a unique word or differential encoding.
• Using phase() after the update as the phase of the sample just processed: it is the phase for the next sample, which biases error measurements by one freq.
• Passing the loop bandwidth in Hz where Bn/fs is expected.
• Modulated data in the detector (data-aided vs non-data-aided): non-data-aided detectors add self-noise at low SNR.

Variants and alternatives

In production

Software-defined radio receivers (GNU Radio, liquid-dsp), GNSS software receivers (carrier and code tracking), digital TV and satellite modems (carrier recovery), PLLs inside every radio and clock chip (analogue and fractional-N). Power electronics: grid-connected converters synchronize to the mains phase with a PLL (Kaura and Blasko).

Implementations

Sources

• Floyd M. Gardner, Phaselock Techniques, 3rd ed., Wiley, 2005 (loop theory, pull-in, noise bandwidth).
• Michael Rice, Digital Communications: A Discrete-Time Approach, Pearson, 2009 (discrete-time loop filter design, the theta, alpha, beta formulas).
• John P. Costas, Synchronous communications, Proc. IRE 44(12), 1956. paper
• Vikram Kaura, Vladimir Blasko, Operation of a phase locked loop system under distorted utility conditions, IEEE Trans. Industry Applications 33(1), 1997. paper
• GNU Radio control_loop.cc (alpha/beta formulas). code