Differences
This shows you the differences between two versions of the page.
|
ps:labs_python:06 [2023/10/19 13:26] constantin.savu1510 created |
ps:labs_python:06 [2023/11/07 01:14] (current) ionut.gorgos |
||
|---|---|---|---|
| Line 1: | Line 1: | ||
| ===== Laboratorul 06. ===== | ===== Laboratorul 06. ===== | ||
| - | <hidden> | + | /*<hidden>*/ |
| ==== SNR, decibeli și DFT ==== | ==== SNR, decibeli și DFT ==== | ||
| + | **<hidden>** | ||
| Prezentarea PowerPoint pentru acest laborator poate fi găsită aici: [[https://docs.google.com/presentation/d/1s85BCzqZ3ziRwyqL1ecqilqaMDBMRxq-/edit?usp=share_link&ouid=110538702824281541719&rtpof=true&sd=true|aici]] | Prezentarea PowerPoint pentru acest laborator poate fi găsită aici: [[https://docs.google.com/presentation/d/1s85BCzqZ3ziRwyqL1ecqilqaMDBMRxq-/edit?usp=share_link&ouid=110538702824281541719&rtpof=true&sd=true|aici]] | ||
| + | **</hidden>** | ||
| În acest laborator vom rezolva câteva exerciții legate de calcularea raportului semnal-zgomot(eng. signal-to-noise ratio - SNR), vom exprima acest raport în decibeli și vom căpăta experiență cu transformata Fourier discretă (DFT). Deși nu am discutat prea mult despre ea, în toate exercițiile care folosesc DFT vom utiliza funcția //fft// (Fast Fourier Transform) și inversa acesteia //ifft//. | În acest laborator vom rezolva câteva exerciții legate de calcularea raportului semnal-zgomot(eng. signal-to-noise ratio - SNR), vom exprima acest raport în decibeli și vom căpăta experiență cu transformata Fourier discretă (DFT). Deși nu am discutat prea mult despre ea, în toate exercițiile care folosesc DFT vom utiliza funcția //fft// (Fast Fourier Transform) și inversa acesteia //ifft//. | ||
| Line 19: | Line 20: | ||
| [<color red>6p</color>] | [<color red>6p</color>] | ||
| - | În acest exercițiu aveți dat un semnal cu zgomot, ({{:ps:labs:noisy_signal.mat|click aici}}), pentru care știți că informația de interes se află în primele 10 componente DFT, corespunzătoare, de fapt cu $k \in \{ -9, ..., 9\}$, adică componenta continuă / directă (DC: $k = 0$) și componentele frecvențelor pozitive $k = 1, ..., 9$ și negative $k = −1, ... , −9$. | + | În acest exercițiu aveți dat un semnal cu zgomot, ({{:ps:labs_python:lab6_noisy_signal.zip|click aici}}), pentru care știți că informația de interes se află în primele 10 componente DFT, corespunzătoare, de fapt cu $k \in \{ -9, ..., 9\}$, adică componenta continuă / directă (DC: $k = 0$) și componentele frecvențelor pozitive $k = 1, ..., 9$ și negative $k = −1, ... , −9$. |
| Mai jos puteți vedea semnalul cu zgomot precum și spectrul pozitiv( folosind doar $k = \{0, 1, ..., N/2 - 1\}$). | Mai jos puteți vedea semnalul cu zgomot precum și spectrul pozitiv( folosind doar $k = \{0, 1, ..., N/2 - 1\}$). | ||
| Line 31: | Line 32: | ||
| <code> | <code> | ||
| - | load('noisy_signal.mat') | + | npz = np.load('noisy_signal.npz') |
| + | noisy_signal = npz['noisy_signal'] | ||
| </code> | </code> | ||
| Acest semnal are $N = 128$ eșantioane și puteți considera că frecvența de eșantionare este $f_s = N = 128$ Hz. | Acest semnal are $N = 128$ eșantioane și puteți considera că frecvența de eșantionare este $f_s = N = 128$ Hz. | ||
| - | 2. Calculați DFT pentru semnalul zgomotos folosind formula de mai sus, apoi comparați rezultatul cu cel al funcției fft din MATLAB. Plotați spectrul, atât în forma dată de //fft//, cât și în forma centrată în zero, formată aplicând funcția //fftshift// pe rezultatul de la //fft//. Ne amintim că funcția //fft// conține coeficienții pentru frecvențele pozitive între $0$ și $N/2 - 1$, iar următorii coeficienți aparțin de fapt spectrului negativ, de la $-N/2$ la $-1$. [<color red>1p</color>] | + | 2. Calculați DFT pentru semnalul zgomotos folosind formula de mai sus, apoi comparați rezultatul cu cel al funcției fft din scipy.fft. Plotați spectrul, atât în forma dată de //fft//, cât și în forma centrată în zero, formată aplicând funcția //fftshift// pe rezultatul de la //fft//. Ne amintim că funcția //fft// conține coeficienții pentru frecvențele pozitive între $0$ și $N/2 - 1$, iar următorii coeficienți aparțin de fapt spectrului negativ, de la $-N/2$ la $-1$. [<color red>1p</color>] |
| 3. Calculați raportul semnal-zgomot în domeniul frecvență, ca raportul dintre puterea semnalului și puterea zgomotului. Pentru aceasta vom considera că semnalul util se află doar în componentele de frecvență $k \in \{-9, ..., 9\}$ pe când zgomotul se află în toate componentele de frecvență. Formula pentru calcularea puterii unui semnal este: $power = \frac{1}{N}\sum^{N-1}_{k = 0}{|S(k)|^2}$, unde N este numărul total de componente considerate în semnal. [<color red>1p</color>] | 3. Calculați raportul semnal-zgomot în domeniul frecvență, ca raportul dintre puterea semnalului și puterea zgomotului. Pentru aceasta vom considera că semnalul util se află doar în componentele de frecvență $k \in \{-9, ..., 9\}$ pe când zgomotul se află în toate componentele de frecvență. Formula pentru calcularea puterii unui semnal este: $power = \frac{1}{N}\sum^{N-1}_{k = 0}{|S(k)|^2}$, unde N este numărul total de componente considerate în semnal. [<color red>1p</color>] | ||
| Line 44: | Line 46: | ||
| 5. Din ieșirea funcției //fft// eliminați frecvențele din afara zonei de interes (faceți-le zero). Aveți grijă la indicii frecvențelor pozitive și negative și luați în considerare că primul element corespunde componentei directe, care nu are corespondent negativ. [<color red>1p</color>] | 5. Din ieșirea funcției //fft// eliminați frecvențele din afara zonei de interes (faceți-le zero). Aveți grijă la indicii frecvențelor pozitive și negative și luați în considerare că primul element corespunde componentei directe, care nu are corespondent negativ. [<color red>1p</color>] | ||
| - | 6. Pentru a reconstrui semnalul corespunzător acestui spectru, care are componente diferite de zero doar pentru frecvențele de interes, calculați transformata Fourier discretă inversă (//IDFT//) folosind formula de mai sus, iar apoi comparați rezultatul cu cel al funcției //ifft// din MATLAB. Plotați acest semnal și comparați-l cu semnalul original zgomotos. Astfel, practic ați implementat un tip foarte simplu de filtru trece-jos. [<color red>2p</color>] | + | 6. Pentru a reconstrui semnalul corespunzător acestui spectru, care are componente diferite de zero doar pentru frecvențele de interes, calculați transformata Fourier discretă inversă (//IDFT//) folosind formula de mai sus, iar apoi comparați rezultatul cu cel al funcției //ifft// din scipy.fft. Plotați acest semnal și comparați-l cu semnalul original zgomotos. Astfel, practic ați implementat un tip foarte simplu de filtru trece-jos. [<color red>2p</color>] |
| Aici aveți o imagine cu semnalul fără zgomot. | Aici aveți o imagine cu semnalul fără zgomot. | ||
| Line 57: | Line 59: | ||
| [<color red>4p</color>] | [<color red>4p</color>] | ||
| - | Acest exercițiu este similar cu cel precedent. Acum aveți dat o înregistrare a unei secvențe vorbite ({{:ps:labs:noisy_sound.mat|click aici}}) care a fost amestecată cu zgomot de frecvență înaltă. Trebuie sa aplicați pașii precedenți pentru a elimina zgomotul și pentru a determina secvența vorbită. | + | Acest exercițiu este similar cu cel precedent. Acum aveți dat o înregistrare a unei secvențe vorbite ({{:ps:labs_python:lab6_noisy_sound.zip|click aici}}) care a fost amestecată cu zgomot de frecvență înaltă. Trebuie sa aplicați pașii precedenți pentru a elimina zgomotul și pentru a determina secvența vorbită. |
| Line 72: | Line 74: | ||
| <code> | <code> | ||
| - | S = load('noisy_sound.mat'); | + | npz = np.load('noisy_sound.npz') |
| - | s = S.noisy_sound; | + | noisy_sound = npz['noisy_sound'] |
| - | fs = S.fs; | + | fs = npz['fs'] |
| </code> | </code> | ||
| Pentru a asculta sunetul stocat în vectorul //s//, folosiți codul următor: | Pentru a asculta sunetul stocat în vectorul //s//, folosiți codul următor: | ||
| <code> | <code> | ||
| - | sound(s) | + | import sounddevice as sd |
| + | import time | ||
| + | |||
| + | sd.play(s, fs) | ||
| + | time.sleep(T) # Ne asigură că putem auzi tot semnalul, unde T = durata semnalului s | ||
| + | sd.stop() | ||
| + | </code> | ||
| + | |||
| + | Pentru a salva într-un fișier sunetul stocat în vectorul //s//, folosiți codul următor: | ||
| + | <code> | ||
| + | from scipy.io.wavfile import write | ||
| + | |||
| + | sound = np.int16(s/np.max(np.abs(s)) * 32767) | ||
| + | write('s.wav', fs, sound) | ||
| </code> | </code> | ||
| Line 92: | Line 107: | ||
| Până acum am descompus semnale 1D în serii de sinusoide cu transformata Fourier. Însă putem aplica același procedeu și pentru semnale 2D, descompunându-le în sinusoide 2D. Putem să ne gândim la o sinusoidă 2D ca la un val. Cel mai comun semnal 2D este o imagine cu nivele de gri. | Până acum am descompus semnale 1D în serii de sinusoide cu transformata Fourier. Însă putem aplica același procedeu și pentru semnale 2D, descompunându-le în sinusoide 2D. Putem să ne gândim la o sinusoidă 2D ca la un val. Cel mai comun semnal 2D este o imagine cu nivele de gri. | ||
| - | Putem aplica transformata Fourier discretă pe imagini cu funcția fft2 din MATLAB. | + | Putem aplica transformata Fourier discretă pe imagini cu funcția fft2 din bibiloteca PYTHON scipy.fft. Pentru a aplica inversa acesteia vom folosi ifft2. Puteți folosi următorul cod pentru a importa fft2, repsectiv ifft2: |
| - | în acest exercițiu veți face următoarele: | + | <code> |
| - | - încărcați o imagine folosind funcția imread din MATLAB: ex: img = imread('peppers.png') | + | from scipy.fft import fft2 |
| - | - convertiți imaginea RGB în imagine gri și normalizați: img = double(rgb2gray(img))/255 | + | from scipy.fft import ifft2 |
| - | - aplicați fft2 pe imagine și obțineți un spectru 2D (o matrice) S(k) | + | </code> |
| - | - tăiați frecvențele joase din spectru obținând un nou spectru S1(k). Frecvențele joase 2D sunt cele de frecvență mică pe ambele axe, adică, pe rezultatul dat de fft sunt în colturile matricii. | + | |
| - | - tăiați frecvențele înalte din spectru obținând un nou spectru S2(k). Frecvențele înalte sunt formate din tot spectrul (toată matricea) mai puțin colțurile. | + | |
| - | - obțineți din fiecare spectru o imagine: reconstructed_img = ifft2(S1) | + | |
| - | - afișați cele 2 imagini cu funcția imshow() | + | În acest exercițiu veți face următoarele: |
| + | - Încărcați și afișați o imagine folosind funcția imread din matplotlib.image: ex:<code> | ||
| + | import matplotlib.image as image | ||
| + | import matplotlib.pyplot as plt | ||
| + | |||
| + | |||
| + | img = image.imread("peppers.png") | ||
| + | plt.figure() | ||
| + | plt.imshow(img, cmap=plt.get_cmap('gray'), vmin=0, vmax=1) | ||
| + | plt.title("Original Image") | ||
| + | plt.show() | ||
| + | </code> | ||
| + | - Convertiți imaginea RGB în imagine gri și normalizați (dacă este nevoie). Pentru a converti imaginea în gri puteți folosi următoarea funcție:<code> | ||
| + | def rgb2gray(rgb): | ||
| + | return np.dot(rgb[..., :3], [0.2989, 0.5870, 0.1140]) | ||
| + | </code> | ||
| + | - Aplicați fft2 pe imagine și obțineți un spectru 2D (o matrice) S(k) | ||
| + | - Tăiați frecvențele joase din spectru obținând un nou spectru S1(k). Frecvențele joase 2D sunt cele de frecvență mică pe ambele axe, adică, pe rezultatul dat de fft sunt în colturile matricii. | ||
| + | - Tăiați frecvențele înalte din spectru obținând un nou spectru S2(k). Frecvențele înalte sunt formate din tot spectrul (toată matricea) mai puțin colțurile. Hint: Vă puteți folosi și de spectrul calculat anterior pentru a determina mai ușor S2(k) | ||
| + | - Obțineți din fiecare spectru o imagine: reconstructed_img = ifft2(S1) | ||
| + | - Afișați cele 2 imagini cu funcția <code> | ||
| + | plt.imshow(img1, cmap=plt.get_cmap('gray'), vmin=0, vmax=1)</code> | ||
| Folosiți această imagine: | Folosiți această imagine: | ||
| Line 113: | Line 149: | ||
| {{:ps:labs:lab07_imagine_high_pass.png?300|}} | {{:ps:labs:lab07_imagine_high_pass.png?300|}} | ||
| - | </hidden> | + | /*</hidden>*/ |
