Reading Amateur Radio Frequencies with RTLSDR device and Python

There are many cheap ($5-20) USB RTL-SDR devices you can find on Newegg or Amazon that have an unexpected extra of being able to pick up broadcast radio. This is a handy feature, and in fact we can even read ham-radio frequencies by changing the bandwidth and frequency to read. There are excellent howtos on the easy installation of GQRX for just listening to the radio and setting your frequency and listening, and that’s a good way to start with testing your device. In this post I’ll show how to read and listen to radio using just Python and the Python RTLSDR library, which will show a lot more detail in how decoding radio works, and lets you run a clean interface to listen or record radio.

NYU Polytechnic School of Engineering has an excellent tutorial on using the RTLSDR connected to a computer to read the data and process it using Python. The blog post is nearly a class in a blog post, but you can put all the code together to simply listen to your favorite local broadcast radio, as seen below. When you browse through their code, consider that

  1. Your local machine with the RTLSDR plugged in will not need the plt.savefig() and plt.close() on local system with a screen… nor the “scp -o …” command to copy off a remote server. Generally you can just plug in your device to the computer you’re using so a plt.show() should pause and let you look at the chart.
  2. You can show multiple of those plots on the same screen with the figure subplot in Matplotlib, rather than clicking “x” on each window you open with plt.show():
    Instead of:
    plt.title(…title…) and other commands, use a subplot figure:
    ax1 = figure.add_subplot(rows,cols,1)
    ax1.set_title(…)
    and use only one plt.show() at the end, like so:
import matplotlib.pyplot as plt
from rtlsdr import RtlSdr
from scipy import signal
from time import sleep
import numpy as np
from scipy.io import wavfile
import sounddevice

sdr = RtlSdr()

# configure device
Freq = 91.7e6 #mhz
Fs = 1140000
F_offset = 250000
Fc = Freq - F_offset
sdr.sample_rate = Fs
sdr.center_freq = Fc
sdr.gain = 'auto'
samples = sdr.read_samples(8192000)
print(samples)
print(samples.shape)
print(np.max(samples))
#continue
#sdr.close()
figure = plt.figure(figsize=(30,10))
rows = 2
cols = 2

x1 = np.array(samples).astype("complex64")
fc1 = np.exp(-1.0j*2.0*np.pi* F_offset/Fs*np.arange(len(x1))) 
x2 = x1 * fc1
ax1 = figure.add_subplot(rows,cols,1)
ax1.specgram(x2, NFFT=2048, Fs=Fs)  
ax1.set_title("x2")  
ax1.set_xlabel("Time (s)")  
ax1.set_ylabel("Frequency (Hz)")  
ax1.set_ylim(-Fs/2, Fs/2)  
ax1.set_xlim(0,len(x2)/Fs)  
ax1.ticklabel_format(style='plain', axis='y' )  

bandwidth = 200000#khz broadcast radio.
n_taps = 64
# Use Remez algorithm to design filter coefficients
lpf = signal.remez(n_taps, [0, bandwidth, bandwidth+(Fs/2-bandwidth)/4, Fs/2], [1,0], Hz=Fs)  
x3 = signal.lfilter(lpf, 1.0, x2)
dec_rate = int(Fs / bandwidth)
x4 = x3[0::dec_rate]
Fs_y = Fs/dec_rate
f_bw = 200000  
dec_rate = int(Fs / f_bw)  
x4 = signal.decimate(x2, dec_rate) 
# Calculate the new sampling rate
# Fs_y = Fs/dec_rate  
ax2 = figure.add_subplot(rows,cols,2)
ax2.scatter(np.real(x4[0:50000]), np.imag(x4[0:50000]), color="red", alpha=0.05)  
ax2.set_title("x4")  
ax2.set_xlabel("Real")
ax2.set_xlim(-1.1,1.1)  
ax2.set_ylabel("Imag")  
ax2.set_ylim(-1.1,1.1)   

y5 = x4[1:] * np.conj(x4[:-1])
x5 = np.angle(y5)
ax3  = figure.add_subplot(rows,cols,3)
ax3.psd(x5, NFFT=2048, Fs=Fs_y, color="blue")  
ax3.set_title("x5")  
ax3.axvspan(0,             15000,         color="red", alpha=0.2)  
ax3.axvspan(19000-500,     19000+500,     color="green", alpha=0.4)  
ax3.axvspan(19000*2-15000, 19000*2+15000, color="orange", alpha=0.2)  
ax3.axvspan(19000*3-1500,  19000*3+1500,  color="blue", alpha=0.2)  
ax3.ticklabel_format(style='plain', axis='y' )  
plt.show()


# The de-emphasis filter
# Given a signal 'x5' (in a numpy array) with sampling rate Fs_y
d = Fs_y * 75e-6   # Calculate the # of samples to hit the -3dB point  
x = np.exp(-1/d)   # Calculate the decay between each sample  
b = [1-x]          # Create the filter coefficients  
a = [1,-x]  
x6 = signal.lfilter(b,a,x5)  
audio_freq = 48100.0  
dec_audio = int(Fs_y/audio_freq)  
Fs_audio = Fs_y / dec_audio
x7 = signal.decimate(x6, dec_audio) 

sounddevice.play(x7,Fs_audio) 
x7 *= 10000 / np.max(np.abs(x7))  

#sounddevice.play(x7,Fs_audio)
x7.astype("int16").tofile("wbfm-mono.raw")  #Raw file.
wavfile.write('wav.wav',int(Fs_audio), x7.astype("int16"))
print('playing...')
sounddevice.play(x7.astype("int16"),Fs_audio,blocking=True)

Before continuing you may want to install the libraries if you haven’t already:

pip3 install sounddevice
pip3 install pyrtlsdr
pip3 install scipy
pip3 install numpy

Installing this way is not always considered the “correct” way as you can install it specific to a directory using virtualenv, but this way will let you use all the above libraries from any command line Python.

One big difference from broadcast FM radio and amateur radio FM frequencies, is that amateur low bandwidth FM audio is 5khz where broadcast radio is 100khz. Make some adjustments and you should be able to hear your local ham radio repeater:

Now you can read a number of samples from the radio device, process it, then send it to the sounddevice library, write to a file and or listen! What if you want to just keep listening though? If you loop this code you will have a long wait while it reads more samples again, processes more, then outputs it to the sound card, one step at a time. I tried it with a thread and it was better, but stops for a sample between playing samples. Use a thread for continuously reading from the radio, a thread for continuously processing the input… and on the main thread, reading the processed data, works much better. Queue is a thread-safe way to make the radio reading thread provide that in an ordered way to the processing thread, and the processing thread provide final sounds to the sounddevice thread:

import matplotlib.pyplot as plt
from rtlsdr import RtlSdr
from scipy import signal
from time import sleep
import numpy as np
from scipy.io import wavfile
import sounddevice
import threading
import queue

sdr = RtlSdr()
samples = queue.Queue()
sounds = queue.Queue()

class readThread (threading.Thread):
   def __init__(self, sdr, samples):
      threading.Thread.__init__(self)
      self.srd = sdr
      self.samples = samples
      self.stopit = False
   def run(self):
        print ("Starting read")
        while not self.stopit:
            self.samples.put(sdr.read_samples(8192000/4))

class processThread (threading.Thread):
    def __init__(self, samples, sounds):
        threading.Thread.__init__(self)
        self.samples = samples
        self.sounds = sounds
        self.stopit = False
    def run(self):
        while not self.stopit:
            print('getting')
            samples = self.samples.get()

            x1 = np.array(samples).astype("complex64")
            fc1 = np.exp(-1.0j*2.0*np.pi* F_offset/Fs*np.arange(len(x1))) 
            x2 = x1 * fc1
            bandwidth = 2500#khz broadcast radio.
            n_taps = 64
            # Use Remez algorithm to design filter coefficients
            lpf = signal.remez(n_taps, [0, bandwidth, bandwidth+(Fs/2-bandwidth)/4, Fs/2], [1,0], Hz=Fs)  
            x3 = signal.lfilter(lpf, 1.0, x2)
            dec_rate = int(Fs / bandwidth)
            x4 = x3[0::dec_rate]
            Fs_y = Fs/dec_rate
            f_bw = 200000  
            dec_rate = int(Fs / f_bw)  
            x4 = signal.decimate(x2, dec_rate) 
            # Calculate the new sampling rate
            Fs_y = Fs/dec_rate  
            y5 = x4[1:] * np.conj(x4[:-1])
            x5 = np.angle(y5)

            # The de-emphasis filter
            # Given a signal 'x5' (in a numpy array) with sampling rate Fs_y
            d = Fs_y * 75e-6   # Calculate the # of samples to hit the -3dB point  
            x = np.exp(-1/d)   # Calculate the decay between each sample  
            b = [1-x]          # Create the filter coefficients  
            a = [1,-x]  
            x6 = signal.lfilter(b,a,x5)  
            audio_freq = 41000.0  
            dec_audio = int(Fs_y/audio_freq)  
            Fs_audio = Fs_y / dec_audio
            self.Fs_audio = Fs_audio
            x7 = signal.decimate(x6, dec_audio) 

            #sounddevice.play(x7,Fs_audio) 
            x7 *= 10000 / np.max(np.abs(x7))  
            self.sounds.put(x7)
            print("processed")

thread1 = readThread(sdr, samples)
thread2 = processThread(samples, sounds)
thread1.start()
thread2.start()

# configure device
Freq = 440.713e6 #mhz
Fs = 1140000
F_offset = 2500
Fc = Freq - F_offset
sdr.sample_rate = Fs
sdr.center_freq = Fc
sdr.gain = 'auto'
Fs_audio = 0
try:
    while True:
        x7 = sounds.get()
        #sounddevice.play(x7,Fs_audio)
        #x7.astype("int16").tofile("wbfm-mono.raw")  #Raw file.
        #wavfile.write('wav.wav',int(thread2.Fs_audio), x7.astype("int16"))
        print('playing...')
        sounddevice.play(x7.astype("int16"),thread2.Fs_audio,blocking=True)
except KeyboardInterrupt:
    print("bye")
    thread1.stopit = True
    thread2.stopit = True

Notice that the KeyboardInterrupt handler makes it so you don’t have to Ctrl+C multiple times for the multiple threads. With this code there is only a slight lag while the play function opens up and plays a sound, then closes down. This apparently can be fixed with sounddevice’s outputstream, which is left as an exercise for the reader.

Leave a Reply

Your email address will not be published. Required fields are marked *

+ seventy one = seventy six