Salta al contenuto
Note per Studenti Esercizio - laboratorio 2, convoluzione numerica

Esercizio - laboratorio 2, convoluzione numerica

In questa pagina 7

Testo (laboratorio 2 del corso Teoria dei Segnali, UniPD, file lab2.m, parti 4-8). (4) Scrivere la funzione calcolatore_convoluzione che calcola la convoluzione di due segnali campionati con lo stesso passo (due cicli annidati e il fattore Δt\Delta t, con l'asse temporale giusto) e usarla su due impulsi rettangolari, rect⁡1\operatorname{rect}_1 non nullo tra −0,8-0{,}8 e 0,80{,}8 e rect⁡2\operatorname{rect}_2 tra −3-3 e 33. (5) Verificare la proprietà commutativa. (6) Verificare che l'area della convoluzione è il prodotto delle aree. (7) Rifare il calcolo con il comando della libreria. (8) Convoluzione ciclica di due rect periodici (periodo Tp=4T_p=4) e confronto con la convoluzione lineare dei segnali periodici troncati.

Teoria usata: ConvoluzioneLa convoluzione $x*y(t)=\int x(u),y(t-u),du$ combina due segnali ribaltando e traslando il secondo, moltiplicando e integrando. È commutativa, associativa, lineare; l'area del risultato è il prodotto delle aree; l'estensione è la somma delle estensioni (estremo con estremo); l'impulso $\delta$ è l'elemento neutro; la convoluzione con il gradino integra. Per due segnali periodici di uguale periodo si usa la convoluzione ciclica (integrale su un periodo). È l'operazione del filtraggio: l'uscita di un filtro è la convoluzione dell'ingresso con la risposta impulsiva.Convoluzione →, Segnali notevoli - gradino, rect, tri, sinc ed esponenzialiI segnali di uso più frequente sono la costante, la sinusoide $A_0\cos(2\pi f_0t+\varphi_0)$ e l'esponenziale complesso $Ae^{i2\pi f_0t}$ (periodici, a potenza finita), il gradino $\mathbf 1(t)$ e il segno, e gli impulsi a energia finita: $\operatorname{rect}$ (area $D$), $\operatorname{tri}$, $\operatorname{sinc}$ (area $1$), la gaussiana $e^{-\pi t^2}$ e gli esponenziali smorzati. Per ciascuno si sanno a memoria forma, area ed energia; gli altri segnali si ottengono da questi con traslazioni, scalature, somme e differenze.Segnali notevoli - gradino, rect, tri, sinc ed esponenziali →, Energia, potenza e valor medio dei segnaliSu un segnale continuo si calcolano quattro numeri riassuntivi: l'area $\int s,dt$, il valor medio (componente continua) $\lim\frac1{2T}\int_{-T}^Ts,dt$, l'energia $\int|s|^2dt$ e la potenza media $\lim\frac1{2T}\int_{-T}^T|s|^2dt$. Un segnale ad energia finita ha potenza nulla e uno a potenza finita non nulla ha energia infinita. Per un segnale periodico di periodo $T_p$ area ed energia si calcolano su un periodo; un segnale periodico è la ripetizione periodica $\sum_ku(t-kT_p)$ di un suo periodo. Valgono per traslazione l'invarianza, per scala $s(at)$ la divisione per $|a|$.Energia, potenza e valor medio dei segnali → (area), Trasformata di Fourier discreta (DFT) e FFTUn segnale discreto periodico di periodo $NT$ ha una trasformata discreta e periodica, la DFT: $S(kF)=\sum_{n=0}^{N-1}T,s(nT)e^{-i2\pi kn/N}$, con $F=1/(NT)$, e $s(nT)=\sum_{k=0}^{N-1}F,S(kF)e^{i2\pi kn/N}$. Dipende da soli $N$ numeri, si calcola senza approssimazioni e con la FFT costa $N\log_2N$ invece di $N^2$. I coefficienti di Fourier del segnale periodico sono $S_k=F,S(kF)$. La DFT dà anche campioni della trasformata di un segnale discreto o continuo a durata limitata (con zero-padding), ma la scalatura $T$ e l'asse delle frequenze vanno gestiti con cura.Trasformata di Fourier discreta (DFT) e FFT → (convoluzione ciclica).

Cosa si vuole calcolare e perché

La convoluzione di due segnali è x∗y(t)=∫x(τ) y(t−τ) dτx*y(t)=\int x(\tau)\,y(t-\tau)\,d\tau (ConvoluzioneLa convoluzione $x*y(t)=\int x(u),y(t-u),du$ combina due segnali ribaltando e traslando il secondo, moltiplicando e integrando. È commutativa, associativa, lineare; l'area del risultato è il prodotto delle aree; l'estensione è la somma delle estensioni (estremo con estremo); l'impulso $\delta$ è l'elemento neutro; la convoluzione con il gradino integra. Per due segnali periodici di uguale periodo si usa la convoluzione ciclica (integrale su un periodo). È l'operazione del filtraggio: l'uscita di un filtro è la convoluzione dell'ingresso con la risposta impulsiva.Convoluzione →): per ogni istante tt si ribalta e si trasla yy, si moltiplica per xx e si prende l'area del prodotto. Con i campioni xi=x(tx0+i Δt)x_i=x(t_{x0}+i\,\Delta t) e yk=y(ty0+k Δt)y_k=y(t_{y0}+k\,\Delta t) l'integrale diventa una somma e dτd\tau diventa Δt\Delta t:

zn≈∑kΔt  xn−k yk.z_n\approx\sum_k\Delta t\;x_{n-k}\,y_k.

Tre cose da ricordare, e da giustificare nel codice:

  1. Il fattore Δt\Delta t moltiplica ogni termine: senza di esso si ottiene la somma dei prodotti, non l'area (stesso discorso di Esercizio - laboratorio 2, area ed energia numeriche).
  2. Il numero di campioni: se xx ha LxL_x campioni e yy ne ha LyL_y, i valori di n=i+kn=i+k vanno da 00 a (Lx−1)+(Ly−1)(L_x-1)+(L_y-1), cioè Lz=Lx+Ly−1L_z=L_x+L_y-1 campioni.
  3. L'asse dei tempi. Il campione znz_n ha tempo tx0+i Δt+ty0+k Δt=(tx0+ty0)+n Δtt_{x0}+i\,\Delta t+t_{y0}+k\,\Delta t=(t_{x0}+t_{y0})+n\,\Delta t: l'asse della convoluzione parte dalla somma degli istanti iniziali e arriva alla somma degli istanti finali. È la regola dell'estensione: se xx vive in [ax,bx][a_x,b_x] e yy in [ay,by][a_y,b_y], il supporto di x∗yx*y è [ax+ay, bx+by][a_x+a_y,\ b_x+b_y].

Parte 4 - la funzione e la convoluzione di due rect

python
import time as orologio
import numpy as np
import matplotlib.pyplot as plt

# funzione dei trapezi (laboratorio 2, parte 1), serve per le aree
def calcolatore_area_trapezi(f, t):
    if len(t) != len(f):
        raise ValueError("Il vettore della funzione e quello del tempo non hanno la stessa lunghezza!")
    I = 0.0
    for a in range(len(f) - 1):
        I = I + (f[a + 1] + f[a]) * (t[a + 1] - t[a]) / 2
    return I

# Convoluzione numerica di due segnali campionati con lo stesso passo
def calcolatore_convoluzione(x, t_x, y, t_y):
    """x, y: campioni; t_x, t_y: i loro istanti. Restituisce (z, t_z) con z = x*y."""
    dt_x = np.diff(t_x)
    dt_y = np.diff(t_y)
    tol = 1e-12

    if np.any(np.abs(dt_x - dt_x[0]) > tol):
        raise ValueError('Il vettore t_x non ha passo costante.')
    if np.any(np.abs(dt_y - dt_y[0]) > tol):
        raise ValueError('Il vettore t_y non ha passo costante.')
    if abs(dt_x[0] - dt_y[0]) > tol:
        raise ValueError('I due vettori temporali non hanno lo stesso passo di campionamento.')

    dt = dt_x[0]                       # passo temporale comune
    Lx = len(x)
    Ly = len(y)

    # asse temporale della convoluzione: parte dalla somma degli istanti iniziali
    Lz = Lx + Ly - 1                   # numero di campioni
    t_min = t_x[0] + t_y[0]            # inizio della convoluzione
    t_z = t_min + np.arange(Lz) * dt   # MATLAB: [t_min:dt:t_max], t_max = t_x(end) + t_y(end)

    z = np.zeros(Lz)
    for n in range(Lz):
        for k in range(Ly):
            if 0 <= n - k < Lx:                         # MATLAB: n-k+1 tra 1 e Lx
                z[n] = z[n] + dt * x[n - k] * y[k]      # z(t_n) = somma dt x(t_n - t_k) y(t_k)
    return z, t_z

L'indice n - k deve restare tra 00 e Lx−1L_x-1: la condizione 0 <= n - k < Lx salta i prodotti che cadono fuori dal vettore di xx (fuori dal suo supporto xx vale zero). In MATLAB gli indici partono da 1 e la condizione diventa n-k+1 >= 1 e <= Lx. Il ciclo interno ha Ly=1001L_y=1001 passi, quello esterno Lz=2001L_z=2001: due milioni di prodotti.

python
# PARTE 4 - convoluzione di due rect (segnali aperiodici a tempo continuo)
t_start = -5
t_end = 5
dt = 0.01
time = t_start + np.arange(round((t_end - t_start) / dt) + 1) * dt      # 1001 campioni

# primo rect: 1 tra -0.8 e 0.8 (estremi compresi)
t11, t12 = -0.8, 0.8
rect1 = np.zeros(len(time))
i_t11 = round((t11 - t_start) / dt)          # indice da 0 (in MATLAB: (t11-t_start)/dt + 1)
i_t12 = round((t12 - t_start) / dt)
rect1[i_t11:i_t12 + 1] = 1                   # +1 perche' l'estremo della fetta e' escluso

# secondo rect: 1 tra -3 e 3
t21, t22 = -3, 3
rect2 = np.zeros(len(time))
i_t21 = round((t21 - t_start) / dt)
i_t22 = round((t22 - t_start) / dt)
rect2[i_t21:i_t22 + 1] = 1

t0 = orologio.perf_counter()
y, t_y = calcolatore_convoluzione(rect1, time, rect2, time)
tempo_cicli = orologio.perf_counter() - t0

print(f"campioni: rect1 = {int(rect1.sum())} uni, rect2 = {int(rect2.sum())} uni, convoluzione = {len(y)} (1001 + 1001 - 1)")
print(f"asse della convoluzione: da {t_y[0]:.2f} a {t_y[-1]:.2f} s")
print(f"massimo della convoluzione: {y.max():.4f}  (tempo {t_y[y.argmax()]:.2f} s)")
print(f"primo e ultimo istante con y > 0: {t_y[y > 0][0]:.2f} e {t_y[y > 0][-1]:.2f} s")
print(f"tempo dei due cicli: circa {tempo_cicli:.1f} s")

plt.figure(figsize=(8, 7))
plt.subplot(3, 1, 1)
plt.plot(time, rect1, linewidth=2, color='r')
plt.xlim(t_y[0], t_y[-1]); plt.ylim(0, 1.3); plt.grid(True); plt.title('Primo segnale')
plt.subplot(3, 1, 2)
plt.plot(time, rect2, linewidth=2, color='b')
plt.xlim(t_y[0], t_y[-1]); plt.ylim(0, 1.3); plt.grid(True); plt.title('Secondo segnale')
plt.subplot(3, 1, 3)
plt.plot(t_y, y, linewidth=2, color='k')
plt.xlim(t_y[0], t_y[-1]); plt.ylim(0, 1.7); plt.grid(True); plt.title('Convoluzione')
for x_ in (t11, t12):
    plt.axvline(x_, linestyle='--', color='r')          # MATLAB: xline
for x_ in (t21, t22):
    plt.axvline(x_, linestyle='--', color='b')
for x_ in (t11 + t21, t12 + t22):                       # estremi della convoluzione = somme degli estremi
    plt.axvline(x_, linestyle='--', color='k')
plt.xlabel('t [s]')
plt.tight_layout()
plt.show()

Risultato:

campioni: rect1 = 161 uni, rect2 = 601 uni, convoluzione = 2001 (1001 + 1001 - 1)
asse della convoluzione: da -10.00 a 10.00 s
massimo della convoluzione: 1.6100  (tempo -2.20 s)
primo e ultimo istante con y > 0: -3.80 e 3.80 s
tempo dei due cicli: circa 0.3 s

Costruire i rect. Il rect è un vettore di zeri con degli uni nelle posizioni giuste. L'istante tt corrisponde all'indice t−tstartΔt\frac{t-t_{start}}{\Delta t} (da 0; in MATLAB +1+1), che per t=−0,8t=-0{,}8 è 4,20,01=420\frac{4{,}2}{0{,}01}=420. Con round si evitano gli errori di arrotondamento dei decimali; nella fetta si somma 1 all'ultimo indice perché l'estremo finale è escluso. Il primo rect ha quindi 161161 uni (da −0,8-0{,}8 a 0,80{,}8, estremi compresi: 160160 intervalli e 161161 campioni), il secondo 601601.

Cosa viene fuori. L'asse della convoluzione va da −5−5=−10-5-5=-10 a 5+5=105+5=10 s con 1001+1001−1=20011001+1001-1=2001 campioni, come stampato. La convoluzione di due rect è un trapezio: nonnulla tra −0,8−3=−3,8-0{,}8-3=-3{,}8 e 0,8+3=3,80{,}8+3=3{,}8 (la somma degli estremi, e le linee tratteggiate nere nel grafico del laboratorio sono proprio lì), sale linearmente fino a t=0,8−3=−2,2t=0{,}8-3=-2{,}2, resta piatta fino a t=−0,8+3=2,2t=-0{,}8+3=2{,}2 e scende fino a 3,83{,}8 (anche i quattro vertici sono somme di estremi). Si vede ragionando sul significato: y(t)y(t) è la lunghezza della parte comune tra [−0,8, 0,8][-0{,}8,\,0{,}8] e l'intervallo [t−3, t+3][t-3,\,t+3] (il secondo rect ribaltato e traslato in tt):

y(t)=max⁡(0, min⁡(t+3, 0,8)−max⁡(t−3, −0,8))={t+3,8−3,8<t<−2,21,6−2,2≤t≤2,23,8−t2,2<t<3,8y(t)=\max\Bigl(0,\ \min(t+3,\,0{,}8)-\max(t-3,\,-0{,}8)\Bigr)=\begin{cases}t+3{,}8&-3{,}8<t<-2{,}2\\1{,}6&-2{,}2\le t\le2{,}2\\3{,}8-t&2{,}2<t<3{,}8\end{cases}

Il plateau vale 1,61{,}6, cioè l'area del rect più piccolo, perché per ∣t∣≤2,2|t|\le2{,}2 il rect corto cade interamente dentro il rect lungo.

Grafico interattivo: I due rect del laboratorio (rect1 di larghezza 1,6 e altezza 1, rect2 di larghezza 6 e altezza 1) e la loro convoluzione, un trapezio con vertici in ±3,8 (valore 0) e ±2,2 (valore 1,6: l'area del rect corto)

Perché 1,61 e non 1,60. Il massimo numerico è 1,611{,}61 e non 1,601{,}60 (l'errore massimo rispetto al trapezio teorico è 0,01=Δt0{,}01=\Delta t). Il motivo sta nei campioni: il vettore rect1 ha 161161 uni, quindi la sua "larghezza" numerica è 161 Δt=1,61161\,\Delta t=1{,}61 s (ogni campione pesa Δt\Delta t), non 1,601{,}60. Nel rect del corso i due estremi ±0,8\pm0{,}8 valgono 12\frac12 (regola dell'emivalore); il vettore del laboratorio li mette a 11. Se si mettesse 0,50{,}5 ai due estremi il vettore pesava 160 Δt=1,60160\,\Delta t=1{,}60. Le altre differenze col trapezio teorico sono dello stesso ordine di Δt\Delta t: y(3)=0,81y(3)=0{,}81 contro 0,800{,}80, ultimo campione non nullo y(3,8)=0,01y(3{,}8)=0{,}01.

python
# Confronto con il risultato "teorico" (trapezio): sovrapposizione dei due rect a distanza t
y_teorico = np.maximum(0, np.minimum(t_y + 3, 0.8) - np.maximum(t_y - 3, -0.8))
print(f"massimo teorico: {y_teorico.max():.2f},  errore massimo |numerico - teorico|: {np.max(np.abs(y - y_teorico)):.4f}")
print("y(0)       =", round(y[np.argmin(np.abs(t_y))], 4))
print("y(2.2)     =", round(y[np.argmin(np.abs(t_y - 2.2))], 4), "(inizio della discesa)")
print("y(3.0)     =", round(y[np.argmin(np.abs(t_y - 3.0))], 4), "(meta' discesa, teorico 0.8)")
print("y(3.8)     =", round(y[np.argmin(np.abs(t_y - 3.8))], 4), "(ultimo campione non nullo)")
massimo teorico: 1.60,  errore massimo |numerico - teorico|: 0.0100
y(0)       = 1.61
y(2.2)     = 1.61 (inizio della discesa)
y(3.0)     = 0.81 (meta' discesa, teorico 0.8)
y(3.8)     = 0.01 (ultimo campione non nullo)

Parte 5 - proprietà commutativa

python
# PARTE 5 - proprieta' commutativa: x*h = h*x
y2, t_y2 = calcolatore_convoluzione(rect2, time, rect1, time)
print("assi uguali:", np.allclose(t_y, t_y2))
print(f"massima differenza tra rect1*rect2 e rect2*rect1: {np.max(np.abs(y - y2)):.2e}")

# (i grafici sono identici a quelli della parte 4, con i primi due riquadri scambiati)
assi uguali: True
massima differenza tra rect1*rect2 e rect2*rect1: 0.00e+00

x∗h=h∗xx*h=h*x (ConvoluzioneLa convoluzione $x*y(t)=\int x(u),y(t-u),du$ combina due segnali ribaltando e traslando il secondo, moltiplicando e integrando. È commutativa, associativa, lineare; l'area del risultato è il prodotto delle aree; l'estensione è la somma delle estensioni (estremo con estremo); l'impulso $\delta$ è l'elemento neutro; la convoluzione con il gradino integra. Per due segnali periodici di uguale periodo si usa la convoluzione ciclica (integrale su un periodo). È l'operazione del filtraggio: l'uscita di un filtro è la convoluzione dell'ingresso con la risposta impulsiva.Convoluzione →): scambiare i due segnali dà esattamente lo stesso vettore (differenza zero) e lo stesso asse. Si capisce guardando la somma: zn=Δt∑i+k=nxi ykz_n=\Delta t\sum_{i+k=n}x_i\,y_k è simmetrica nello scambio dei nomi i↔ki\leftrightarrow k, perché la condizione è i+k=ni+k=n.

Parte 6 - l'area della convoluzione è il prodotto delle aree

python
# PARTE 6 - l'area della convoluzione e' il prodotto delle aree
area_rect1 = calcolatore_area_trapezi(rect1, time)
area_rect2 = calcolatore_area_trapezi(rect2, time)

area_trapezio = calcolatore_area_trapezi(y, t_y)
aree_moltiplicate = area_rect1 * area_rect2

print(f"area di rect1 (trapezi) = {area_rect1:.4f},  area di rect2 = {area_rect2:.4f}")
print(f"area della convoluzione (trapezi) = {area_trapezio:.4f}")
print(f"prodotto delle aree               = {aree_moltiplicate:.4f}")
print(f"differenza: {abs(area_trapezio - aree_moltiplicate):.1e}")

# i valori "teorici" dei due rect continui sarebbero 1.6 e 6 (prodotto 9.6): vedi il testo
# stessa verifica con la regola "area = dt * somma dei campioni"
A1 = dt * rect1.sum()
A2 = dt * rect2.sum()
Ay = dt * y.sum()
print(f"dt*somma: {A1:.4f}, {A2:.4f}, prodotto {A1 * A2:.4f}, area di y {Ay:.4f}")
area di rect1 (trapezi) = 1.6100,  area di rect2 = 6.0100
area della convoluzione (trapezi) = 9.6761
prodotto delle aree               = 9.6761
differenza: 5.0e-13
dt*somma: 1.6100, 6.0100, prodotto 9.6761, area di y 9.6761

La regola ∫(x∗y) dt=(∫x dt)(∫y dt)\int(x*y)\,dt=\bigl(\int x\,dt\bigr)\bigl(\int y\,dt\bigr) (ConvoluzioneLa convoluzione $x*y(t)=\int x(u),y(t-u),du$ combina due segnali ribaltando e traslando il secondo, moltiplicando e integrando. È commutativa, associativa, lineare; l'area del risultato è il prodotto delle aree; l'estensione è la somma delle estensioni (estremo con estremo); l'impulso $\delta$ è l'elemento neutro; la convoluzione con il gradino integra. Per due segnali periodici di uguale periodo si usa la convoluzione ciclica (integrale su un periodo). È l'operazione del filtraggio: l'uscita di un filtro è la convoluzione dell'ingresso con la risposta impulsiva.Convoluzione →) è verificata con l'accuratezza del calcolatore (5⋅10−135\cdot10^{-13}): l'area della convoluzione è 9,67619{,}6761 e il prodotto delle aree dei due rect è 1,61×6,01=9,67611{,}61\times6{,}01=9{,}6761. Le aree dei vettori sono 1,611{,}61 e 6,016{,}01, e non i 1,61{,}6 e 66 dei rect ideali (con prodotto 9,69{,}6), per il motivo visto nella parte 4: i campioni a 11 sono 161161 e 601601, ciascuno pesa Δt\Delta t, e i trapezi applicati a un vettore che vale zero agli estremi dell'asse coincidono con la somma Δt∑nxn\Delta t\sum_nx_n. Su campioni la regola vale esattamente, perché

Δt∑nzn=Δt⋅Δt∑i∑kxiyk=(Δt∑ixi)(Δt∑kyk),\Delta t\sum_nz_n=\Delta t\cdot\Delta t\sum_i\sum_kx_iy_k=\Bigl(\Delta t\sum_ix_i\Bigr)\Bigl(\Delta t\sum_ky_k\Bigr),

(ogni coppia (i,k)(i,k) compare in un solo n=i+kn=i+k). Il fattore Δt\Delta t nella definizione della convoluzione è proprio ciò che dà questa proprietà; senza di esso l'area di zz sarebbe 1Δt\frac1{\Delta t} volte il prodotto.

Parte 7 - il comando della libreria

python
# PARTE 7 - il comando di numpy: np.convolve (somma dei prodotti, senza il passo)
t0 = orologio.perf_counter()
y_conv = dt * np.convolve(rect2, rect1)
tempo_np = orologio.perf_counter() - t0
print("stessa lunghezza dei cicli:", len(y_conv) == len(y))
print(f"massima differenza dalla versione a cicli: {np.max(np.abs(y_conv - y)):.2e}")
print(f"tempo di np.convolve: {tempo_np * 1000:.1f} ms (contro circa {tempo_cicli:.1f} s dei cicli)")
stessa lunghezza dei cicli: True
massima differenza dalla versione a cicli: 3.44e-14
tempo di np.convolve: 0.1 ms (contro circa 0.3 s dei cicli)

np.convolve(a, b) (in MATLAB conv) calcola ∑kan−kbk\sum_ka_{n-k}b_k per n=0,…,La+Lb−2n=0,\dots,L_a+L_b-2: la stessa doppia somma dei cicli, ma senza Δt\Delta t e senza asse dei tempi. Perciò si moltiplica per dt e si usa l'asse tx0+ty0+nΔtt_{x0}+t_{y0}+n\Delta t costruito sopra. Il risultato coincide con la versione a cicli fino a 10−1410^{-14} (arrotondamenti) e il tempo scende da circa 0,30{,}3 s a 0,10{,}1 ms, perché il ciclo gira in C invece che in Python (i tempi cambiano da una macchina all'altra). Per vettori molto lunghi si può fare ancora meglio con la FFT (scipy.signal.fftconvolve), vedi i laboratori 3 e 4.

Parte 8 - convoluzione ciclica contro convoluzione lineare

Cos'è. Due segnali periodici di periodo TpT_p hanno energia infinita e la loro convoluzione lineare ∫−∞∞\int_{-\infty}^{\infty} non esiste. Per i segnali periodici si definisce la convoluzione ciclica (o periodica) integrando su un solo periodo: (x⊛y)(t)=∫Tpx(τ) y(t−τ) dτ(x\circledast y)(t)=\int_{T_p}x(\tau)\,y(t-\tau)\,d\tau, che è a sua volta periodica di periodo TpT_p. Sui campioni di un periodo (LL campioni, LΔt=TpL\Delta t=T_p)

zn=Δt∑k=0L−1x(n−k) mod L  yk.z_n=\Delta t\sum_{k=0}^{L-1}x_{(n-k)\bmod L}\;y_k.

L'indice (n−k) mod L(n-k)\bmod L "si riavvolge": se n−k<0n-k<0 si riparte dalla fine del vettore, perché il segnale si ripete. In Python l'operatore % dà già un risultato tra 00 e L−1L-1, mentre MATLAB, con indici da 1, scrive mod(n-k, Lx) + 1 (e il controllo > 0 && <= Lx del file è sempre vero).

python
# PARTE 8 - convoluzione ciclica
def calcolatore_convoluzione_ciclica(x, t_x, y, t_y):
    dt_x = np.diff(t_x)
    dt_y = np.diff(t_y)
    tol = 1e-12
    if np.any(np.abs(dt_x - dt_x[0]) > tol) or np.any(np.abs(dt_y - dt_y[0]) > tol):
        raise ValueError('Passo non costante.')
    if abs(dt_x[0] - dt_y[0]) > tol:
        raise ValueError('I due vettori temporali non hanno lo stesso passo di campionamento.')
    dt = dt_x[0]

    Lx = len(x)
    Ly = len(y)
    if Lx != Ly:
        raise ValueError('x e y devono avere la stessa lunghezza per la convoluzione ciclica.')

    Lz = Lx                      # nella ciclica la lunghezza resta la stessa
    t_z = t_x                    # stesso asse temporale di un periodo dei due segnali
    z = np.zeros(Lz)
    for n in range(Lz):
        for k in range(Ly):
            z[n] = z[n] + dt * x[(n - k) % Lx] * y[k]     # indice (n-k) "modulo" Lx: si riavvolge
    return z, t_z


Tp = 4
N = 800
dt = Tp / N                                         # 0.005 s
t_one = -Tp / 2 + np.arange(N) * dt                 # un periodo: da -2 a 1.995

x_one = (np.abs(t_one) < 0.7) * 1.0                 # rect di semi-larghezza 0.7
y_one = (np.abs(t_one) < 0.7) * 1.0

nPeriods = 7
t_rep = -nPeriods * Tp / 2 + np.arange(N * nPeriods) * dt      # 5600 campioni, da -14 a 13.995
x_rep = np.tile(x_one, nPeriods)                    # MATLAB repmat(x_one, 1, nPeriods)
y_rep = np.tile(y_one, nPeriods)

# convoluzione lineare dei segnali troncati: con i cicli sarebbero 11199*5600 = 62 milioni di passi, si usa np.convolve
z_lin = dt * np.convolve(x_rep, y_rep)
t_ylin = t_rep[0] + t_rep[0] + np.arange(len(z_lin)) * dt
z_circ, t_ycirc = calcolatore_convoluzione_ciclica(x_one, t_one, y_one, t_one)
z_circ_rep = np.tile(z_circ, nPeriods)

print(f"campioni a 1 in un periodo: {int(x_one.sum())}  (larghezza {dt * x_one.sum():.3f} s)")
print(f"lineare: {len(z_lin)} campioni da {t_ylin[0]:.0f} a {t_ylin[-1]:.3f} s, massimo {z_lin.max():.3f} in t = {t_ylin[z_lin.argmax()]:.2f}")
print(f"ciclica: {len(z_circ)} campioni, massimo {z_circ.max():.3f} nel campione di indice {z_circ.argmax()} (t = {t_ycirc[z_circ.argmax()]:.2f} secondo l'asse del file)")

Risultato:

campioni a 1 in un periodo: 279  (larghezza 1.395 s)
lineare: 11199 campioni da -28 a 27.990 s, massimo 9.765 in t = 0.00
ciclica: 800 campioni, massimo 1.395 nel campione di indice 0 (t = -2.00 secondo l'asse del file)

Con Tp=4T_p=4, N=800N=800 campioni e Δt=0,005\Delta t=0{,}005 s, il rect |t| < 0.7 ha 279279 campioni a 11 (i punti con ∣t∣=0,7|t|=0{,}7 esatto sono esclusi dalla disuguaglianza stretta), cioè larghezza numerica 1,395≈1,41{,}395\approx1{,}4. La convoluzione ciclica di due rect periodici di larghezza 1,41{,}4 in un periodo di 44 è una ripetizione periodica di triangoli, 1,4tri⁡t1,41{,}4\operatorname{tri}\frac{t}{1{,}4}: non nulla per ∣t∣<1,4|t|<1{,}4, con picco 1,41{,}4 (l'area di un rect); i triangoli di periodi adiacenti non si toccano perché 2⋅1,4=2,8<42\cdot1{,}4=2{,}8<4. Il picco numerico è 1,3951{,}395.

Grafico interattivo: Convoluzione ciclica di due rect di larghezza 1,4 con periodo 4: triangoli di base 2,8 e altezza 1,4 ripetuti ogni 4 secondi (tratteggiato: il rect periodico)

Lineare contro ciclica. Per la convoluzione lineare dei segnali troncati a 7 periodi non si possono usare i due cicli (sarebbero 11199×5600≈6211199\times5600\approx62 milioni di passi): si usa np.convolve, che dà lo stesso risultato (verificato nella parte 7). Il confronto:

python
# Rapporto tra lineare e ciclica negli istanti 0, 4, 8, 12 (picchi del lineare)
for m in range(4):
    i = np.argmin(np.abs(t_ylin - Tp * m))
    print(f"t = {Tp * m:>2}: lineare = {z_lin[i]:.3f},  lineare / ciclica(picco) = {z_lin[i] / z_circ.max():.1f}   (periodi sovrapposti: {nPeriods - m})")
t =  0: lineare = 9.765,  lineare / ciclica(picco) = 7.0   (periodi sovrapposti: 7)
t =  4: lineare = 8.370,  lineare / ciclica(picco) = 6.0   (periodi sovrapposti: 6)
t =  8: lineare = 6.975,  lineare / ciclica(picco) = 5.0   (periodi sovrapposti: 5)
t = 12: lineare = 5.580,  lineare / ciclica(picco) = 4.0   (periodi sovrapposti: 4)

La convoluzione lineare ha la stessa forma a triangoli, ma con un inviluppo triangolare: al tempo t=4mt=4m vale (7−∣m∣)(7-|m|) volte il picco della ciclica, perché lì si sovrappongono solo 7−∣m∣7-|m| periodi dei due segnali troncati, e ogni periodo sovrapposto aggiunge un contributo uguale. La ciclica integra su un solo periodo e vale 11 volta. Quindi le due curve hanno le stesse forme locali ma ampiezze diverse: la ciclica è la risposta "giusta" per segnali periodici (è la convoluzione del segnale periodico, e la sua serie di Fourier ha coefficienti prodotto, Sn(x⊛y)=TpXnYnS_n(x\circledast y)=T_pX_nY_n); la lineare di segnali troncati cresce con il numero di periodi.

Attenzione all'asse della convoluzione ciclica. Nel file MATLAB l'asse della ciclica è t_z = t_x, cioè va da −2-2 a 22. Ma, come per la convoluzione lineare, il campione nn corrisponde al tempo tx0+ty0+nΔt=−4+nΔtt_{x0}+t_{y0}+n\Delta t=-4+n\Delta t, che modulo il periodo 44 vale nΔtn\Delta t: il campione di indice 00 è il tempo 00 (dove i due rect si sovrappongono, cioè il picco), non −2-2. Eseguendo il codice, infatti, il massimo cade nel campione di indice 0 (ultima riga dell'output della parte 8), che l'asse del file chiama t=−2t=-2: in un confronto sovrapposto alla lineare, i picchi della ciclica sembrano spostati di mezzo periodo (2 s) rispetto a quelli della lineare. Per mettere il picco a t=0t=0 basta far scorrere il vettore di mezzo periodo con np.roll(z, N // 2):

python
# Correzione dell'asse della ciclica: il campione n corrisponde al tempo t_x[0] + t_y[0] + n dt = -4 + n dt,
# cioe' a n dt modulo 4. Portando il picco a t = 0 (meta' vettore) si ha l'asse di un periodo, da -2 a 2.
z_circ_centrata = np.roll(z_circ, N // 2)
print(f"dopo np.roll: massimo {z_circ_centrata.max():.3f} in t = {t_one[z_circ_centrata.argmax()]:.3f}")

# la ciclica e' anche ifft(fft(x) fft(y)) * dt (vedi laboratorio 3-4): stessi valori
z_fft = dt * np.real(np.fft.ifft(np.fft.fft(x_one) * np.fft.fft(y_one)))
print(f"differenza tra ciclica a cicli e ciclica con fft: {np.max(np.abs(z_circ - z_fft)):.2e}")

# confronto con il triangolo teorico 1.4 tri(t/1.4)
tri_teorico = np.maximum(0, 1.4 - np.abs(t_one))
print(f"errore massimo rispetto a 1.4 tri(t/1.4): {np.max(np.abs(z_circ_centrata - tri_teorico)):.4f}")
dopo np.roll: massimo 1.395 in t = 0.000
differenza tra ciclica a cicli e ciclica con fft: 2.98e-14
errore massimo rispetto a 1.4 tri(t/1.4): 0.0050

Il triangolo teorico 1,4tri⁡(t/1,4)1{,}4\operatorname{tri}(t/1{,}4) si discosta dal risultato numerico al più di 0,005=Δt0{,}005=\Delta t. L'ultimo controllo mostra anche che la convoluzione ciclica è la trasformata inversa del prodotto delle DFT, z=Δt⋅ifft(fft(x) fft(y))z=\Delta t\cdot\text{ifft}\bigl(\text{fft}(x)\,\text{fft}(y)\bigr) (Trasformata di Fourier discreta (DFT) e FFTUn segnale discreto periodico di periodo $NT$ ha una trasformata discreta e periodica, la DFT: $S(kF)=\sum_{n=0}^{N-1}T,s(nT)e^{-i2\pi kn/N}$, con $F=1/(NT)$, e $s(nT)=\sum_{k=0}^{N-1}F,S(kF)e^{i2\pi kn/N}$. Dipende da soli $N$ numeri, si calcola senza approssimazioni e con la FFT costa $N\log_2N$ invece di $N^2$. I coefficienti di Fourier del segnale periodico sono $S_k=F,S(kF)$. La DFT dà anche campioni della trasformata di un segnale discreto o continuo a durata limitata (con zero-padding), ma la scalatura $T$ e l'asse delle frequenze vanno gestiti con cura.Trasformata di Fourier discreta (DFT) e FFT →): è il motivo per cui la DFT si usa per calcolare le convoluzioni in tempo Nlog⁡NN\log N invece di N2N^2.

python
# grafici della parte 8
plt.figure(figsize=(8, 6))
plt.subplot(3, 1, 1)
plt.plot(t_rep, x_rep, linewidth=1.5); plt.ylim(-0.5, 1.5); plt.grid(True); plt.ylabel('x(t)')
plt.subplot(3, 1, 2)
plt.plot(t_ylin, z_lin, linewidth=1.5, label='lineare'); plt.grid(True); plt.ylabel('z_lin(t)')
plt.plot(t_rep, z_circ_rep, linewidth=1.5, label='ciclica (asse del file)')
plt.plot(t_rep, np.tile(z_circ_centrata, nPeriods), '--', linewidth=1.5, label='ciclica centrata')
plt.legend()
plt.xlabel('t [s]')
plt.tight_layout()
plt.show()

Controllo

  • Supporto della convoluzione: [−3,8, 3,8][-3{,}8,\,3{,}8], la somma degli estremi dei due supporti, con 20012001 campioni nell'asse [−10,10][-10,10].
  • Plateau 1,611{,}61 (campioni a 1 pesati con Δt\Delta t) contro 1,601{,}60 teorico, errore Δt=0,01\Delta t=0{,}01.
  • Area: 1,61×6,01=9,67611{,}61\times6{,}01=9{,}6761 uguale all'area di yy.
  • Ciclica contro triangolo teorico: errore massimo 0,0050{,}005; lineare = (7−∣m∣)×(7-|m|)\times ciclica nei picchi; ciclica = ifft del prodotto delle fft, errore 3⋅10−143\cdot10^{-14}.

Versione ripasso

  • Formula. zn=Δt∑kxn−kykz_n=\Delta t\sum_kx_{n-k}y_k (due cicli, o dt * np.convolve). Il fattore Δt\Delta t dà l'area e fa valere l'area di zz = prodotto delle aree.
  • Asse. Parte da tx0+ty0t_{x0}+t_{y0} e finisce a tx[−1]+ty[−1]t_{x}[-1]+t_{y}[-1], con Lx+Ly−1L_x+L_y-1 campioni (supporto = somma dei supporti).
  • Due rect (larghezza 1,61{,}6 e 66): trapezio con vertici in ±3,8\pm3{,}8 e ±2,2\pm2{,}2, plateau 1,61{,}6 (numerico 1,611{,}61: estremi contati come campioni interi). Commutativa verificata.
  • Ciclica. zn=Δt∑kx(n−k) mod Lykz_n=\Delta t\sum_kx_{(n-k)\bmod L}y_k: rect periodici di larghezza 1,41{,}4 e periodo 44 danno triangoli 1,4tri⁡t1,41{,}4\operatorname{tri}\frac t{1{,}4}; il campione 0 è il tempo 0 (con t_z = t_x i picchi si vedono spostati di mezzo periodo, serve np.roll).
  • Lineare troncata a 7 periodi: picco in 4m4m pari a (7−∣m∣)(7-\lvert m\rvert) volte la ciclica. Ciclica = ifft(fft·fft).

Teoria collegata