PAM(t)$.Thefollowingparagraphsillustrateaprocedurethatisbasedonsquaring$r(t)$.$

It is more effective to recover the carrier from $r^2(t)$ than $r(t)$ because, while $r(t)$ may have a DC value $𝔼[r(t)]\; =\; 𝔼[x_{}$

PAM(t)] $\approx \; 0$ close to zero, the new signal $r^2(t)$ has discrete spectrum components²⁰ imposed by $𝔼[x_{}^{}$

$PAM(t)]\ne 0$. Using trigonometry again:

= x^2_PAM(t)4[(1+cos (2ω ct))(1+cos(2(ωc+ωΔ)t+2ϕ))]

= x^2_PAM(t)4 [1+cos (2ω ct)+cos(2(ωc+ωΔ)t+2ϕ)+

0.5 (cos((4ωc + ωΔ)t + 2ϕ) + cos(2ωΔt + 2ϕ))](4.70)

= term at DC + terms at high frequencies + x^2_PAM(t) 8 cos(2ωΔt+2ϕ).

Eq. (4.70) indicates that, depending on the values of ${\omega }_c$ and ${\omega }_{\mathrm{\Delta }}$, it may be possible to filter out all terms other than $\frac{x_{}^{}}{}$

$PAM(t)$$8cos(2{\omega }_{\mathrm{\Delta }}t+2\varphi )$ using a bandpass filter. A numerical example is provided in the sequel.

Example 4.21. Frequencies originated by squaring a PAM signal. A baseband signal with bandwidth BW when raised to the power of $n$ has its bandwidth expanded to $n\mathrm{\text{BW}}$. This can be seen by considering that multiplication in time-domain corresponds to convolution in frequency-domain, which expands the support of the original spectrum. Listing 4.17 provides an example of the frequencies obtained via Eq. (4.70) assuming a carrier frequency $f_c=\; 900$ MHz (fc in the code) and a frequency offset of 27 kHz (fdelta), where these frequencies correspond to ${\omega }_c=\; 2\pi f_c$ and ${\omega }_{\mathrm{\Delta }}=\; 2\pi$ fo, respectively. The PAM signal is assumed to have $\mathrm{\text{BW}}=\; 200$ kHz.

1BWpam = 0.2; %signal bandwidth in MHz (200 kHz) 
2fc = 900; %carrier frequency in MHz (900 MHz) 
3fdelta = 0.027; %frequency offset in MHz (27 kHz) 
4squaredBW = 2*BWpam; %BW expands by 2 after squaring 
5%positive (negative freqs. not shown) center frequency of each term: 
6f_DC=0+[-squaredBW, squaredBW] %DC frequency term bandwidth 
7f_2wc=2*fc+[-squaredBW, squaredBW] %two times carrier frequency term 
8f_2wcwoffset = 2*(fc+fdelta)+[-squaredBW, squaredBW] %3rd term 
9f_4wcwoffset = 4*fc+fdelta+[-squaredBW, squaredBW] %4th term 
10f_2woffset = 2*fdelta+[-squaredBW, squaredBW] %term of interest

Listing 4.17 outputs:

     f_DC = [-0.4000    0.4000]
     f_2w0 = 1.0e+03 * [1.7996    1.8004]
     f_2w0woffset =  1.0e+03 * [1.7997    1.8005]
     f_4w0woffset =  1.0e+03 * [3.5996    3.6004]
     f_2woffset = [-0.3460    0.4540]

This indicates that the DC term $\frac{x_{}^{}}{}$

$PAM(t)$$4$ in Eq. (4.70) occupies a band from $-\; 400$ to 400 kHz, overlapping with the term of interest $\frac{x_{}^{}}{}$

$PAM(t)$$8cos(2{\omega }_{\mathrm{\Delta }}t+2\varphi )$, which is located in the band from $-\; 346$ to $454$ kHz.

The overlap can be a problem for a routine that estimates offsets. Hence, the receiver may impose a value distinct from ${\omega }_c$ to purposely have ${\omega }_{\mathrm{\Delta }}>\mathrm{\text{BW}}$, where $\mathrm{\text{BW}}$ is the original PAM bandwidth, as discussed in the sequel.    $\square$

To avoid overlap, ${\omega }_{\mathrm{\Delta }}$ will be assumed to be

where ${\omega }_o$ is the actual (and unknown) offset to be estimated and ${\omega }_{\mathrm{\text{IF}}}$ is an intermediate frequency used by the receiver that can be chosen by the designer to avoid the mentioned overlap.

The FFT-based estimation of frequency and phase offsets is simplified when three values are known (or can be estimated):

f 1.

carrier frequency ${\omega }_c$;

f 2.

PAM bandwidth BW;

f 3.

maximum of the absolute offset value ${\omega }_{\mathrm{\text{max}}}=\; max\; |{\omega }_o|$, which may be determined by a standard (see Application 1.5), obtained from an oscillator datasheet, etc.

In order to choose ${\omega }_{\mathrm{\text{IF}}}$, it is often possible to assume that ${\omega }_c»{\omega }_{\mathrm{\text{max}}}$. But choosing a large value for ${\omega }_{\mathrm{\text{IF}}}$ may be inconvenient, for example, if an ADC conversion should be used afterward. Fortunately, knowing ${\omega }_c$, BW (assumed to be in Hz) and ${\omega }_{\mathrm{\text{max}}}$, allows to conveniently choose ${\omega }_{\mathrm{\text{IF}}}$ as follows.

The angular frequency that will be used for the offsets estimation is within the range $2{\omega }_{\mathrm{\text{IF}}}\pm \; 2{\omega }_{\mathrm{\text{max}}}$, i. e., the minimum search range is $4{\omega }_{\mathrm{\text{max}}}$. Its minimum value $2{\omega }_{\mathrm{\text{IF}}}-\; 2{\omega }_{\mathrm{\text{max}}}$ should be higher than $2\pi (2\mathrm{\text{BW}})$, to avoid overlapping with the “DC term” in Eq. (4.70). Besides, its highest frequency $2{\omega }_{\mathrm{\text{IF}}}+\; 2{\omega }_{\mathrm{\text{max}}}$ should be lower than $2({\omega }_c-{\omega }_{\mathrm{\text{IF}}}-{\omega }_{\mathrm{\text{max}}}-\; 2\pi \mathrm{\text{BW}})$. Hence, ${\omega }_{\mathrm{\text{IF}}}$ (in rad/s) can be chosen within the range

The following example illustrates the usage of Eq. (4.72).

Example 4.22. Choosing the intermediate frequency for IF-based estimation of frequency and phase offsets. Listing 4.18 implements Eq. (4.72) to exemplify a possible range for ${\omega }_{\mathrm{\text{IF}}}$.

1BWpam = 0.2; %signal bandwidth in MHz (200 kHz) 
2fc = 900; %carrier frequency in MHz (900 MHz) 
3fmax = 40/1e6*fc; %maximum absolute frequency offset in MHz (40 ppm) 
4fIFmin = fmax + BWpam %minimum intermediate frequency 
5fIFmax = 0.5*(fc - BWpam) - fmax %maximum intermediate frequency

In this case the values of fIFmin and fIFmax were 236 kHz and 449.9 MHz, respectively. Choosing such extreme values is avoided because they could impose severe restrictions to filtering operations.   $\square$

The following example illustrates a complete FFT-based estimation procedure.

Example 4.23. Complete example of frequency and phase offsets correction for PAM. Listing 4.19 illustrates the complete procedure for estimating the frequency and phase offsets corresponding to a PAM transmission In this example, the carrier frequency is $f_c=\; 1$ GHz, the frequency offset is $-\; 500$ kHz (${\omega }_o=\; -2\pi 500$ rad/s) and the phase offset is $\varphi \; =\; -0.3$ rad. The intermediate frequency (fIF in the code) is chosen as 130 MHz.

6%% 0) Define simulation parameters and also check some choices 
7fc=1e9; %nominal (transmitter's) carrier frequency in Hz (1 GHz) 
8freqOffset=-500/1e6*fc; %frequency offset (-500 ppm) 
9worstFreqOffset=(fc/1e6)*800; %worst accuracy is 800 ppm 
10desiredResolutionHz=1/1e6*fc; %desired max error/FFT resolution,1 ppm 
11%force freqOffset to coincide with a FFT bin via Rsym => Fs 
12Rsym=desiredResolutionHz*1e4; %symbol rate in bauds (10 Mbauds) 
13rollOff=0.5; %roll-off factor for raised cosine filter 
14BWpam=(Rsym/2)*(1+rollOff); %PAM signal BW (~7.5 MHz if raised cos.) 
15fIF = 130e6; %chosen intermediate frequency in Hz 
16P=2; %power the signal will be raised to: x^P(t) 
17phaseOffset=-0.3; %phase offset (to be imposed to the signal) in rad 
18%choose a high enough sampling frequency to represent x^2(t): 
19maxSquaredSignalFreq = 4*fc+fIF+worstFreqOffset+2*BWpam; 
20L=ceil(3*maxSquaredSignalFreq/Rsym); % oversampling factor 
21Fs=L*Rsym; %sampling frequency in Hz 
22M=4; % modulation order (cardinality of the set of symbols) 
23N=5000; % number of symbols to be generated 
24fIFmin = worstFreqOffset + BWpam; %minimum intermediate frequency 
25fIFmax = 0.5*(fc-BWpam)-worstFreqOffset; %max intermediate frequency 
26if fIF > fIFmax || fIF < fIFmin 
27    error('fIF is out of suggested range!') 
28end 
29if freqOffset > worstFreqOffset %based on accuracy of oscillator 
30    error('freqOffset > worstFreqOffset') 
31end 
32%% 1) Transmitter processing 
33%PAM generation: 
34indicesTx = floor(M*rand(1,N))+1; % random indices from [1, M] 
35constellation = -(M-1):2:M-1; % symbols ..., -3, -1, 1, 3, 5, 7,... 
36m=constellation(indicesTx); % obtain N random symbols 
37m_upsampled = zeros(1,N*L); %pre-allocate space with zeros 
38m_upsampled(1:L:end)=m; % symbols with L -1 zeros in - between 
39Dsym=3; % group delay, at rate Rsym (input symbols rate) 
40p=ak_rcosine(1,L,'fir/sqrt',rollOff,Dsym); %square-root raised cosine 
41%p=ak_rcosine(1,L,'fir/normal',rollOff,Dsym); %normal raised cosine 
42gdelayRC=(length(p)-1)/2; %group delay of this even order RC filter 
43sbb=conv(m_upsampled,p);%baseband: convolve symbols and shaping pulse 
44n=0:length(sbb)-1; %discrete-time index 
45wc=(2*pi*fc)/Fs; %convert fc from Hz to radians 
46carrierTx=cos(wc*n); %transmitter carrier (assume phase is 0) 
47s=sbb.*carrierTx; %upconvert to carrier frequency 
48%% 2) First part of the receiver processing 
49wOffset=(2*pi*freqOffset)/Fs; %convert freqOffset from Hz to radians 
50wIF=(2*pi*fIF)/Fs; %convert fIF from Hz to radians 
51carrierRx=cos((wc+wIF+wOffset)*n+phaseOffset);%Rx carrier w/ offsets 
52r=2*s.*carrierRx; %downconvert to approximately baseband 
53%% 3) Estimate and correct offsets: 
54r_carrierRecovery = r.^P; %raise to power P 
55%Will use at least desired resolution: 
56Nfft = max(ceil(Fs/desiredResolutionHz),length(r_carrierRecovery)); 
57if Nfft > 2^24 %avoid too large FFT size. 2^24 is an arbitrary number 
58    Nfft = 2^24; %Modify according to your computer's RAM 
59end 
60if rem(Nfft,2) == 1 %force Nfft to be an even number 
61    Nfft=Nfft+1; 
62end 
63minFreqOffset=P*(fIF-worstFreqOffset); %range for search using ... 
64maxFreqOffset=P*(fIF+worstFreqOffset); %the FFT (in Hz) 
65resolutionHz = Fs/Nfft; %FFT analysis resolution 
66minIndexOfInterestInFFT = floor(minFreqOffset/resolutionHz);%min ind. 
67maxIndexOfInterestInFFT = ceil(maxFreqOffset/resolutionHz);%max ind. 
68R=fft(r_carrierRecovery,Nfft); %calculate FFT of the squared signal 
69R(1:minIndexOfInterestInFFT-1)=0; %eliminate values at the left (DC) 
70%strongest peak within the range of interest (start from index 1): 
71[maxPeak,indexMaxPeak]=max(abs(R(1:maxIndexOfInterestInFFT))); 
72estPhaseOffset=angle(R(indexMaxPeak))/P; %obtain phase in rad 
73%R=[]; %maybe useful, to discard R (invite for freeing memory) 
74estDigitalFreqOffset=(2*pi/Nfft*(indexMaxPeak-1))/P; %frequen. in rad 
75estDigitalFreqOffset=estDigitalFreqOffset-wIF; %deduct IF 
76estFrequencyOffset = estDigitalFreqOffset*Fs/(2*pi); %from rad to Hz 
77%% 4) Demodulate and estimate symbol error rate (SER) 
78%correct carrier offsets (taking in account the imposed wIF): 
79carOffSet=exp(-1j*((estDigitalFreqOffset+wIF)*n+estPhaseOffset)); 
80rc=r.*carOffSet; %generates complex signal with offsets subtracted 
81rc2=2*real(rc); %convert to real signal, take factor of 2 in account 
82rf=conv(p,rc2); %uses matched filter (note p is symmetric and real) 
83startSample=2*gdelayRC+1; %take in account the filtering processes 
84symbolsRx = rf(startSample:L:startSample+L*(N-1)); %extract symbols 
85gain=sqrt(mean(abs(constellation).^2)/mean(abs(symbolsRx).^2)); 
86symbolsRx = gain*symbolsRx; %normalize average energy of Rx symbols 
87indicesRx=1+ak_pamdemod(symbolsRx,M); %add one to make 1,...,M 
88disp(['SER = ' num2str(100*sum(indicesRx~=indicesTx)/N) '%']) %SER 
89disp(['EVM = ' num2str(ak_evm(m, symbolsRx, 1)) '%'])

Figure 4.46 and Figure 4.47 illustrate the baseband PAM signal and its frequency-upconverted version, respectively.

Figure 4.48 depicts the PAM signal after frequency-downconversion to an IF of 130 MHz for the estimation of the offsets.

Figure 4.49 depicts the PSD of the signal used to estimate the offsets, which is obtained by squaring the received signal r. Note that the specific frequency in this case is 2*(fIF+freqOffset)/1e6=259 MHz (corresponding to $2({\omega }_{\mathrm{\text{IF}}}+{\omega }_o)$ rad/s), which is indicated in Figure 4.49 by a datatip. Note that a relatively high sampling frequency of Fs=12.44 GHz was adopted in this simulation in order to properly represent the squared signal. Care should be exercised when using smaller values of Fs in order to avoid aliasing.

The squared received signal (whose PSD is in Figure 4.49) is converted by an FFT to the frequency domain. The values of range=minFreqOffset:maxFreqOffset indicate a search range of 258.4 to 261.6 MHz. As discussed, this range is centered at twice fIF (in this case 260 MHz) and has a band of 4*worstFreqOffset (3.2 MHz). Because in this example the FFT resolution is 1 kHz, the search range corresponds to the FFT indices minIndexOfInterestInFFT=258,400 to minIndexOfInterestInFFT=261,600. Within this range, the FFT index indexMaxPeak=259001 corresponds to the maximum absolute value. Figure 4.50 depicts 20*log(abs(R(range)))), which clearly shows the peak. The FFT value at this frequency (i. e., R(indexMaxPeak)) has a magnitude of approximately 6,200 and an angle of -0.6 rad. The magnitude value is not used but this angle allows to estimate the phase offset as $\varphi \; =\; -0.6∕2\; =\; -0.3$ rad. The frequency value corresponding to ${\omega }_{\mathrm{\Delta }}=\; 2\pi 259∕2$ Mrad/s is obtained from indexMaxPeak as indicated in Listing 4.19.

Having estimates of ${\omega }_{\mathrm{\Delta }}$ and $\varphi$, a complex exponential $e^{-j({\omega }_{\mathrm{\Delta }}n+\varphi )}$ can be created and used to correct the offsets. The squared signal can be eventually discarded, because the next steps will depend on r. Figure 4.51 shows the PSD of rc, which is obtained by multiplying r by the mentioned complex exponential.

Figure 4.51 indicates that rc is complex-valued because its PSD does not have Hermitian symmetry. Discarding its imaginary part leads to the signal with PSD depicted in the top plot of Figure 4.52. This signal is then matched-filtered using the same square-root raised cosine (SRC) adopted by the transmitter, and the final signal is rf, with PSD showed in the bottom plot of Figure 4.52.

The adopted SRC has a group delay of gdelayRC=3,732 samples. Taking in account the two filtering stages (at transmitter and receiver), the first symbol is extracted at sample startSample=2*gdelayRC+1.

The final received symbols constellation is depicted in Figure 4.53.

In this simulation, the SER was 0% (no symbol errors) and the EVM = 1.94%.    $\square$

Example 4.23 illustrated a “perfect” estimation, but there are several difficulties related to this FFT-based procedure. For example:

• Phase ambiguity: Because the phase offset $\varphi$ shows off as a phase of $2\varphi$ at the frequency of interest, there is an ambiguity²¹ in this estimation, unless $\varphi$ is restricted to be within a range of $\pi$ rad (e. g., $[0,\pi [$ or $[-\pi ∕2,\pi ∕2[$). For example, using phaseOffset=2 in Listing 4.19 leads to an erroneous estimated phase offset of $2\; -\; \pi \; \approx -1.1416$ rad.

• FFT leakage: The FFT leakage can decrease the accuracy in case the frequency of interest does not coincide with a FFT bin. For example, using fc=pi/4*1e9 in Listing 4.19 leads to errors in frequency and phase estimation of approximately 110.6 Hz and -0.22 rad, respectively. In this specific case, the SER remains zero and the EVM increases to 2.1%.

The latter issue and the high-computational cost of the FFT-based method motivate the adoption of adaptive algorithms, such as the ones based on the PLL. Before discussing them, the FFT-based method is applied to QAM in the sequel.

QAM carrier recovery using FFT

The offsets estimation procedure used in Listing 4.19 is based on the discrete spectrum component $\frac{x_{}^{}}{}$

$PAM(t)$$8cos(2{\omega }_{\mathrm{\Delta }}t+2\varphi )$ that appears in Eq. (4.70). The implicit condition for having this component available for the estimation is that $𝔼[x_{}^{}$

$PAM(t)]\ne 0$, which is a consequence of having $𝔼[m^2]\ne 0$, where $m$ is the PAM constellation symbol.

As discussed in Appendix A.19.2, some QAM constellations (called here “balanced”) lead to $𝔼[x^2(t)]\; =\; 0$, where $x(t)$ is the transmit signal, and a discrete spectrum component would not be created. In such cases, a power of $P\ne 2$ is used to create $x^P(t)$ as the input to the FFT. The reader is invited to vary the modulation order $M$ and $P$ to evaluate, for example:

1M=16; m=ak_qamSquareConstellation(M); mean(m.^2) %returns 0 
2M=16; m=ak_qamSquareConstellation(M); mean(m.^3) %returns 0 
3M=16; m=ak_qamSquareConstellation(M); mean(m.^4) %it's not zero

With a similar test, it can also be seen that, for a modulation order $M\; =\; 8$ (which corresponds to a “rectangular”, not “square” constellation), $𝔼[m^2]\ne 0$ and, in this specific case, $P\; =\; 2$ would work. But in other cases $P\; =\; 4$ is required. This choice makes a little more complicated the FFT-based offsets estimation, as discussed in the sequel.

A QAM receiver uses both a cosine and sine to recover the in-phase $x_i(t)$ and quadrature $x_q(t)$ components, but the offsets estimation procedure is illustrated here with just the cosine. The transmit QAM signal is $s(t)\; =x_i(t)cos({\omega }_{c}t)\; -x_q(t)sin({\omega }_{c}t)$ and (again considering an ideal channel) the receiver then generates $r(t)\; =\; s(t)cos(({\omega }_c+{\omega }_{\mathrm{\Delta }})t\; +\; \varphi )$ (and it is assumed here that ${\omega }_c>{\omega }_{\mathrm{\Delta }}$). It is the analysis of $r^4(t)$ that will lead to the estimation of ${\omega }_{\mathrm{\Delta }}$ and $\varphi$. The Matlab Symbolic Math Toolbox can be used as follows:

1syms xi xq wt Dt p %define symbolic variables and calculate: 
2simplify(expand(((xi*cos(wt)-xq*sin(wt))*cos(wt+Dt+p))^4))

where Dt corresponds to ${\omega }_{\mathrm{\Delta }}t$, p to $\varphi$, etc. This allows to decompose $r^4(t)$ in 54 parcels that after careful inspection can be written as:

The (real-valued) signal centered at DC is given by $x_0(t)\; =\; 9∕64[2x_i^2(t)x_q^2(t)\; +x_i^4(t)\; +x_q^4(t)]$, and has a bandwidth expanded by approximately $P\; =\; 4$ times the one of the original QAM signal. The other signals $x_{k,i}(t)$ and $x_{k,q}(t)$ are described in Table 4.4.

Some frequency values of interest in Eq. (4.73) are: ${\Theta }_1(t)\; =\; 2{\omega }_{\mathrm{\Delta }}t\; +\; 2\varphi$, ${\Theta }_2(t)\; =\; 4{\omega }_{\mathrm{\Delta }}t\; +\; 4\varphi$ and ${\Theta }_3(t)\; =\; 2{\omega }_{c}t\; -\; 2{\omega }_{\mathrm{\Delta }}t\; -\; 2\varphi$, which correspond to the three lowest frequencies above DC. The term of interest to be used to estimate the offsets is ${\Theta }_2(t)$, which centers the signals $x_{2,i}(t)$ and $x_{2,q}(t)$ at frequency $4{\omega }_{\mathrm{\Delta }}$ rad/s. The highest center frequency $8{\omega }_c+\; 4{\omega }_{\mathrm{\Delta }}$ is associated to the term ${\Theta }_{12}(t)\; =\; 8{\omega }_c+\; 4{\omega }_{\mathrm{\Delta }}+\; 4\varphi$. The other frequencies can be found in snip_channel_qam_squared_freqs.m.

As expected, the signals in Table 4.4 depend on the PAM signals $x_i(t)$ and $x_q(t)$. To better understand why some of these signals do not present a discrete spectral line (as indicated in the last column of Table 4.4), it is useful to simplify the discussion assuming that the continuous-time signals are perfectly sampled at the symbol rate $R_{\mathrm{\text{sym}}}$ such that, for example, instead of statistical moments of continuous-time signals such as $𝔼[x_i^4(t)]$, one can deal directly with the associated constellation symbols and their moments.

Because we are interested in balanced QAM constellations, assume that the two underlying PAMs have the same order $\sqrt{M}$. In this case, the signal $x_{1,i}(t)\; =\; (x_i^4(t)\; -x_q^4(t))∕8$ sampled at $R_{\mathrm{\text{sym}}}$ has an average value of $(𝔼[m_i^4]\; -\; 𝔼[m_q^4])∕8\; =\; 0$ because $𝔼[m_i^4]\; =\; 𝔼[m_q^4]$, where $m_i$ and $m_q$ are the PAM transmit symbols. Similarly, the average of $x_{1,q}(t)$ at the baud rate is $(𝔼[m_i{m}_q^3]\; +\; 𝔼[m_i^3{m}_q])∕4\; =\; 0$ because $m_i$ and $m_q$ are independent, such that e. g. $𝔼[m_i{m}_q^3]\; =\; 𝔼[m_i]𝔼[m_q^3]\; =\; 0$ because $𝔼[m_i]\; =\; 0$. In contrast, $x_{r,1}(t)$ has a non-zero average dictated by $(3∕16)(𝔼[m_i^4]\; +\; 𝔼[m_q^4]\; +\; 2𝔼[m_i^2]𝔼[m_q^2])$. Along this reasoning one can explain the presence of a spectrum line for $k\; =\; 5$ and the absence for $k\; =\; 1$ as indicated in Table 4.4, for example.

A numerical example that further discusses these relations is provided in the sequel.

Example 4.24. Complete example of frequency and phase offsets correction for QAM. The companion script MatlabOctaveBookExamples/ex_fftBasedQAMCarrierRecovery.m is similar to Listing 4.19 but incorporates some modifications to support offsets estimation for balanced QAM constellations, restricted to a phase offset $\varphi \; \in \; [-\pi ∕4,\pi ∕4[$ rad due to the inherent phase ambiguity created by the process of raising the signal to the fourth-power.

Figure 4.54 depicts the PSD of $r^4(t)$, obtained via Eq. (4.73). In this case, snip_channel_qam_squared_freqs.m informs that the 12 frequencies related to ${\Theta }_i(t),i\; =\; 1,\dots ,12$ are: 0.259, 0.518, 1.741, 2.0, 2.259, 2.518, 4.0, 4.259, 4.518, 6.259, 6.518 and 8.518. respectively, all in GHz, as can be identified in Figure 4.54.

As indicated in Table 4.4, there is a factor of 4 between the signals associated to ${\Theta }_3(t)$ and ${\Theta }_{12}(t)$ (located in Figure 4.54 at frequencies 1.741 and 8.518 GHz, respectively), which corresponds to $20{log}_{10}4\; \approx \; 12$ dB. This is confirmed by the difference of the respective PSD values (see the datatips in Figure 4.54): $143.4\; -\; 131.2\; \approx \; 12.2$ dB.

Datatips for the DC and eight of the center frequencies are shown in Figure 4.54. With the exception of the datatip at 259 MHz (from ${\Theta }_1(t)$), all other eight parcels have a discrete line. In other words, among the 13 parcels in Figure 4.54, there are 5 that do not have a discrete line, and the one corresponding to 259 MHz is the only one with a datatip. This is better visualized by zooming Figure 4.54 as done in Figure 4.55.

(a) Lowest band

(b) Second band

(c) Third band

(d) Highest band

The results generated by the script ex_fftBasedQAMCarrierRecovery.m are SER equals to zero and EVM=1.9 %. Figure 4.56 shows the transmit and received constellations after carrier recovery.

Different from the PAM case of Listing 4.19, where the estimated offset was obtained via a normalization by $P\; =\; 2$, the relation between the actual phase offset $\varphi$ and the phase $\Phi$ at the chosen FFT (of the $r^4(t)$ signal) peak is more involved. One indication of this fact is that the signal of interest, related to ${\Theta }_2(t)$ in Table 4.4, is composed by both in-phase $x_{2,i}(t)$ and quadrature $x_{2,q}(t)$ components. To illustrate the issue, the script MatlabBookFigures/figs_synchronism_carrier_recovery.m generated Figure 4.57, which shows in its top plot the relation between $\Phi$ (called fftPeakPhase in the script) versus $\varphi$.

Listing 4.20 shows the lines of the mentioned script that specifically take care of obtaining the estimated offset $\widehat{\varphi }$ (estPhaseOffset in the code) from $\Phi$, according to the top plot in Figure 4.57.

76[maxPeak,indexMaxPeak]=max(abs(R(1:maxIndexOfInterestInFFT))); 
77fftPeakPhase=angle(R(indexMaxPeak)); %extract the phase at FFT peak 
78if fftPeakPhase > 0 %deal with phase ambiguity assume offset is ... 
79    estPhaseOffset = (fftPeakPhase-pi)/4; %in range [-pi/4, pi/4[ 
80elseif fftPeakPhase < 0 
81    estPhaseOffset = (fftPeakPhase+pi)/4; 
82else 
83    estPhaseOffset = 0; %avoid discontinuity at zero

The bottom plot of Figure 4.57 shows the estimation error $\varphi \; -\widehat{\varphi }$. Note that the abscissa range is $[-\pi ∕4,\pi ∕4[$.   $\square$

The FFT-based method used in Listing 4.19 (and ex_fftBasedQAMCarrierRecovery.m) is pedagogical but it is not used in practice due to its relatively high computational cost and the discussed limitations. The next section details the PLL, which is widely used in practice for synchronization purposes.

4.11.3  Phase-locked loop (PLL)

As the name suggests, a phase-locked loop (PLL) is a closed-loop control system that generates an output signal with a phase ${𝜃}_o$ that depends on the input signal phase ${𝜃}_i$. In other words, ${𝜃}_o$ is said to be “locked” to ${𝜃}_i$ because it is a function of the input phase. This relation typically corresponds to a time-delay such that ${𝜃}_o={𝜃}_i-\mathrm{\Delta }𝜃$ or, for modeling purposes and without concerns with causality, a simple equality ${𝜃}_o={𝜃}_i$. Notice that keeping the phases in lock, implies that the frequencies are also locked.

Some of the PLL applications are:

• deriving a stable signal with frequency that is a multiple or submultiple of signal with a reference frequency;

• regenerating the carrier frequency for coherent demodulation;

• generate (synthesize) periodic signals with fundamental frequencies that are multiples of the PLL input frequency, such as a clock tree synchronized with a master clock;

• demodulate signals such as a frequency-modulated signal;

• symbol-timing recovery from a line code such as AMI.

A typical PLL is composed of three main blocks: phase detector, filter and voltage-controlled oscillator (VCO). The PLL can be implemented in the analog or digital domain. A substantial fraction of the literature is devoted to the analog PLL. The digital PLL is, of course, closely related to its analog counterpart, but there are differences on their operations and even terminology. For example, the VCO in the digital PLL is often called numerically controlled oscillator (NCO). This discussion starts by the analog PLL and then concentrates on the digital version.

As pointed out, the PLL is a versatile block and, therefore, used in many applications. There is a wide range of distinct PLL circuits and the name PLL denotes the whole family not a specific circuit. Besides, the PLL output signal depends on the application. For example, the PLL can output a sinusoid $y(t)$ to be used on frequency downconversion in synchronous reception or can eventually already output the signal of interest (e. g., a demodulated FM signal).

AK-TODO: finish this brief introduction on different types of PLL.

In many cases the PLL generates a pulse train (square waveform) or sinusoid.

The PLL can work with sinusoids at its input or digitally modulated signals at the baud rate.

The digital PLL (DPLL) can use NDA (non-data aided).

The digital PLL (DPLL) can use decision-directed schemes.

The strategy is to discuss continuous-time for sinusoids. Then...

4.11.4  Continuous-time PLL for tracking sinusoids

The behavior of a continuous-time PLL is briefly described in the sequel with help from Figure 4.58.

When the PLL is in “lock”, the VCO outputs a signal $y(t)$ that is compared, by the phase detector with $r(t)$. The error signal is then used to drive the VCO to adjust $y(t)$ to track $r(t)$.

More specifically, the phase detector estimates the phase difference $𝜖(t)\; ={\varphi }_i(t)\; -{\varphi }_o(t)$, where ${\varphi }_o(t)$ is the phase of the sinusoid generated by the PLL. This difference $𝜖(t)$ passes through the loop filter, which outputs a control signal $c(t)$ that drives the voltage-controlled oscillator (VCO).

To keep the PLL in lock, the error must be within a range, such as , which is called the lock range. Hence, the PLL operates in one of the modes: acquisition or tracking mode. There are special techniques to bring the PLL to its lock range, i. e., from the acquisition to the tracking mode.

To make the discussion concrete and simple, it is assumed hereafter that the PLL outputs a sinusoid and that its input is $r(t)\; =\; cos(2\pi f_{i}t\; +{\varphi }_i(t))\; +\; \nu (t)$, where $\nu (t)$ is the noise. The PLL goal it to deliver a “clean” version

of $r(t)$.

As mentioned, keeping the phases in lock, implies that the frequencies are also locked. Assume that $f_i\ne f_o$, which can be modeled as $f_o=f_i+{\mathrm{\Delta }}_f$. It is convenient to rewrite

where ${\varphi }_o(t)\; =\; 2\pi {\mathrm{\Delta }}_{f}t\; +{\widehat{\varphi }}_o(t)$. This way, without loss of generality, the PLL treatment can assume that if a frequency offset exists, it is then incorporated in ${\varphi }_o(t)$ as a linear function of $t$ with slope $2\pi {\mathrm{\Delta }}_f$.

The next sections provide more details about the PLL components, assuming a sinusoidal phase detector.

VCO

Before discussing the VCO, it is useful to recall that the instantaneous angular frequency $\omega (t)$ of a sinusoid $sin(𝜃(t))$ is simply the derivative

Inspired in Eq. (4.74) or Eq. (4.75), $𝜃(t)\; =\; 2\pi f_{n}t\; +\; \varphi (t)$ is assumed to be composed by a nominal frequency $f_n$ (in Hz) and time-varying phase $\varphi (t)$ and, from Eq. (4.76):

Now, it should be considered that the VCO has its oscillation frequency controlled by the input voltage $c(t)$. The actual relation between input and output can be nonlinear but it is often modeled as the linear relation

where $f_n$ is the nominal (or quiescent) linear frequency of the oscillator and $g$ is a gain (also called oscillator sensitivity). For simplicity, assume that $g\; =\; 1$, and compare Eq. (4.77) and Eq. (4.78) to obtain the control signal as

Hence, when used in a PLL, the VCO output phase is

such that the VCO output $y(t)\; =\; cos(2\pi f_{i}t\; +{\varphi }_o(t))$, as in Eq. (4.75), is controlled by $c(t)$.

Phase detector

Among many options, a phase detector can be implemented in continuous-time domain with a mixer and lowpass filter. One alternative is to use the trigonometric identity

and, assuming $f_o=f_i$, obtain at the mixer output

A lowpass filter is used to get rid of the term with angular frequency $\omega \; =\; 4\pi f_i$ such that, ignoring the scaling factor $1∕2$, the output of the phase detector composed by the mixer and filter is

when $𝜖(t)$ is sufficiently small such that $sin(x)\; \approx \; x$.

Similar to the variety of PLL circuits, there are many other phase detectors.

Loop filter

From a simplistic point of view, $L(s)$ can be seen as a linear filter. However, because $L(s)$ controls several important PLL features such as stability, settling time, etc., its study and design is based on control systems theory as discussed in the sequel.

A linear model for the continuous-time PLL in the phase domain

It is convenient to perform a PLL analysis that does not depend on the signals amplitudes but only on their phases. This can be done in the so-called “phase domain”, where ${\varphi }_i(t)$ and ${\varphi }_o(t)$ are the input and output phases, respectively. Figure 4.59 represents a linear model that is suitable to design and analyze a PLL when operating in its linear range.

All signals in Figure 4.59 are represented using Laplace, with ${\Phi }_i(s),{\Phi }_o(s),E(s)$ and $C(s)$ being the transforms of ${\varphi }_i(t),{\varphi }_o(t),𝜖(t)$ and $c(t)$, respectively. $N_{\varphi }(s)$ is the Laplace transform of the signal ${\nu }_{\varphi }(t)$, which is the PLL phase noise. In Figure 4.59, the integral of Eq. (4.80) is represented as $1∕s$.

The PLL phase transfer function is defined as $H_{\varphi }(s)\; ={\Phi }_o(s)∕{\Phi }_i(s)$, with the assumption that $N_{\varphi }(s)\; =\; 0$. From Figure 4.59, $E(s)\; ={\Phi }_i(s)\; -{\Phi }_o(s)$ and ${\Phi }_o(s)\; =\; E(s)L(s)(1∕s)$, such that

Similarly, it is possible to obtain the phase error transfer function as

It is useful to recall from control theory that a closed-loop with unitary feedback as in Figure 4.59 has a transfer function

where $G(s)$ is the open-loop transfer function. Eq. (4.82) can be obtained from Eq. (4.84) because in this case $G(s)\; =\; L(s)∕s$.

One alternative for the loop filter is

which corresponds to a PI (proportional-integral) controller. This option is called a second-order PLL given that Eq. (4.82) becomes a second-order system when $L(s)$ is given by Eq. (4.85).

Alternatively, using $L(s)\; =\; \alpha$ would correspond to a first-order PLL. The first-order PLL has a non-zero steady-state phase error when there is a frequency offset (difference between $f_i$ and $f_o$). Hence, due to its improved robustness, the second-order PLL is discussed in the sequel.

Substituting Eq. (4.85) in Eq. (4.82) leads to

which corresponds to Eq. (E.31) when ${\omega }_n=\sqrt{\alpha \beta }$ and $\zeta \; =\; 0.5\sqrt{\alpha ∕\beta }$. As indicated by Figure E.20b, the parameters ${\omega }_n$ and $\zeta$ (or, equivalently, $\alpha$ and $\beta$) can be tuned to obtain the desired behavior of the PLL. For example, a larger bandwidth (primarily controlled by ${\omega }_n$) implies in more phase noise at the output, given that ${\nu }_{\varphi }(t)$ has a wide bandwidth, but allows a faster tracking interval.

The transient behavior of a PLL can be observed by assuming that the input ${\varphi }_i(t)\; =\; u(t)$ is a step function and obtaining the step response associated to Eq. (4.86) (or, equivalently, Eq. (E.31)). Figure 4.60 illustrates the step responses for ${\omega }_n=\; 20$ rad/s and four values of $\zeta$.

From Figure 4.60 it can be seen that a larger $\zeta$ implies an improved tracking performance (faster and without overshoot) but a larger bandwidth, and consequently, output phase noise, as illustrated in Figure E.20b. Due to this tradeoff, $\zeta \; =\; 0.707$ is often used in practice.

Example 4.25. PLL output when there is a frequency offset. Complementing the step response, the ramp response allows to evaluate the PLL when there is a frequency offset ${\mathrm{\Delta }}_f$ such that $f_o=f_i+{\mathrm{\Delta }}_f$ and ${\varphi }_i(t)\; =\; 2\pi {\mathrm{\Delta }}_{f}tu(t)$. Because the assumed model is linear, the output will scale with any input scaling factor, and it suffices to adopt $2\pi {\mathrm{\Delta }}_f=\; 1$ and use the ramp ${\varphi }_i(t)\; =\; tu(t)$ as input.

Figure 4.61 depicts ramp responses for ${\omega }_n=\; 20$ rad/s and two values of $\zeta$. It can be seen that the output presents a transient but then follows the input, which can be seen especially for $\zeta \; =\; 0.5$.

In complement to graphs of transfer functions as Eq. (4.82), it is illustrative to observe the phase error transfer function as in Eq. (4.83), which for the second-order PLL with $L(s)$ given by Eq. (4.85), corresponds to Eq. (E.30).

Figure 4.62 illustrates the PLL error for ${\omega }_n=\; 20$ rad/s and three values of $\zeta$ when the input is the same ramp as in Figure 4.61.

Recall that, when obtaining Figure 4.62, it was assumed $2\pi {\mathrm{\Delta }}_f=\; 1$. The values in Figure 4.62 should be scaled by $2\pi {\mathrm{\Delta }}_f$ when the goal is to determine the lock range, for example. To illustrate, consider the maximum error for $\zeta \; =\; 0.5$ in Figure 4.62 is 0.027 and that the error is restricted to . In this case, , which leads to a lock range given by  Hz.   $\square$

4.11.5  Discrete-time PLL (DPLL) for tracking sinusoids

AK-TODO: include phase detector from Miao, page 397, before or later here.

The block diagram of a discrete-time PLL (also called digital PLL or DPLL) is depicted in Figure 4.63, where the VCO is sometimes called NCO, as previously mentioned.

The PLL input and output in discrete-time are $r[n]\; =\; cos({\mathrm{\Omega }}_{i}n\; +{\varphi }_i[n])\; +\; \nu [n]$ and $y[n]\; =\; cos({\mathrm{\Omega }}_{i}n\; +{\varphi }_o[n])$, respectively. Similarly, the discrete-time version of Eq. (4.79) is $c[n]\; =\frac{\varphi [n]-\varphi [n-1]}{T_s}$, where the sampling period $T_s$ can be substituted by another value, denoted here as $\gamma$, such that Eq. (4.87) can be written as

and $\gamma$ is a gain (or “learning rate”) that can be 1 or another value.

AK-TODO: Livro [FB10] (e parece que [Mia07]) usa(m):

The Z-transform of Eq. (4.87) gives the transfer function

which plays the role of the integrator in continuous-time.

A linear model for the discrete-time PLL in the phase domain

Figure 4.64 describes a linear model that is a discrete-time version of the model in Figure 4.59.

Assuming there is no phase noise (i. e. $N_{\varphi }(z)\; =\; 0$) and from Eq. (4.84) with the open loop transfer function as $G(z)\; =\; (\gamma L(z))∕(1\; -z^{-1})$, the transfer function is

A second-order DPLL is obtained with the adoption of

where $\alpha$ and $\beta$ are the proportional and integral gains, respectively. Eq. (4.91) can be implemented according to the signal flow graph depicted in Figure 4.65.

Using $L(z)$ given by Eq. (4.91) and the integrator of Eq. (4.89) leads to the open loop transfer function²²

Figure 4.65 can be obtained from Eq. (4.91) by adopting a parallel realization of $L(z)\; =H_1(z)\; +H_2(z)$, where $H_1(z)\; =\; \alpha$ and $H_2(z)\; =\; \beta ∕(1\; -z^{-1})$. The corresponding difference equation is $c[n]\; =c_1[n]\; +c_2[n]$, where $c_1[n]\; =\; \alpha 𝜖[n]$ and $c_2[n]\; =\; \beta 𝜖[n]\; +c_2[n\; -\; 1]$, which leads to

Other realizations can be developed by rewriting Eq. (4.91) as

which naturally leads to

The value $c[n\; -\; 1]$ in Eq. (4.95) is the previous output value²³ and allows Eq. (4.95) to be implemented, e. g., as the transposed direct form II illustrated in Figure E.56 or any other standard realization.

Incorporating Figure 4.65 into Figure 4.64 leads to Figure 4.66.

Substituting Eq. (4.94) into Eq. (4.90) leads to

4.11.6  Design of PLL and DPLL

There are several distinct options for implementing the blocks that compose a PLL. Hence, PLLs differ in many details. The reader is warned that AQUI and AQUI, for example, are only options among many alternatives.²⁴

Design of continuous-time PLL

A second-order PLL is assumed here. The loop filter is given by Eq. (4.85) and the VCO is modeled as $1∕s$. Hence, the PLL phase transfer function is assumed to be given by Eq. (4.86), which corresponds to Eq. (E.31). AK-TODO TO-DO - Hence, tabl sosParameters4 provides the relations among parameters of interest.

Design of discrete-time PLL

4.11.7  Discrete-time PLL (DPLL) for tracking digitally modulated... decision directed...