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.

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.

gsm
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.

gsm
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

gsm
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
#data = sdr.read_samples(25600*3)
#sdr.close()
In [4]:
#numpy.save('gsm.npy',data)
data = numpy.load('gsm.npy')
# 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>

Task I:

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.

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.

Task II:

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.

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.

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.

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.

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

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 []: