In [1]:
from rtlsdr import RtlSdr
In [2]:
import numpy as np
import scipy.signal as signal
In [3]:
import matplotlib.pyplot as plt
In [4]:
F_rcv = 100.6e6 # 100.6 MHz = Radio City
F_offset = 250e3 # 250 kHz offset to avoid noise
SAMPLE_RATE = 1140000
SECONDS = 10
In [5]:
F_center = F_rcv - F_offset # Tune to (desired freq) - offset
N_SAMPLES = 8192000 #SECONDS * SAMPLE_RATE
Open RTL-SDR
In [6]:
sdr = RtlSdr()
sdr.sample_rate = SAMPLE_RATE
sdr.center_freq = F_center
sdr.gain = "auto"
Capture N_SAMPLES samples and close the device
In [28]:
%%time
raw_samples = sdr.read_samples(N_SAMPLES)
sdr.close()
Convert to complex number array
In [8]:
samples = np.array(raw_samples).astype("complex64")
Raw data
In [9]:
plt.specgram(samples, NFFT=2048, Fs=SAMPLE_RATE)
plt.show()
Undo the frequency shift (center on F_rcv)
In [10]:
shift_func = np.exp(-1.0j*2.0*np.pi* F_offset/SAMPLE_RATE*np.arange(len(samples)))
In [11]:
centered = samples * shift_func
In [12]:
plt.specgram(centered, NFFT=2048, Fs=SAMPLE_RATE)
plt.show()
Crop signal to width BANDWIDTH
In [13]:
BANDWIDTH = 200e3 # Wideband FM has a bandwidth of 200 kHz
In [14]:
decimate_rate = int(SAMPLE_RATE / BANDWIDTH)
focused = signal.decimate(centered, decimate_rate)
decimated_sample_rate = SAMPLE_RATE / decimate_rate
In [15]:
plt.specgram(focused, NFFT=2048, Fs=decimated_sample_rate)
plt.show()
And a polar plot, because why not?
In [16]:
N = 5000
plt.scatter(np.real(focused[0:N]), np.imag(focused[0:N]), alpha=0.05)
plt.show()
Polar discriminator
In [17]:
shifted = focused[1:]
conjugated = np.conj(focused[:-1])
demodulated = np.angle(shifted * conjugated)
In [18]:
plt.specgram(demodulated, NFFT=2048, Fs=decimated_sample_rate)
plt.show()
De-emphasis filter
In [19]:
d = decimated_sample_rate * 75e-6 # Calculate the # of samples to hit the -3dB point
decay = np.exp(-1 / d) # Calculate the decay between each sample
de_emphasised = signal.lfilter([1 - decay], [1, -decay], demodulated)
In [20]:
plt.specgram(de_emphasised, NFFT=2048, Fs=decimated_sample_rate)
plt.show()
Power spectral density plot (i.e. a look from the side)
In [21]:
plt.psd(de_emphasised, NFFT=2048, Fs=decimated_sample_rate)
plt.show()
Finally time to make audio out of this!
In [22]:
AUDIO_FREQ = 44100.0
audio_decimate_rate = int(decimated_sample_rate / AUDIO_FREQ)
audio_sample_rate = decimated_sample_rate / audio_decimate_rate
raw_audio = signal.decimate(de_emphasised, audio_decimate_rate)
In [23]:
plt.specgram(raw_audio, NFFT=2048, Fs=audio_sample_rate)
plt.show()
Increase volume
In [24]:
GAIN = 10000
final_audio = raw_audio * GAIN / np.max(np.abs(raw_audio))
In [25]:
audio_sample_rate
Out[25]:
Output raw PCM to file
In [26]:
final_audio.astype("int16").tofile("audio.pcm")
Play it!
In [27]:
!aplay audio.pcm -r 45600 -f S16_LE -t raw -c 1