Î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). Modul de funcționare al acestui dispozitiv poate fi observat în imaginea următoare preluată din cartea menționată:

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

Veți putea folosi un dispozitiv RTL-SDR dacă aveți/puteți împrumuta de la laborant, sau două fișiere pre-generate de noi, pe care le puteți descărca de 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. 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 și la începutul fișierului.

Dacă nu utilizați dispozitivul RTL-SDR, puteți sări peste pașii 1, 2, 3, 4 și 6 prezentați mai jos.

Pentru a realiza laboratorul local folosind dispozitivul RTL-SDR, urmați acești pași:

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 scoateți și reinserați dispozitivul RTL în portul USB.

Pe partea de cod veți esantiona un semnal de pe o anumită bandă FM, de exemplu Trinitas FM, (puteți schimba cu altă bandă) și apoi veți folosi diverse tehnici/metode discutate la cursurile și laboratoarele precedente: decimare, FFT, spectrograma, etc.

Urmăriți codul de mai jos cu atentie. Scheletul cuprinde o mare parte din rezolvare cu scopul de a ajunge mai ușor la rezultate și pentru a le analiza (diagramele, spectrele, semnalele, sunetele etc). Rezolvarea TODO-urilor reprezintă task-ul vostru pentru acest laborator pentru a ajunge la rezultate.

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
import math
import cmath
 
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 după downsampling intern)
N = 8192000             # numarul esantioanelor (t = N / Fs)
 
if use_sdr == 0:
    x1 = np.load('x1.npy')
else:
    # initializati device-ul RTL-SDR
    # inlocuiti linia cu sdr = 0 in cazul in care nu folositi dispozitivul, iar pachetul rtlsdr nu se instaleaza cum trebuie
    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 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
    # inlocuiti linia cu sdr = 0 in cazul in care nu folositi dispozitivul, iar pachetul rtlsdr nu se instaleaza cum trebuie
    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(x2, NFFT=2**10, Fs=Fs)
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 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)
# 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 = 8 # incercati mai multe variante: 8, 16, 32, vedeti apoi diferentele, prin spectrograma urmatoare
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
# 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)
 
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) - important: trebuie să fie un întreg (int(...))
xd = 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
 
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 (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 (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.
xd_int = xd[1:]
xd_int_conj = 0 # calculati conjugatul cu numpy.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)
# 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
y3 = 0 # calculati aici spectrul lui x3 folosind scipy.fftpack.fft
y3 = fft(x3)
nf3 = len(y3)
nf3_2 = int(nf3/2)
xx3 = np.linspace(0, Fs, nf3)
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
# 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 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)
 
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 (tot un întreg)
dec_audio = int(Fs_new / 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 (ar trebui să se înțeleagă ceva)
 
# TODO Bonus: încercați să îmbunătățiți calitatea semnalului!
 
print("Ați terminat! Felicitări!\n")
input("Press Enter to continue...")

În acest laborator vom face câteva exerciții pentru a ne familiariza cu operația de convoluție precum și cu filtre cu răspunsul finit la impuls (finite impulse response - FIR) și cu metode de proiectare a acestora prin metode folosind ferestre.

Materiale utile:

• Vedeți slide-urile de la R. Lyons aici

[4p]

Precum am spus la curs, putem defini operația de convoluție dintre două secvențe h(k) și x(n) după cum urmează:

$y(n) = h(k) \ast x(n) = \sum_{k=0}^{M-1} h(k) \cdot x(n-k)$,

unde $M$ este lungimea secvenței h(k). De reținut că această operație definește un singur element de ieșire y(n). Pentru următorul element de ieșire, y(n+1), trebuie să shiftăm secvența h(k) astfel încât să se potrivească cu elementele [x(n+1), x(n), …, x(n-M+2)]. Vedeți slide-urile 5 și 6 aici.

În general presupunem că secvența h(k) are un număr $M$ finit și, în general, mic de elemente.

Rezolvați următoarele exerciții:

1. Dacă x(n) are N elemente și h(k) are M elemente, câte elemente are secvența obținută prin convoluție $h(k) \ast x(n)$? (presupunând că efectuăm convoluția doar pe elementele diferite de zero, adică atunci când elementele celei mai scurte secvențe se suprapun în totalitate cu elementele celeilalte)
2. Fie x(n) secvența [1, 3, 5, 7, 5, 4, 2] și h(k) secvența [0.1, 0.3, 0.1]. Secvența obținută prin convoluție $y(n) = h(k) \ast x(n), n \in \{1,2,3,4,5\}$ ar trebui să aibă 5 elemente. Scrieți fiecare dintre aceste elemente ca un produs scalar. Observați că trebuie să inversați ordinea secvenței x(n) (doar partea care se înmulțeste cu filtrul) înainte să o înmulțiți cu h(k) (alternativ, puteți inversa o singură dată elementele filtrului).
3. Fie x(n) o secvență de $N=64$ elemente corespunzătoare unei sinusoide de frecvență $f=3$ kHz, eșantionată cu $f_s=64$ kHz. Fie h(k) secvența [0.1, 0.2, 0.2, 0.2, 0.1]. Generați aceste secvențe în Python și implementați convoluția pentru a obține elementele $y(n) = h(k) \ast x(n)$. Plotați inputul x(n) și ieșirea y(n); folosiți funcția stem în locul funcției plot pentru acest exercițiu.
4. Înlocuiți x(n) cu secvența de impuls cu 9 elemente [0, 0, 0, 0, 1, 0, 0, 0, 0] și efectuați convoluția cu aceeași secvență h(k) ca mai sus. Ce obțineți ca y(n)? Cum se numește aceasta?
5. Încercați operațiile de mai sus folosind funcția np.convolve din NumPy sau signal.convolve din scipy.signal. Obțineți aceleași rezultate? Care sunt diferențele?
6. [BONUS] Teorema Convoluției. Convoluția în timp [ $h(k) \ast x(n)$ ] este echivalentă cu înmulțirea element cu element a spectrelor în frecvență [ $H(m) \cdot X(m)$ ] și invers. Pentru a testa asta faceți următorii pași:

• Folosiți secvența h(k) = [0.1, 0.2, 0.2, 0.2, 0.1], dar adăugați încă 59 de zerouri la finalul ei, precum am făcut în laboratoarele precedente.
• Aplicați FFT peste această secvență și afișați spectrul ei (cu stem).
• Aplicați FFT peste sinusoida de la subpunctul 3 și afișați spectrul ei (cu stem).
• Înmulțiți element cu element cele două spectre și afișați rezultatul cu stem. S-a modificat ceva?
• Reveniți în domeniul timp folosind ifft și np.real.
• Comparați rezultatul de mai sus cu cel al funcției np.convolve din NumPy sau signal.convolve, setând parametrul mode='full' și păstrând doar primele N elemente. Obțineți aceleași rezultate? Care sunt diferențele?

[4p]

În acest exercițiu vom crea secvența x(n) pentru un filtru trece-jos plecând de la un filtru ideal în domeniul frecvență și, folosind IFFT, vom obține secvența x(n). Vom folosi de asemenea diferite ferestre pentru a le compara performanța în proiectarea de filtre trece-jos.

Încercați să urmați acești pași:

1. Generați o secvență de filtru ideal trece-jos HK având N = 256 elemente, reprezentând spectrul de frecvență al unui filtru trece-jos. Folosiți o frecvență de cut-off de $f_s/16$. Adică totul înainte de $f_s/16$ trebuie sa treacă, pe când totul după $f_s/16$ trebuie sa fie oprit (folosiți un dreptunghi care se oprește la $f_s/16$). Observați că trebuie să generați un spectru simetric pentru a obține o secvență reală la următorul pas. Plotați această secvență (folosind plot). Notați axa frecvențelor (axa x) ca o funcție de $F_s$, adică de la 0 la 1. Ar trebui să obțineți ceva precum:

• Țineți minte că acest spectru poate fi văzut ca ieșirea din DFT(FFT), adică primul element corespunde frecvenței 0, pe când următoarele N/2-1 corespund frecvențelor pozitive, iar ultimele N/2 componente reprezintă frecvențele negative.

1. Acum aplicați inversa DFT(în practică inversa FFT, ifft în Python) pentru a obține secvența corespunzătoare în domeniul timp hk(n). Rețineți: trebuie să aplicați funcția ifftshift pe rezultat pentru a obține o funcție sinc simetrică, adică folosiți ceva precum:

   hk = ifftshift(ifft(HK));

2. Trunchiați secvența hk(n) prin selectarea a doar L=65 de eșantioane din centru(32 din stânga maximului funcției sinc, maximul funcției, și 32 de eșantioane din dreapta). Aceasta corespunde multiplicării secvenței hk(n) cu o fereastră dreptunghiulară centrată în punctul maxim al funcției sinc. Plotați secvența.
3. Aplicați DFT(fft) pe secvența trunchiată, adică cea înmulțită cu fereastra dreptunghiulară (care conține doar 1) și plotați spectrul (cu plot). Rețineți: este important aici, precum și la primul plot pentru filtru trece-jos ideal, să notăm axa frecvențelor (axa x) ca o funcție de $F_s$, adică de la 0 la 1. Vedeți diferențe față de filtrul ideal trece-jos? Acestea sunt efectele ferestrei dreptunghiulare.
4. Folosiți aceeași secvență trunchiată (hk(n)), dar înmulțiți-o cu o fereastră precum Blackman (np.blackman în Python). Efectuați din nou DFT și plotați spectrul (cu plot). Arată mai bine?.
5. În final, folosiți ca intrare sinusoida din Exercițiul 1 ca x(n) și filtrați-o printr-o convoluție cu secvența obținută mai sus după folosirea ferestrei Blackman(folosiți funcția np.convolve din NumPy sau signal.convolve din scipy.signal). Plotați intrarea și ieșirea în aceeași figură folosind stem pentru a observa efectele filtrului.

[2p]

Putem descrie un filtru (sistem liniar) cu feedback (IIR având termenii $a_i$ mai jos) sau fără feedback (FIR) folosind o ecuație cu diferențe precum:

$y(n) = b_0 \cdot x(n) + b_1 \cdot x(n-1) + \ldots + b_q \cdot x(n-q) + a_1 \cdot y(n-1) + \ldots + a_p \cdot y(n-p)$

Putem reprezenta întârzierile $x(n-1)$ ca $z^{-1} \cdot x(n)$, unde $z=e^{j2\pi}$. Apoi, obținem o ecuație care depinde doar de $x(n)$ și $y(n)$ și obținem funcția de transfer a filtrului $H(z) = \frac{Y(z)}{X(z)}$ precum:

$H(z) = \frac{\sum_{k=0}^q b_k \cdot z^{-k}}{1 - \sum_{k=1}^p a_k \cdot z^{-k}}$

În Python, puteți folosi funcția signal.firwin din SciPy, pentru a obține coeficienții ($b_i$) ai unui filtru FIR trece-jos, trece-bandă sau trece-sus(ignorați coeficienții $a_i$ deocamdată). Apoi puteți folosi funcția signal.lfilter din SciPy, pentru a filtra orice secvență folosind coeficienții $b_i$ dați de signal.firwin, iar $a_0 = 1$.

Pentru acest exercițiu, trebuie să proiectați filtre FIR trece-jos (folosiți cutoff = 0.1), trece-bandă (folosiți cutoff = [0.2, 0.5]) și trece-sus (folosiți cutoff = 0.75). Puteți alege numtaps = 65. Apoi folosiți funcția signal.lfilter pentru a testa filtrele cu niște sinusoide ca în Exercițiul 1, cu $f=3$ kHz, $15$ kHz, $30$ kHz. Afișați cu stem, în subplot-uri sinusoidele inițiale și pe cele filtrate. Pentru răspunsul în frecvență, vă puteți folosi de următorul cod:

freq, H = signal.freqz(b_low, 1, fs=fs)
plt.figure()
plt.plot(freq, 20 * np.log10(abs(H)))
plt.xlabel('Frequency [Hz], from 0 to fs/2')
plt.ylabel('Amplitude [dB]')
plt.title('Digital Filter Frequency Response')
plt.grid()
plt.show()

, unde b_low sunt coeficienții $b_i$ ai unui filtru FIR trece-jos.

[BONUS]

În acest exercițiu, încercăm să urmăm pașii de la exercițiul 3, dar în MATLAB. Pentru a accesa MATLAB, puteți folosi MATLAB instalat pe PC-urile din laborator sau să folosiți MATLAB Online, din browser. Pentru MATLAB Online, trebuie să vă creați un cont cu adresa de e-mail de la facultate pe site-ul acesta.

În MATLAB, puteți folosi funcția fir1 pentru a obține rapid elementele $b_i$ ale unui filtru FIR trece-jos, trece-bandă sau trece-sus. Apoi puteți folosi funcția filter pentru a filtra orice secvență folosind coeficienții $b_i$ dați de fir1, iar $a_0 = 1$.

Pentru acest exercițiu, folosiți funcția fir1 pentru a proiecta filtre FIR trece-jos, trece-bandă și trece-sus. Apoi folosiți funcția filter pentru a testa filtrele cu aceleași secvențe ca în exercițiul precedent. Puteți verifica designul filtrelor folosind tool-ul fvtool.

De asemenea puteți folosi tool-ul MATLAB fdatool pentru a proiecta și analiza rapid performanțele filtrelor FIR și IIR:

1. Generați un filtru FIR trece-jos folosind o fereastră Kaiser, având Fs = 48000 Hz, Fc = 10000 Hz, cu N=10 coeficienți(ordin)
2. Generați un filtru IIR trece-jos (Butterworth sau Chebyshev type I) cu aceiași parametri.
3. Comparați amplitudinea răspunsului.
4. Încercați să creați un filtru FIR cu un răspuns similar cu al unui filtru IIR, prin creșterea numărului de coeficienți.

După proiectarea filtrelor precum ați dorit, le puteți salva (File→Generate MATLAB Code) și să le folosiți direct în alte scripturi MATLAB pentru a filtra diverse semnale.

[8p]

Pentru acest exercițiu vom utiliza metoda proiectării cu fereastră pentru a crea filtre FIR trece-bandă și trece-sus. Precum am văzut la curs, putem folosi același principiu pentru a crea filtre trece-jos, trece-bandă sau trece-sus. Tot ce trebuie să facem este să înmulțim coeficienții filtrelor (adică secvența $h_{ideal}$) cu valorile unei sinusoide de o anumită frecvență (centrul frecvențelor pentru filtrul trece-bandă).

Pentru a crea un filtru trece-bandă cu frecvența centrală $f_B$ ar trebui să procedați în felul următor:

• Generați filtrul în timp $h_{ideal}$ precum în laboratorul 9, ex.2 (creați filtrul în frecvență, treceți în domeniul timp, înmulțiți cu o fereastră precum Blackman sau alta)
• Înmulțiți $h_{ideal}$ element cu element cu secvența $\cos(\frac{2\pi f_B n}{f_s})$, unde $f_B$ este frecvența din centrul benzii dorite, iar $f_s$ este frecvența de eșantionare.

Câteva cazuri particulare:

• $f_B = \frac{f_s}{4}$, în acest caz secvența cosinus devine [1, 0, -1, 0, 1, 0, -1, 0, …]. Acesta este un tip de filtru eficient trece-bandă centrat în frecvența $f_s / 4$.
• $f_B = \frac{f_s}{2}$, în acest caz secvența cosinus devine [1, -1, 1, -1, …] și obținem un filtru trece-sus.

Acum că știți toate acestea (sperăm că ați reținut și de la curs), aveți de făcut următoarele:

1. Generați o secvență de filtru ideal trece-jos $H_{ideal}$ având N = 256 elemente, reprezentând spectrul de frecvență al unui filtru trece-jos. Folosiți o frecvență de cut-off de fs/16. Adică totul înainte de fs/16 trebuie să treacă, pe când totul mai sus trebuie sa fie oprit (folosiți un dreptunghi care se oprește la fs/16). Observați că trebuie să generați un spectru simetric pentru a obține o secvență reală la următorul pas. Plotați această secvență (folosind plot). Notați axa frecvențelor (axa x) ca o funcție de $F_s$, adică de la 0 la 1. Ar trebui să obțineți ceva precum:
2. Acum aplicați inversa DFT (în practică inversa FFT) pentru a obține secvența corespunzătoare în domeniul timp $h_{ideal}$.
3. Trunchiați secvența $h_{ideal}$ prin selectarea a doar $L = 65$ de eșantioane din centru (32 din stânga maximului funcției sinc, maximul funcției, și 32 de eșantioane din dreapta). Aceasta corespunde multiplicării secvenței $h_{ideal}$ cu o fereastră dreptunghiulară centrată în punctul maxim al funcției sinc. Plotați secvența.
4. Aplicați DFT($fft$) pe secvența trunchiată înmulțită cu fereastra dreptunghiulară (care conține doar 1) și plotați spectrul. Rețineți: este important aici, precum și la primul plot pentru filtru trece-jos ideal, să notăm axa frecvențelor (axa x) ca o funcție de Fs, adică de la 0 la 1.
5. Folosiți aceeași secvență trunchiată ca mai sus, dar înmulțiți-o cu o fereastră precum $Blackman$. Aplicați din nou $fft$ și plotați spectrul.
6. Din această secvență puteți obține o secvență corespunzătoare unui filtru trece-bandă cu $f_B = \frac{f_s}{4}$. Afișați spectrul.
7. Obțineți o secvență corespunzătoare unui filtru trece-sus cu $f_B = \frac{f_s}{2}$. Afișați spectrul.
8. Generați trei sinusoide cu frecvențe diferite (ex: $ f = 3 kHz$, $15 kHz$, $30 kHz$, cu $f_s = 64000$ și $N = 64$) și filtrați-le (folosind funcția np.convolve din NumPy sau signal.convolve din scipy.signal) cu filtrele trece-sus și trece-bandă obținute mai sus. Plotați atât input-ul cât și output-ul în același plot folosind stem pentru a observa efectele filtrelor.

[2p]

Putem descrie un filtru (sistem liniar) cu feedback (IIR având termenii $a_i$ mai jos) sau fără feedback (FIR) folosind o ecuație cu diferențe precum:

$y(n) = b_0 \cdot x(n) + b_1 \cdot x(n-1) + \ldots + b_q \cdot x(n-q) + a_1 \cdot y(n-1) + \ldots + a_p \cdot y(n-p)$

Putem reprezenta întârzierile $x(n-1)$ ca $z^{-1} \cdot x(n)$, unde $z=e^{j2\pi}$. Apoi, obținem o ecuație care depinde doar de $x(n)$ și $y(n)$ și obținem funcția de transfer a filtrului $H(z) = \frac{Y(z)}{X(z)}$ precum:

$H(z) = \frac{\sum_{k=0}^q b_k \cdot z^{-k}}{1 - \sum_{k=1}^p a_k \cdot z^{-k}}$

În Python, puteți folosi funcția signal.butter din SciPy, pentru a obține coeficienții ($b_i$ și $a_i$) ai unui filtru IIR trece-jos, trece-bandă sau trece-sus. Apoi puteți folosi funcția signal.lfilter din SciPy, pentru a filtra orice secvență folosind coeficienții $b_i$ și $a_i$ dați de signal.butter.

Pentru acest exercițiu, trebuie să proiectați filtre IIR trece-jos (folosiți cutoff = 0.1), trece-bandă (folosiți cutoff = [0.2, 0.5]) și trece-sus (folosiți cutoff = 0.75). Puteți alege numtaps = 5. Apoi folosiți funcția signal.lfilter pentru a testa filtrele cu niște sinusoide ca în Exercițiul 1, cu $f=3$ kHz, $15$ kHz, $30$ kHz. Afișați cu stem, în subplot-uri sinusoidele inițiale și pe cele filtrate. Pentru răspunsul în frecvență, vă puteți folosi de următorul cod:

freq, H = signal.freqz(b_low, a_low, fs=fs)
plt.figure()
plt.plot(freq, 20 * np.log10(abs(H)))
plt.xlabel('Frequency [Hz], from 0 to fs/2')
plt.ylabel('Amplitude [dB]')
plt.title('Digital Filter Frequency Response')
plt.grid()
plt.show()

, unde b_low sunt coeficienții $b_i$ și a_low coeficienții $a_i$ ai unui filtru IIR trece-jos. Observați diferențe față de cele FIR?