Laboratorul 08.

În acest laborator veți face diverse procesări pe un semnal obținut de dispozitivul RTL-SDR.

Găsiți aici o carte informativă despre acest dispozitiv și despre cum funcționează. Pe scurt, acest dispozitiv captează semnale dintr-o gamă mare de frecvențe (25MHz – 1.75GHz) pe care apoi le poate transforma prin IQ downconversion în semnale cu o frecvență de eșantionare mult mai mică (selectând bineînțeles frecvența de interes pe care o dorim). O schemă din acea carte care arată funcționalitatea este următoarea:

Spectrul unui semnal FM arată în felul următor:

Veți putea folosi un dispozitiv RTL-SDR dacă aveți, sau un fișier pre-generat de noi, pe care îl puteți descărca de aici (câteva secunde de pe Radio Trinitas – 88.5MHz – trebuie dezarhivat înainte de folosire).

Pentru a urma acest laborator, folosiți scheletul de cod de mai jos, din Python. Acolo veți vedea 7 TODO-uri, fiecare este punctat 2 puncte, deci aveti un total de 14 puncte la acest laborator (4 sunt considerate bonus).

Pentru a rula scriptul aveți nevoie de Python 3 și de câteva biblioteci instalate, așa cum este specificat la începutul fișierului.

1. Instalați Python3 și pip pentru Python3

sudo apt-get install pip3

2. Instalați librtlsdr (pentru a putea folosi dispozitivul RTL-SDR)

Pe MAC OS:

 sudo port install rtl-sdr 

Sau pentru brew: https://macappstore.org/librtlsdr/

Pe Linux:

 sudo apt-get install  librtlsdr-dev 

3. Pe linux, instalati si libportaudio2

sudo apt-get install libportaudio2

4. Dacă vreți să folosiți un mediu izolat pentru Python3 puteți instala virtualenv:

 pip3 install virtualenv 

Apoi în folderul în care doriți să creați mediul izolat dați comanda:

 virtualenv . 

Pentru activarea mediului virtual dați comanda din folderul respectiv:

 source bin/activate 

5. Instalați matplotlib, scipy, sounddevice, ipython pentru Python3:

pip3 install pyrtlsdr matplotlib scipy ipython sounddevice

6. Tot pentru Linux, posibil sa trebuiasca sa adugati dispozitivul in udev, asa cum este descris aici. Pe scurt, ca root faceti un fisier in /etc/udev/rules.d/20.rtlsdr.rules care contine urmatoarea linie:

SUBSYSTEM=="usb", ATTRS{idVendor}=="0bda", ATTRS{idProduct}=="2838", GROUP="adm", MODE="0666", SYMLINK+="rtl_sdr"

Apoi scoateti si reinserati dispozitivul RTL in portul USB.

Acum laboratorul efectiv… Aici veti esantiona un semnal de pe o anumita banda FM (e.g. Trinitas FM, dar puteti schimba cu alta daca vreti) si apoi veti folosi diverse tehnici/metode discutate la cursurile si laboratoarele precedente: decimare, FFT, spectrograma, etc.

Urmariti codul de mai jos cu atentie, pe rand, fara graba. Multe lucruri sunt deja facute pentru voi, ca sa puteti mai degraba sa vizualizati/auziti rezultatele. Mai sunt cateva locuri de completat, care ar trebui sa fie simple.

Dacă folosiți un dispozitiv RTL-SDR, atunci setați use_sdr = 1, altfel puneți use_sdr = 0. Dacă nu folosiți dispozitivul RTL-SDR, de asemenea comentați linia care încarcă biblioteca: #from rtlsdr import RtlSdr

lab_rtl_sdr.py
# -*- coding: utf-8 -*-
"""
    Laborator PS cu RTL-SDR
    Cod initial de Matei Simtinica
    Unele modificari de Marios Choudary
 
    Nota: trebuie instalate cateva lucruri pentru ca acest script sa functioneze:
 
    1. Python3 și pip pentru Python3 (acest script se bazeaza pe folosirea Python3)
 
    2. Install librtlsdr (pentru a putea folosi dispozitivul RTL-SDR)
 
    Pe MAC OS:
    sudo port install rtl-sdr
 
    Sau pentru brew:
    https://macappstore.org/librtlsdr/
 
    Pe Linux
    sudo apt-get install  librtlsdr-dev
 
    3. Install matplotlib, scipy, sounddevice pentru Python3:
    pip3 install pyrtlsdr matplotlib scipy sounddevice
 
"""
 
from rtlsdr import RtlSdr # comentati asta daca nu folositi dispozitiv RTL
import matplotlib.pyplot as plt
import numpy as np 
import scipy.signal as signal
#from IPython.display import Audio
from scipy.fftpack import fft
import sounddevice as sd
from scipy.io.wavfile import write
 
plt.rcParams['figure.dpi'] = 170
plt.rcParams['figure.figsize'] = (8, 4)
 
# Setati aici 1 daca folositi un RTL-SDR, 0 altfel
use_sdr = 0
 
F_station = int(88.5e6) # Trinitas FM
Fs = 2280000            # frecventa de esantionare (samplerate)
N = 8192000             # numarul esantioanelor
 
if use_sdr == 0:
    x1 = np.load('x1.npy')
else:
    # initializati device-ul RTL-SDR
    sdr = RtlSdr()
 
    # setati parametri device-ului
    sdr.sample_rate = Fs
    sdr.center_freq = F_station
    sdr.gain = 'auto'
 
    # cititi semnalul brut
    x1 = sdr.read_samples(N)
    np.save('x1.npy', x1)
 
    # eliberati resursele device-ului
    sdr.close()
 
# plotati spectrograma semnalului
# Pentru detalii, vedeti aici: https://matplotlib.org/stable/api/_as_gen/matplotlib.pyplot.specgram.html
plt.figure()
plt.specgram(x1, NFFT=2**10, Fs=Fs)
plt.title('Spectrograma semnalului (X1)')
plt.ylim(-Fs/2, Fs/2)
plt.ticklabel_format(style='plain')
plt.show(block=False)
# sau salvati imaginea ca png, daca nu va merge plt.show
# plt.savefig('fig1.png')
 
# semnalele FM sunt difuzate pe o latime de banda de 200kHz
f_bw = 200000
 
# filtrarea semnalului
# TODO 1: creati un filtru trece-jos digital folosind scipy.signal.firwin
# 1. Selectati un numar de coeficienti (taps) ai filtrului (incercati 8, 16, 32)
# 2. Calculati frecventa de taiere (cutoff), ca raport intre
#    frecventa maxima de interes (e.g. f_bw) si fs/2
# 3. alegeti o fereastra pentru atenuarea DFT leakage rezultat din filtru (e.g. 'hamming')
# 4. folositi functia firwin din scipy.signal pentru a obtine coeficientii unui filtru FIR
# 5. filtrati semnalul initial (x1) cu coeficientii obtinuti anterior, folosind functia lfilter din scipy.signal
ntaps = 0 # incercati mai multe variante: 8, 16, 32, vedeti apoi diferentele, prin spectrograma urmatoare
fcut = 0 # calculati cat este frecventa de cutoff necesara, ca o valoare intre 0 si 1, unde 1 corespunde lui fs/2
b = 0 # folositi signal.firwin pentru a obtine coeficientii (b) ai unui filtru trece-jos:
      # https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.firwin.html#scipy.signal.firwin
x1f = signal.lfilter(b, 1, x1)
 
# plotati spectrograma semnalului dupa filtrare
# TODO 2: faceti spectrograma, ca mai sus, dar pentru semnalul x1f
# Observati diferentele intre spectrograma semnalului initial si cel filtrat
# pentru diferite valori ale numarului de coeficienti ai filtrului (ntaps)
 
 
# decimare
# TODO 3: decimati (faceti subsampling) semnalul, astfel incat sa aveti noua frecventa
# de esantionare f_bw   (i.e. fs' = f_bw)
dec_rate = 0 # calculati factorul de decimare (raportul intre frecvente)
x2 = signal.decimate(x1f, dec_rate)
 
# noua frecventa de esantionare
Fs_new = Fs / dec_rate
 
# plotati spectrograma semnalului decimat
# TODO 4: afisati spectrograma pentru semnalul obtinut dupa decimare
 
# demodulati semnalul FM
# TODO 5: demodulati semnalul FM prin folosirea unui discriminator simplu de frecvente
# 1. Folosind semnalul decimat (x2), obtineti o copie a acestuia dar intarziat cu un element (x2_int = x2[1:])
# 2. Calculati conjugatul acestui semnal intarziat (puteti folosi numpy.conj)
# 3. Inmultiti semnalul decimat (x2) cu semnalul intarziat conjugat
# 4. calculati faza (unghiul) semnalului rezultat din inmultirea precedenta (puteti folosi numpy.angle)
# Aceasta faza va fi direct proportionala cu semnalul transmis, asa ca poate fi folosita direct.
x2_int = x2[1:]
x2_int_conj = 0 # calculati conjugatul cu numpy.conj
xx = x2[:-1] * x2_int_conj
x3 = 0 # calculati faza semnalului inmultit (xx) folosind numpy.angle
 
 
# vizualizati semnalul FM obtinut prin afisarea puterii spectrale (power spectral density)
# zonele colorate corespund componentelor relevante semnalului
# rosu -> audio mono
# portocaliu -> audio stereo
plt.figure()
plt.psd(x3, NFFT=2048, Fs=Fs_new, color='blue')
plt.title('Semnal FM decimat, apoi demodulat')
plt.axvspan(0,             15000,         color='red',    alpha=0.2)
plt.axvspan(19000-500,     19000+500,     color='green',  alpha=0.2)
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.show(block=False)
# sau salvati imaginea ca png, daca nu va merge plt.show
# plt.savefig('fig4.png')
 
# Calculati si afisati FFT pentru semnalului demodulat
# TODO 6:
# 1. Calculati FFT din x3 prin functia fft din scipy.fftpack
# 2. afisati valoarea absoluta a spectrului (np.abs) intre 0 si fs/2
# Nota: puteti folosi np.linspace pentru a crea un vector de frecvente intre 0 si fs, util la plot
yf = 0 # calculati aici spectrul lui x3 folosind scipy.fftpack.fft
nf = len(yf)
xf = np.linspace(0, Fs_new, nf)
plt.figure()
# plt.plot(...)
plt.grid()
plt.show(block=False)
# sau salvati imaginea ca png, daca nu va merge plt.show
# plt.savefig('fig5.png')
 
# extragerea semnalului audio (mono) prin filtrare si decimare
# TODO 7:
# 1. filtrati semnalul precedent (x3) pentru a retine doar componentele intre 0 si 15 kHz (corespunde semnalului mono)
# 2. Calculati si afisati FFT pentru semnalul filtrat si comparati cu FFT-ul semnalului precedent (nefiltrat)
# 3. decimati semnalul astfel incat sa aveti frecventa finala de esantionare f_audio = 44100 Hz pentru a-l putea asculta
ntaps2 = 32
fmaxa = 15000
fcut2 = 0 # calculati frecventa de cutoff necesara pentru a opri doar semnalele intre 0 si 15 kHz (tot in raport cu fs/2)
b2 = signal.firwin(ntaps2, fcut2, window='hamming')
x3f = signal.lfilter(b2, 1, x3)
 
yff = fft(x3f)
nff = len(yff)
nff2 = int(nff/2)
xff = np.linspace(0, Fs_new, nff)
plt.figure()
plt.plot(xff[0:nff2], 2.0/nff * np.abs(yff[0:nff2]))
plt.title('spectrul (FFT) semnalului FM demodulat dupa filtrare pe canalul mono')
plt.grid()
plt.show(block=False)
# sau salvati imaginea ca png, daca nu va merge plt.show
# plt.savefig('fig6.png')
 
f_audio = 44100
dec_audio = 0 # calculati factorul de decimare pentru a face subsampling de la Fs_new la f_audio
Fs_audio = Fs_new / dec_audio
xa = signal.decimate(x3f, dec_audio)
 
input("Press Enter to play audio from signal...")
 
# Redarea semnalului si salvarea lui ca wav
xs = np.int16(xa/np.max(np.abs(xa)) * 32767)
sd.play(xs, f_audio) # daca nu va merge, comentati linia
write('x1.wav', f_audio, xs) # si ascultati intr-un media player fisierul acesta
 
print("Am terminat!\n")
input("Press Enter to continue...")
ps/labs/08.txt · Last modified: 2023/12/05 22:16 by darius.necula
CC Attribution-Share Alike 3.0 Unported
www.chimeric.de Valid CSS Driven by DokuWiki do yourself a favour and use a real browser - get firefox!! Recent changes RSS feed Valid XHTML 1.0