Esercizio - processi aleatori, media d'insieme, correlazione e potenza temporale in Python (laboratorio 6)

Questa pagina non ha ancora la versione ripasso: qui sotto c'è il testo completo.

Testo (laboratorio 6, parti 5-7, corso Teoria dei Segnali, UniPD). Con $f_0=5$ Hz si considerano i tre processi aleatori a tempo continuo $$x_1(t)=A\cos(2\pi f_0t),\qquad x_2(t)=\cos(2\pi f_0t+\varphi),\qquad x_3(t)=A\cos(2\pi f_0t+\varphi),$$ con $A\sim U[-1,1]$ e $\varphi\sim U[0,2\pi)$ indipendenti. Si osserva ciascun processo su $T_{obs}=1$ s, campionato a $F_s=500$ Hz ($500$ istanti, cioè esattamente $5$ periodi), con $M=5000$ realizzazioni. Per ciascun processo:

1. disegnare alcune realizzazioni;
2. calcolare la media d'insieme $m_x(t)$ (media sulle realizzazioni a $t$ fissato);
3. calcolare la correlazione d'insieme $r_x(t_1,t_1+\tau)$ per due istanti $t_1$ diversi ($0{,}02$ s e $0{,}16$ s);
4. calcolare la potenza temporale di ogni realizzazione e la correlazione temporale di alcune realizzazioni;
5. concludere su stazionarietà in media, stazionarietà in correlazione ed ergodicità.

Teoria usata: Processi aleatori - definizioni, media e autocorrelazioneUn processo aleatorio x(t), t in I (R o Z(T)), è una famiglia di variabili aleatorie sullo stesso spazio di probabilità; fissato l'esito si ottiene una realizzazione (un segnale). Si descrive con le densità di ordine N (complete), in particolare del primo e del secondo ordine, oppure solo con media m_x(t) e correlazione r_x(t,s) = E[x(t)x*(s)] (descrizione di potenza). Un processo gaussiano è determinato da media e correlazione.Processi aleatori - definizioni, media e autocorrelazione →, Processi stazionari e densità spettrale di potenzaUn processo è stazionario in media, in potenza, nel primo ordine o in correlazione se la rispettiva grandezza non cambia traslando il tempo; stazionario in senso lato (sl) = media costante e r_x(t,s) = r_x(t-s); in senso stretto (ss) = tutte le densità invarianti. La densità spettrale di potenza R_x(f) = F[r_x(τ)] è non negativa e ha integrale = potenza statistica. Rumore bianco: R costante. Ciclostazionario: statistiche periodiche. Ergodico in media: la media temporale di una realizzazione converge a m_x.Processi stazionari e densità spettrale di potenza →, Valore attesoIl valore atteso E[X] = Σ x p_X(x) è la media dei valori di X pesata con le loro probabilità (esiste se la serie converge assolutamente); per una funzione g vale E[g(X)] = Σ g(x) p_X(x) senza trovare la legge di g(X), ed E è lineare: E[aX + bY + c] = aE[X] + bE[Y] + c.Valore atteso →, Distribuzioni uniforme continua ed esponenzialeU(a, b) ha densità costante 1/(b − a) su [a, b], media (a + b)/2 e varianza (b − a)²/12; Exp(λ) ha densità λe^(−λx) per x ≥ 0, FdD 1 − e^(−λx), P(X > t) = e^(−λt), media 1/λ, varianza 1/λ², ed è l'unica legge continua senza memoria (versione continua della geometrica).Distribuzioni uniforme continua ed esponenziale →, 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 →.

(1) Due tipi di media

Un processo è un insieme di realizzazioni (segnali) $x(t;\omega)$: ogni colonna della matrice del codice è una realizzazione, ogni riga è un istante. Si possono fare due medie diverse:

• media d'insieme (statistica): fissato $t$, si media sulle realizzazioni, cioè si stima $E[x(t)]$ e $r_x(t_1,t_2)=E[x(t_1)x(t_2)]$ (qui reali, quindi senza coniugio);
• media temporale: fissata una realizzazione, si media su $t$: la potenza $P=\lim\frac1{T}\int x^2\,dt$ e la correlazione temporale $\overline{x(t)x(t+\tau)}$ (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 →).

Se le due medie coincidono il processo si dice ergodico: allora una sola realizzazione, osservata a lungo, basta a conoscere le statistiche. Per questo servono prima le statistiche d'insieme (calcolo analitico) e poi il confronto con le medie temporali (simulazione).

(2) Calcolo analitico

Si usa $E[A]=0$, $E[A^2]=\int_{-1}^1\frac{a^2}2da=\frac13$, $E[\cos(c+\varphi)]=0$ e $E[\cos(c+2\varphi)]=0$ per ogni $c$ (media di un coseno su un periodo intero, $\varphi$ uniforme). Si pone $\omega_0=2\pi f_0$ e $\tau=t_1-t_2$.

Processo 1, $x_1=A\cos\omega_0t$ (casuale solo l'ampiezza):

• media: $m_1(t)=E[A]\cos\omega_0t=0$ ($t$ qualunque);
• correlazione: $r_1(t_1,t_2)=E[A^2]\cos\omega_0t_1\cos\omega_0t_2=\frac13\cos\omega_0t_1\cos\omega_0t_2=\frac16\left[\cos\omega_0(t_1-t_2)+\cos\omega_0(t_1+t_2)\right]$;
• potenza statistica: $M_1(t)=r_1(t,t)=\frac13\cos^2\omega_0t$, che dipende da $t$.

La media è costante ($0$), ma la correlazione dipende da $t_1+t_2$ e non solo da $t_1-t_2$: il processo è stazionario in media ma non in correlazione, quindi non è nemmeno stazionario in senso lato.

Processo 2, $x_2=\cos(\omega_0t+\varphi)$ (casuale solo la fase):

• $m_2=0$;
• $r_2(t_1,t_2)=E[\cos(\omega_0t_1+\varphi)\cos(\omega_0t_2+\varphi)]=\frac12\cos\omega_0(t_1-t_2)+\frac12E[\cos(\omega_0(t_1+t_2)+2\varphi)]=\frac12\cos\omega_0\tau$;
• $M_2=r_2(0)=\frac12$ costante.

Media costante e correlazione funzione di $\tau$: stazionario in senso lato. È anche stazionario in senso stretto: traslando il tempo di $t_0$ si ha $\cos(\omega_0t+\omega_0t_0+\varphi)$, e $\varphi+\omega_0t_0$ (modulo $2\pi$) è ancora uniforme su $[0,2\pi)$, quindi tutte le densità congiunte restano le stesse.

Processo 3, $x_3=A\cos(\omega_0t+\varphi)$ ($A$ e $\varphi$ indipendenti):

• $m_3=E[A]\,E[\cos(\omega_0t+\varphi)]=0$;
• $r_3(t_1,t_2)=E[A^2]\cdot\frac12\cos\omega_0\tau=\frac16\cos\omega_0\tau$;
• $M_3=\frac16$.

Anche questo è stazionario in senso lato, e in senso stretto per lo stesso motivo (la fase uniforme "assorbe" la traslazione).

Potenze temporali e correlazioni temporali. Ogni realizzazione è un coseno di ampiezza $a$ e fase $\theta$ fissate; la media di $\cos^2$ su un periodo intero è $\frac12$ e quella di $\cos(\omega_0t+\theta)\cos(\omega_0(t+\tau)+\theta)$ è $\frac12\cos\omega_0\tau$, quindi:

Ergodicità.

• Processo 2: potenza temporale $\frac12=M_2$ e correlazione temporale $\frac12\cos\omega_0\tau=r_2(\tau)$ per ogni realizzazione: ergodico in media e in correlazione.
• Processo 3: ergodico nella media (ogni media temporale è $0=m_3$), ma non nella correlazione: la correlazione temporale $\frac{A_k^2}2\cos\omega_0\tau$ cambia con $k$ (ampiezza da $0$ a $\frac12$), mentre l'insieme ha $\frac16\cos\omega_0\tau$. Si ritrova il valore d'insieme solo mediando sulle realizzazioni: $E[A^2/2]=\frac16$. Il processo è stazionario ma non ergodico.
• Processo 1: non stazionario in correlazione; la media temporale è $0$ come quella d'insieme, ma la potenza $P_k=A_k^2/2$ è casuale e la potenza statistica $M_1(t)=\frac13\cos^2\omega_0t$ oscilla fra $0$ e $\frac13$ (valore medio sul periodo $\frac16$): l'ergodicità non ha senso senza stazionarietà.

Per i processi 2 e 3 la densità spettrale di potenza (Processi stazionari e densità spettrale di potenzaUn processo è stazionario in media, in potenza, nel primo ordine o in correlazione se la rispettiva grandezza non cambia traslando il tempo; stazionario in senso lato (sl) = media costante e r_x(t,s) = r_x(t-s); in senso stretto (ss) = tutte le densità invarianti. La densità spettrale di potenza R_x(f) = F[r_x(τ)] è non negativa e ha integrale = potenza statistica. Rumore bianco: R costante. Ciclostazionario: statistiche periodiche. Ergodico in media: la media temporale di una realizzazione converge a m_x.Processi stazionari e densità spettrale di potenza →) è la trasformata di $r(\tau)$: due righe in $\pm f_0$, $\frac14\left[\delta(f-f_0)+\delta(f+f_0)\right]$ per il processo 2 (potenza $\frac12$) e $\frac1{12}\left[\delta(f-f_0)+\delta(f+f_0)\right]$ per il processo 3 (potenza $\frac16$).

(3) Simulazione in Python

La matrice X[i, k] contiene il valore all'istante t[i] della realizzazione k. Le due medie sono semplicemente due assi di mean: axis=1 (sulle realizzazioni, per ogni t: media d'insieme) e axis=0 (sul tempo, per ogni realizzazione: media temporale). Il laboratorio sceglie gli istanti $t_1$ con gli indici $10$ e $80$ (Matlab, base $1$: $t(10)=0{,}018$ s, $t(80)=0{,}158$ s); qui con gli indici Python $10$ e $80$ si hanno $0{,}02$ e $0{,}16$ s, differenza di un campione senza effetto sulle conclusioni.

import numpy as np
import matplotlib.pyplot as plt

rng = np.random.default_rng(5)
f0, Fs, Tobs, M = 5, 500, 1.0, 5000              # M = numero di realizzazioni
t = np.arange(0, Tobs, 1 / Fs)                   # 500 istanti da 0 a 0.998 s: esattamente 5 periodi di f0
tc = t[:, None]                                  # colonna: tempo (righe); le realizzazioni sono le colonne

A = rng.uniform(-1, 1, M)                        # A ~ U[-1,1]
phi = rng.uniform(0, 2 * np.pi, M)               # phi ~ U[0,2pi)
X = {
    1: A * np.cos(2 * np.pi * f0 * tc),                    # A cos(2 pi f0 t)
    2: np.cos(2 * np.pi * f0 * tc + phi),                  # cos(2 pi f0 t + phi)
    3: A * np.cos(2 * np.pi * f0 * tc + phi),              # A cos(2 pi f0 t + phi)
}

def r_teo(caso, t1, t2):                         # correlazione d'insieme teorica r(t1, t2)
    w = 2 * np.pi * f0
    if caso == 1:
        return np.cos(w * t1) * np.cos(w * t2) / 3         # E[A^2] = 1/3
    if caso == 2:
        return 0.5 * np.cos(w * (t1 - t2))
    return np.cos(w * (t1 - t2)) / 6                       # E[A^2]/2 = 1/6

i1s, tau_max = (10, 80), 150                     # t1 = 0.02 s e 0.16 s; lag fino a 150 campioni = 0.3 s
for k, Xk in X.items():
    m_t = Xk.mean(axis=1)                        # media d'insieme: media sulle realizzazioni, per ogni t
    P = (Xk**2).mean(axis=0)                     # potenza temporale: media nel tempo, per ogni realizzazione
    print(f"--- caso {k}: max|m(t)| = {np.abs(m_t).max():.4f}; potenza temporale: "
          f"min {P.min():.4f}, max {P.max():.4f}, media {P.mean():.4f}, dev.st. {P.std():.4f}")
    for i1 in i1s:
        r = np.array([np.mean(Xk[i1] * Xk[i1 + lag]) for lag in range(tau_max + 1)])   # media sulle realizzazioni
        teo = r_teo(k, t[i1], t[i1 : i1 + tau_max + 1])
        print(f"   t1={t[i1]:.2f} s: r(t1,t1)={r[0]:.4f} (teoria {teo[0]:.4f}); "
              f"r(t1,t1+0.1 s)={r[50]:.4f} (teoria {teo[50]:.4f}); errore max sui lag {np.abs(r - teo).max():.4f}")
    for j in range(3):                           # correlazione temporale (divisa per N) di tre realizzazioni
        x = Xk[:, j]
        rt = np.array([x[: len(x) - l] @ x[l:] / len(x) for l in range(tau_max + 1)])
        att = (A[j]**2 / 2 if k != 2 else 0.5) * (1 - np.arange(tau_max + 1) / Fs / Tobs) * np.cos(2 * np.pi * f0 * np.arange(tau_max + 1) / Fs)
        print(f"   realiz. {j}: r_t(0)={rt[0]:.4f}, r_t(0.1 s)={rt[50]:.4f}, r_t(0.2 s)={rt[100]:.4f}; "
              + (f"A^2/2={A[j]**2 / 2:.4f}; " if k != 2 else "") + f"scarto da P(1-tau/T)cos: {np.abs(rt - att).max():.4f}")
plt.show()

Nota sulla correlazione temporale. Sommando solo sui campioni disponibili e dividendo per $N$ (come xcorr(..., 'biased')) si ottiene $\hat r_t(\tau)\approx P_k\left(1-\frac{|\tau|}{T_{obs}}\right)\cos\omega_0\tau$: la finestra di $1$ s accorcia la somma di $\tau$ secondi, e il fattore $1-\frac{|\tau|}{T_{obs}}$ è l'effetto della finestra (triangolare). Per $\tau=0$ è esatta; per $\tau=0{,}1$ s e $\tau=0{,}2$ s (multipli di mezzo periodo) il termine residuo si annulla e il fattore è esatto, in altri lag c'è una piccola correzione dipendente dalla fase (fino a $0{,}016$ in questa simulazione).

(4) Risultati e conclusioni

Cosa si vede.

• Realizzazioni: nel caso 1 sono coseni in fase tra loro (tutti passano per zero negli stessi istanti) con ampiezza casuale in $[-1,1]$; nel caso 2 sono coseni di ampiezza $1$ con fasi diverse; nel caso 3 coseni con ampiezza e fase diverse. Le tre nuvole di curve oscillano attorno allo zero a $5$ Hz.
• Media d'insieme: nei tre casi è una curva praticamente nulla; il massimo scarto da zero su $500$ istanti è $0{,}013$ (caso 1), $0{,}013$ (caso 2), $0{,}010$ (caso 3), dello stesso ordine dell'errore standard $\sigma_x/\sqrt M$, che vale $0{,}006$-$0{,}010$ per $M=5000$.
• Potenza temporale: caso 2, tutte le $5000$ realizzazioni hanno $P=0{,}5000$ (punti allineati su una retta orizzontale); casi 1 e 3, i punti riempiono la fascia fra $0$ e $0{,}5$ (le potenze sono $A_k^2/2$, quindi molte vicine a $0$: la densità di $A^2/2$ è più alta vicino a $0$), con media $0{,}168$ (teoria $\frac16=0{,}1667$) e deviazione standard $0{,}151$.
• Correlazione d'insieme ($\tau$ da $0$ a $0{,}3$ s):

Nel caso 1 le due curve (per $t_1=0{,}02$ e $0{,}16$ s) sono diverse (ampiezze $0{,}22$ e $0{,}03$): la correlazione dipende da $t_1$. Nei casi 2 e 3 le due curve coincidono e sono un coseno di ampiezza $0{,}5$ e $\frac16$. L'errore massimo sui $151$ lag è $\le0{,}0026$.

• Correlazione temporale: caso 2, tutte le realizzazioni hanno $\hat r_t(0)=0{,}5000$, $\hat r_t(0{,}1)=-0{,}4500$, $\hat r_t(0{,}2)=0{,}4000$ (cioè $\frac12\cos\omega_0\tau\,(1-\tau/T_{obs})$) e sono sovrapposte alla correlazione d'insieme; casi 1 e 3, le curve hanno ampiezze diverse: per esempio le prime tre realizzazioni hanno $\hat r_t(0)=0{,}1861;\ 0{,}1897;\ 0{,}0005$, uguali a $A_k^2/2$ ($A_0=0{,}610$, $A_1=0{,}616$, $A_2=0{,}031$) e diverse dal valore d'insieme $0{,}1667$.

Conclusioni.

Il punto del laboratorio: per valutare la potenza di un processo da una sola realizzazione serve l'ergodicità, non solo la stazionarietà. Un processo come il 3, stazionario ma non ergodico, ha la stessa $r(\tau)$ del 2 a parte l'ampiezza, ma ogni realizzazione ha una sua potenza.

Errori tipici.

• Confondere le due medie: mean(axis=1) (d'insieme) e mean(axis=0) (temporale); la stessa matrice, risultati con significato diverso.
• Dedurre l'ergodicità dalla sola stazionarietà (processo 3) o dalla media nulla (processo 1).
• Dimenticare che la correlazione temporale stimata su una finestra finita contiene il fattore $1-|\tau|/T_{obs}$ e confrontarla senza correzione con $\frac12\cos\omega_0\tau$.

Vedi anche: Esercizio - rumore bianco gaussiano e filtro passa-basso in Python (laboratorio 6), Processi aleatori stazionari e densità spettrale di potenzaUn processo aleatorio è un segnale i cui valori a ogni istante sono variabili aleatorie. Se è stazionario in senso lato (WSS) la media è costante e l'autocorrelazione $r_x(\tau)$ dipende solo dalla differenza dei tempi; la sua trasformata è la densità spettrale di potenza $\mathcal P_x(f)$, il cui integrale è la potenza statistica $r_x(0)$. Un filtro LTI dà $m_y=m_xH(0)$ e $\mathcal P_y=\mathcal P_x\lvert H\rvert^2$; se l'ingresso è gaussiano anche l'uscita lo è. Il rumore bianco ha $\mathcal P(f)=\frac{N_0}2$.Processi aleatori stazionari e densità spettrale di potenza →.