We hope you enjoy it!
# Import functions and libraries
from __future__ import division
import numpy as np, matplotlib.pyplot as plt, time, IPython.display
from numpy import *
from numpy.fft import *
from matplotlib.pyplot import *
from rtlsdr import RtlSdr
from scipy import signal
import pyaudio
import time
%matplotlib inline
# 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, fig = None) :
eps = 10.0**(-dbf/20.0) # minimum signal
# 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)*(1-eps) + eps )
# rescale image intensity to 256
img = 256*(y_log + dbf)/dbf - 1
fig=figure(figsize=(16,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, fig = None):
# Plot the spectrogram of x.
# First take the original signal x and split it into blocks of length m
# check if x is real
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 = fft(xmw,len(xmw),axis=0)
fig = sg_plot(t_range, f_range, xmf[0:m//2,:], dbf, fig);
else:
f_range = [-fs / 2.0 + fc, fs / 2.0 + fc]
xmf = abs(fftshift( fft( xmw ,len(xmw),axis=0), axes=0 ))
fig = sg_plot(t_range, f_range, xmf, dbf, fig);
return fig
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 therefore 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, so the absolute frequency offset 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.
Before we start talking about how to use the GSM network to calibrate our rtl-sdr, let's learn a bit about the system
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), which means that several phones use the same channel by time interleaving. 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 have the 820T tuner based rtl-sdrs, and the 1900MHz band is outside of its range (28-1700MHz).
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 phones and base stations 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.
GSM uses Gaussian Minimum Shift Keying (GMSK in short) as a modulation method. We will cover this in part II of this lab. But, in short.... 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)
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 the so called "frequency correction" burst, which is a burst of a pure frequency tone at 1/4th the bitrate of GSM or (1625000 / 6) / 4 = 67708.3 Hz.
GSM has hirarchial structure. At the bottom is the burst, which was discussed above, the next is a TDMA frame, which consists of several bursts and the next is a multi-frame, which consists of several frames.
A TDMA frame consists of 8 time slots and is 4.615ms long. A control channel multi-frame is composed of 51 TDMA frames (numbered 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.
If we receive a signal from a base station, we should see these tone bursts at 67.7083KHz offset with respect to the center frequency of the channel. If the frequency offset of the bursts is different than 67.7083KHz, then we can compute its offset and find the correction that is needed to calibrate our software defined radio.
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. It is posible to develop routines in order to do so automatically, but for the sake of simplicity (and time) we will manually search for GSM base station in our area.
Use either Gqrx or SDR# to look for GSM base station. Here's a screen shot of GSM and CDMA
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 the right frequency.
I've searched the spectrum around Cory hall and found the following frequencies to have GSM base stations: 869.2 (med), 869.4 (hi) , 869.6 (low), 874.4 (med), 874.6 (hi), 879.2 (low), 879.4 (med), 879.6 (hi), 879.8 (med), 890.2 (med), 890.4 (med), 890.6 (low), 890.8 (med), 891.0 (hi), 891.4 (med), where (low), (med), (high) are the signal intensity I measured at my office at 5th floor of cory. NOTE: These are good options to try, but reception varies across the building and campus.
In case you are in a place without GSM basestations, then we recorder some data for you to practice and develop your code on. It can be downloaded from: https://inst.eecs.berkeley.edu/~ee123/sp14/lab/lab4/gsm.npy
sdr.set_freq_correction(ppm)
to correct for some of the offset. Setting sdr.set_freq_correction(1)
will shift the spectrum by 880Hz around 880MHzfs = 240000
fc =
gain =
ppm =
# your code here
sdr = RtlSdr()
sdr.sample_rate = fs # sampling rate
sdr.gain = gain
sdr.center_freq = fc
sdr.set_freq_correction(ppm)
# collect data and throw away the first 2000 samples, since they are no good!
data = sdr.read_samples(25600*3)[2048:]
sdr.close()
# In case you don't have a base station, download gsm.npy and uncomment the next line
#data = np.load('gsm.npy')
# throw away the first 2000 samples, since they are no good!
data = data[2048:]
# 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:1200.0]/fs*1000,abs(data[:1200]))
title('Magnitude GSM signal, showing TDMA frames')
xlabel('t [ms]')
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.
# Your code here
There are many methods for robust frequency correction for GSM. The approach we will take here is using an envelope detection over a frequency band to detect the position in time of the FCCH bursts. Once the the samples of a burst are recognized, we can use the DFT to calculate the frequency offset. Magnitude averaging over several bursts will give us a more 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, and the rest of the GSM signal, we would like to use a narrow-band filter. However, because of the frequency offset the FCCH burst might be outside of our filter bandwwidth and we will miss it! We can solve this in general by using a filter-bank to cover all the possible range of frequency drifts. In addition, narrow-band filters require many filter taps for selectivity. We can save some computation by 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.
For simplicity, in this lab, we will use one narrow-band filter, not a filter bank. To make sure we are looking at the right range, look at the spectrogram above. You should have a good idea what is the offset of the FCCH bursts.
# Your code here
The next step is to detect the position of the FCCH bursts in time and compute their frequency.
There are many ways to detect bursts, ... we leave that as a design challenge for you. Some approaches might be more robust then others.
def locateBurst(x, fs):
# fcch_idx_pos = locateBurst ( x, fs)
# function to find the starting and ending FCCH bursts in a given signal segment
# Inputs:
# x - signal array
# fs - sampling frequency
#
# Outputs:
# fcch_idx_pos and N x 2 array in which each row represend starting and ending indexes of bursts
#
# Your code here
return pulses
# Your test 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.
Of course there are several bursts in our segment. We can combine the estimates by computing the mean power spectrum of all the bursts and then extract the frequency by finding the peak in the mean power spectrum.
def findOffset(data,fcch_pos, fcc, nbw, disp=0):
# freq = findFreq ( data, fcch_pos, fcc, nbw, disp=0)
#
# 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
# fcc - estimated fcch frequency, which you used to demodulate the data
# nbw - sampling rate of the data
# disp - optinally display the average power spectrum
#
# Outputs:
# offset - The estimated frequency offset with respect to ideal fcch frequency.
# Your code here
return (fcch + fcc) - (1625000.0 / 6.0) / 4.0
# Your code here
offset = findOffset( )
print('offset=',offset)
# Your code here
Fantastic! now, you have found the offset in PPM, you can set the driver of the rtlsdr to correct for this offset by using the
sdr.set_freq_correction(-ppm)
instruction. The negative sign is because of your offset was -1ppm, you should tell the sdr to compensate by setting its offset to +1ppm. Unfortunately, the driver only accepts integer ppm corrections, so you can only set it to the closes integer value.
Task:
Write a new function, Kalibrate(Nmframes, fc, fs,gain, fcc, ppm, h_bw)
, which automates all the above.
Nmframe * 235.4ms
, which is effectively collecting Nmframe
number of multi-frames worth of data. It collects them with sdr settings of center frequency fc
, sampling rate of fs
, frequency correction offset in ppm of ppm
, and gain of gain
. The function also takes the estimate FCCH frequency fcc
and the filter bandwidth h_bw
. Debug your function well!
def Kalibrate(Nmframe, fc, fs,ppm,gain, fcc, h_bw, disp=0):
# (offset, ppm) = Kalibrate(Nmframe, fc, fs,ppm,gain, fcc, h_bw, disp=0)
#
# Inputs:
# Nmfram - number of multi-frames length of data to acquire
# fc - center frequency of the sdr driver
# fs - the sampling rate of the sdr driver
# ppm - the frequency correction in ppm for the sdr driver
# gain - gain of the sdr driver
# fcc - estimated fcch frequency, which you used to demodulate the data
# h_bw - Bandwidth of the low-pass filter
# disp - optinally display the average power spectrum
#
# Outputs:
# offset - The estimated frequency offset with respect to ideal fcch frequency.
# ppm - frequency offset in ppm
# Your code here
return (offset, offset / (fc + (1625000.0 / 6.0) / 4.0)*1e6)
# Your debug code here
fs =
fc =
Ndata = 1
gain =
ppm =
offset = Kalibrate(Ndata, fc, fs,ppm,gain, (1625000.0 / 6.0) / 4.0, 15000,disp=1)
print(offset)
Track the frequency drift by calling Kalibrate consequitively for about 4 minutes. Plot the frequency drift in Hz as a function of time. To record the time, record the value of time.time()
before each call to Kalibrate. Can you see the effect of the heating of the crystal oscilator as the chip is being used? Try blowing air on the sdr, to see the effect on the frequency drift.
%matplotlib inline
import time
import pylab as pl
from IPython import display
fs = 240000
fc =
Ndata = 1
gain =
ppm = ;
t = np.zeros(300)
fset=np.zeros(300)
fig = figure(figsize=(16,4))
for n in r_[0:300]:
t[n]=time.time()
fset[n] = Kalibrate(Ndata, fc, fs,ppm, gain, (1625000.0 / 6.0) / 4.0, 15000,disp=0)[0]
if n % 10 == 0:
pl.plot((t[:n]-t[0])/60.0,fset[:n],'b')
display.clear_output(wait=True)
display.display(pl.gcf())
fig = figure(figsize=(16,4))
plt.plot((t[:n]-t[0])/60.0,fset[:n])