# PI Loop Filter

Est. read time: 11 minutes | Last updated: August 12, 2024 by John Gentile

*Contents*

```
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.ticker import EngFormatter
import control as con
from rfproto import pi_filter, nco
```

## Overview

A **Proportional + Integrator (PI) Loop Filter** is used in closed-loop control systems to track input offsets/errors (e.g. frequency or phase offsets in the case of Frequency Error Detectors (FED) or Phase Error Detectors (PEDs), respectively) and provide a corresponding output correction value to drive the loop/system to zero error. Other closed-loop examples using PI filters are Phase Locked Loops (PLLs) and Frequency Locked Loops (FLLs).

It’s basic structure involves an error integrator (scaled by integral gain constant $K_{i}$) summed with a portion of the current error (scaled by proportional gain constant $K_{p}$). For example, when used in a phase correction loop (phase error detector (PED) and phase-rotating NCO), the structure looks like:

As we will see later in deriving PI filter time constants, there are two main knobs for the PI Loop filter: **damping factor** and **loop bandwidth**.

**Damping Factor ($\zeta$):** a damping factor, $\zeta$, is used to control the oscillatory behavior of the impulse response. When $\zeta < 1$, a system is *underdamped* and can exhibit overshoots and undershoots, but fast loop convergence time. When $\zeta > 1$, a system is *overdamped* and oscillatory behavior dissapears, however results in long convergence time. A good balance between convergence time and oscillatory behavior is $\zeta = \frac{1}{\sqrt{2}} \approx 0.707$

```
zeta = 1/np.sqrt(2)
print(f"Damping factor: {zeta}")
```

Damping factor: 0.7071067811865475

**Loop Bandwidth ($B_{n}$):** the frequency response of a PI filter is a *lowpass filter*, and therefore has a an effective bandwidth of frequencies it can respond to (e.g. DC to some 3dB point $B_{n}$). In the case of a PLL/FLL, this is the range of phase and frequency variations the loop can track and adjust to.

Given symbol rate, $R$, the intended loop bandwidth (also known as the *Noise Bandwidth* to refer to the range of frequencies the loop should respond to), $B_{n}$, is suggested to be between:

A higher loop bandwidth allows faster acquisition and can track a larger range of variations, while a lower loop bandwidth can allow for low-noise tracking after acquisition. In this case, some PLL/FLLs use a different set of coefficients to lower the loop bandwidth, and attenuate more noise/variation, when a carrier has been successfully acquired and locked (tracking stage).

```
R = 15e6 # Example symbol rate: 15 MSymbol/s
Bn = 1/100 # Fraction of symbol rate to set loop bandwidth (here 1% of R)
w3db = 2 * np.pi * R * Bn # Target 3dB loop bandwidth (rad/s)
print(f"PI Loop filter bandwidth: {w3db} (rad/s) [{w3db/(2*np.pi*1000)} kHz]")
```

PI Loop filter bandwidth: 942477.7960769379 (rad/s) [150.0 kHz]

## Loop Model

In communications systems, the PI loop filter forms a *Type 2 loop* in order to track offsets to zero error, which means it has two poles at $z = 1$: one pole is the integrator of the PI filter, and the other is the pole formed by the phase accumulator in the NCO. Here we can model a notional Carrier Recovery (CR) loop in the $z$-domain consisting of a Phase Error Detector (PED), PI loop filter and phase compensation NCO.

### PED Model

A PED can take many forms, however here we will model a simple cross product phase detector. The product of one signal and the complex conjugate of the other gives a vector with an angle to the x-axis equivalent to the phase difference between the two signals:

Based on a simplified (e.g. floating-point) PED, the output of this PED ranges from $0 \rightarrow \pi$ radians, as the exact phase difference/measurement. For practical PEDs, the output range (and associated *detector gain*, $K_{d}$) may be different (e.x. fixed-point integer math). Given the aforementioned range, we can say the detector gain (aka discriminator gain) is $K_{d} = \pi$.

```
Kd = np.pi
## PED as a cross-product phase detector:
# NOTE: cross-product detector has gain of A^2 so in FXP, to normalize error signal to
# +/- 1.0, one could adjust Kp/Ki values, or multiply error signal by: Kd = 1/(amp**2)
# However, this error signal directly gives the cross-product angle measurement in a
# range of 0..pi radians.
# See Dan B.'s class notes for more info on how vector product's angle == phase difference and other PEDs
def PED(y0: complex, y1: complex) -> float:
return np.angle(np.conj(y0) * y1)
```

As well, we are assuming this PED has a directly linear relationship between phase error and detected output/error signal. In practical PEDs, this can also be assumed that they are approximately linear when operating within a nominal lock range (e.g. when tracking in vicinity of zero-error). However, different detectors like Decision Directed PEDs (ones that give phase error to a nearest constellation value or other symbol-aided acquisition techniques) may operate in non-linear regions and is beyond the scope of this doc. But for this model, we can assume a linear PED, yielding a transfer function of:

$H_{PED}(z) = K_{d} = \pi$```
H_ped = Kd
```

### NCO Model

An NCO’s basic operation is to accumulate a *Frequency Control Word* (FCW) into a phase accumulator, and use the most-significant-bits as an index into a *look-up table* (LUT) holding sinusoidal sample values.

The scaling of the FCW by the symbol period, $T = 1/R$, and the associated phase accumulator gain, $K_{v}$ is convienent for modeling the loop gain with normalized error (e.g. +/- 1.0). In implementation, the NCO gain is dependent on NCO construction (e.g. phase accumulator width, LUT table size, NCO clock frequency, etc.). For simplicity in demonstrating phase to frequency conversion of the NCO, we’ll model the gain as the inverse of symbol period:

$K_{v} = \frac{1}{T}$The associated transfer function of the NCO as an integrator with parasitic delay in control loop (e.g. a simple $z^{-1}$ delay stage in loop implementation) can be shown as a simple integrator with no zeros:

$H_{NCO}(z) = K_{v}\left(\frac{Tz}{z-1}\right)z^{-1} = \frac{1}{T}\left(\frac{T}{z-1}\right) = \boldsymbol{\frac{1}{z-1}}$```
T = 1/R # Symbol duration (seconds)
Kv = 1/T # NCO gain
# the numerator z in NCO cancels with denominator z from parasitic delay in loop
H_nco = con.tf(Kv*T, [1, -1], T) # NCO transfer function (phase accumulator) w/o numerator z
print(H_nco)
```

Natural frequency: 942477.7960769379 (rad/s) Loop bandwidth: 942477.7960769379 (rad/s)

Now the time constants $$\tau_{1}$$ and $$\tau_{2}$$ can be calculated as: $$ \tau_{1} = \frac{K_{v}K_{d}}{\omega_{n}^{2}},\quad \tau_{2} = \frac{2\zeta}{\omega_{n}} $$ ```python tau1 = Kv*Kd/(wn**2) tau2 = 2*zeta/wn H_lf = con.tf([(T+tau2)/tau1, -tau2/tau1], [1, -1], T) # PI loop filter # Show the combined, open-loop gain response as a function with two-poles and one zero Gol = H_ped * H_nco * H_lf print(Gol) print(f"Tau1: {tau1}") print(f"Tau2: {tau2}") ```

Loop estimate to respond to impulse: 35.0 symbols

We can also show the location of poles (again, at the origin given the two poles at $$z=1$$ from the PI filter and NCO) and zeros of the system: ```python plt.figure() con.pzmap(1/(1+Gol), grid = True) plt.axis([0.95, 1.02, -0.05, 0.05]); plt.show() ``` ![png](PI_filter_files/PI_filter_23_0.png) From here, we can finally calculate the discrete-time, PI loop filter gain constants. **Proportional Gain Constant ($$K_{p}$$):** scalar value directly multiplied to input error value. Therefore it direclty adjusts the amount the current input error contributes to system output response (e.g higher $$K_{p}$$ means a larger change in output for a given change in input error). This proportional term also provides stability control (e.g. phase margin) of the loop: $$ K_{p} = \frac{\tau_{2}}{\tau_{1}} $$ ```python Kp = tau2/tau1 print(f"Kp: {Kp}") ```Kp: 0.0282842712474619

**Integral Gain Constant ($$K_{i}$$):** scalar value directly multiplied to the integrated error value. The integral term is simply an accumulator of error values over time and gives the accumulated offset that should have been corrected previously. This term moves the system output towards a steady-state, setpoint: $$ K_{i} = \frac{T}{\tau_{1}} $$ ```python Ki = T/tau1 print(f"Ki: {Ki}") ```Ki: 0.0012566370614359175

#### $$K_{i}$$ & $$K_{p}$$ using Normalized Loop Bandwidth When a synchronous (every thing runs at same clock/sample rate) 2nd order PI loop is used, and the 2nd pole ($$z = 1$$) stage which uses the output of the PI filter (e.x. an NCO in PLL/FLL), the $$K_{v}$$ term which scales by sample frequency can be eliminated. In this case, $$K_{i}$$ and $$K_{p}$$ derivation can be drastically simplified and only depend on damping factor, $$\zeta$$, and normalized, fractional (e.g. $$0 \rightarrow 2\pi$$) loop bandwidth $$\omega_{frac}$$. Given previous equations: $$ T = 1 / \omega_{f_{s}},\quad K_{v} = 1/T \rightarrow K_{v} = \omega_{f_{s}} $$ $$ \omega_{n} = \omega_{f_{s}} * \omega_{frac} $$ $$ \tau_{1} = \frac{K_{v}K_{d}}{\omega_{n}^{2}},\quad \tau_{2} = \frac{2\zeta}{\omega_{n}} $$ $$ K_{p} = \frac{\tau_{2}}{\tau_{1}},\quad K_{i} = \frac{T}{\tau_{1}} $$ **NOTE:** sample/symbol frequency $$\omega_{f_{s}}$$ assumed pre-scaled by $$\alpha$$ equation to find natural frequency $$\omega_{n}$$, and in cases where damping factor $$\zeta = 1/\sqrt{2}$$, the sample/symbol and natural frequencies are equivalent. Since $$\omega_{f_{s}}$$ cancels out, we can still scale $$\omega_{frac}$$ by that $$\alpha$$ factor equation for other damping factors, thus assume: $$ \omega_{frac} \equiv \frac{\omega_{frac}}{\sqrt{\alpha + \sqrt{\alpha^{2} + 1}}},\quad \alpha = 1 - 2\zeta^{2} $$ **NOTE:** that the $$\omega_{frac}$$ is normalized to radians around a unit circle ($$0 \rightarrow 2\pi$$), so if the fractional bandwidth is given as a percentage of total bandwidth (e.g. sample/symbol frequency), it should be multiplied by $$2\pi$$ before use in these simplified equations. Plug-in constants: $$ K_{p} = \frac{\frac{2\zeta}{\omega_{f_{s}}\omega_{frac}}}{\frac{\omega_{f_{s}}K_{d}}{(\omega_{f_{s}}\omega_{frac})^{2}}} \rightarrow K_{p} = \frac{\frac{2\zeta}{\omega_{f_{s}}\omega_{frac}}}{\frac{K_{d}}{\omega_{f_{s}}\omega_{frac}^{2}}} \rightarrow K_{p} = \frac{2\zeta}{\cancel{\omega_{f_{s}}\omega_{frac}}} * \frac{\cancel{\omega_{f_{s}}}\omega_{frac}^{\cancel{2}}}{K_{d}} \rightarrow \boldsymbol{ K_{p} = \frac{2\zeta\omega_{frac}}{K_{d}} } $$ $$ K_{i} = \frac{\frac{1}{\omega_{f_{s}}}} {\frac{\omega_{f_{s}} K_{d}}{(\omega_{f_{s}} \omega_{frac})^{2}}} \rightarrow K_{i} = \frac{\frac{1}{\omega_{f_{s}}}} {\frac{K_{d}}{\omega_{f_{s}} \omega_{frac}^{2}}} \rightarrow K_{i} = \frac{1}{\cancel{\omega_{f_{s}}}} * \frac{\cancel{\omega_{f_{s}}} \omega_{frac}^{2}}{K_{d}} \rightarrow \boldsymbol{ K_{i} = \frac{\omega_{frac}^{2}} {K_{d}} } $$ ```python w_frac = 2 * np.pi * Bn Kp_new = 2 * zeta * w_frac / Kd Ki_new = w_frac * w_frac / Kd print(f"Kp: {Kp_new}") print(f"Ki: {Ki_new}") assert(Kp_new == Kp) assert(Ki_new == Ki) test_pi_filt = pi_filter.PiFilter(Bn, Kd) assert(Kp_new == test_pi_filt.Kp) assert(Ki_new == test_pi_filt.Ki) ```Kp: 0.0282842712474619 Ki: 0.0012566370614359175

## PI Loop Filter Simulation in Carrier Recovery (CR) Loop