Tutorial: Phase-Locked Loop
Keywords: tutorial phase locked loop PLL
This tutorial demonstrates the functionality of a carrier phase-locked loop and introduces the iirfilt object.
You will need on your local machine:
- the liquid DSP libraries built and installed (see[ref:section-installation] )
- a text editor such as Vim
- a C compiler such as gcc , and
- a terminal
The problem statement and a brief theoretical description of phase-locked loops is given in the next section. A walk-through of the source code follows.
Problem Statement∞
Wireless communications systems up-convert the data signal with a high-frequency carrier before transmitting over the air. This transmitted signal is orthogonal to other signals so long as their bandwidths don't overlap and can be recovered at the receiver by mixing it back down to baseband. Many digital communications systems modulate information in the phase of the carrier requiring the receiver to demodulate the signal coherently in order to recover the original data message. In this regard the receiver must synchronize its carrier oscillator to that of the transmitter. To put it simply, the receiver must lock on to the phase of transmitter's carrier. One of the key advantages to performing signal processing in software is that the radio can operate at complex baseband.
In this simulation, the received signal is simply a complex sinusoid with an unknown initial carrier phase and frequency. The carrier holds no information-bearing symbols and is simply a tone whose frequency and phase represent the residual mismatch between the transmitter and receiver. The received signal \(x\) at time step \(k\) can be described as
$$ x_k = \exp\bigl\{ j(\theta + k\omega) \bigr\} $$where \(j \triangleq \sqrt{-1}\) and\(\theta\) and \(\omega\) represent the unknown initial carrier phase and frequency offsets, respectively. The receiver generates a complex sinusoid with a phase \(\phi_k\) as the phase difference between \(x_k\) and \(y_k\) and can be computed as
$$ y_k = \exp\bigl\{j\phi_k\bigr\} $$The phase error at time step \(k\) is expressed as
$$ \Delta\phi_k = \arg\bigl\{ x_k y_k^* \bigr\} $$where \((^*)\) denotes complex conjugation.
.. footnote
Those who are savvy with communications techniques will
appreciate that we are dealing in complex baseband and can easily
compute the phase error estimate simply as the argument of the
product of $x_k$ and $y_k$.
Conventional PLLs which have operated strictly in the real domain
multiply only the real components of $x_k$ and $y_k$ for a phase
error estimate, assume that the loop filter rejects the
high-frequency component, and make the approximation
$\Delta\phi\approx \sin(\Delta\phi) = \sin(\phi-\hat{\phi})$
for small phase errors.
The goal of the receiver is to control \(\phi_k\) (the phase of the output signal \(y\) at time \(k\) ) to lock onto the input phase of \(x\) , hence the name "phase-locked loop." If the phase of the output sample \(y_k\) is behind that of the input (\(\Delta\phi \gt 0\) ) then \(\phi\) needs to be advanced appropriately for the next time step. Conversely, if the phase of \(y_k\) is ahead of the phase of \(x_k\) (\(\Delta\phi \lt 0\) ) then the receiver need to retard \(\phi\) .
Without going into a great amount of detail, this control is accomplished using a special filter within the loop. This filter, known as a "loop filter," is designed to reject high-frequency noise and is described with the transfer function \(H(z)\) . Specifically \(H(z)\) is a 2\(^{nd}\) -order integrating low-pass recursive filter with a natural frequency \(\omega_n\) , a damping factor \(\zeta\) , and a loop gain \(K\) . The natural frequency is the resonant frequency of \(H(z)\) and for all practical purposes is the filter's bandwidth. Increasing \(\omega_n\) permits the loop to track to the input signal faster (reduces lock time), but also increases the amount of noise passed through the loop. Decreasing \(\omega_n\) reduces this noise but also increases the loop's acquisition time. The damping factor \(\zeta\) controls the stability of the filter and is typically set to a value near \(1/\sqrt{2} \approx 0.707\) . The loop gain \(K\) is typically very large (on the order of \(1000\) or so). For more detailed information on loop filter design the interested reader is referred to [ref:section-nco-pll] .
The estimated phase error \(\Delta\phi_k\) is filtered using \(H(z)\) resulting in an output phase estimate \(\phi_{k+1}\) which is used for the subsequent output sample \(y_{k+1}\) as
$$ y_{k+1} = \exp\bigl\{ j\phi_{k+1} \bigr\} $$.. algorithm [alg-tutoriall_pll] Phase-locked Loop Control
\algsetup{indent=2em}
\begin{algorithmic}[1]
\STATE $\vec{x} \leftarrow \{x_0,x_1,x_2,\ldots\}$ \COMMENT{input array}
\STATE $\hat{\phi}_0 \leftarrow 0$ \COMMENT{initial output phase}
\FOR{$k=0,\,1,\,2,\,\ldots$}
\STATE $y_k \leftarrow \exp\bigl\{ j\hat{\phi}_k \bigr\}$ \COMMENT{compute output sample}
\STATE $\Delta\phi_k \leftarrow \arg\bigl\{ x_k y_k^* \bigr\}$ \COMMENT{phase detector}
\STATE $\hat{\phi}_{k+1} \leftarrow \text{filter}(\Delta\phi_k)$ \COMMENT{update output phase estimate}
\ENDFOR
\end{algorithmic}
A summary of the algorithm is given in [ref:alg-tutoriall_pll] . In the next section we will create a simple C program to simulate a phase-locked loop with liquid .
Setting up the Environment∞
For this tutorial and others, I assume that you are using the GNU compiler collection for compiling source and linking objects and that you have a familiarity with the C (or C++) programming language. Create a new file pll.c and open it with your favorite editor. Include the headers stdio.h , complex.h , math.h , and liquid/liquid.h and add the int main() definition so that your program looks like this:
#include <stdio.h>
#include <complex.h>
#include <math.h>
#include <liquid/liquid.h>
int main() {
printf("done.\n");
return 0;
}
Compile and link the program using gcc :
$ gcc -Wall -o pll pll.c -lm -lc -lliquid
The flag " -Wall " tells the compiler to print all warnings (unused and uninitialized variables, etc.), " -o pll " specifies the name of the output program is " pll ", and " -lm -lc -lliquid " tells the linker to link the binary against the math, standard C, and liquid DSP libraries, respectively. Notice that the above command invokes both the compiler and the linker collectively. If the compiler did not give any errors, the output executable pll is created which can be run as
$ ./pll
and should simply print " done. " to the screen. You are now ready to add functionality to your program.
We will now edit the file to set up the basic simulation but without controlling the phase of the output sinusoid. As such the output won't track to the input resulting in a significant amount of phase error. This simulation will operate one sample at a time and is organized into three sections. First, set up the simulation parameters: the initial phase and frequency offsets ( float ), and number of samples to run ( unsigned int ). Next, initialize the complex input and output variables ( x and y ) to zero, as well as the state of the phase error ( phase_error ) and output phase ( phi_hat ) estimates. Finally, set up the computational loop which generates the input and output samples, computes the phase error between them, and then prints the results to the screen. Edit pll.c to set up the basic simulation:
#include <stdio.h>
#include <complex.h>
#include <math.h>
#include <liquid/liquid.h>
int main() {
// simulation parameters
float phase_offset = 0.8f; // initial phase offset
float frequency_offset = 0.01f; // initial frequency offset
unsigned int n = 40; // number of iterations
float complex x = 0; // input sample
float phase_error = 0; // phase error estimate
float phi_hat = 0; // output sample phase
float complex y = 0; // output sample
unsigned int i;
for (i=0; i<n; i++) {
// generate input sample
x = cexpf(_Complex_I*(phase_offset + i*frequency_offset));
// generate output sample
y = cexpf(_Complex_I*phi_hat);
// compute phase error
phase_error = cargf(x*conjf(y));
// print results to screen
printf("%3u : phase = %12.8f, error = %12.8f\n", i, phi_hat, phase_error);
}
printf("done.\n");
return 0;
}
The variables x and y are of type float complex which contains both real and imaginary components of type float . The function cexpf() computes the complex exponential of its argument which for a purely imaginary input \(j\alpha\) is simply\(e^{j\alpha} = \cos\alpha + j\sin\alpha\) .
Compile and run the program as before. The program should now output something like this:
0 : phase = 0.00000000, error = 0.80000001
1 : phase = 0.00000000, error = 0.81000000
2 : phase = 0.00000000, error = 0.81999999
3 : phase = 0.00000000, error = 0.82999998
4 : phase = 0.00000000, error = 0.84000003
...
35 : phase = 0.00000000, error = 1.14999998
36 : phase = 0.00000000, error = 1.15999997
37 : phase = 0.00000000, error = 1.17000008
38 : phase = 0.00000000, error = 1.18000007
39 : phase = 0.00000000, error = 1.19000006
done.
Notice that because we aren't controlling the output phase yet the error increases with the input phase. In the next section we will design the loop filter to adjust the output phase to lock onto the input signal given the phase error.
Designing the Loop Filter∞
Our program so far has not used any of the liquid DSP libraries for computation and has only relied on the standard C libraries for dealing with complex math operations. In this section we will introduce liquid 's iirfilt_rrrf object to realize a recursive (infinite impulse response) filter with real inputs, coefficients, and outputs. Additionally we will use the function iirdes_pll_active_lag() to design the coefficients for the PLL's filter, specifically an "active lag" design. While the explanation in this section is fairly long, relax! We will only need to add about 15 lines of code to our program. If you are eager to edit your program you may skip to[ref:section-tutoriall_pll-completed] .
Digital representations of infinite impulse response (IIR) filters have two sets of coefficients: feedback and feedforward. In the digital domain the transfer function is a ratio of the polynomials in \(z^{-1}\) where the feedforward coefficients\(\vec{b} = \{b_0, b_1, b_2, \ldots, b_{N-1}\}\) are in the numerator and the feedback coefficients\(\vec{a} = \{a_0, a_1, a_2, \ldots, a_{M-1}\}\) are in the denominator. Specifically, the transfer function is
$$ H(z) = \frac{ b_0 + b_1 z^{-1} + b_2 z^{-2} + \ldots + b_{N-1}z^{-(N-1)} }{ a_0 + a_1 z^{-1} + a_2 z^{-2} + \ldots + a_{M-1}z^{-(M-1)} } $$This transfer function means that the output of the filter is the linear combination of the \(N\) previous filter inputs (\(\vec{x}\) ) and \(M-1\) previous filter outputs (\(\vec{y}\) ), viz
eqnarray:
y[k] =
\frac{1}{a_0}
\Bigl(
b_0 x[k] &+& b_1 x[k-1] + \cdots + b_{N-1} x[k-N]\\
&-& a_1 y[k-1] - \cdots - a_{M-1} y[k-M]
\Bigr)
Typically the number of feedback and feedforward coefficients are equal (\(M=N\) ), and the coefficients themselves are normalized so that \(a_0=1\) .
liquid implements IIR filters with the iirfilt_xxxt family of objects where " xxxt " denotes the type definition (see [ref:section-datastructures] for details). In our example we will be using the iirfilt_rrrf object which indicates that this is an IIR filter with real inputs, outputs, and coefficients with precision of type float . The IIR filter objects in liquid maintain their state internally, storing the previous inputs and outputs in its internal buffers. Nearly every object in liquid (filter or otherwise) has at least four basic methods: create() , print() , execute() , and destroy() . For our program we will need to create the filter object by passing to it a vector of each the feedback and feedforward coefficients. The infinite impulse response (IIR) filter we are designing is of order two which means that \(\vec{a}\) and \(\vec{b}\) have three coefficients each.
Generating the loop filter coefficients is fairly straightforward. As stated before, the loop filter has parameters for natural frequency \(\omega_n\) , damping factor \(\zeta\) , and loop gain \(K\) . Furthermore the filter is 2\(^{nd}\) -order which means that it has three coefficients each for \(\vec{a}\) and \(\vec{b}\) . liquid provides a method for computing such a filter with the iirdes_pll_active_lag() function which accepts \(\omega_n\) , \(\zeta\) , and \(K\) as inputs and generates the coefficients in two output arrays. The coefficients can be computed as follows:
float wn = 0.1f; // pll bandwidth
float zeta = 0.707f; // pll damping factor
float K = 1000.0f; // pll loop gain
float b[3]; // feedforward coefficients array
float a[3]; // feedback coefficients array
iirdes_pll_active_lag(wn, zeta, K, b, a);
The life cycle of the IIR filter can be summarized as follows
iirfilt_rrrf loopfilter = iirfilt_rrrf_create(b,3,a,3);
float sample_in = 0.0f;
float sample_out;
{
// repeat as necessary
iirfilt_rrrf_execute(loopfilter, sample_in, &sample_out);
}
iirfilt_rrrf_destroy(loopfilter);
noting that the execute() method can be repeated as many times as necessary before the object is destroyed.
Using the code snippets above, modify your program to include the loop filter to adjust the output signal's phase. The input to the filter will be the phase_error variable, and its output will be phi_hat . Don't forget to destroy your filter object once the loop has finished running.
Final Program∞
The final program is listed below, and a copy of the source is located in the doc/tutorials/ subdirectory.
#include <stdio.h>
#include <complex.h>
#include <math.h>
#include <liquid/liquid.h>
int main() {
// simulation parameters
float phase_offset = 0.8f; // initial phase offset
float frequency_offset = 0.01f; // initial frequency offset
float wn = 0.10f; // pll bandwidth
float zeta = 0.707f; // pll damping factor
float K = 1000; // pll loop gain
unsigned int n = 40; // number of iterations
// generate IIR loop filter coefficients
float b[3]; // feedforward coefficients
float a[3]; // feedback coefficients
iirdes_pll_active_lag(wn, zeta, K, b, a);
// create and print the loop filter object
iirfilt_rrrf loopfilter = iirfilt_rrrf_create(b,3,a,3);
iirfilt_rrrf_print(loopfilter);
float complex x = 0; // input sample
float phase_error = 0; // phase error estimate
float phi_hat = 0; // output sample phase
float complex y = 0; // output sample
unsigned int i;
for (i=0; i<n; i++) {
// generate input sample
x = cexpf(_Complex_I*(phase_offset + i*frequency_offset));
// generate output sample
y = cexpf(_Complex_I*phi_hat);
// compute phase error
phase_error = cargf(x*conjf(y));
// run error through loop filter
iirfilt_rrrf_execute(loopfilter, phase_error, &phi_hat);
// print results to screen
printf("%3u : phase = %12.8f, error = %12.8f\n", i, phi_hat, phase_error);
}
// destroy IIR filter object
iirfilt_rrrf_destroy(loopfilter);
printf("done.\n");
return 0;
}
Compile the program as before, creating the executable " pll ." Running the program should produce an output similar to this:
iir filter [normal]:
b : 0.32277358 0.07999840 -0.24277516
a : 1.00000000 -1.99995995 0.99996001
0 : phase = 0.25821885, error = 0.80000001
1 : phase = 0.75852644, error = 0.55178112
2 : phase = 1.12857747, error = 0.06147351
3 : phase = 1.27319980, error = -0.29857749
4 : phase = 1.23918116, error = -0.43319979
...
35 : phase = 1.15999877, error = 0.00000751
36 : phase = 1.17000139, error = 0.00000122
37 : phase = 1.18000150, error = -0.00000131
38 : phase = 1.19000030, error = -0.00000140
39 : phase = 1.19999886, error = -0.00000024
done.
Notice that the phase error at the end of the output is very small. The initial error (at \(k=0\) ) is 0.8 which is the value of the phase_offset parameter at the beginning of the program. Notice also that the difference in phase of the last several samples (i.e. the difference between the phase at steps 38 and 39 ) is approximately 0.01 which is the initial frequency offset that was given in the beginning. Play around with the input parameters, particularly the frequency offset and the phase-locked loop bandwidth. Increasing the PLL bandwidth ( wn ) should reduce the resulting phase error more quickly. The downside of having a PLL with a large bandwidth is that when the input signal has been corrupted by noise then the phase error estimate is also noisy. In this tutorial no noise term was introduced.