Differences
This shows you the differences between two versions of the page.
|
ps:labs_python:05 [2023/10/23 01:01] ionut.gorgos |
ps:labs_python:05 [2023/10/31 10:14] (current) ionut.gorgos |
||
|---|---|---|---|
| Line 1: | Line 1: | ||
| ===== Laboratorul 05. ===== | ===== Laboratorul 05. ===== | ||
| - | <hidden> | + | /*<hidden>*/ |
| - | ==== Semnale digitale - procesarea simplă, eşantionarea și recuperarea semnalului ==== | + | ==== Semnale digitale - procesarea simplă, eşantionarea și modularea în amplitudine ==== |
| În acest laborator vom începe să experimentăm efectele semnalelor eşantionate şi procesarea semnalelor digitale. | În acest laborator vom începe să experimentăm efectele semnalelor eşantionate şi procesarea semnalelor digitale. | ||
| Line 11: | Line 11: | ||
| === Exerciţiul 1 -- Procesarea simplă === | === Exerciţiul 1 -- Procesarea simplă === | ||
| - | [<color red>4p</color>] | + | [<color red>5p</color>] |
| Vom încerca să implementăm un sistem de procesare digitală simplă, care funcţionează după cum urmează: | Vom încerca să implementăm un sistem de procesare digitală simplă, care funcţionează după cum urmează: | ||
| - | * luaţi 5 eşantioane (samples) din semnalul de intrare, $d_1, d_2, \ldots, d_5$. | + | * luaţi 5 eşantioane (samples) din semnalul de intrare, $x_1, x_2, \ldots, x_5$. |
| - | * faceţi o medie a celor 5 eşantioane (samples), $a = 0.2 \cdot (d_1 + d_2 + d_3 + d_4 + d_5)$ | + | * faceţi o medie a celor 5 eşantioane (samples), $a = 0.2 \cdot (x_1 + x_2 + x_3 + x_4 + x_5)$ |
| * puneţi la ieşire valoarea curentă $y(n) = a$ | * puneţi la ieşire valoarea curentă $y(n) = a$ | ||
| Urmăriţi aceşti paşi: | Urmăriţi aceşti paşi: | ||
| - | - Creaţi câteva secvenţe digitale care reprezintă sinusoide de diferite frecvenţe (1, 2, 10, 20, 100 Hz), la aceeaşi frecvenţă de eşantionare (folosiţi acelaşi număr de eşantioane, e.g. $N=128$). [<color red>1p</color>] | + | - Creaţi câteva secvenţe digitale care reprezintă sinusoide de diferite frecvenţe (1, 2, 10, 20, 100 Hz) și amplitudine 1, la aceeaşi frecvenţă de eşantionare (folosiţi acelaşi număr de eşantioane, e.g. $N=128$). [<color red>1p</color>] |
| - Implementaţi sistemul de procesare menţionat mai sus. [<color red>1p</color>] | - Implementaţi sistemul de procesare menţionat mai sus. [<color red>1p</color>] | ||
| - Creaţi secvenţele corespunzătoare de ieşire. [<color red>1p</color>] | - Creaţi secvenţele corespunzătoare de ieşire. [<color red>1p</color>] | ||
| - Plotaţi toate secvenţele de intrare/ieşire. [<color red>1p</color>] | - Plotaţi toate secvenţele de intrare/ieşire. [<color red>1p</color>] | ||
| - Ce fel de sistem de procesare este acesta ? | - Ce fel de sistem de procesare este acesta ? | ||
| - | - Cum ați putea implementa sistemul de procesare menționat mai sus, folosind funcția 'filter'? (opțional) | + | - Cum ați putea implementa sistemul de procesare menționat mai sus, folosind funcția [[https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.lfilter.html|lfilter]] din biblioteca //scipy//? (opțional) |
| + | - Adăugați la sinusoidele generate mai sus, un zgomot alb de deviație standard 0.1 și medie 0, folosind următorul cod: | ||
| + | <code python> | ||
| + | import numpy as np | ||
| + | noise = np.random.normal(0, 0.1, x.shape) | ||
| + | x = x + noise | ||
| + | </code>, unde x este semnalul sinusoidal. Plotați din nou toate secvenţele de intrare/ieşire.[<color red>1p</color>] | ||
| <note tip> | <note tip> | ||
| Pentru fiecare eşantion de intrare ($x(n)$) trebuie să produceţi o valoare de ieşire ($y(n)$). Pentru primele patru eşantioane nu aveţi suficiente valori de intrare diferite de zero, şi pentru simplitate ori setaţi ieşirea pe zero, ori faceţi media dintre 1, 2, 3, 4 eşantioane, în timp ce consideraţi ca celelalte să fie zero. Dar de la al 5-lea eşantion veţi avea câte 5 eşantioane la îndemână (actualul $x(n)$ şi anterioarele patru eşantioane). | Pentru fiecare eşantion de intrare ($x(n)$) trebuie să produceţi o valoare de ieşire ($y(n)$). Pentru primele patru eşantioane nu aveţi suficiente valori de intrare diferite de zero, şi pentru simplitate ori setaţi ieşirea pe zero, ori faceţi media dintre 1, 2, 3, 4 eşantioane, în timp ce consideraţi ca celelalte să fie zero. Dar de la al 5-lea eşantion veţi avea câte 5 eşantioane la îndemână (actualul $x(n)$ şi anterioarele patru eşantioane). | ||
| </note> | </note> | ||
| - | <note> Pentru crearea timpilor de eșantionare puteți folosi 0:perioada_esantionare:timp_maxim sau funcția linspace(0,timp_maxim,numar_esantioane)</note> | + | <note> Pentru crearea timpilor de eșantionare puteți folosi funcția //np.linspace(0,timp_maxim,numar_esantioane)// sau //np.arange//, cu pasul //1/frecventa_de_esantionare//. Pentru acest exercițiu puteți considera //timp_maxim = 1//.</note> |
| + | |||
| + | === Exercițiul 2 – modularea în amplitudine === | ||
| + | [<color red>5p</color>] | ||
| + | |||
| + | În acest exercițiu va trebui să încercați să efectuați modularea în amplitudine asupra următorului semnal exponențial (eng. exponential decay signal) : $s(t) = e^{-at}u(t)$, unde $a>0$ și $u(t)$ este treapta unitară (i.e. egală cu 1 pentru $t\ge0$, 0 altfel). | ||
| + | |||
| + | Obiectivul vostru: | ||
| + | |||
| + | Folosiți o frecvență purtătoare $f_c = \frac{20}{T}$, unde $T$ este numărul eșantioanelor (samples) și calculați cu ajutorul DFT (Discrete Fourier Transform) și FFT (Fast Fourier Transform) spectrul semnalelor. | ||
| + | |||
| + | Vom considera că DFT are același număr de componente ca și semnalul de intrare ($K = N$). Formulele pentru DFT și pentru inversa ei (IDFT) sunt următoarele: | ||
| + | \begin{equation} | ||
| + | DFT: S(k) = \sum^{N-1}_{n = 0}{s(n)e^{\frac{-j 2 \pi n k}{K}}}, k \in \{0, ..., K-1\} \\ | ||
| + | |||
| + | IDFT\ (inversa\ DFT): s(n) = \frac{1}{K}\sum^{K-1}_{k = 0}{S(k)e^{\frac{j 2 \pi n k}{K}}}, n \in \{0, ..., N-1\} | ||
| + | \end{equation} | ||
| + | |||
| + | Pentru asta va trebui să urmăriți acești pași: | ||
| + | * Creați semnalul $s(t)$ pentru $t\in\{1,\ldots,T=128\}$, folosind $a=0.05$. [<color red>1p</color>] | ||
| + | * Calculați și plotați (cu //stem//) spectrul folosind FFT (Transformata Fourier Rapidă), pe care nu am făcut-o încă la curs, dar o vom face în următoarele cursuri. Pentru moment, puteți folosi acest cod: [<color red>1p</color>] | ||
| + | <code python> | ||
| + | import numpy as np | ||
| + | from scipy.fft import fft, fftshift | ||
| + | import matplotlib.pyplot as plt | ||
| + | |||
| + | fx = np.linspace(-T//2, T//2 - 1, T) # pentru T par | ||
| + | spectru = fft(s) | ||
| + | plt.figure() | ||
| + | plt.plot(fx, np.abs(fftshift(spectru))) | ||
| + | plt.stem(fx, np.abs(fftshift(spectru)), basefmt=" ") | ||
| + | plt.xlabel('Frequency component (k)') | ||
| + | plt.ylabel('Magnitude of component') | ||
| + | plt.title('Fourier coefficients before amplitude modulation') | ||
| + | </code> | ||
| + | * Calculați și plotați (cu //stem//) spectrul cu ajutorul DFT, folosind formula de mai sus, apoi comparați rezultatul cu cel al funcției //fft//. [<color red>1p</color>] | ||
| + | * Modulați semnalul în amplitudine folosind frecvența purtătoare $f_c = \frac{20}{T}$, i.e. face 20 de perioade complete în $T=128$ eșantioane ale semnalului $s(t)$. O variantă simplă de modulare este să calculați: $x(t) = (1+s(t)) \cdot \cos(2\pi f_c t)$. [<color red>1p</color>] | ||
| + | * Calculați și plotați (cum am făcut mai devreme, cu funcția //fft//) spectrul semnalului modulat în amplitudine. Comparați-l cu spectrul semnalului original. Este ceea ce v-ați așteptat? [<color red>1p</color>] | ||
| + | |||
| + | <note> Atenție: Vom învăța la curs că spectrul obținut prin transformata Fourier Discretă este periodic. Funcțiile //fft// / //ifft// consideră primii jumătate plus unu coeficienți pentru frecvențele pozitive, apoi următoarea jumătate corespunzătoare coeficienților negativi. Pentru a obține un spectru centrat în zero (doar ca să îl vizualizăm precum ne-am obișnuit) vom folosi funcția fftshift. </note> | ||
| + | |||
| + | <hidden> | ||
| + | === Exercițiul 3 – modularea în amplitudine [<color red>3p</color>] === | ||
| + | |||
| + | În acest exercițiu veți încerca să refaceți un semnal modulat în amplitudine ({{ps:labs:lab05_modulated_sound.mat|click aici}}). | ||
| + | Urmați următorii pași: | ||
| + | * Descărcați semnalul modulat și încărcați-l în MATLAB folosind funcția: <code matlab>load('lab05_modulated_sound.mat')</code> Se va încărca semnalul, precum și Fs, frecvența de eșantionare. | ||
| + | * Semnalul modulat provine dintr-un semnal original $s(t)$ care a fost înmulțit cu un semnal cosinus de frecvență $fc$, cos($2 \pi \cdot fc \cdot t$), unde t=0:1/Fs:timp_maxim. Precum ați observat și în exercițiul anterior spectrul semnalului modulat conține spectrul semnalului original shiftat la frecvența fc și -fc. Plotați spectrul pentru a vedea acest lucru. | ||
| + | * Dacă înmulțim (din nou) semnalul primit cu același semnal cosinus ar trebui să putem recupera semnalul original după aplicarea unui filtru trece-jos. Folosiți fc = Fs / 8. | ||
| + | * Pentru a realiza efectul unui filtru trece-jos vom egala cu zero coeficienții corespunzători frecvențelor peste un prag. | ||
| + | * Având noii coeficienți putem reface semnalul în timp folosind inversa Fourier. Realizați acest lucru folosind funcția //ifft//. | ||
| + | * Puteți asculta semnalul rezultat folosind funcția MATLAB //sound(semnal)//. Pentru a opri sunetul folosiți comanda <code matlab> clear sound </code> | ||
| + | |||
| + | <note> Atenție: Vom învăța la curs că spectrul obținut prin transformata Fourier Discretă este periodic. Funcțiile //fft// / //ifft// consideră primii jumătate plus unu coeficienți pentru frecvențele pozitive, apoi următoarea jumătate corespunzătoare coeficienților negativi. Pentru a obține un spectru centrat în zero (doar ca să îl vizualizăm precum ne-am obișnuit) vom folosi funcția fftshift. </note> | ||
| + | </hidden> | ||
| + | |||
| + | <hidden> | ||
| === Exercițiul 2 -- Subeşantionare (Bandpass sampling) === | === Exercițiul 2 -- Subeşantionare (Bandpass sampling) === | ||
| [<color red>3p</color>] | [<color red>3p</color>] | ||
| Line 90: | Line 152: | ||
| Pentru semnalele continue (analogice), puteți folosi pentru afișare //plot//, iar pentru cele digitale (eșantionate), //stem//. | Pentru semnalele continue (analogice), puteți folosi pentru afișare //plot//, iar pentru cele digitale (eșantionate), //stem//. | ||
| </note> | </note> | ||
| - | |||
| - | === Exercițiul 3 – modularea în amplitudine === | ||
| - | [<color red>2p</color>] | ||
| - | |||
| - | În acest exercițiu va trebui să încercați să efectuați modularea în amplitudine asupra următorului semnal exponențial (eng. exponential decay signal) : $s(t) = e^{-at}u(t)$, unde $a>0$ și $u(t)$ este treapta unitară (i.e. egală cu 1 pentru $t\ge0$, 0 altfel). | ||
| - | |||
| - | Obiectivul vostru: | ||
| - | |||
| - | Folosiți o frecvență purtătoare $f_c = \frac{20}{T}$, unde $T$ este numărul eșantioanelor (samples). | ||
| - | |||
| - | /*<hidden>*/ | ||
| - | Vom considera că DFT are același număr de componente ca și semnalul de intrare ($K = N$). Formulele pentru DFT și pentru inversa ei (IDFT) sunt următoarele: | ||
| - | \begin{equation} | ||
| - | DFT: S(k) = \sum^{N-1}_{n = 0}{s(n)e^{\frac{-j 2 \pi n k}{K}}}, k \in \{0, ..., K-1\} \\ | ||
| - | |||
| - | IDFT\ (inversa\ DFT): s(n) = \frac{1}{K}\sum^{K-1}_{k = 0}{S(k)e^{\frac{j 2 \pi n k}{K}}}, n \in \{0, ..., N-1\} | ||
| - | \end{equation} | ||
| - | De asemenea, calculați cu ajutorul FFT (FFT este varianta rapidă a DFT) și plotați spectrul folosind acest cod: | ||
| - | /*</hidden>*/ | ||
| - | |||
| - | Pentru asta va trebui să urmăriți acești pași: | ||
| - | * Creați semnalul $s(t)$ pentru $t\in\{1,\ldots,T=128\}$, folosind $a=0.05$. [<color red>0.5p</color>] | ||
| - | * Calculați și plotați (cu //stem//) spectrul folosind FFT (Transformata Fourier Rapidă), pe care nu am făcut-o încă la curs, dar o vom face în următoarele cursuri. Pentru moment, puteți folosi acest cod: [<color red>0.5p</color>] | ||
| - | <code python> | ||
| - | import numpy as np | ||
| - | import scipy as sp | ||
| - | import matplotlib.pyplot as plt | ||
| - | |||
| - | fx = np.zeros(T) | ||
| - | findex = T/2 * np.linspace(0, 1, T // 2) | ||
| - | fx[T // 2:] = findex | ||
| - | fx[1:(T//2)] = -findex[:0:-1] | ||
| - | fx[0] = -T // 2 | ||
| - | fs = sp.fft.fft(s) | ||
| - | plt.figure() | ||
| - | plt.stem(fx, abs(sp.fft.fftshift(fs))) | ||
| - | plt.xlabel('Frequency component (k)') | ||
| - | plt.ylabel('Magnitude of component') | ||
| - | plt.title('Fourier coefficients before amplitude modulation') | ||
| - | </code> | ||
| - | * Modulați semnalul în amplitudine folosind frecvența purtătoare $f_c = \frac{20}{T}$, i.e. face 20 de perioade complete în $T=128$ eșantioane ale semnalului $s(t)$. O variantă simplă de modulare este să calculați: $x(t) = (1+s(t)) \cdot \cos(2\pi f_c t)$. [<color red>0.5p</color>] | ||
| - | * Calculați și plotați (cum am făcut mai devreme, cu funcția //fft//) spectrul semnalului modulat în amplitudine. Comparați-l cu spectrul semnalului original. Este ceea ce v-ați așteptat? [<color red>0.5p</color>] | ||
| - | |||
| - | /*<hidden>*/ | ||
| - | === Exercițiul 4 – modularea în amplitudine [<color red>Bonus</color>] === | ||
| - | |||
| - | În acest exercițiu veți încerca să refaceți un semnal modulat în amplitudine ({{ps:labs:lab05_modulated_sound.mat|click aici}}). | ||
| - | Urmați următorii pași: | ||
| - | * Descărcați semnalul modulat și încărcați-l în MATLAB folosind funcția: <code matlab>load('lab05_modulated_sound.mat')</code> Se va încărca semnalul, precum și Fs, frecvența de eșantionare. | ||
| - | * Semnalul modulat provine dintr-un semnal original $s(t)$ care a fost înmulțit cu un semnal cosinus de frecvență $fc$, cos($2 \pi \cdot fc \cdot t$), unde t=0:1/Fs:timp_maxim. Precum ați observat și în exercițiul anterior spectrul semnalului modulat conține spectrul semnalului original shiftat la frecvența fc și -fc. Plotați spectrul pentru a vedea acest lucru. | ||
| - | * Dacă înmulțim (din nou) semnalul primit cu același semnal cosinus ar trebui să putem recupera semnalul original după aplicarea unui filtru trece-jos. Folosiți fc = Fs / 8. | ||
| - | * Pentru a realiza efectul unui filtru trece-jos vom egala cu zero coeficienții corespunzători frecvențelor peste un prag. | ||
| - | * Având noii coeficienți putem reface semnalul în timp folosind inversa Fourier. Realizați acest lucru folosind funcția //ifft//. | ||
| - | * Puteți asculta semnalul rezultat folosind funcția MATLAB //sound(semnal)//. Pentru a opri sunetul folosiți comanda <code matlab> clear sound </code> | ||
| - | |||
| - | |||
| - | /*</hidden>*/ | ||
| - | |||
| - | <note> Atenție: Vom învăța la curs că spectrul obținut prin transformata Fourier Discretă este periodic. Funcțiile //fft// / //ifft// consideră primii jumătate plus unu coeficienți pentru frecvențele pozitive, apoi următoarea jumătate corespunzătoare coeficienților negativi. Pentru a obține un spectru centrat în zero (doar ca să îl vizualizăm precum ne-am obișnuit) vom folosi funcția fftshift. </note> | ||
| </hidden> | </hidden> | ||
