Show page

Differences

This shows you the differences between two versions of the page.

--- ps:labs:08 [2021/12/05 23:40]
ionut.gorgos
+++ ps:labs:08 [2024/11/20 16:55] (current)
marios.choudary
@@ Line 1: / Line 1: @@
 ===== Laboratorul 08. =====
+/*<hidden>*/
 În acest laborator veți face diverse procesări pe un semnal obținut de dispozitivul
@@ Line 7: / Line 8: @@
 {{:ps:labs:screenshot_2021-12-04_at_23.24.06.png?direct&800|}}
-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  {{:ps:labs:x1.npy.zip|aici}} (câteva secunde de pe Radio Trinitas -- 88.5MHz -- trebuie dezarhivat înainte de folosire).
+Spectrul unui semnal FM arată în felul următor:
+{{:ps:labs:fm-spectrum.png?800|}}
+Veți putea folosi un dispozitiv RTL-SDR dacă aveți, sau două fișiere pre-generate de noi, pe care le puteți descărca de
+{{:ps:labs:x1_2.npy.zip|aici}} (câteva secunde de pe Radio Trinitas, bandă centrată pe frecvența bună de 88.5MHz cât și
+un semnal centrat pe o frecvență alăturată -- trebuie dezarhivat înainte de folosire).
 Pentru a urma acest laborator, folosiți scheletul de cod de mai jos, din Python.
@@ Line 15: / Line 22: @@
 . Instalați Python3 și pip pentru Python3
+<code>
+sudo apt-get install pip3
+</code>
 . Instalați librtlsdr (pentru a putea folosi dispozitivul RTL-SDR)
 Pe MAC OS:
-sudo port install rtl-sdr
+<code bash> sudo port install rtl-sdr </code>
 Sau pentru brew:
@@ Line 25: / Line 36: @@
 Pe Linux:
-sudo apt-get install  librtlsdr-dev
+<code bash> sudo apt-get install  librtlsdr-dev </code>
+. Pe linux, instalati si libportaudio2
+<code>
+sudo apt-get install libportaudio2
+</code>
+. Dacă vreți să folosiți un mediu izolat pentru Python3 puteți instala virtualenv:
+<code bash> pip3 install virtualenv </code>
+Apoi în folderul în care doriți să creați mediul izolat dați comanda:
+<code bash> virtualenv . </code>
+Pentru activarea mediului virtual dați comanda din folderul respectiv:
+<code bash> source bin/activate </code>
+. Instalați matplotlib, scipy, sounddevice, ipython pentru Python3:
+<code bash>
+pip3 install pyrtlsdr matplotlib scipy ipython sounddevice
+</code>
+. Tot pentru Linux, posibil sa trebuiasca sa adugati dispozitivul in udev,
+asa cum este descris [[https://www.instructables.com/rtl-sdr-on-Ubuntu/|aici]].
+Pe scurt, ca root faceti un fisier in /etc/udev/rules.d/20.rtlsdr.rules care contine urmatoarea linie:
+<code>
+SUBSYSTEM=="usb", ATTRS{idVendor}=="0bda", ATTRS{idProduct}=="2838", GROUP="adm", MODE="0666", SYMLINK+="rtl_sdr"
+</code>
+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.
-. Instalați matplotlib, scipy, sounddevice pentru Python3:
+Urmariti codul de mai jos cu atentie, pe rand, fara graba. Multe lucruri sunt deja facute pentru voi,
-pip3 install matplotlib scipy sounddevice
+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.
@@ Line 58: / Line 99: @@
 . Install matplotlib, scipy, sounddevice pentru Python3:
-    pip3 install matplotlib scipy sounddevice
+    pip3 install pyrtlsdr matplotlib scipy sounddevice
 """
@@ Line 64: / Line 105: @@
 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
@@ Line 70: / Line 111: @@
 import sounddevice as sd
 from scipy.io.wavfile import write
+import math
+import cmath
 plt.rcParams['figure.dpi'] = 170
@@ Line 78: / Line 121: @@
 F_station = int(88.5e6) # Trinitas FM
-Fs = 2280000            # frecventa de esantionare (samplerate)
+Fs = 2280000            # frecventa de esantionare (samplerate după downsampling intern)
-N = 8192000             # numarul esantioanelor
+N = 8192000             # numarul esantioanelor (t = N / Fs)
 if use_sdr == 0:
@@ Line 98: / Line 141: @@
     # eliberati resursele device-ului
     sdr.close()
+# plotati spectrograma semnalului cu specgram
+# Pentru detalii, vedeti aici: https://matplotlib.org/stable/api/_as_gen/matplotlib.pyplot.specgram.html
+fig, ax = plt.subplots()
+ax.specgram(x1, NFFT=2**10, Fs=Fs)
+ax.set_title('Spectrograma semnalului (X1)')
+ax.set_ylim(-Fs/2, Fs/2)
+ax.ticklabel_format(style='plain')
+ax.set_xlabel('Timp (s)')
+ax.set_ylabel('Frecvență')
+plt.show(block=False)
+# sau salvati imaginea ca png, daca nu va merge plt.show
+# plt.savefig('fig1.png')
+print('Apasati o tastă pentru a continua...\n')
+input()
+# Acum să luăm eșantioane de la o frecvență puțin depărtată de centrul frecvenței anterioare
+# Acest sceariu seamănă cu cazul în care facem întâi o eșantionare intermediară (IF=intermediate frequency)
+# și apoi aplicăm downconvert ca să ajungem în baseband (ca să putem face filtrare și downsampling ușor)
+F_station_bad = int(88.3e6) # La 200 kHz de Trinitas FM (sau la o distanță față o altă stație)
+if use_sdr == 0:
+    x2 = np.load('x2.npy')
+else:
+    # initializati device-ul RTL-SDR
+    sdr = RtlSdr()
+    # setati parametri device-ului
+    sdr.sample_rate = Fs
+    sdr.center_freq = F_station_bad
+    sdr.gain = 'auto'
+    # cititi semnalul brut
+    x2 = sdr.read_samples(N)
+    np.save('x2.npy', x2)
+    # 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.specgram(x2, NFFT=2**10, Fs=Fs)
-plt.title('Spectrograma semnalului (X1)')
+plt.title('Spectrograma semnalului (X2)')
 plt.ylim(-Fs/2, Fs/2)
 plt.ticklabel_format(style='plain')
 plt.show(block=False)
+plt.xlabel('Timp')
+plt.ylabel('Frecvență')
+# sau salvati imaginea ca png, daca nu va merge plt.show
+# plt.savefig('fig2.png')
+# Q1: ce observați între cele două spectrograme (pt x1 și pt x2) ? Care este diferența ? De ce ?
+print('Apasati o tastă pentru a continua...\n')
+input()
+# Acum să vedem spectrul celor două semnale cu FFT (varianta rapidă pentru DFT):
+y1 = fft(x1)
+nf1 = len(y1)
+nf1_2 = int(nf1/2)
+xf1 = np.linspace(0, Fs, nf1)
+y2 = fft(x2)
+nf2 = len(y2)
+nf2_2 = int(nf2/2)
+xf2 = np.linspace(0, Fs, nf2)
+fig, (ax1, ax2) = plt.subplots(2)
+fig.suptitle('Spectrul (FFT) pentru x1 (sus) și x2 (jos)')
+ax1.plot(xf1[0:nf1_2], np.abs(y1[0:nf1_2]))
+ax1.grid()
+ax2.plot(xf2[0:nf2_2], np.abs(y2[0:nf2_2]))
+ax2.set_xlabel('Frecvența')
+ax2.grid()
+plt.show(block=False)
+# sau salvati imaginea ca png, daca nu va merge plt.show
+# plt.savefig('fig6.png')
+print('Apasati o tastă pentru a continua...\n')
+input()
+# Să folosim IQ downconvert pentru a aduce semnalul util din x2 în origine (baseband)
+# Ar trebui să obținem în cele din urmă un semnal/spectru asemănător cu cel pentru x1
+# Reminder: pentru a coborâ spectrul din f_c în 0 trebuie să înmulțim în timp cu o
+# exponențială de frecvență f_c (în domeniul digital este raportul f_c/fs)
+F_c = F_station - F_station_bad # sau puteți să puneți manual uitându-vă la spectru
+lenx2 = len(x2)
+x2n = np.zeros(lenx2, dtype=np.complex128)
+fr = F_c / Fs
+# TODO: înmulțiți aici fiecare element  x2[i] cu exponențiala e^(-j*2*pi*fr*i)
+# pentru exponențială complexă folosiți cmath.exp
+for i in range(lenx2):
+    x2n[i] = 0
+# Acum să vedem spectrul celor două semnale cu FFT (x2 și x2 după downconvert):
+y2n = fft(x2n)
+nf2n = len(y2n)
+nf2n_2 = int(nf2n/2)
+xx2n = np.linspace(0, Fs, nf2n)
+fig, (ax1, ax2) = plt.subplots(2)
+fig.suptitle('Spectrul (FFT) pentru x2 (sus) și x2 după downconvert (jos)')
+ax1.plot(xf2[0:nf2_2], np.abs(y2[0:nf2_2]))
+ax1.grid()
+ax2.plot(xx2n[0:nf2n_2], np.abs(y2n[0:nf2n_2]))
+ax2.set_xlabel('Frecvența')
+ax2.grid()
+plt.show(block=False)
+# sau salvati imaginea ca png, daca nu va merge plt.show
+# plt.savefig('fig6.png')
+print('Apasati o tastă pentru a continua...\n')
+input()
+# TODO: Să facem și spectrograma pentru acest nou semnal (x2n) și să comparăm cu cea pentru x1
+plt.figure()
+plt.specgram(x2n, NFFT=2**10, Fs=Fs)
+plt.title('Spectrograma semnalului (X2) după downconvert')
+plt.ylim(-Fs/2, Fs/2)
+plt.ticklabel_format(style='plain')
+plt.show(block=False)
+plt.xlabel('Timp')
+plt.ylabel('Frecvență')
+# sau salvati imaginea ca png, daca nu va merge plt.show
+# plt.savefig('fig2.png')
+print('Apasati o tastă pentru a continua...\n')
+input()
+# Folosiți acest nou semnal în loc de x1 (BONUS: încercați lab cu ambele variante: cu x1 original și cu acesta):
+# TODO: scoați comentariul după ce ați implementat îmnulțirea de la IQ downconvert
+#x1 = x2n
 # semnalele FM sunt difuzate pe o latime de banda de 200kHz
 f_bw = 200000
-# filtrarea semnalului
+# filtrarea semnalului x1
 # 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)
@@ Line 119: / Line 288: @@
 # 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
+ntaps = 8 # 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
+fcut = 2*f_bw/Fs  # 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)
+# Să vedem spectrul semnalului filtrat
+y1f = fft(x1f)
+nf1f = len(y1f)
+nf1f_2 = int(nf1f/2)
+xx1f = np.linspace(0, Fs, nf1f)
+fig, ax = plt.subplots()
+fig.suptitle('Spectrul (FFT) pentru x1 filtrat')
+ax.plot(xx1f[0:nf1f_2], np.abs(y1f[0:nf1f_2]))
+ax.set_xlabel('Frecvența')
+ax.grid()
+plt.show(block=False)
 # plotati spectrograma semnalului dupa filtrare
@@ Line 130: / Line 313: @@
 # pentru diferite valori ale numarului de coeficienti ai filtrului (ntaps)
+fig, ax = plt.subplots()
+ax.specgram(x1f, NFFT=2**10, Fs=Fs)
+ax.set_title('Spectrograma semnalului filtrat (x1f)')
+ax.set_ylim(-Fs/2, Fs/2)
+ax.ticklabel_format(style='plain')
+ax.set_xlabel('Timp (s)')
+ax.set_ylabel('Frecvență')
+plt.show(block=False)
+print('Apasati o tastă pentru a continua...\n')
+input()
 # 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)
+dec_rate = 0 # calculati factorul de decimare (raportul intre frecvente) - important: trebuie să fie un întreg (int(...))
-x2 = signal.decimate(x1f, dec_rate)
+xd = signal.decimate(x1f, dec_rate)
 # noua frecventa de esantionare
@@ Line 142: / Line 337: @@
 # plotati spectrograma semnalului decimat
 # TODO 4: afisati spectrograma pentru semnalul obtinut dupa decimare
+fig, ax = plt.subplots()
+ax.specgram(xd, NFFT=2**10, Fs=Fs_new)
+ax.set_title('Spectrograma semnalului filtrat decimat (xd)')
+ax.set_ylim(-Fs_new/2, Fs_new/2)
+ax.ticklabel_format(style='plain')
+ax.set_xlabel('Timp (s)')
+ax.set_ylabel('Frecvență')
+plt.show(block=False)
+print('Apasati o tastă pentru a continua...\n')
+input()
 # 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:])
+# 1. Folosind semnalul decimat (xd), obtineti o copie a acestuia dar intarziat cu un element (xd_int = xd[1:])
 # 2. Calculati conjugatul acestui semnal intarziat (puteti folosi numpy.conj)
-# 3. Inmultiti semnalul decimat (x2) cu semnalul intarziat conjugat
+# 3. Inmultiti semnalul decimat (xd) 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:]
+xd_int = xd[1:]
-x2_int_conj = 0 # calculati conjugatul cu numpy.conj
+xd_int_conj = 0 # calculati conjugatul cu numpy.conj
-xx = x2[:-1] * x2_int_conj
+xx = xd[:-1] * xd_int_conj
 x3 = 0 # calculati faza semnalului inmultit (xx) folosind numpy.angle
 # vizualizati semnalul FM obtinut prin afisarea puterii spectrale (power spectral density)
@@ Line 168: / Line 374: @@
 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
@@ Line 174: / Line 382: @@
 # 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
+y3 = 0 # calculati aici spectrul lui x3 folosind scipy.fftpack.fft
-nf = len(yf)
+y3 = fft(x3)
-xf = np.linspace(0, Fs_new, nf)
+nf3 = len(y3)
-plt.figure()
+nf3_2 = int(nf3/2)
-# plt.plot(...)
+xx3 = np.linspace(0, Fs, nf3)
-plt.grid()
+fig, ax = plt.subplots()
+fig.suptitle('Spectrul (FFT) pentru semnalul demodulat')
+ax.plot(xx3[0:nf3_2], np.abs(y3[0:nf3_2]))
+ax.set_xlabel('Frecvența')
+ax.grid()
 plt.show(block=False)
+# sau salvati imaginea ca png, daca nu va merge plt.show
+# plt.savefig('fig5.png')
+print('Apasati o tastă pentru a continua...\n')
+input()
 # extragerea semnalului audio (mono) prin filtrare si decimare
@@ Line 189: / Line 406: @@
 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)
+fcut2 = 0 # calculati frecventa de cutoff pentru semnalele intre 0 si 15 kHz (tot in raport cu fs/2, atentie acum Fs_new)
 b2 = signal.firwin(ntaps2, fcut2, window='hamming')
 x3f = signal.lfilter(b2, 1, x3)
@@ Line 199: / Line 416: @@
 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.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
+dec_audio = 0 # calculati factorul de decimare pentru a face subsampling de la Fs_new la f_audio (tot un întreg)
+dec_audio = int(Fs_new / f_audio)
 Fs_audio = Fs_new / dec_audio
 xa = signal.decimate(x3f, dec_audio)
@@ Line 210: / Line 430: @@
 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)
+sd.play(xs, f_audio) # daca nu va merge, comentati linia
-write('x1.wav', f_audio, xs)
+write('x1.wav', f_audio, xs) # si ascultati intr-un media player fisierul acesta (ar trebui să se înțeleagă ceva)
-print("Am terminat!\n")
+# TODO Bonus: încercați să îmbunătățiți calitatea semnalului!
+print("Ați terminat! Felicitări!\n")
 input("Press Enter to continue...")
-</file>
+</file>
+/*</hidden>*/

Laboratoare

Resurse

Tema_2024

ps/labs/08.1638740437.txt.gz · Last modified: 2021/12/05 23:40 by ionut.gorgos

Show page Old revisions

Media Manager Back to top