Show page

Differences

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

--- ps:labs:08 [2020/10/28 18:58]
ionut.gorgos
+++ ps:labs:08 [2023/12/05 22:16] (current)
darius.necula
@@ Line 1: / Line 1: @@
 ===== Laboratorul 08. =====
 /*<hidden>*/
-==== DFT în detaliu: DFT leakage, zero-padding ====
+În acest laborator veți face diverse procesări pe un semnal obținut de dispozitivul
+[[https://www.rtl-sdr.com/|RTL-SDR]].
-În acest laborator vom continua să explorăm Transformata Fourier Discretă (DFT), urmărind efectul eșantionării în domeniul frecvență (apariția sinc-ului din cauza fenomenului de leakage) și metode de rezolvare a acestuia (zero-padding, ferestre, creșterea numărului de eșantioane).
+Găsiți [[https://www.desktopsdr.com/dl/book-pdf/|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:
+{{:ps:labs:screenshot_2021-12-04_at_23.24.06.png?direct&800|}}
-=== Exercițiul 1 -- DFT leakage și zero-padding [4p] ===
+Spectrul unui semnal FM arată în felul următor:
-În acest exercițiu veți reface exemplul pe care l-am făcut la curs cu 2 sinusoide pentru a vedea efectul fenomenului de leakage și a experimenta zero-padding. Pentru asta vom folosi următorul semnal:
+{{:ps:labs:fm-spectrum.png?800|}}
-$s(n) = A_1 \sin(2\pi f_1 n t_s) + A_2 \sin(2\pi f_2 n t_s)$,
+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).
+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).
-unde $f_1$ si $f_2$ sunt frecvențele celor două sinusoide care compun semnalul, și $t_s = 1/f_s$ e perioada de eșantionare ($f_s = 1/t_s$).
+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.
-Pentru a face asta urmăriți următorii pași:
+. Instalați Python3 și pip pentru Python3
-  - Creați si plotați semnalul, folosind $A_1=1$, $A_2=0.5$, $f_s = 8000$ Hz, $f_1 = 1000$ Hz, $f_2 = 2000$ Hz pentru $N=8$ eșantioane.
-  - Calculați DFT pentru acest semnal și plotați magnitudinea acesteia, ca în laboratoarele anterioare. Ar trebui să obțineti ceva de genul următor:
-{{:ps:labs:lab7_sinewaves_a.png?300|}} {{:ps:labs:lab7_sinewaves_a_fft.png?300|}}
-  - Apoi eliminați prima sinusoidă (ex.: făcând $f_1=0$) și verificați dacă aveți semnal doar la 2kHz.
-  - Schimbați $f_2$ de la 2kHz la $f_2=2500$ Hz. Plotați spectrul. Ce putem observa? Ar trebui să vedeti că toată energia de la frecvența de 2.5kHz a fost distribuită pe frecvențele adiacente. După cum am învățat la curs, acesta este fenomenul cunoscut ca "DFT leakage", și apare din cauza faptului că folosirea unui număr finit de eșantioane poate fi modelată ca înmulțirea unui semnal infinit eșantionat cu o funcție rectangulară (al cărui spectru este un sinc). Înainte, nu vedeam acest efect pentru că eșantioanele sinc-ului erau exact în punctele unde sinc-ul era 0.
-  - Pentru a vedea mai bine efectul de leakage trebuie să creștem numărul de eșantioane folosite pentru DFT. Pentru asta adăugați zerouri semnalului vostru. De exemplu adăugati 56 de zerouri ca să obțineți un total de $K=64$ eșantioane (din care doar $N=8$ sunt diferite de 0). Apoi calculați DFT pentru acest semnal. Ar trebui să vedeti forma sinc-ului mult mai clară și de asemenea că e centrată în jurul frecvenței semnalului (2.5kHZ).
-  - Acum schimbati din nou frecvența la $f2=2000$ Hz, dar folosind în continuare zero-padding și plotați DFT. Ar trebui să vedeți că într-adevăr sinc-ul era acolo, dar eșantioanele de la $1000, 3000, 4000 \ldots$ erau 0.
-=== Exercițiul 2 -- DFT leakage și ferestre [5p] ===
+<code>
+sudo apt-get install pip3
+</code>
-În acest exercițiu aveți dat un semnal({{:ps:labs:notes_signal.mat|click aici}}) care conține două note (două unde sinusoidale). Însă, una dintre ele este mult mai puternică decât cealaltă, așa că a doua, cea mai slabă, nu e ușor de detectat din spectrul semnalului. În acest exercițiu vom încerca să folosim o funcție fereastră pentru a determina cele două note.
+. Instalați librtlsdr (pentru a putea folosi dispozitivul RTL-SDR)
-Să facem următoarele:
+Pe MAC OS:
-  - Incărcați și plotați semnalul dat. Ar trebui să observați că se vor încărca variabilele 'notes_signal' și 'fs', unde fs este frecvența de eșantionare (amintiți-vă ca aveți nevoie de ea pentru a înțelege rezultatul dat de DFT).
+<code bash> sudo port install rtl-sdr </code>
-  - Calculați DFT pentru semnal și plotați magnitudinea, ca în laboratoarele precedente. Ar trebui să obțineți ceva precum aceasta: {{:ps:labs:notes_signal_fftpos.png?300|}}. În acest moment probabil nu puteți spune care sunt cele două frecvențe ale semnalului, din cauza faptului că funcția sinc a primei sinusoide acoperă componenta celei de-a doua sinusoide.
-  - Folosind zero-padding în acest caz nu va ajuta prea mult (încercați). Așa că vom aplica semnalului o funcție fereastră (ex. 'Hanning' sau 'Hamming'; căutați aceste funcții în Matlab folosind help), precum am discutat la curs. Ideea este să generăm o funcție fereastră pe care o vom înmulții cu semnalul original. Plotați semnalul după aplicarea funcției fereastră.
-  - Calculați DFT pentru semnalul rezultat după aplicarea funcției fereastră. Puteți spune, cel puțin aproximativ, care sunt cele două frecvente conținute de semnal?
-Semnalul, fereastra hanning și semnalul atenuat ar trebui să arate așa:
+Sau pentru brew:
+https://macappstore.org/librtlsdr/
-{{:ps:labs:lab08_notes_signal.png?300|}} {{:ps:labs:lab08_hanning_window.png?300|}}
+Pe Linux:
-{{:ps:labs:lab08_notes_signal_window.png?300|}}
+<code bash> sudo apt-get install  librtlsdr-dev </code>
+. Pe linux, instalati si libportaudio2
+<code>
+sudo apt-get install libportaudio2
+</code>
-=== Exercițiul 3 -- Tratarea DFT leakage prin creșterea numărului de eșantione [1p] ===
+. 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>
-Precum am discutat la curs, frecvențele date de sinc pot fi reduse prin creșterea numărului de eșantioane diferite de zero ale semnalului nostru. Așa că, atunci când este posibil, aceasta ne va ajuta să vizualizăm semnale foarte aproape în frecvență.
+. Instalați matplotlib, scipy, sounddevice, ipython pentru Python3:
+<code bash>
+pip3 install pyrtlsdr matplotlib scipy ipython sounddevice
+</code>
-Pentru a vedea acest efect, să utilizăm aceleași note, dar cu un semnal mult mai lung {{:ps:labs:notes_signal_long.mat|click aici}}.
+. 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.
-Procedați la fel ca înainte:
+Acum laboratorul efectiv... Aici veti esantiona un semnal de pe o anumita banda FM (e.g. Trinitas FM,
-  - Plotați semnalul și spectrul său. Verificați dacă puteți distinge cele două frecvențe(ar trebui).
+dar puteti schimba cu alta daca vreti) si apoi veti folosi diverse tehnici/metode discutate la
-  - Aplicați funcția fereastră și verificați spectrul. Ar trebui sa fie mult mai clar.
+cursurile si laboratoarele precedente: decimare, FFT, spectrograma, etc.
-  - Ce note muzicale reprezintă aceste frecvențe? Puteți să redați acest sunet folosind funcția Matlab 'sound'.
+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
+<file python 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:
+. Python3 și pip pentru Python3 (acest script se bazeaza pe folosirea Python3)
+. 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
+. 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...")
+</file>
+/*</hidden>*/

Laboratoare

Resurse

Tema

ps/labs/08.1603904329.txt.gz · Last modified: 2020/10/28 18:58 by ionut.gorgos

Show page Old revisions

Media Manager Back to top