In the previous notebook, we have seen that a small carrier frequency offset (CFO) can have a large impact on the quality of the constellation diagram. Hence, we concluded that it's necessary to estimate and compensate any CFO that occurs during the transmission.

In this notebook, we will use the Schmidl&Cox metric to estimate and subsequently compensate the CFO with the aim of improving the constellation diagram in the audio channel. Before we start, let us quickly see mathematically how CFO influences a signal.

Assume a baseband signal $x(t)$ with bandwith $B$ is upconverted using carrier frequency $f_c$. At the receiver side, it is downconverted using carrier frequency $f_c-\Delta_f$. Here, $\Delta_f$ describes the difference between the oscillators at the transmitter and receiver and hence mixing frequencies. Then, the downconverted signal $y(t)$ can be written as

$$ y(t) = \exp\left(-j2\pi (f_c-\Delta_f)t\right)\exp\left(j2\pi f_c t\right) x(t)=\exp\left(j2\pi \Delta_f t\right)x(t) $$Hence, the received signal is the transmitted signal multiplied by a complex exponetial. In the frequency domain we have

$$Y(f)=X(f-\Delta_f),$$i.e. the received signal is a frequency-shifted version of the transmitted signal.

In their 1997 paper "Robust Frequency and Timing Synchronization for OFDM" Schmidl&Cox describe a technique for both time- and frequency offset estimation for OFDM. Before, we have already used their technique for the estimation of the start of the OFDM symbol in the time domain. In this notebook, we will use their metric to estimate the CFO. To recap, the metric $P(d)$ is given by

$$P(d)=\sum_{m=0}^{L-1} r[d+m]^*r[d+m+L].$$Where $r[n]=x[n]\exp(j2\pi \Delta_f n/B)$ is the received signal and the transmitted signal $x[n]$ obeys the periodicity $x[n+L]=x[n]$ where $L$ is half of an OFDM symbol length. Hence we can write $r[d+m+L]=r[d+m]\exp(j2\pi \Delta_f \frac{T}{2})$, since $\frac{T}{2}=L/B$. So, the expression for $P(d)$ becomes

$$ P(d)=\sum_{m=0}^{L-1} r[d+m]^*r[d+m]\exp\left(j\pi \Delta_f T\right)=\exp\left(j\pi \Delta_f T\right) \sum_{m=0}^{L-1}r[d+m]^*r[d+m]$$Now, the sum above is a real number. Hence, the phase of $P(d)$ tells us something about the frequency offset:

$$ \phi=\angle P(d)=\pi \Delta_f T$$and hence we can estimate the frequency offset by

$$\hat{\Delta}_f=\frac{\phi}{\pi T}. $$Let us now gow ahead and check if above derivations numerically hold. First, some common imports:

In [5]:

```
from audioComms.components import Component, TX, Environment, Recorder
from audioComms.plotting import PlotWaveform, PlotSpectrum, PlotConstellation, PlotBlock
from audioComms.channels import AudioChannel, SimulatedChannel
from audioComms.channels import IdentityChannelEffect, SimulatedChannelEffect, NoiseChannelEffect, InitialLossChannelEffect, ChainedChannelEffects, LoseSamplesChannelEffect
from audioComms.utils import ZadoffChu
```

In [6]:

```
from audioComms.passband import PassbandChannel
```

We also import specific OFDM blocks that have been developed in the previous notebooks:

In [7]:

```
from audioComms.ofdm import OFDM, random_qam, ofdm_modulate, MostBasicOFDMReceiver, createSchmidlCoxPreamble
from audioComms.ofdm import BlockPilotRandomDataOFDMTransmitter, BlockPilotLSChannelEstimator, BlockPilotOFDMReceiver
from audioComms.synchronization import SchmidlCoxDetectPeaks, SchmidlCoxMetric
```

In [8]:

```
# Source code has been redacted in online version
# Omitting 17 lines of source code
```

Let's run this function with some CFO of e.g. 5Hz:

In [9]:

```
B = 441*5*4 # baseband bandwidth
cfo = 5 # the experienced CFO
S = recordPMetric(SimulatedChannel, B=B, cfo=cfo)
M_metric = S[1]
P_metric = S[2]
```

Once we have recorded the metrics, let's see the overall metric and the angle of $P$.

In [10]:

```
# Source code has been redacted in online version
# Omitting 4 lines of source code
```

In [11]:

```
plt.figure(figsize=(8,2))
useful_angle = np.angle(P_metric) * (abs(M_metric) > 0.9).astype(int) # filter out parts where the metric is below 0.9
plt.subplot(121)
plt.plot(M_metric)
plt.plot(useful_angle, lw=2, color='r')
plt.grid(True); plt.ylim((-4,4))
plt.subplot(122) # plot a zoomed version of above plot
plt.plot(M_metric)
plt.plot(useful_angle, lw=2, color='r')
plt.grid(True); plt.xlim((1000, 1400)); plt.ylim((-4, 4));
```

In [12]:

```
plt.figure(figsize=(8,2))
plt.plot(M_metric)
K = 256 # OFDM K (needs to match the value in recordPMetric(...))
T = K/B
plt.plot(useful_angle / (np.pi*T), color='r');
plt.grid(True); plt.ylim((-abs(cfo)-1, abs(cfo)+1));
plt.axhline(cfo, ls='--', color='k', lw=3); plt.xlabel('sample index'); plt.ylabel('Measured CFO');
```

Great! Using the Schmidl-Cox formula, we can clearly infer the frequency offset from the metric!

Let us now in the following write a component that measures the current CFO given the S&C P-metric:

In [13]:

```
class CFOExtractor(Component):
"""This component expects to be called with the P-metric over the full OFDM frame. It then extracts the CFO phase
from the position of the plateau of the metric"""
def __init__(self, environment, ofdm, rate):
super().__init__(environment, None)
self._ofdm = ofdm
self._rate = rate
def receive(self, P_metric):
# The end of the plateau is at the end of the first OFDM symbol, i.e. at sample K+CP. Hence, we choose the center of the
# plateau for extraction of the CFO (i.e. K+CP/2 since the plateau consists of CP samples)
extractionPoint = self._ofdm.K + self._ofdm.CP // 2
c = P_metric[extractionPoint + np.arange(1)]
# convert the phase into an actual CFO value and forward it to the next components
self._send(np.angle(c)/(np.pi*self._ofdm.K)*self._rate, stream='CFO')
```

We can now include this component into an actual signal processing chain:

In [14]:

```
def runCFOMeasurement(channelFunc, B=441*5*4, Fc=8000, cfo=0):
ofdm = OFDM(K=256,Kon=200,CP=64,Npayload=1, qam=16)
frameLen = (ofdm.K+ofdm.CP) * (1+ofdm.Npayload)
# Create the components
env = Environment(samplerate=44100)
# transmitter
tx = BlockPilotRandomDataOFDMTransmitter(env, ofdm, usePilotBlock=False)
# channel
chan = PassbandChannel(env, channelFunc=channelFunc, Fc=Fc, B=B, cfo=cfo)
# synchronization
scMetric = SchmidlCoxMetric(env, ofdm=ofdm, minEnergy=1)
scDetect = SchmidlCoxDetectPeaks(env, K=ofdm.K, CP=ofdm.CP, frameLen=frameLen) # extract the frames based on the metric magnitude
cfoExtract = CFOExtractor(env, ofdm, B)
# visualization
fig = env.figure(figsize=(10,6))
showSpectrum = PlotSpectrum(env, windowDuration=0.3, logScale=1, axes=fig.add_subplot(211), xlim=(0, 44100/2), ylim=(-100, -40))
showMetric = PlotWaveform(env, axes=fig.add_subplot(223), signalTransform=abs, ylim=(-0.1, 2), title='Sync metric', numLines=2)
showCFOstream = PlotWaveform(env, axes=fig.add_subplot(224), integerXaxis=True, windowDuration=100, ylim=(-5,5), title='Measured CFO [Hz]')
# set up the connections
tx.transmitsTo(chan)
chan.transmitsTo(showSpectrum, stream='RF')
chan.transmitsTo(scMetric)
scMetric.transmitsTo(scDetect)
scMetric.transmitsTo(showMetric)
scDetect.transmitsTo(cfoExtract, stream='P_metric')
cfoExtract.transmitsTo(showCFOstream, stream='CFO')
env.run()
```

Let's run this chain, first in a simulated channel with a deterministic CFO of 2Hz:

In [15]:

```
runCFOMeasurement(SimulatedChannel, cfo=2, B=441*5*4*1)
```

In [16]:

```
runCFOMeasurement(AudioChannel, cfo=2, B=441*5*4)
```

Now, that we could correctly estimate the CFO, we are required to compensate or correct it before attempting demodulation. The CFO correction simpley applies the inverse operation

$$\hat{y}(t)=y(t)\exp\left(-j2\pi \hat{\Delta}_ft\right)=x(t)\exp\left(j2\pi (\Delta_f-\hat{\Delta}_f)t\right) $$Hence, if we estimated correctly, the CFO is completely removed. Let's write an according component that takes one OFDM frame and corrects its CFO based on the last measurement:

In [17]:

```
class CFOCorrection(Component):
def __init__(self, environment, rate, bypass=False):
super().__init__(environment, None)
self._rate = rate
self._bypass = bypass
self._currentCFO = 0
def receive(self, signal, stream):
if stream == 'CFO':
# update with a new CFO measure
self._currentCFO = signal[0]
elif stream == 'frame':
# correct a frame signal
if self._bypass:
self._send(signal, stream='frame')
else:
t = np.arange(len(signal))/self._rate
factor = np.exp(2j*np.pi*self._currentCFO*t)
self._send(signal*factor, stream='frame')
else:
raise RuntimeError("Unknown stream \"%s\"!" % stream)
```

In [18]:

```
def runTransmissionWithCFOCorrection(channelFunc, B=441*5*4, Fc=8000, cfo=0, bypassCFO=False):
ofdm = OFDM(K=256,Kon=200,CP=64,Npayload=5, qam=16)
frameLen = (ofdm.K+ofdm.CP) * (1+1+ofdm.Npayload)
# Create the components
env = Environment(samplerate=44100)
# transmitter
tx = BlockPilotRandomDataOFDMTransmitter(env, ofdm, usePilotBlock=True)
# channel
chan = PassbandChannel(env, channelFunc=channelFunc, Fc=Fc, B=B, cfo=cfo)
# synchronization
scMetric = SchmidlCoxMetric(env, ofdm=ofdm, minEnergy=1)
scDetect = SchmidlCoxDetectPeaks(env, K=ofdm.K, CP=ofdm.CP, frameLen=frameLen)
cfoExtract = CFOExtractor(env, ofdm, rate=B)
cfoCorrect = CFOCorrection(env, rate=B, bypass=bypassCFO)
# Receiver, use the previously defined MostBasicChannelEstimation
rx = BlockPilotOFDMReceiver(env, ofdm, BlockPilotLSChannelEstimator(ofdm))
# visualization
fig = env.figure(figsize=(10,6))
showSpectrum = PlotSpectrum(env, windowDuration=0.3, logScale=1, axes=fig.add_subplot(311), xlim=(0, 44100/2), ylim=(-100, -40))
showMetric = PlotWaveform(env, axes=fig.add_subplot(323), signalTransform=abs, ylim=(-0.1, 2), title='Sync metric', numLines=2)
showCFO = PlotWaveform(env, axes=fig.add_subplot(324), integerXaxis=True, windowDuration=100, ylim=(-5,5), title='Measured CFO')
showHest = PlotBlock(env, axes=fig.add_subplot(325), numLines=2, ylim=(-4,4), title='Channel Estimation')
showConstellation = PlotConstellation(env, axes=fig.add_subplot(326), title='Rx Constellation', numLines=5)
# set up the connections
tx.transmitsTo(chan)
chan.transmitsTo(showSpectrum, stream='RF')
chan.transmitsTo(scMetric)
scMetric.transmitsTo(scDetect)
scMetric.transmitsTo(showMetric)
scDetect.transmitsTo(cfoExtract, stream='P_metric')
scDetect.transmitsTo(cfoCorrect, stream='frame')
cfoExtract.transmitsTo(showCFO, stream='CFO')
cfoExtract.transmitsTo(cfoCorrect, stream='CFO')
cfoCorrect.transmitsTo(rx, stream='frame')
rx.transmitsTo(showHest, stream='Hest')
rx.transmitsTo(showConstellation, stream='constellation')
env.run()
```

Let's first run this with a CFO of 2Hz and see if the constellation is sharp:

In [19]:

```
runTransmissionWithCFOCorrection(SimulatedChannel, cfo=2)
```

In [20]:

```
runTransmissionWithCFOCorrection(SimulatedChannel, cfo=2, bypassCFO=True)
```

In [23]:

```
runTransmissionWithCFOCorrection(AudioChannel, cfo=0)
```

In [22]:

```
runTransmissionWithCFOCorrection(AudioChannel, cfo=0, bypassCFO=True)
```

Hm, the constellation points have a similar behaviour with or without CFO correction. Hence, we need another method to improve the constellation diagram. Nevertheless, the correction of CFO is an important part in any OFDM system, as we cannot ensure that both transmitter and receiver exhibit no CFO. For example, transmitting from one soundcard to another might have a bit of CFO, when the sound cards are not synchronized enough.

Unfortunately, the successful CFO correction did not really have the expected positive effect on the constellation diagram. Also, we could not measure a significant CFO when using the soundcard. So, the journey for a clearer constellation diagram goes on.

In the next notebook, we will again focus on channel estimation to improve the constellation diagram. We introduce a phase-correction pilot and comb-type pilots estimation procedure, which allows a better equalization in time-variant channels, as the audio channel is one.

Copyright (C) 2018 - dspillustrations.com

DSPIllustrations.com is a participant in the Amazon Services LLC Associates Program, an affiliate advertising program designed to provide a means for sites to earn advertising fees by advertising and linking to amazon.com, amazon.de, amazon.co.uk, amazon.it.