Segnali e sistemi in Python - campioni, convoluzione e filtri

Il programma del corso comprende le "attività pratiche e dimostrazioni al calcolatore" (Matlab o Python): rappresentare segnali, convolvere, calcolare trasformate con la FFT. Gli esami contengono una domanda di calcolo numerico (FFT e zero-padding - TF, TFtd e asse delle pulsazioniLa FFT dei campioni di un segnale a durata finita, moltiplicata per il passo $T_c$, approssima la trasformata di Fourier: $X(\omega_k)\approx T_c,\mathtt{fft}(x,M)[k]$ con $\omega_k=\frac{2\pi k}{MT_c}$ (e fattore di fase $e^{-j\omega t_0}$ se l'asse parte da $t_0$). Lo zero-padding ($M>N$) infittisce i punti della stessa TFtd senza aggiungere informazione; la risoluzione dipende dalla durata osservata. Per un segnale reale $|X|$ è simmetrico: il picco in $k$ ha un gemello in $M-k$.FFT e zero-padding - TF, TFtd e asse delle pulsazioni →). Questa nota mostra gli strumenti di base con numpy; il codice è stato eseguito e i risultati riportati sono quelli ottenuti. (I grafici si ottengono con matplotlib, per esempio plt.plot(t, y); non sono riprodotti qui.)

Rappresentare un segnale

Un segnale continuo $x(t)$ si campiona con un passo piccolo $\Delta$ (qui $10^{-3}$ s): si usa un vettore di tempi t e il vettore dei valori x = f(t). L'asse dei tempi serve sempre: senza di esso non si sa a quale istante appartiene ogni campione.

import numpy as np
from scipy.signal import lfilter

rect = lambda t: (np.abs(t) < 0.5).astype(float)   # finestra rettangolare
u    = lambda t: (t >= 0).astype(float)            # gradino

dt = 1e-3
t = np.arange(-5, 5, dt)                           # asse dei tempi

np.sinc(x) è già $\frac{\sin\pi x}{\pi x}$, cioè la definizione di $\operatorname{sinc}$ usata nel corso (in Matlab: sinc, con la stessa convenzione; non è $\frac{\sin x}x$). Gli operatori *, +, ** agiscono elemento per elemento.

Area ed energia. Gli integrali diventano somme moltiplicate per il passo (approssimazione rettangolare):

x = np.exp(-3*t) * u(t)
print(np.sum(np.abs(x)**2) * dt)     # energia -> 0.16717  (esatta 1/6 = 0.16667)

L'errore è dell'ordine del passo: qui il campione in $t=0$ vale $1$ (il gradino numerico usa $u(0)=1$), mentre la convenzione dell'emivalore darebbe $\frac12$, e la differenza è circa $\frac{\Delta}{2}=5\cdot10^{-4}$, proprio lo scarto osservato (Segnali - supporto, area, valor medio, energia e potenzaUn segnale è una funzione del tempo (continuo $t$ o discreto $n$). Si descrive con pochi numeri: estensione, area, valor medio, energia $\int|x|^2$ e potenza (energia media). Energia finita implica potenza nulla; potenza finita non nulla implica energia infinita; per i segnali periodici tutto si calcola su un periodo.Segnali - supporto, area, valor medio, energia e potenza →).

Convoluzione di segnali continui

np.convolve(x, y) calcola la convoluzione discreta $\sum_kx(k)y(n-k)$. Per approssimare l'integrale di convoluzione (Sistemi LTI, risposta impulsiva e convoluzioneUn sistema lineare tempo-invariante (LTI) è completamente descritto dalla sua risposta impulsiva $h=\Sigma[\delta]$: l'uscita è la convoluzione $y=x*h$, cioè $y(t)=\int x(u)h(t-u)du$ (somma $\sum_k x(k)h(n-k)$ nel discreto). Il teorema discende da linearità e tempo-invarianza applicate alla scomposizione del segnale in impulsi.Sistemi LTI, risposta impulsiva e convoluzione →) servono due accorgimenti:

1. fattore $\Delta$: $\int x(u)y(t-u)du\approx\Delta\sum_kx(k\Delta)y(n\Delta-k\Delta)$;
2. asse dei tempi: l'uscita comincia all'istante $t_x+t_y$ (somma degli istanti iniziali, per la regola del supporto in Calcolo della convoluzione e sue proprietàIl supporto della convoluzione è la somma dei supporti, $\operatorname{rect}*\operatorname{rect}=\Lambda$, e due esponenziali causali danno $(e^{-bt}-e^{-at})/(a-b)$. Si calcola con il metodo grafico a casi (ribaltare, traslare, individuare gli intervalli di sovrapposizione). Proprietà: lineare, commutativa, associativa, $\delta$ è l'elemento neutro, la traslazione si somma, l'area è il prodotto delle aree.Calcolo della convoluzione e sue proprietà →) e ha $L_x+L_y-1$ campioni.

x = rect(t)
z  = np.convolve(x, x) * dt
tz = t[0] + t[0] + np.arange(len(z)) * dt
print(z.max(), tz[z.argmax()], z.sum()*dt)   # 1.0  -0.001  1.0

Il risultato è il triangolo $\Lambda$: picco $1$ in $t\approx0$ e area $1\cdot1=1$ (regola dell'area).

Uscita di un filtro RC. $h(t)=e^{-t}u(t)$ con ingresso $\operatorname{rect}\left(\frac{t-1}2\right)$ (esempio 5 di Calcolo della convoluzione e sue proprietàIl supporto della convoluzione è la somma dei supporti, $\operatorname{rect}*\operatorname{rect}=\Lambda$, e due esponenziali causali danno $(e^{-bt}-e^{-at})/(a-b)$. Si calcola con il metodo grafico a casi (ribaltare, traslare, individuare gli intervalli di sovrapposizione). Proprietà: lineare, commutativa, associativa, $\delta$ è l'elemento neutro, la traslazione si somma, l'area è il prodotto delle aree.Calcolo della convoluzione e sue proprietà →):

h = np.exp(-t) * u(t)
x = rect((t - 1) / 2)
y  = np.convolve(x, h) * dt
ty = t[0] + t[0] + np.arange(len(y)) * dt

Confrontando in alcuni istanti con la formula esatta $1-e^{-t}$ per $0<t<2$ e $e^{-t}(e^2-1)$ per $t>2$:

(errori dell'ordine di $\Delta$).

Convoluzione discreta esatta. Per le sequenze non serve nessun fattore:

print(np.convolve([1, 2, 1], [1, -1, 2]))     # [1 1 1 3 2]

che è il risultato calcolato a mano in Sistemi LTI, risposta impulsiva e convoluzioneUn sistema lineare tempo-invariante (LTI) è completamente descritto dalla sua risposta impulsiva $h=\Sigma[\delta]$: l'uscita è la convoluzione $y=x*h$, cioè $y(t)=\int x(u)h(t-u)du$ (somma $\sum_k x(k)h(n-k)$ nel discreto). Il teorema discende da linearità e tempo-invarianza applicate alla scomposizione del segnale in impulsi.Sistemi LTI, risposta impulsiva e convoluzione →.

Filtri con equazioni alle differenze

Un sistema $y(n)-\frac12y(n-1)=x(n)$ (Trasformata Z ed equazioni alle differenzeLa trasformata Z, $X(z)=\sum x(n)z^{-n}$, è la versione discreta di Laplace (la TFtd è $X(e^{j\omega})$). Un sistema descritto da un'equazione alle differenze $\sum a_ky(n-k)=\sum b_kx(n-k)$ ha $H(z)=\frac{\sum b_kz^{-k}}{\sum a_kz^{-k}}$ ed è BIBO stabile (se causale) quando tutti i poli stanno dentro il cerchio unitario, $|p|<1$. Cenni.Trasformata Z ed equazioni alle differenze →) si simula con lfilter(b, a, x), dove b e a sono i coefficienti $b_k$ e $a_k$ di $\sum a_ky(n-k)=\sum b_kx(n-k)$:

d = np.zeros(8); d[0] = 1                      # impulso
print(lfilter([1], [1, -0.5], d))              # 1, 0.5, 0.25, ...

L'uscita è $\left(\frac12\right)^n$ per $n=0,\ldots,7$, cioè la risposta impulsiva attesa $h(n)=\left(\frac12\right)^nu(n)$. In Matlab la funzione omonima è filter(b, a, x). Nel caso FIR (a = [1]) lfilter coincide con la convoluzione con b.

Aliasing

Si campiona $\cos(2\pi\cdot7\,t)$ a $F_c=10$ Hz e $\cos(2\pi\cdot3\,t)$ a $10$ Hz:

n = np.arange(10)
print(np.allclose(np.cos(2*np.pi*7*n/10), np.cos(2*np.pi*3*n/10)))   # True

I due toni producono gli stessi campioni: dopo il campionamento a $10$ Hz il tono a $7$ Hz non si distingue da quello a $3$ Hz (Teorema del campionamento, interpolazione e aliasingTeorema di Shannon: un segnale a banda limitata $\omega_M$ si ricostruisce esattamente dai campioni se $T_c<\pi/\omega_M$ (frequenza di campionamento maggiore di quella di Nyquist $2f_{\max}$), con la formula di interpolazione ideale $x(t)=\sum_nx(nT_c)\operatorname{sinc}\left(\frac{t-nT_c}{T_c}\right)$. Sotto Nyquist c'è aliasing: le frequenze alte si confondono con quelle basse e l'informazione è persa.Teorema del campionamento, interpolazione e aliasing →).

Somme parziali della serie di Fourier e fenomeno di Gibbs

def S(t, K):
    s = 0.5 * np.ones_like(t)
    for k in range(1, K + 1, 2):
        s += 2/(k*np.pi) * np.sin(k*np.pi/2) * np.cos(k*np.pi*t)
    return s

tt = np.linspace(0.3, 0.7, 200001)
for K in (9, 49, 199):
    print(K, S(tt, K).max())       # 1.0912  1.0896  1.0895

La somma parziale dell'onda quadra (Serie di Fourier - analisi e sintesiUn segnale periodico di periodo $T$ si scrive come somma di esponenziali in relazione armonica, $x(t)=\sum_ka_ke^{jk\omega_0t}$ con $\omega_0=2\pi/T$, e i coefficienti si ottengono per proiezione, $a_k=\frac1T\int_Tx(t)e^{-jk\omega_0t}dt$. L'ortogonalità degli esponenziali dà la formula; la convergenza è in media quadratica (Riesz-Fischer), con il fenomeno di Gibbs nei salti.Serie di Fourier - analisi e sintesi →) ha sempre un picco di circa $1{,}09$ vicino al salto: il fenomeno di Gibbs, che non scompare aumentando $K$.

Verificare numericamente le proprietà dei sistemi

Per un sistema (una funzione che trasforma un vettore) si possono eseguire i test di linearità e tempo-invarianza di Sistemi e loro proprietàUn sistema trasforma un segnale di ingresso in uno di uscita, $y=\Sigma[x]$. Le proprietà da saper verificare sono: memoria (statico/dinamico), causalità, linearità (additività + omogeneità), tempo-invarianza, stabilità BIBO, invertibilità, realtà. Per negare una proprietà basta un controesempio; per affermarla serve una dimostrazione con segnali generici.Sistemi e loro proprietà → su segnali casuali. Un risultato True non è una dimostrazione, ma un risultato False è un controesempio (e quindi una dimostrazione di non linearità).

def lineare(sis, x1, x2, a=2.0, b=-3.0):
    return np.allclose(sis(a*x1 + b*x2), a*sis(x1) + b*sis(x2))
def tempo_inv(sis, x, n0=3):
    xr = np.roll(x, n0)                          # ritardo ciclico: i segnali finiscono a zero
    return np.allclose(sis(xr), np.roll(sis(x), n0))

Con segnali che hanno 20 campioni casuali e poi zeri: $y=x^2$ $\to$ (lineare: False, tempo-invariante: True); media mobile $\frac{x(n)+x(n-1)}2$ $\to$ (True, True); $y(n)=n\,x(n)$ $\to$ (True, False). Sono le tre classificazioni della tabella di Sistemi e loro proprietàUn sistema trasforma un segnale di ingresso in uno di uscita, $y=\Sigma[x]$. Le proprietà da saper verificare sono: memoria (statico/dinamico), causalità, linearità (additività + omogeneità), tempo-invarianza, stabilità BIBO, invertibilità, realtà. Per negare una proprietà basta un controesempio; per affermarla serve una dimostrazione con segnali generici.Sistemi e loro proprietà →.

Tabella Python $\leftrightarrow$ Matlab

Le differenze da ricordare: in Python gli indici partono da $0$ e a:b esclude b; in Matlab si parte da $1$ e a:b include b; il prodotto elemento per elemento è * in numpy e .* in Matlab.

Errori comuni

• Dimenticare il fattore $\Delta$ nella convoluzione e nell'integrale approssimato (i risultati sono sbagliati di un fattore $1/\Delta$).
• Dimenticare l'asse dei tempi della convoluzione (la somma degli istanti iniziali).
• Usare np.sinc pensando che sia $\frac{\sin x}x$ (è $\frac{\sin\pi x}{\pi x}$).
• Prendere un passo $\Delta$ troppo grande rispetto alla scala di variazione del segnale.