Phase-locked loops¶

import matplotlib.pyplot as plt
import numpy as np

import sdr

%config InlineBackend.print_figure_kwargs = {"facecolor" : "w"}
%matplotlib inline
# %matplotlib widget

loop_filter = sdr.LoopFilter(0.05, 1)
print(loop_filter)

plt.figure(figsize=(10, 5))
sdr.plot.magnitude_response(loop_filter.iir, x_axis="log", decades=6)
plt.show()

N = 75  # samples
theta_0 = 2 * np.pi / 10  # radians/sample
x = theta_0 * np.arange(N) + np.pi  # Input phase signal, radians
y = np.zeros(x.size)  # Output phase signal, radians
phase_error = np.zeros(x.size)  # Measured phase error, radians
freq_estimate = np.zeros(x.size + 1)  # Estimated frequency, radians/sample

# Create a proportional-plus-integrator (PPI) loop filter with a normalized
# noise bandwidth of 0.05 and a damping factor of 1 (critically damped)
loop_filter = sdr.LoopFilter(0.05, 1)

# Create a numerically controlled oscillator (NCO) with NCO gain of 1,
# constant phase accumulation of theta_0 radians/sample, and initial phase
# offset of theta_0 radians.
nco = sdr.NCO(1, increment=theta_0)

for i in range(N):
    # Process the variable phase increment through the NCO
    y[i] = nco(freq_estimate[i], output="phase")

    # Phase error detector (PED)
    phase_error[i] = x[i] - y[i]

    # Compute the frequency estimate
    freq_estimate[i + 1] = loop_filter(phase_error[i])

plt.figure(figsize=(10, 6))

plt.subplot(2, 1, 1)
sdr.plot.time_domain(x, label="Input phase")
sdr.plot.time_domain(y, label="Output phase")
plt.grid(which="both", linestyle="--")
plt.ylabel("Phase (radians)")
plt.title("Input and output phase signals")

plt.subplot(2, 1, 2)
sdr.plot.time_domain(phase_error)
plt.grid(which="both", linestyle="--")
plt.ylabel("Phase (radians)")
plt.title("Phase error between input and output phase signals")

plt.suptitle("Figure C.2.5 from Digital Communications: A Discrete-Time Approach")
plt.tight_layout()
plt.show()

N = 75  # samples
theta_0 = 2 * np.pi / 10  # radians/sample
x = np.exp(1j * (theta_0 * np.arange(N) + np.pi))  # Input signal
y = np.ones(x.size, dtype=complex)  # Output signal
phase_error = np.zeros(x.size)
freq_estimate = np.zeros(x.size + 1)

# Create a numerically controlled oscillator (NCO) with NCO gain of 1,
# constant phase accumulation of theta_0 radians/sample, and initial phase
# offset of theta_0 radians.
nco = sdr.NCO(1, increment=theta_0)

# Create a heuristic phase error detector (PED).
ped = sdr.PED()

# Create a proportional-plus-integrator (PPI) loop filter with a normalized
# noise bandwidth of 0.05 and a damping factor of 1 (critically damped)
loop_filter = sdr.LoopFilter(0.05, 1, K0=nco.gain, Kp=ped.gain)

for i in range(N):
    # Process the variable phase increment through the DDS
    y[i] = nco(freq_estimate[i])

    # Heuristic phase error detector (PED)
    phase_error[i] = ped(x[i], y[i])

    # Compute the frequency estimate
    freq_estimate[i + 1] = loop_filter(phase_error[i])

plt.figure(figsize=(10, 8))

plt.subplot(3, 1, 1)
sdr.plot.time_domain(x.real, label="Input")
sdr.plot.time_domain(y.real, label="Output")
plt.grid(which="both", linestyle="--")
plt.ylabel("Amplitude")
plt.title("Input and output signals (real part only)")

plt.subplot(3, 1, 2)
sdr.plot.time_domain(phase_error)
plt.grid(which="both", linestyle="--")
plt.ylabel("Phase (radians)")
plt.title("Phase error between input and output signals")

plt.subplot(3, 1, 3)
sdr.plot.time_domain(np.unwrap(np.angle(x)), label="Input")
sdr.plot.time_domain(np.unwrap(np.angle(y)), label="Output")
plt.grid(which="both", linestyle="--")
plt.ylabel("Phase (radians)")
plt.title("Unwrapped phase of input and output signals")

plt.suptitle("Figure C.2.5 from Digital Communications: A Discrete-Time Approach")
plt.tight_layout()
plt.show()

N = 250  # symbols
qpsk = sdr.PSK(4, phase_offset=45)
a = qpsk.map_symbols(np.random.randint(0, 4, N))  # QPSK complex symbols
a *= np.exp(1j * np.pi / 4)  # Rotate symbols by pi/4 radians
lo = np.ones(a.size, dtype=complex)  # Local oscillator
a_tilde = np.ones(a.size, dtype=complex)  # De-rotated symbols

phase_error = np.zeros(a.size)
freq_estimate = np.zeros(a.size + 1)

# Create a numerically controlled oscillator (NCO) with NCO gain of 1.
nco = sdr.NCO(1)

# Create a maximum-likelihood phase error detector (ML-PED).
ped = sdr.MLPED()

# Create a proportional-plus-integrator (PPI) loop filter with a normalized
# noise bandwidth of 0.05 and a damping factor of 1 (critically damped)
loop_filter = sdr.LoopFilter(0.02, 1 / np.sqrt(2), K0=nco.gain, Kp=ped.gain)

for i in range(N):
    # Process the variable phase increment through the DDS
    lo[i] = nco(freq_estimate[i])

    # De-rotate the input symbols
    a_tilde[i] = a[i] * lo[i].conj()

    # Make a symbol decision
    s_hat, a_hat = qpsk.decide_symbols(a_tilde[i])

    # Maximum-likelihood phase error detector (ML-PED)
    phase_error[i] = ped(a_tilde[i], a_hat)

    # Compute the frequency estimate
    freq_estimate[i + 1] = loop_filter(phase_error[i])

plt.figure(figsize=(10, 10))

plt.subplot(4, 1, 1)
sdr.plot.time_domain(a)
plt.title("Input symbols")

plt.subplot(4, 1, 2)
sdr.plot.time_domain(a_tilde)
plt.title("Output symbols (de-rotated)")

plt.subplot(4, 1, 3)
sdr.plot.time_domain(np.angle(lo, deg=True))
plt.grid(which="both", linestyle="--")
plt.ylabel("Phase (degrees)")
plt.title("Phase estimate of local oscillator")

plt.subplot(4, 1, 4)
sdr.plot.time_domain(np.rad2deg(phase_error))
plt.grid(which="both", linestyle="--")
plt.ylabel("Phase (degrees)")
plt.title("Phase error between input and output phase signals")

plt.suptitle("Figure 7.2.6 from Digital Communications: A Discrete-Time Approach")
plt.tight_layout()
plt.show()

plt.figure(figsize=(10, 5))
for zeta in [1 / 2, 1 / np.sqrt(2), 1, np.sqrt(2), 2]:
    pll = sdr.ClosedLoopPLL(0.01, zeta)
    sdr.plot.step_response(pll.iir, N=500, label=rf"$\zeta = {zeta}$")
plt.title("Step response of the closed-loop PLL across damping factor")
plt.show()

plt.figure(figsize=(10, 5))
for zeta in [1 / 2, 1 / np.sqrt(2), 1, np.sqrt(2), 2]:
    pll = sdr.ClosedLoopPLL(0.01, zeta)
    sdr.plot.magnitude_response(pll.iir, sample_rate=2 * np.pi / pll.omega_n, x_axis="log", label=rf"$\zeta = {zeta}$")
plt.xlim([10**-2, 10**2])
plt.ylim([-25, 5])
plt.xlabel(r"Normalized Frequency ($\omega / \omega_n$)")
plt.title("Frequency response of the closed-loop PLL across damping factor")
plt.show()

plt.figure(figsize=(10, 5))
for BnT in [0.001, 0.005, 0.01, 0.05, 0.1]:
    pll = sdr.ClosedLoopPLL(BnT, 1)
    sdr.plot.step_response(pll.iir, N=500, label=f"$B_nT = {BnT}$")
plt.title("Step response of the closed-loop PLL across normalized noise bandwidth")
plt.show()

plt.figure(figsize=(10, 5))
for BnT in [0.001, 0.005, 0.01, 0.05, 0.1]:
    pll = sdr.ClosedLoopPLL(BnT, 1)
    sdr.plot.magnitude_response(pll.iir, x_axis="log", label=f"$B_nT = {BnT}$")
plt.ylim([-25, 5])
plt.title("Frequency response of the closed-loop PLL across normalized noise bandwidth")
plt.show()

freq = np.linspace(0, 0.0005, 100)

plt.figure(figsize=(10, 5))
for BnT in [0.01, 0.0125, 0.015, 0.0175, 0.02]:
    pll = sdr.ClosedLoopPLL(BnT, 1)
    t_lock = pll.lock_time(freq)
    plt.plot(freq, t_lock, label=f"$B_nT = {BnT}$")
plt.grid(which="both", linestyle="--")
plt.legend()
plt.xlabel("Normalized frequency offset ($f / f_s$)")
plt.ylabel("Lock time (samples)")
plt.title("Lock time of the closed-loop PLL across input signal frequency offset")
plt.show()