[2p]
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, $x_1, x_2, \ldots, x_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$
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) și amplitudine 1, la aceeaşi frecvenţă de eşantionare (folosiţi acelaşi număr de eşantioane, e.g. $N=128$).
Implementaţi sistemul de procesare menţionat mai sus.[1p]
Creaţi secvenţele corespunzătoare de ieşire.
Plotaţi toate secvenţele de intrare/ieşire. [0.5p]
Ce fel de sistem de procesare este acesta ?
Cum ați putea implementa sistemul de procesare menționat mai sus, folosind funcția
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. Plotați din nou toate secvenţele de intrare/ieşire.[0.5p]
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 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.
[3p]
Î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$.
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:
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')
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. [1p]
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)$. [1p]
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? [1p]
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.
[3p]
În acest exercițiu vom experimenta subeșantionarea, metodă care reduce frecvența de eșantionare necesară pentru a obține întreg spectru al semnalului.
Un semnal real (un cosinus sau o sinusoidă) cu o frecvență purtătoare (centrată) la 3kHz este eşantionată la 8 kHz. Lăţimea de bandă a semnalului este 1 kHz, cum puteţi vedea mai jos:
Task-uri/întrebări:
Generați semnalul $signal = cos(2 \cdot \pi \cdot f_c \cdot t)$, unde $f_c = 3000$ Hz, $t=0:T_s:100 \cdot T_s$ și $T_s = \frac{1}{f_s} = \frac{1}{8000}$
De ce avem energie în ambele frecvenţe $f_c$ şi $-f_c$ ?
Generaţi cod Python pentru a calcula transformata Fourier Discretă a unei sinusoide de $3kHz$ eşantionată la $8kHz$. Folosiți $N=64$ pentru fft. Folosiți acest cod pentru DFT: [
1p]
figure;
fx = zeros(1, N);
fidx = (fs/N) * linspace(0,N-1,N);
spectrum = fft(signal, N);
stem(fidx, abs(spectrum));
xlabel('Frequency (Hz)');
ylabel('Amplitude');
title('Spectrum of signal');
Plotaţi (cu 'stem') rezultatele şi observaţi frecvenţele. Acum folosiţi “fftshift” pe rezultatul obţinut de la “fft” şi plotaţi din nou rezultatele. Folosiți acest cod: [
1p]
figure;
fidx = (fs/N)*linspace(-N/2, N/2-1, N);
stem(fidx, abs(fftshift(spectrum)));
xlabel('Frequency (Hz)');
ylabel('Amplitude');
title('Zero-centred frequency spectrum of signal');
De ce într-un caz avem semnal la $3 kHz$ şi $5 kHz$ în timp ce, în celălalt caz pare să avem la $-3 kHz$ şi $3 kHz$ ?
Care este semnificaţia semnalului de $5 kHz$ ? Când plotaţi rezultatele de la “fft” luaţi în considerare că, în mod implicit output-ul transformatei Discrete Fourier, deci şi “fft”, arată frecvenţele de la $0$ la $f_s$ (frecvenţa de eşantionare).
Putem să reducem frecvenţa de eşantionare în timp ce încă detectăm spectrul dorit de frecvenţe ? (vedeţi bandpass sampling
Fig 2-7). Cât de mult putem reduce frecvenţa de eşantionare? Ce se întâmplă dacă reducem frecvenţa de eşantionare sub 2B? [
1p]