Как я могу пересмотреть скрипт Python для чтения радиочастот?
Я нашел отличную и, возможно, слишком детальную публикацию о том, как получить вещание FM-радио из RTLSDR с использованием Python и библиотеки rtlsdr. https://witestlab.poly.edu/blog/capture-and-decode-fm-radio/
С несколькими изменениями я смог напрямую воспроизвести клип радио-аудио:
#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()
Теперь все вышеперечисленное работает так, как они были спроектированы, но когда я внес несколько изменений, чтобы он читал ретранслятор радиолюбителя (ширина 5 кГц, а не 200 кГц), я получил его на работу, с парой проблем:
- Звук имеет немного более высокую статичность, чем прослушивание в GQRX
- Не бывает шумоподавителя - он не нужен для постоянной Broadcast FM, как бы я обработал шумоподавитель, чтобы не было громкого звука, когда никто не разговаривает?
Есть ли какие-либо шаги или обработка, которую я забыл изменить здесь, пытаясь декодировать любительские радиочастоты?
#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()