Esercizio - laboratorio 2, area ed energia numeriche

Testo (laboratorio 2 del corso Teoria dei Segnali, UniPD, file lab2.m, parti 1-3). (1) Calcolare l'area sottesa da $f(t)=k\,e^{-\alpha t}$ con $k=3$, $\alpha=0{,}1$ tra $0$ e $10$ s, con il metodo dei trapezi, con la somma dei campioni e con la formula analitica, e confrontarle. (2) Calcolare l'energia del segnale $s(t)=\sin(2\pi t)\,e^{-7(t-0{,}5)^2}$ (in volt) per $t\in[0,3]$ s con tre comandi: sum, trapezi e norm. (3) Calcolare la potenza su un periodo del segnale $s(t)=5\sin(2\pi f_1t+\frac\pi2)+3\cos(2\pi f_2t)+2\sin(2\pi f_3t)$ con $f_1=5$ Hz, $f_2=2f_1$, $f_3=3f_1$ e confrontarla con la potenza teorica.

Teoria usata: 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 →, Segnali a tempo continuo - definizioni e trasformazioniUn segnale a tempo continuo è una funzione complessa di variabile reale $s:\mathbb R\to\mathbb C$ (reale se $s=s^$). Si descrive con parte reale e immaginaria oppure con modulo e fase (valore principale in $(-\pi,\pi]$). Sui segnali si fanno tre operazioni elementari: ribaltamento $s(-t)$, traslazione $s(t-t_0)$ e cambio di scala $s(at)$. Ogni segnale si scompone in una parte pari e una dispari, in una reale e una immaginaria, in una hermitiana e una antihermitiana. Estensione e durata dicono dove il segnale è diverso da zero; i segnali con estensione nel semiasse positivo sono causali.Segnali a tempo continuo - definizioni e trasformazioni →, Serie di FourierUn segnale periodico di periodo $T_p$ si scrive come somma di esponenziali alle frequenze multiple della fondamentale $F=1/T_p$: $s(t)=\sum_nS_ne^{i2\pi nFt}$, con $S_n=\frac1{T_p}\int_{T_p}s(t)e^{-i2\pi nFt}dt$. Si basa sull'ortogonalità degli esponenziali su un periodo. $S_0$ è il valor medio; per segnali reali $S_{-n}=S_n^$ e si passa alla forma con coseni e seni; vale il teorema di Parseval $P=\sum|S_n|^2$. La convoluzione ciclica diventa il prodotto $T_pX_nY_n$ dei coefficienti. Le somme troncate presentano il fenomeno di Gibbs vicino ai salti.Serie di Fourier → (per la potenza delle armoniche), Segnali a tempo discretoUn segnale a tempo discreto è una funzione complessa $s(nT)$ definita sui multipli interi del quanto temporale $T$ (insieme $\mathbb Z(T)$, velocità $F_p=1/T$). Le definizioni sono quelle dei segnali continui con la somma al posto dell'integrale e il quanto $T$ al posto di $dt$: area $\sum T,s(nT)$, energia $\sum T|s(nT)|^2$, convoluzione $\sum T,x(kT)y(nT-kT)$. L'impulso ideale discreto vale $1/T$ nell'origine. Esponenziali e sinusoidi discreti sono periodici solo se $f_0/F_p$ è razionale e hanno frequenza ambigua a meno di multipli di $F_p$. I segnali periodici con periodo $NT$ sono descritti da $N$ valori e si trattano al calcolatore.Segnali a tempo discreto → (per l'area e l'energia di un vettore di campioni).

Cosa si vuole calcolare e perché

L'area di un segnale è $\int s(t)\,dt$, l'energia è $\int|s(t)|^2dt$, la potenza su un periodo $T_p$ è $\frac1{T_p}\int_{T_p}|s(t)|^2dt$ (definizioni in 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 →). Al calcolatore un integrale non si fa: si ha un vettore di campioni $s_n=s(t_n)$ presi con passo $\Delta t$, e l'integrale si approssima con una somma. Tutta la questione è il fattore di scala giusto: la somma $\sum_ns_n$ da sola dipende da quanti campioni si sono presi; l'area si ottiene moltiplicando per il passo $\Delta t$, perché ogni campione rappresenta una striscia di larghezza $\Delta t$. È lo stesso fattore $T$ che compare nelle definizioni dei segnali discreti del corso (area $\sum_nT\,s(nT)$, energia $\sum_nT|s(nT)|^2$).

Il codice usa la funzione dei trapezi, definita una volta e riusata nelle tre parti:

import numpy as np
import matplotlib.pyplot as plt

# Funzione "calcolatore_area_trapezi": integrale di f(t) con il metodo dei trapezi
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                                       # inizializzazione dell'integrale
    for a in range(len(f) - 1):                   # MATLAB: for a = 1:length(f)-1
        I = I + (f[a + 1] + f[a]) * (t[a + 1] - t[a]) / 2   # area del trapezio tra due campioni
    return I

Il ciclo scorre le $N-1$ coppie di campioni consecutivi: ciascuna delimita un trapezio di basi $f_a$ e $f_{a+1}$ e altezza $t_{a+1}-t_a$, di area $\frac{f_a+f_{a+1}}2\,(t_{a+1}-t_a)$; l'integrale è la somma di queste aree. In Python gli indici partono da 0, quindi range(len(f) - 1) percorre $0,\dots,N-2$ (in MATLAB 1:length(f)-1).

Parte 1 - area sottesa da $k\,e^{-\alpha t}$

# PARTE 1 - area sottesa da f(t) = k e^(-alpha t) tra 0 e 10
t_inizio = 0
t_fine = 10
dt = 0.05
t = t_inizio + np.arange(round((t_fine - t_inizio) / dt) + 1) * dt     # 201 campioni
alpha = 0.1
k = 3
f = k * np.exp(-alpha * t)

plt.figure()
plt.plot(t, f, color='red')
plt.xlabel('Tempo [s]')
plt.ylabel('Ampiezza')
plt.title('f(t)')
plt.xlim(t_inizio - 0.5, t_fine + 0.5)
plt.ylim(f.min() - 0.5, f.max() + 0.5)
plt.grid(True)
plt.show()

# metodo dei trapezi
A_trapezi = calcolatore_area_trapezi(f, t)
print(f"L'area calcolata con il metodo dei trapezi è: {A_trapezi:.12f}")

# integrale analitico: int_0^10 k e^(-alpha t) dt = (-k/alpha)(e^(-alpha*10) - e^0)
A = (-k / alpha) * (np.exp(-alpha * t_fine) - np.exp(-alpha * t_inizio))
print(f"L'area calcolata analiticamente è:           {A:.12f}")
print(f"La differenza (analitica - trapezi) è:       {A - A_trapezi:.12f}")

# metodo dei rettangoli: dt * somma dei campioni, tolto l'ultimo (somma sinistra)
A_sum = dt * np.sum(f) - f[-1] * dt
print(f"L'area calcolata con la somma dei rettangoli è: {A_sum:.12f}")
print(f"La differenza (analitica - rettangoli) è:    {A - A_sum:.12f}")

Risultato:

L'area calcolata con il metodo dei trapezi è: 18.963656272375
L'area calcolata analiticamente è:           18.963616764857
La differenza (analitica - trapezi) è:       -0.000039507518
L'area calcolata con la somma dei rettangoli è: 19.011065314287
La differenza (analitica - rettangoli) è:    -0.047448549431

Il valore esatto. $\displaystyle\int_0^{10}3\,e^{-0{,}1t}dt=\frac{3}{0{,}1}\bigl(1-e^{-1}\bigr)=30\,(1-0{,}36788)=18{,}96362$ (Integrali immediatiTabella degli integrali immediati, ottenuti leggendo al contrario le derivate delle funzioni elementari: potenze, 1/x, esponenziali, seno e coseno, tangente, arcoseno, arcotangente, funzioni iperboliche e loro inverse. Per gli integrali non ci sono regole come per le derivate, solo metodi: linearità, sostituzione, per parti.Integrali immediati →). Il codice scrive la stessa cosa come $(-k/\alpha)(e^{-\alpha t_{fine}}-e^{-\alpha t_{inizio}})$.

Quanto sbagliano i due metodi. Con $\Delta t=0{,}05$ (201 campioni):

• i trapezi differiscono dal valore esatto di $4\cdot10^{-5}$ (relativamente, due parti su un milione);
• la somma dei rettangoli, dt * sum(f) - f(end) * dt, cioè $\Delta t\sum_{n=0}^{N-1}f_n$ (si toglie l'ultimo campione perché il rettangolo dell'ultimo campione sporgerebbe oltre $t=10$), differisce di $0{,}047$ (quasi $0{,}25\%$): mille volte peggio.

Il perché è geometrico. La funzione è decrescente: il rettangolo costruito sul valore a sinistra sta sempre sopra la curva e in ogni intervallo sovrastima di circa $\frac{\Delta t}2\,(f_n-f_{n+1})$ (metà della variazione di $f$ nell'intervallo, per la larghezza); sommando le variazioni di tutti gli intervalli si ottiene in totale $\frac{\Delta t}2\,(f(0)-f(10))=0{,}025\,(3-1{,}104)=0{,}0474$, che è l'errore trovato. Il trapezio invece sostituisce la curva con la corda, e l'errore si riduce a $\frac{\Delta t^2}{12}\bigl(f'(10)-f'(0)\bigr)=\frac{0{,}0025}{12}\cdot0{,}1896=3{,}95\cdot10^{-5}$, che è proprio il numero stampato: dimezzando il passo l'errore dei rettangoli si dimezza, quello dei trapezi diventa un quarto.

Una scrittura equivalente. Il metodo dei trapezi è la somma dei campioni con gli estremi contati a metà:

$$\int f\,dt\approx\Delta t\Bigl(\sum_{n=0}^{N-1}f_n-\frac{f_0+f_{N-1}}2\Bigr).$$

Infatti ogni campione interno appartiene a due trapezi, ciascuno gli assegna metà della striscia, quindi in totale $\Delta t\,f_n$; i due campioni estremi appartengono a un trapezio solo. Il blocco seguente confronta tre scritture (senza ciclo, np.trapezoid di numpy, formula con gli estremi a metà): danno lo stesso numero fino all'ultima cifra.

# Altre due scritture dei trapezi
A_vett = np.sum((f[1:] + f[:-1]) * np.diff(t)) / 2           # senza ciclo
A_numpy = np.trapezoid(f, t)                                  # funzione di numpy (np.trapz nelle versioni vecchie)
A_pesi = dt * (np.sum(f) - (f[0] + f[-1]) / 2)               # somma dei campioni con i due estremi a meta'
print(A_vett, A_numpy, A_pesi)
print(np.allclose([A_vett, A_numpy, A_pesi], A_trapezi))

18.963656272375196 18.963656272375196 18.963656272375193
True

Parte 2 - energia di $s(t)=\sin(2\pi t)\,e^{-7(t-0{,}5)^2}$

# PARTE 2 - energia di s(t) = sin(2 pi t) e^(-7 (t-0.5)^2), t tra 0 e 3 secondi
dt = 0.001
t = np.arange(round(3 / dt) + 1) * dt            # 3001 campioni, da 0 a 3
s = np.sin(2 * np.pi * t) * np.exp(-7 * (t - 0.5) ** 2)       # [V]

plt.figure()
plt.plot(t, s, 'r--')
plt.grid(True)
plt.show()

energia_s_sum = dt * np.sum(np.abs(s) ** 2)                       # [V^2 s]
energia_s_trapezi = calcolatore_area_trapezi(np.abs(s) ** 2, t)
energia_s_norm = dt * np.linalg.norm(s) ** 2

print(f"L'energia del segnale s(t) (con la somma)  è: {energia_s_sum:.4f} V^2*s")
print(f"L'energia del segnale s(t) (con i trapezi) è: {energia_s_trapezi:.4f} V^2*s")
print(f"L'energia del segnale s(t) (con norm)      è: {energia_s_norm:.4f} V^2*s")

Risultato:

L'energia del segnale s(t) (con la somma)  è: 0.2224 V^2*s
L'energia del segnale s(t) (con i trapezi) è: 0.2224 V^2*s
L'energia del segnale s(t) (con norm)      è: 0.2224 V^2*s

L'energia è $E=\int|s|^2dt\approx\Delta t\sum_n|s_n|^2$, in $\text{V}^2\cdot\text{s}$ perché $s$ è in volt. Il modulo quadro serve per i segnali complessi (per un segnale reale come questo è semplicemente $s_n^2$). I tre comandi sono tre scritture della stessa somma:

• dt * np.sum(np.abs(s)**2): somma dei quadrati per $\Delta t$;
• calcolatore_area_trapezi(np.abs(s)**2, t): trapezi applicati al segnale $|s|^2$;
• dt * np.linalg.norm(s)**2: np.linalg.norm(s) è la norma euclidea $\sqrt{\sum_ns_n^2}$ del vettore, il suo quadrato è la somma dei quadrati, da moltiplicare per $\Delta t$.

Danno lo stesso valore, $0{,}2224$, perché agli estremi il segnale vale (quasi) zero: $s(0)=\sin0\cdot e^{-7/4}=0$ e $s(3)\approx0$, e quando gli estremi contano zero la somma coincide coi trapezi (vedi l'identità della parte 1). Il segnale ha energia finita ed è concentrato intorno a $t=0{,}5$: è un seno "ritagliato" da una campana gaussiana larga circa $\pm0{,}3$ s.

Controllo con la formula chiusa. Si usa $\sin^2x=\frac12(1-\cos2x)$ e si pone $u=t-\frac12$: $\cos(4\pi t)=\cos(4\pi u+2\pi)=\cos4\pi u$. Allora, su tutto l'asse reale, usando $\int e^{-14u^2}du=\sqrt{\pi/14}$ e $\int\cos(4\pi u)e^{-14u^2}du=\sqrt{\pi/14}\,e^{-(4\pi)^2/56}$,

$$E=\frac12\sqrt{\frac\pi{14}}\Bigl(1-e^{-2\pi^2/7}\Bigr)\approx0{,}2227.$$

# Riscontro con l'integrale "vero" (scipy) e con la formula chiusa su tutto l'asse reale
from scipy.integrate import quad
E_quad, _ = quad(lambda u: (np.sin(2 * np.pi * u) * np.exp(-7 * (u - 0.5) ** 2)) ** 2, 0, 3, points=[0.5, 1, 2])
E_chiusa = 0.5 * np.sqrt(np.pi / 14) * (1 - np.exp(-2 * np.pi ** 2 / 7))
print(f"somma dei campioni con 6 cifre:   {energia_s_sum:.6f}")
print(f"integrale numerico preciso su [0,3]: {E_quad:.6f}")
print(f"formula chiusa su tutto l'asse reale: {E_chiusa:.6f}")
print(f"energia fuori da [0,3] (quasi tutta per t < 0): {E_chiusa - E_quad:.6f}")

somma dei campioni con 6 cifre:   0.222391
integrale numerico preciso su [0,3]: 0.222391
formula chiusa su tutto l'asse reale: 0.222735
energia fuori da [0,3] (quasi tutta per t < 0): 0.000344

I $0{,}2224$ del laboratorio e i $0{,}2227$ della formula chiusa differiscono di $3{,}4\cdot10^{-4}$: è l'energia della parte di segnale per $t<0$, che il vettore (da $0$ a $3$) non contiene. L'integrale numerico preciso su $[0,3]$ (scipy.integrate.quad) riproduce a sei cifre la somma dei campioni ($0{,}222391$), quindi con $\Delta t=1$ ms la somma è affidabile: il segnale è liscio e va a zero agli estremi.

Parte 3 - potenza su un periodo di una somma di sinusoidi

# PARTE 3 - potenza su un periodo di una somma di tre sinusoidi
dt = 0.001
t = np.arange(round(3 / dt) + 1) * dt

f1 = 5
f2 = 2 * f1
f3 = 3 * f1

s = 5 * np.sin(2 * np.pi * f1 * t + np.pi / 2) + 3 * np.cos(2 * np.pi * f2 * t) + 2 * np.sin(2 * np.pi * f3 * t)   # [V]

plt.figure()
plt.plot(t, s)
plt.show()

periodo = 1 / f1                                   # 0.2 s: il piu' lungo dei tre periodi
n_punti_periodo = round(periodo / dt)              # 200 campioni (round e non floor: 0.2/0.001 e' decimale)
potenza_periodo = np.sum(np.abs(s[:n_punti_periodo]) ** 2) / n_punti_periodo

print(f"campioni per periodo: {n_punti_periodo}")
print(f"La potenza media sul periodo del segnale è: {potenza_periodo:.4f} V^2")
print(f"potenza teorica 1/2 (5^2 + 3^2 + 2^2) = {0.5 * (5**2 + 3**2 + 2**2):.4f}")

Risultato:

campioni per periodo: 200
La potenza media sul periodo del segnale è: 19.0000 V^2
potenza teorica 1/2 (5^2 + 3^2 + 2^2) = 19.0000

Quale periodo. Le tre sinusoidi hanno frequenze $5$, $10$, $15$ Hz, cioè $f_1$, $2f_1$, $3f_1$: sono la fondamentale e due armoniche. Il segnale somma si ripete quando si ripete la fondamentale, quindi ha periodo $T_p=\frac1{f_1}=0{,}2$ s, che a $\Delta t=1$ ms vale $N_p=\frac{T_p}{\Delta t}=200$ campioni. Nel codice si usa round(periodo/dt) al posto del floor del file MATLAB: il rapporto di due numeri decimali può venire $199{,}999\ldots$ e floor darebbe 199; con round si è al sicuro.

Perché si divide per il numero di campioni. Per definizione $P=\frac1{T_p}\int_{T_p}|s|^2dt\approx\frac1{T_p}\,\Delta t\sum_{n=0}^{N_p-1}|s_n|^2$, ma $T_p=N_p\Delta t$, quindi $\Delta t$ si semplifica e resta $P\approx\frac1{N_p}\sum_n|s_n|^2$: la potenza è la media dei campioni al quadrato (l'energia di un periodo meno la lunghezza del periodo).

Il valore teorico. Una sinusoide di ampiezza $A$ ha potenza $\frac{A^2}2$ (il quadrato oscilla tra $0$ e $A^2$ con media $\frac{A^2}2$); la fase non conta ($\sin(x+\frac\pi2)=\cos x$ ha lo stesso quadrato medio). Le tre componenti sono armoniche diverse della stessa fondamentale, quindi sono ortogonali sul periodo: nel quadrato della somma i termini misti ($2\cdot5\cdot3\sin\cdot\cos$ e simili) hanno media zero (Serie di FourierUn segnale periodico di periodo $T_p$ si scrive come somma di esponenziali alle frequenze multiple della fondamentale $F=1/T_p$: $s(t)=\sum_nS_ne^{i2\pi nFt}$, con $S_n=\frac1{T_p}\int_{T_p}s(t)e^{-i2\pi nFt}dt$. Si basa sull'ortogonalità degli esponenziali su un periodo. $S_0$ è il valor medio; per segnali reali $S_{-n}=S_n^*$ e si passa alla forma con coseni e seni; vale il teorema di Parseval $P=\sum|S_n|^2$. La convoluzione ciclica diventa il prodotto $T_pX_nY_n$ dei coefficienti. Le somme troncate presentano il fenomeno di Gibbs vicino ai salti.Serie di Fourier →). Resta solo la somma delle potenze:

$$P=\frac{5^2}2+\frac{3^2}2+\frac{2^2}2=12{,}5+4{,}5+2=19\ \text{V}^2.$$

Il calcolo numerico dà $19{,}0000$: l'ortogonalità vale esattamente anche sui campioni, perché 200 campioni contengono un numero intero di periodi di tutte e tre le sinusoidi.

# Cosa cambia se si media su tutto l'intervallo, o su un periodo non intero
print(f"media su tutti i 3001 campioni: {np.mean(s ** 2):.4f}")
print(f"media su 150 campioni (3/4 di periodo): {np.mean(s[:150] ** 2):.4f}")
print(f"media del quadrato di UNA sola sinusoide, 5 sin(2 pi 5 t + pi/2): {np.mean((5 * np.sin(2 * np.pi * f1 * t + np.pi / 2))[:200] ** 2):.4f}")

media su tutti i 3001 campioni: 19.0150
media su 150 campioni (3/4 di periodo): 20.7360
media del quadrato di UNA sola sinusoide, 5 sin(2 pi 5 t + pi/2): 12.5000

Se invece si media su una finestra che non è un periodo intero i termini misti non si annullano e il risultato è sbagliato: $20{,}74$ sui primi $150$ campioni (tre quarti di periodo). Mediando su tutti i 3001 campioni (15 periodi esatti più un campione in più) si ottiene $19{,}015$, quasi giusto; e la media di una sola sinusoide, $5\sin(\cdot)$, dà $12{,}5$, il primo addendo.

Controllo

I tre metodi di integrazione concordano con la formula chiusa (parte 1: $18{,}9636$ contro $18{,}9636$; parte 2: $0{,}2224$ contro l'integrale numerico $0{,}222391$; parte 3: $19{,}0000$ contro $\frac12(25+9+4)=19$). Gli errori trovati sono quelli previsti dalla teoria: $\frac{\Delta t}2\Delta f$ per i rettangoli e $\frac{\Delta t^2}{12}\Delta f'$ per i trapezi.