# Lab 2: Time-Frequency part II, Calibration of the SDR frequency using GSM signals

### Michael Lustig

Note: The instruction in this lab are now less detailed and require though and common sense with respect to choosing parameters and implementations.

• The lab was inspired by the Kalibrate project by Joshua Lackey

We hope you enjoy it!

In [1]:
# Import general functions and libraries
import numpy as np
import matplotlib.pyplot as plt
import pyaudio
from numpy import pi
from numpy import sin
from numpy import zeros
from numpy import ones
from numpy import r_
from scipy import signal

from numpy import hanning
from numpy import outer

from rtlsdr import RtlSdr
from numpy import mean
from numpy import power
from numpy.fft import fft
from numpy.fft import fftshift
from numpy.fft import ifft
from numpy.fft import ifftshift

from numpy import isreal

%matplotlib inline

In [2]:
# function to compute average power spectrum
def avgPS( x, N=256, fs=1):
M = floor(len(x)/N)
x_ = reshape(x[:M*N],(M,N)) * np.hamming(N)[None,:]
X = np.fft.fftshift(np.fft.fft(x_,axis=1),axes=1)
return r_[-N/2.0:N/2.0]/N*fs, mean(abs(X**2),axis=0)

# Plot an image of the spectrogram y, with the axis labeled with time tl,
# and frequency fl
#
# t_range -- time axis label, nt samples
# f_range -- frequency axis label, nf samples
# y -- spectrogram, nf by nt array
# dbf -- Dynamic range of the spect

def sg_plot( t_range, f_range, y, dbf = 60) :
eps = 1e-3

# find maximum
y_max = abs(y).max()

# compute 20*log magnitude, scaled to the max
y_log = 20.0 * np.log10( abs( y ) / y_max + eps )

fig=figure(figsize=(15,6))

plt.imshow( np.flipud( 64.0*(y_log + dbf)/dbf ), extent= t_range  + f_range ,cmap=plt.cm.gray, aspect='auto')
plt.xlabel('Time, s')
plt.ylabel('Frequency, Hz')
plt.tight_layout()

def myspectrogram_hann_ovlp(x, m, fs, fc,dbf = 60):
# Plot the spectrogram of x.
# First take the original signal x and split it into blocks of length m
# This corresponds to using a rectangular window %

isreal_bool = isreal(x).all()

# pad x up to a multiple of m
lx = len(x);
nt = (lx + m - 1) // m
x = append(x,zeros(-lx+nt*m))
x = x.reshape((m/2,nt*2), order='F')
x = concatenate((x,x),axis=0)
x = x.reshape((m*nt*2,1),order='F')
x = x[r_[m//2:len(x),ones(m//2)*(len(x)-1)].astype(int)].reshape((m,nt*2),order='F')

xmw = x * hanning(m)[:,None];

# frequency index
t_range = [0.0, lx / fs]

if isreal_bool:
f_range = [ fc, fs / 2.0 + fc]
xmf = np.fft.fft(xmw,len(xmw),axis=0)
sg_plot(t_range, f_range, xmf[0:m/2,:],dbf=dbf)
print 1
else:
f_range = [-fs / 2.0 + fc, fs / 2.0 + fc]
xmf = np.fft.fftshift( np.fft.fft( xmw ,len(xmw),axis=0), axes=0 )
sg_plot(t_range, f_range, xmf,dbf = dbf)

return t_range, f_range, xmf


## RTL-SDR frequency Drift

The rtl-sdr has has a 28.8MHz crystal oscillator, which sets the reference frequency used for demodulation on the tuner. The crystal is known to have poor quallity and therefor its frequency can drift considerably with temperature. The drift is measured in parts per millions (PPM) and can go up to $\pm50$ PPM. The result of such a drift is that when you tune to a certain frequency, you will actually be tunning to an offset frequency. This offset is proportional to the center frequency and gets worse for higher frequencies.

This sort of oscillator drifts are very common in other devices, such as your cellphone. The wireless industry has come up with many different techniques in which the clock, or the absolute frequency can be corrected for. In this lab we will use the GSM cellphone network to correct for the frequency drift of the rtl-sdr. GSM base stations are required to have an accurate clock within 0.05 ppm, so they are an excellent source of "accurate" signals. This is very much the first step that a GSM based phone would do in order to connect to the network.

## GSM : Global System for Mobile Communications

This section is heavily based on "GSM for dummies" (http://www.pennula.de/datenarchiv/gsm-for-dummies.pdf)

GSM is a popular digital cellular network. It is based on time-division-multiple-access (TDMA). GSM operates in several different carrier frequencies. In the US and Canada and number of other countries, the frequency bands are 850MHz and 1900MHz. In this lab we will use the 850MHz frequency band, since for the majority of you with 820T tuners the 1900MHz band is outside of the rtl-sdr range.

### Frequencies

For the GSM-850 the band is separated into Uplink frequencies 824-849MHz and downlink frequencies 869-894MHz. The base station transmits on the downlink frequencies and mobile phones transmit on the uplink frequencies. GSM operates in duplex mode in which the phone and base station receive and transmit at the same time. This is the reason that the uplink and the downlink are separated by 50MHz!. A channel number, also known as Absolute Radio Frequency Channel Number (ARFCN) is assigned a pair of uplink and downlink frequencies. There are 124 channels numbered 128-251 in the GSM-850 system. Each channel is allocated a bandwidth of 200KHz. For the downlink frequencies $f_N = 869.2+0.2*(ARFCN-128)$ MHz.

### Modulation

GSM uses Gaussian Minimum Shift Keying (GMSK in short) as a modulation method. It is a binary digital Frequency Modulated scheme in which the phase between each bit period is continuouse. Bits are encoded as different frequencies separated by one-half the bit-rate. The modulation rate in GSM is $1625/6 \approx 270.833$ kb/s. For more information read (http://en.wikipedia.org/wiki/Minimum-shift_keying)

### Time-Division-Multiple-Acess

GSM uses time-division to share a frequency channel among users. Each frequency is divided into blocks of time that are known as time-slots. There are 8 time-slots that are numberes TS0 - TS7. Each time-slot lasts about 576.9 $\mu s$. Given the bit-rate above, a total of 156.25 bits can be transmitted in each slot. Each slot allocates 8.25 of it "bits time" as guard-time, split between the beginning and the end of each time-slot. This time is necessary to prevent overlapping between the different time slots. In addition, 3 bits on both sides of the time-slot do not contain any data and are there to allow for ramping the amplifiers up and down.

Data transmitted within each time-slot is called a burst. There are several type of burts, with different purposes. We are interested in a frequency correction burst, which is a burst of pure frequency at 1/4th the bitrate of GSM or (1625000 / 6) / 4 = 67708.3 Hz.

Figure 1: TDMA frames and bursts

### Frame Structure, and Frequency Correction Channel (FCCH)

A TDMA frame consists of 8 time slots and is 4.615ms long. A control channel multi-frame is composed of 51 TDMA frames (0-50) and is 235.4ms long. A base station transmits a frequency correction burst in regular positions. The frequency correction channel repeats every 51 frames and the burst occurs at the TS0 time-slot in frames 0,10,20,30 and 40 in the control channel multi-frame.

Figure 2: Frequency correction bursts at 67708.3Hz are transmitted at time slot TS0 on frames 0,10,20,30 and 40 in a control hannel multi-frame. That sequence repeates every 51 TDMA frames. So the interval between bursts within a frame is 46.15ms, and the interval between the last burst in a frame to the first burst in the next is 50.765ms

If we receive a signal from a base station, it should have these tone bursts at 67.7083KHz with respect to the center frequency of the channel. If the frequency is different than that, then we can compute its offset and find the correction offset that is needed to calibrate our software defined radio.

### Finding a GSM base station

The first thing to do is to find a GSM base station in your area, that has a good steady signal. If you have poor cell reception, you might want to try this near cory hall -- where we were able to pick up a few. Unfortunately, the downlink frequency band is shared with spread-spectrum LTE and 3G stations as well, so detection based on the energy of a signal in a channel is not enough. It is easy to find a base station based on the spectrum shape in the channel and also based on the frequency correction bursts. But, we need to develop some routines in order to do so automatically.

What we will do first is to acquire samples from each channel, compute the average power spectrum and manually check if the spectrum looks like a GSM spectrum. In addition, we collected ~300ms of data from a GSM base station in our area so you can see what you should be looking for! The data file is located on the class website at: https://inst.eecs.berkeley.edu/~ee123/sp14/lab/lab2/gsm.npy

Another way is to use one of your favorite SDR software, like Gqrx or SDR# to look for GSM base station. Here's a screen shot of GSM and CDMA

Figure 3: Spectrum of GSM and spread spectrum when sampled at 2MS/s

Here's the average power spectrum as well as the magnitude of the time signal. Note that for GSM, the energy is concentrated towards the center of the spectrum and drops at the edges smoothly. You can also see that we have a considerable amount of frequency offset, because the spectrum is not centered around 0. In the time domain signal you can also see the guard time between TDMA frames!

I've searched the spectrum around Cory hall and found the following frequencies to have GSM base stations: 869.2, 869.6, 874.6, 879.8, 890.2, 890.4, 890.6, 890.8, 891.0, 891.4,

In [3]:
fs = 240000
fc = 869.4e6
#sdr = RtlSdr()
#sdr.sample_rate = fs    # sampling rate
#sdr.gain = 36
#sdr.center_freq = fc
#sdr.close()

In [4]:
#numpy.save('gsm.npy',data)
# throw away the first 2000 samples, since they are no good!
data = data[2000:]

In [5]:
# compute average power spectrum
f, sp = avgPS(data,N=256,fs=fs/1000)
fig = figure(figsize=(8,4))
plot(f,10*log10(sp))
title('average power spectrum of GSM')
xlabel('frequency offset [KHz]')

# plot
fig = figure(figsize=(16,4))
plot(r_[0:1220.0]/fs*1000,abs(data[:1220]))
title('Magnitude GSM signal, showing TDMA frames')
xlabel('t [ms]')

Out[5]:
<matplotlib.text.Text at 0x107c25250>


• Find a GSM base-station. Collect 25600*3 samples of data. Discard of the first 2000 samples as they are generally not valid.
• Compute the signals's average power spectrum.
• Display both the power spectrum and a magnitude-time domain section that corresponds to 1.1 TDMA frames (1220 samples).
• Can you see 8 TDMA time-slots?

One way to see the frequency correction bursts is to compute a spectrogram of the signal. This is very effective, though not very computation efficient for the purpose of detecting the bursts.

• Use the function myspectrogram_hann_ovlp from the previous section of the lab.
• For the length of the burst, what length should the spectrogram window be in order to detect the bursts? Play with different parameters.
• You should be able to see about 6 frequency correction bursts. Also note that in many cases the frequency correction bursts can be 30-50KHz away from where they are supposed to be!
In []:
# plot the spectrogram of the GSM signal


### Finding the FCCH bursts and computing the frequency offset

There are many methods for robust frequency correction for GSM. The approach we will take here is using filter-banks to detect the position of the FCCH bursts. Once the position is determined, we can use the DFT to calculate the frequency offset. Averaging over several bursts will give us a robust estimation of the offset.

Of course, there are going to be many tradeoffs in the design. For example, to be able to distinguish the burst from noise, we would like to use narrow-band filters. However this will result in a large filter bank array to cover all the possible range of frequency drifts. At the same time, narrow-band filters require many filter taps for selectivity. We can save some computation be using decimating filters -- so the detection can be performed at a lower rate. A polyphase implementation of these filters would also be efficient, though you are probably not going to observe much improvements in the python environment.

Before going for a filter bank, lets look at the result of a single filter. From the spectrogram, you should have a good idea what is the offset of the FCCH bursts.

• Calculate the length of a FCCH burst in samples?
• To create a bandpass filer, we will first design a low-pass filter with a pass-band bandwidth of $\pm3.75$KHz. The number of taps should be odd and correspond roughly to the length of the FCCHburst. What is the TBW product of the filter in that case??
• Modulate the filter such that its center frequency is centered around the offset FCCH burst frequency.
• Filter the data and plot the magnitude of the result, can you see the FCCH burts?
• plot the spectrum of the filtered signal as well
In []:
# calculate the number of samples in a FCCH burst (including guard time)
# calculate the TBW of the filter

In []:
# your code here


Because the signal is filtered, we can decimate it without aliasing and detect the bursts at a much lower rate.

• Decimate the filtered signal appropriately and display the result.
• Plot the spectrum of the decimated signal as well. Given the decimation rate and the center frequency of the filter, in what aliased frequency the center frequency is going to end up at after decimation? Keep this value! it will be necessary in order to figure out what is the real frequency of the FCCH bursts in the decimated segment.
In []:
In []:
# Calculate where the center frequency is aliasing to:


The next step is to detect the position of the FCCH bursts in time and compute their frequency.

• Write a function that will accept the signal and a threshold (ratio to the max peak) as an input. the function will find the location of the FCCH bursts, and return the starting and ending indexes of each burst.

Optional: The function could also prune falsely detected bursts. For example, burst that are too short or too long with respect to what they are supposed to be.

• Apply the function on the decimated signal to obtain the FCCH bursts positions.
In []:
def locateBurst(x, thresh):
#  fcch_idx_pos = locateBurst ( x, thresh)
#  function to find the starting and ending FCCH bursts in a siven signal segment
#  Inputs:
#           x  - signal array
#           thresh - threshold with respect to the maximum of magnitude of x in which a burst is considered detected
#
#  Outputs:
#           fcch_idx_pos and N x 2 array in which each row represend starting and ending indexes of bursts
#


In []:
# Your code here


Now that we detected the bursts, we can extract just bursts signals and compute their frequency by finding the peak of a zero-padded DFT. We could do this either by extracting the segments from the decimated signal or from its original non-decimated. The former is simpler and faster whereas the latter may be slightly more accurate. We will first use the decimated signal. If you wish, you can later use the original and compare the results.

Of course there are several bursts in our segment. We can combine the estimate in several ways: 1) find the frequency in each segment, and then average the result. If we compute the standard deviation we will also get some information on how reliable is the estimate. 2) combine the result of the DFT's by computing the mean power spectrum and then extract the frequency by finding the peak in the spectrum.

• Write a function that will take the signal the FCCH burst positions, the center frequency of the filter and the bandwidth of the decimated signal. The function will extract the segments of the FCCH bursts, zero-pad them appropriately (why zero-pad?), find the actual frequency of the FCCH bursts within the decimated signal, and calculate the actual frequency.

Optional: The function can also return the standard deviation

In []:
def findFreq(data,fcch_pos, h_fc,nbw):
# freq = findFreq ( data, fcch_pos, h_fc, nbw)
#
# The function takes the signal the FCCH burst positions, the center frequency of the filter and
# the bandwidth of the decimated signal. The function extracts the segments of the FCCH bursts,
# zero-pad them appropriately, finds the actual frequency of the FCCH bursts within the decimated signal,
# and calculate the actual frequency.
# Inputs:
#        data - an array of data
#        fcch_pos - positions of fcch bursts generated by locateBurst
#        h_fc -  The original undecimated center frequency of the band
#        nbw  -   Bandwidth of the decimated band
#
# Outputs:
#       freq - The estimated frequency of the the fcch bursts.


In []:
# your code here


Now that you have found the frequency of the FCCH burst, you can caluclate the frequency correction offset to be f_fcch-(1625000.0 / 6.0) / 4.0

• what is the shift in Hz and in PPM?
In []:

Fantastic! now, you have calibrated your SDR to have a very small frequency offset. Apply the calculated frequency offsrt. Choose a different GSM station from the above list and run the code again to see if your requency is adjusted. Don't forget to use a filter centered around 67708Hz!

### Task III: (Mini-project) Automatic Frequency calibration

Write a function (or a set of functions) to find the frequency offset of your SDR based on a GSM base station signal. The function input would be the frequency of a known GSM station. Your function must be capable of detecting up to $\pm50$ppm offsets. Of course in this case a single filter will not be enough and you will have to use a filter bank.

Optional: Write a function that will search all GSM channels, find a GSM base-station automatically and find the frequency offset from it.

In []: