I found an excellent and perhaps-too-detailed posting on how to get broadcast FM Radio from a RTLSDR using Python and the rtlsdr library. https://witestlab.poly.edu/blog/capture-and-decode-fm-radio/
With a few modifications I was able to directly play a clip of radio audio:
#from pylab import psd, xlabel, ylabel, show
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()
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
# plt.specgram(x2, NFFT=2048, Fs=Fs)
# plt.title("x2")
# plt.xlabel("Time (s)")
# plt.ylabel("Frequency (Hz)")
# plt.ylim(-Fs/2, Fs/2)
# plt.xlim(0,len(x2)/Fs)
# plt.ticklabel_format(style='plain', axis='y' )
# plt.show()
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
# plt.scatter(np.real(x4[0:50000]), np.imag(x4[0:50000]), color="red", alpha=0.05)
# plt.title("x4")
# plt.xlabel("Real")
# plt.xlim(-1.1,1.1)
# plt.ylabel("Imag")
# plt.ylim(-1.1,1.1)
# plt.show()
y5 = x4[1:] * np.conj(x4[:-1])
x5 = np.angle(y5)
# plt.psd(x5, NFFT=2048, Fs=Fs_y, color="blue")
# plt.title("x5")
# plt.axvspan(0, 15000, color="red", alpha=0.2)
# plt.axvspan(19000-500, 19000+500, color="green", alpha=0.4)
# plt.axvspan(19000*2-15000, 19000*2+15000, color="orange", alpha=0.2)
# plt.axvspan(19000*3-1500, 19000*3+1500, color="blue", alpha=0.2)
# plt.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)
#print('try2')
#Sounds no good
#sounddevice.play(x7, blocking=True)
# use matplotlib to estimate and plot the PSD
#psd(samples, NFFT=1024, Fs=sdr.sample_rate/1e6, Fc=sdr.center_freq/1e6)
#xlabel('Frequency (MHz)')
#ylabel('Relative power (dB)')
#print("showing")
#show()
Now the above works as they designed it, but when I make a few changes to make it read a ham radio repeater (5khz width rather than 200khz), I got it to work, with a couple problems:
- The audio has a bit more highpitched static than listening in GQRX
- There isn't a squelch - unneeded with the constant Broadcast FM it appears, how would I process a squelch to not have loud static whenever no one is talking?
Are there any steps or processing I am forgetting to change here trying to decode amateur radio frequencies?
#from pylab import psd, xlabel, ylabel, show
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 = 440.712e6 #mhz
Fs = 1140000
F_offset = 2500
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()
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
plt.specgram(x2, NFFT=2048, Fs=Fs)
plt.title("x2")
plt.xlabel("Time (s)")
plt.ylabel("Frequency (Hz)")
plt.ylim(-Fs/2, Fs/2)
plt.xlim(0,len(x2)/Fs)
plt.ticklabel_format(style='plain', axis='y' )
plt.show()
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
# plt.scatter(np.real(x4[0:50000]), np.imag(x4[0:50000]), color="red", alpha=0.05)
# plt.title("x4")
# plt.xlabel("Real")
# plt.xlim(-1.1,1.1)
# plt.ylabel("Imag")
# plt.ylim(-1.1,1.1)
# plt.show()
y5 = x4[1:] * np.conj(x4[:-1])
x5 = np.angle(y5)
# plt.psd(x5, NFFT=2048, Fs=Fs_y, color="blue")
# plt.title("x5")
# plt.axvspan(0, 15000, color="red", alpha=0.2)
# plt.axvspan(19000-500, 19000+500, color="green", alpha=0.4)
# plt.axvspan(19000*2-15000, 19000*2+15000, color="orange", alpha=0.2)
# plt.axvspan(19000*3-1500, 19000*3+1500, color="blue", alpha=0.2)
# plt.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 = 41000.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)
#print('try2')
#Sounds no good
#sounddevice.play(x7, blocking=True)
# use matplotlib to estimate and plot the PSD
#psd(samples, NFFT=1024, Fs=sdr.sample_rate/1e6, Fc=sdr.center_freq/1e6)
#xlabel('Frequency (MHz)')
#ylabel('Relative power (dB)')
#print("showing")
#show()
