# Import functions and libraries
from __future__ import division
import numpy as np, matplotlib.pyplot as plt
from numpy import *
from numpy.fft import *
import scipy.signal as signal
from matplotlib.pyplot import *
from rtlsdr import RtlSdr
import scipy
import scipy.io
from scipy.io.wavfile import write
%matplotlib inline
Below are some functions from Part I that we will need in this part
# 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 spectrum
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()
return fig
import pyaudio
import wave
# function that reads wav file to array
def read_wav( wavname ):
wf = wave.open(wavname, 'rb')
CHUNK = 1024
frames = []
data_str = wf.readframes(CHUNK) #read a chunk
while data_str != '':
data_int = np.fromstring( data_str, 'int16') # convert from string to int (assumes .wav in int16)
data_flt = data_int.astype( np.float32 ) / 32767.0 # convert from int to float32
frames.append( data_flt ) #append to list
data_str = wf.readframes(CHUNK) #read a chunk
return np.concatenate( frames )
def play_audio( data, p, fs):
# data - audio data array
# p - pyAudio object
# fs - sampling rate
# open output stream
ostream = p.open(format=pyaudio.paFloat32, channels=1, rate=fs,output=True)
# play audio
ostream.write( data.astype(np.float32).tostring() )
The samples that are obtained by the SDR represent a bandwidth of the spectrum around a center frequency. Hence, when demodulating to base-band (i.e. zero frequency) the signal must be imaginary since it has a non symmetric Fourier transform.
In this case, we would like to display both sides of the spectrum.
myspectrogram_hann_ovlp(x,m,fs,fc)
such that it detects whether the input signal x
is complex. In that case, it will compute a double sided spectrum with frequency range centered around fc (center frequency). For this, it would be useful to use the commands: isreal
and fftshift
.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
# Your code here:
We will first look at radio FM spectrum. In the US the broadcast FM radio band is 88.1-107.9Mhz. It is split into 200KHz slots. This is relatively a large bandwidth and therefore it is also called wideband FM as opposed to narrowband FM which can be as low as 5 Khz. In FM radio the information is encoded by modulating the frequency of the carrier,
$$y_c(t) = A\cos \left (2\pi f_c t + 2\pi \Delta f \int_0^t x(\tau) d\tau \right ).$$
Here, $f_c$ is the carrier frequency, $\Delta f$ is the frequency deviation and $x(t)$ is a normalized baseband signal, which contains all the information the station wants to broadcast. This includes stereo audio, as well as digital information on the station and sometimes additional narrow bandwidth subcarrier channels.
More specifically, the broadcast FM baseband signal, $x(t)$, consists of mono (Left+Right) channels from 30Hz to 15 KHz, a pilot signal at $f_p=19$ KHz, amplitude modulated Stereo (Left - Right) channels around $2\cdot f_p = 38$KHz and digital RBDS, which encodes the station information, song name etc. at $3\cdot f_p =57$KHz. (See http://en.wikipedia.org/wiki/FM_broadcasting for more information).
The baseband signal is:
$ \qquad \qquad x(t) = \underbrace{(L+R)}_{\mathrm{mono}} + \underbrace{0.1 \cdot \cos(2\pi f_p t)}_\mathrm{pilot} + \underbrace{(L-R)\cos(2\pi (2f_p) t)}_\mathrm{stereo} + \underbrace{0.05\cdot \mathrm{RBDS}(t)\cos(2\pi (3f_p) t)}_\mathrm{digital~ RBDS} + \mathrm{subcarriers~signals}. $
The signal $\mathrm{RBDS}(t)$ is a $m(t)\cos(2\pi(3(f_p))$ where $m(t)$ is a binary signal consists of $\pm1$ at constant intervals which encode 0, and 1 bits. The subcarriers are often FM modulated signals of foreign stations or carry other information. This is the spectrum of $x(t)$:
We will listen to our Berkeley own KPFA 94.1MHz station. KPFA transmitter is located close to Grizzly Peak Road, close to where it meets Claremont Ave.
Recall that our SDR does IQ demodulation around a center frequency of choice, $f_d$, followed by a low-pass filter and sampling. This results in a complex digital signal: $$\qquad \qquad y_b(t) = Ae^{j2\pi (f_c-f_d) t + j2\pi \Delta f \int_0^t x(\tau) d\tau }.$$
When $f_d=f_c$ we get: $$\qquad \qquad y_b(t) = Ae^{ j2\pi \Delta f \int_0^t x(\tau) d\tau }.$$
We will implement the following system to demodulate and listen to KPFA:
Like always, you will get the best results if you collect samples outside.
Task:
Tips:
fs = 2400000
fc = 94.1e6
sdr = RtlSdr()
sdr.sample_rate = fs # sampling rate
sdr.center_freq = fc # center frequency
sdr.gain = 24.0
N_Samples = 2560000
data = sdr.read_samples(N_Samples) # get samples
sdr.close()
# Plot. Plotting may take a while!
fig = myspectrogram_hann_ovlp(data, 1024, fs,fc,dbf=30)
The bandwidth of an FM station is 200KHz. Since we sampled at 2.4MHz, you should be seeing several stations in the spectrogram: at 93.3MHz, 94.1MHz and 94.9MHz. KPFA, the station we are interested in is in the middle around 94.1MHz. We will need to filter out and select the desired station from the rest of the signals within the spectrum.
Task:
Design an FIR filter low-pass filter. The filter should be of length 513 and the cutoff frequency be 100KHz. This will select the frequencies between -100KHz to 100KHz. To design the filter, we will use the function signal.firwin
(We will talk more about this filter design technique later in class -- but in essence the function creates a windowed sinc function with the appropriate bandwidth for the given length).
Type signal.firwin?
and execute to get help.
h = signal.firwin(513,100000.0,nyq=2400000.0/2,window='hanning')
Plot the double-sided magnitude frequency response (log scale using plt.semilogy
) by computing a zero-padded FFT to length 1024 and using fftshift. Use KHz as the x-axis units.
# Your code here:
Now, filter the acquired SDR signal with the low-pass filter. Use the function signal.fftconvolve
, which uses the FFT to perform fast convolution. Plot the spectrogram of the filtered signal.
# Your code here:
Our signal is now bandlimited with a bandwidth of 200KHz out of 2.4MHz. We can now downsample the signal to reduce computation without aliasing.
Decimate the signal by a factor of 10 by selecting every 10th sample, to get a signal representing a rate of 240KHz. After that, you have successfully implemented a downsampler through low-pass filtering and decimation!
Plot the spectrogram with a window size of 512. Remember to use the new sampling factor fs/10.0
for the spectrogram.
# Your code here:
Now that you have gone through this processing with 1 second worth of data. We will acquire a longer segment, which will allow you to appreciate the audio that comes out of it once we demodulate it
# Your code here:
The plot does not resemble the broadcast FM spectrum figure above, since it is not yet FM demodulated. We can easily see that the signal is frequency modulated -- because its frequency looks like the time-signal of speech or music.
The next step we are going to perform is to demodulate the FM signal and look at its spectrogram. For this we need to find the instantaneous frequency as a function of time, which is the time derivative of the phase. Computing the phase and then taking the derivative will be sensitive to phase wraps, which we would like to avoid.
Instead, we will take the digital implementation version of the classical FM demodulatation using a limiter, followed by a discriminator. The limiter makes the input have constant amplitude, and the discriminator converts frequency deviations into amplitude modulation. Just as a comment, implementation of an accurate analog disciminator is very difficult whereas implementing a digital one is ridiculously easy!
Recall that $y_b(t) = A(t)e^{ j2\pi (f_c-f_d) t + j2\pi \Delta f \int_0^t x(\tau) d\tau }.$ The leading coefficient $A(t)$ is some unwanted amplitude modulation, which can be a result of signal fading or other sources of signal variations. The role of the limiter component in an FM radio is to remove this amplitude modulation so the discriminator will only be sensitive to frequency variations.
In the digital domain, implementing a limiter is done by simply dividing our signal by its amplitude.
# Your code here:
Assuming that $f_c=f_d$, after the limiter our signal is $ y_b(t)=e^{ j2\pi \Delta f \int_0^t x(\tau) d\tau }.$ To get $x(t)$ we can compute:
$$\left(\frac{d}{dt}y_b(t)\right)y_b^*(t) = j2\pi\Delta f \cdot x(t)\cdot e^{ j2\pi \Delta f \int_0^t x(\tau) d\tau }\cdot e^{ -j2\pi \Delta f \int_0^t x(\tau) d\tau } = j2\pi\Delta f \cdot x(t)$$Which gives us our desired signal when we take the imaginary part of the result.
We will need to design a differentiator filter, which its frequency response approximates the ideal frequency response, $H_{diff}(e^{jw}) = w$, of a differentiator. Since we are only going to demodulate up to 105KHz, we can have the differentiator roll of to zero after that.
signal.remez
which implements an equi-ripple min-max optimal based filter design technique. It requires prescribing frequency bands and their corresponding frequency responses. We need to tell the function that it's a differentiator so it knows to match a linear frequency response within the band of interest. We will prescribe a linear frequency extending from 0-110KHz and then tapering to zero at the Nyquist frequency 120KHz. As we will learn later in class, since the filter is of even order (odd number of coefficients) and is antisymmetric, its Nyquis frequency must be zero!Type signal.remez?
for more information.
Specifically:
h_diff = signal.remez(31,[0.0,105000.0,120000.0,120000.0],[1.05/1.2,0],Hz = 240000.0, type='differentiator')
will try to get a linear frequency response from 0-105000 and zero amplitude at 120000.
# Your code here:
Demodulate the FM signal by first:
Note, that the default implementation of signal.fftconvolve
will have a delay of 16 samples with respect to the original signal. To avoid that, use the option mode='same'
.
Note, that after FM demodulating the signal should be real (by taking only the imaginary component) and hence only half the spectrum should be displayed.
# Your code here:
The mono signal covers the frequency range of 30Hz-16KHz. However, there are many other signals present. There's another problem. Our sampled signal is at a rate of 240KHz. The soundcard on most modern computers can only deal with a sampling rate of 48KHz. Similarly to the downsampling operation we did before, filter our signal and decimate it before being able to play it.
Design a 129 length FIR Bandpass filter with a cuttoff frequency of 16KHz using the command:
h = signal.firwin(129,16000.0,nyq=240000.0/2, window='hanning')
Filter the signal and decimate by a factor of 5 to reduce the rate to 48KHz. Store the result in the a variable called LPR
(Left + Right). Display the spectrogram of the LPR
. Use a window length of 256.
# Your code here:
LPR
so it is in the range between -1 and 1. Do you hear radio?
# Your code here:
# Write to wav file
write('fm_LR.wav', 48000, LPR_sc)
There are two other signals in the spectrum, centered at 67 kHz, and 92 kHz. These are subsidiary communications authorization (SCA) services. The FCC specifically doesn't regulate the FM band from 57 kHz up to 100 kHz, only requiring that whatever you transmit there doesn't interfere with the principle FM signal. For a long time this extra spectrum was used for Muzak (a generic name for elevator music), and other targeted FM signals. This has become less common with the advent of the internet, but there are still stations that use these channels. KPFA transmits a Hatian-French subchannel at 67 kHz, and a Punjabi subchannel at 92 kHz.
Oddly enough, the subchannels themselves are FM encoded. This is like the Russian Matrushka dolls, with one inside another. Fortunately, you have all the tools you need to decode these signals!
First, demodulate the subcarrier to baseband. We will then low-pass to separate it from the other signals. Then, we will then downsample it to 48KHz, and then frequency demodulate again using the technique described above.
The entire system is described in the following diagram:
To demodulate the subcarrier to baseband -- or zero frequency we will need to:
t
representing the samples of the signalTask: Demodulate the frequency demodulated signal in Task 5 by $f_0=67.65$KHz and plot the spectrogram. Since the signal is now complex, you should plot both sides of the spectrum. The subcarrier should be placed at the DC frequency.
# Your code here:
do you see aliasing?
# Your code here:
To frequency demodulate, you will need to apply a limiter as before. The discriminator should be implemented in the same way as before. However, you need to design a new filter which approximates a differentiator between $\pm$8KHz and then tapers to zero. Here's an example:
h_diff = signal.remez(31,[0.0,8000.0,12000.0,24000],[8.0/24.0,0],Hz = 48000, type='differentiator')
Can you hear French Hatian?
# Your code here:
Repeat the procedure to play the subcarrier (Punjabi channel) at $f_0=92$ kHz
# Your code here: