Esercizio - rumore bianco gaussiano e filtro passa-basso in Python (laboratorio 6)

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

Testo (laboratorio 6, parte 8, corso Teoria dei Segnali, UniPD). Un processo a tempo discreto $x(nT)$ con $F_s=1/T=1000$ Hz è un rumore bianco gaussiano: campioni indipendenti $N(0,\sigma^2)$, $\sigma=1$, $N=2^{17}$ campioni. Passa in un filtro FIR passa-basso con frequenza di taglio $f_c=80$ Hz e $129$ coefficienti (ordine $L=128$), progettato con il metodo della finestra. Si chiede di:

1. generare il rumore e rappresentare $x$ e l'uscita $y$ nel tempo;
2. stimare la correlazione empirica di $x$ e $y$;
3. stimare la densità spettrale di potenza (PSD) bilatera con il metodo di Welch;
4. confrontare con la teoria: $r_x$, $R_x$, $R_y=\lvert H\rvert^2R_x$, potenza dell'uscita.

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 →, Processi aleatori attraverso sistemi LTIUna trasformazione applicata a ogni realizzazione dà un nuovo processo, ma dalla descrizione statistica di x non sempre si ricava quella di y (controesempi); si può se ogni vettore di y dipende da un vettore finito di x. Per una tf lineare m_y e r_y si calcolano col nucleo; un filtro LTI conserva la stazionarietà in senso lato: m_y = G(0) m_x, R_y = |G|² R_x. Il campionamento di un processo sl dà un processo sl con PSD ripetuta periodicamente. L'interpolazione LTI di un processo discreto sl è ciclostazionaria. Teorema del campionamento per processi: R_x nulla fuori da (-B,B) e Fc ≥ 2B danno ricostruzione esatta in media quadratica; altrimenti la potenza dell'errore è in banda + fuori banda.Processi aleatori attraverso sistemi LTI →, Risposta in frequenza e filtriLa risposta in frequenza G(f) = F[g] di un filtro LTI dà Y = G·X. Gli esponenziali complessi sono autofunzioni (autovalore G(f)), quindi un ingresso sinusoidale esce sinusoidale con ampiezza moltiplicata per |G(f0)| e fase aumentata di arg G(f0). Il filtro è reale se G è hermitiana, invertibile se G non si annulla. I filtri ideali sono rect in frequenza; non distorsione secondo Heaviside: |G| costante e fase lineare.Risposta in frequenza e filtri →, Trasformata di Fourier a tempo discretoLa trasformata di Fourier di un segnale discreto è $S(f)=\sum_nT,s(nT),e^{-i2\pi fnT}$: una funzione continua e periodica in $f$ di periodo $F_p=1/T$. L'antitrasformata è l'integrale su un periodo, $s(nT)=\int_0^{F_p}S(f)e^{i2\pi fnT}df$. Le regole sono quelle del caso continuo (traslazione, convoluzione $\leftrightarrow$ prodotto, Parseval $\sum T|s|^2=\int_0^{F_p}|S|^2$), con incremento e somma corrente al posto di derivata e integrale. Per segnali reali $S(f)=S^*(-f)$, quindi basta $[0,F_p/2]$. Esempi: $\delta\to1$, $1\to\delta_{F_p}$, rect $\to$ sinc periodico, $a^n\mathbf 1_0\to\frac T{1-ae^{-i2\pi fT}}$.Trasformata di Fourier a tempo discreto →, Vettori gaussianiX = (X₁, ..., Xₙ) è un vettore gaussiano N(m, Σ) se ogni combinazione lineare a·X è gaussiana (equivalentemente X = m + AZ con Z gaussiane standard indipendenti); se Σ è invertibile ha densità exp(−½(x−m)ᵀΣ⁻¹(x−m)) / √((2π)ⁿ det Σ). Proprietà chiave: AX + b ~ N(Am + b, AΣAᵀ), le marginali sono gaussiane e componenti non correlate sono indipendenti.Vettori gaussiani →.

(1) Teoria: correlazione e PSD con la convenzione del corso

Ingresso. Campioni indipendenti a media nulla e potenza $\sigma^2$ hanno $$r_x(kT)=E[x(nT+kT)\,x(nT)]=\begin{cases}\sigma^2 & k=0\\0&k\ne0.\end{cases}$$ La densità spettrale di potenza di un processo su $\mathbb Z(T)$ è la trasformata di Fourier discreta (con il peso $T$) della correlazione: $$R_x(f)=\sum_kT\,r_x(kT)\,e^{-i2\pi fkT}=\sigma^2T=\frac{\sigma^2}{F_s}=10^{-3}\ \text{(per }\sigma=1,\ F_s=1000\text{)},$$ costante (da qui "bianco") e periodica di periodo $F_s$; sull'intervallo $[-F_s/2,F_s/2]$ ha integrale $\sigma^2T\cdot F_s=\sigma^2$, la potenza statistica (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 →). Il valore $\sigma^2/F_s$ è proprio quello che restituisce scipy.signal.welch con scaling='density' (densità in V$^2$/Hz).

Uscita. Per un filtro LTI (Processi aleatori attraverso sistemi LTIUna trasformazione applicata a ogni realizzazione dà un nuovo processo, ma dalla descrizione statistica di x non sempre si ricava quella di y (controesempi); si può se ogni vettore di y dipende da un vettore finito di x. Per una tf lineare m_y e r_y si calcolano col nucleo; un filtro LTI conserva la stazionarietà in senso lato: m_y = G(0) m_x, R_y = |G|² R_x. Il campionamento di un processo sl dà un processo sl con PSD ripetuta periodicamente. L'interpolazione LTI di un processo discreto sl è ciclostazionaria. Teorema del campionamento per processi: R_x nulla fuori da (-B,B) e Fc ≥ 2B danno ricostruzione esatta in media quadratica; altrimenti la potenza dell'errore è in banda + fuori banda.Processi aleatori attraverso sistemi LTI →) l'uscita di un processo stazionario è stazionaria con $$m_y=H(0)\,m_x=0,\qquad R_y(f)=\lvert H(f)\rvert^2R_x(f)=\lvert H(f)\rvert^2\,\frac{\sigma^2}{F_s}.$$ Il FIR $y[n]=\sum_kh_kx[n-k]$ ha, nel corso, risposta impulsiva $g(kT)=h_k/T$ e quindi $H(f)=\sum_kT\,g(kT)\,e^{-i2\pi fkT}=\sum_kh_ke^{-i2\pi fkT}$, con $H(0)=\sum_kh_k=1$ per un passa-basso a guadagno unitario. La potenza dell'uscita è l'area della PSD (Parseval sul periodo): $$M_y=\int_{-F_s/2}^{F_s/2}R_y\,df=\frac{\sigma^2}{F_s}\int\lvert H\rvert^2df=\sigma^2\sum_kh_k^2.$$ Per il filtro ideale ($\lvert H\rvert=1$ per $\lvert f\rvert<f_c$, zero altrove) l'area è un rettangolo di altezza $\sigma^2/F_s$ e larghezza $2f_c$: $$M_y=\sigma^2\frac{2f_c}{F_s}=\sigma^2\frac{f_c}{F_s/2}=1\cdot\frac{80}{500}=0{,}16.$$ È la frazione $\frac{f_c}{F_s/2}$ della banda di Nyquist lasciata passare: il rumore bianco ha la stessa potenza per Hz a tutte le frequenze e il filtro ne tiene solo il $16\%$. Anche la correlazione dell'uscita si ottiene da $R_y$: per l'ideale $$r_y(\tau)=\sigma^2\frac{2f_c}{F_s}\,\mathrm{sinc}(2f_c\tau)=0{,}16\ \mathrm{sinc}(160\,\tau)\quad(\tau\text{ in s}),$$ nulla ai lag multipli di $\frac1{2f_c}=6{,}25$ ms; per il FIR reale $r_y(kT)=\sigma^2\sum_mh_mh_{m+k}$. I campioni dell'uscita non sono più indipendenti: sono correlati su circa $\frac1{2f_c}\approx6$ ms, cioè su $\approx6$ campioni a $1$ ms.

(2) Il codice, passo per passo

• Generazione. rng.standard_normal(N) dà campioni indipendenti $N(0,1)$, moltiplicati per $\sigma$ (Matlab randn). Campioni indipendenti implica incorrelati (il viceversa vale perché sono gaussiani).
• Filtro. firwin(L+1, fc, fs=Fs) progetta con il metodo della finestra la risposta all'impulso di un passa-basso ideale (un seno cardinale troncato a $129$ coefficienti) con finestra di Hamming, che è anche la finestra predefinita di fir1 di Matlab (il commento del file dice Hann, ma fir1 usa Hamming se non specificato). Il guadagno in continua è $1$ e, per costruzione, $\lvert H(f_c)\rvert\approx0{,}5$ ($-6$ dB nel taglio, non $-3$ dB come per Butterworth).
• Transitorio. lfilter parte con stato nullo; si scartano i primi $L$ campioni, x[L:] e y[L:], per lavorare a regime.
• Correlazione. Stima biased: $\hat r(k)=\frac1N\sum_nz[n]z[n+k]$ (Matlab xcorr(z, tau_max, 'biased')); è la stima della media d'insieme $E[z(nT+kT)z(nT)]$ fatta con la media temporale, lecita perché il rumore bianco gaussiano stazionario è ergodico.
• PSD. welch divide il segnale in segmenti di $2048$ campioni sovrapposti al $50\%$, moltiplica ciascuno per una finestra di Hann, ne calcola il periodogramma e fa la media sui $126$ segmenti. return_onesided=False dà la PSD bilatera su $[0,F_s)$ (frequenze negative nella seconda metà): fftshift la riordina in $[-F_s/2,F_s/2)$. Risoluzione in frequenza $F_s/2048=0{,}49$ Hz. La media su $K$ periodogrammi riduce lo scarto relativo di una stima di PSD a circa $1/\sqrt K$: per $126$ segmenti $\approx9\%$ (un singolo periodogramma ha scarto del $100\%$, non converge al crescere di $N$).
• Confronto. freqz(h, 1, worN=f, fs=Fs) calcola $H(f)$ esatta sulla stessa griglia, e np.correlate(h, h, 'full') dà $\sum_mh_mh_{m+k}$ per la correlazione teorica dell'uscita.

import numpy as np
from scipy import signal, stats
import matplotlib.pyplot as plt

rng = np.random.default_rng(8)
Fs, N, sigma, fc, L = 1000, 2**17, 1.0, 80, 128
x = sigma * rng.standard_normal(N)                     # rumore bianco gaussiano: campioni indipendenti N(0, sigma^2)

h = signal.firwin(L + 1, fc, fs=Fs)                    # MATLAB: fir1(L, fc/(Fs/2)); L+1 coefficienti, finestra di Hamming
y = signal.lfilter(h, 1, x)                            # uscita: y[n] = sum_k h[k] x[n-k]
x, y = x[L:], y[L:]                                    # scarto il transitorio iniziale (L campioni)

def corr(z, tau_max=150):                              # correlazione empirica 'biased' (divisa per N), lag 0..tau_max
    return np.array([z[: len(z) - k] @ z[k:] / len(z) for k in range(tau_max + 1)])
rx, ry = corr(x), corr(y)

# PSD bilatera con Welch; fftshift per mettere f = 0 al centro
nper = 2048
f, Rx = signal.welch(x, Fs, window="hann", nperseg=nper, noverlap=nper // 2, nfft=nper, return_onesided=False)
_, Ry = signal.welch(y, Fs, window="hann", nperseg=nper, noverlap=nper // 2, nfft=nper, return_onesided=False)
f, Rx, Ry = np.fft.fftshift(f), np.fft.fftshift(Rx), np.fft.fftshift(Ry)
df = f[1] - f[0]

# ---- confronto con la teoria ----
_, Hf = signal.freqz(h, 1, worN=f, fs=Fs)              # H(f) = sum_k h[k] e^{-i2 pi f k T}  (= T * fft(h/T), convenzione del corso)
Ry_teo = np.abs(Hf)**2 * sigma**2 / Fs                 # R_y = |H|^2 R_x, con R_x = sigma^2 T = sigma^2/Fs
print(f"H(0) = sum h = {h.sum():.4f}; |H(fc)| = {abs(signal.freqz(h, 1, worN=[fc], fs=Fs)[1][0]):.4f}")
print(f"r_x(0) = {rx[0]:.4f}; max |r_x(k)|, k >= 1: {np.abs(rx[1:]).max():.4f}")
print(f"R_x medio (Welch) = {Rx.mean():.3e}; teoria sigma^2/Fs = {sigma**2 / Fs:.3e}; scarto relativo tipico {Rx.std() / Rx.mean():.3f}")
print(f"R_y per |f|<50: {Ry[np.abs(f) < 50].mean():.3e}; in f=+-fc: {Ry[np.abs(np.abs(f) - fc).argmin()]:.3e} "
      f"(teoria {Ry_teo[np.abs(np.abs(f) - fc).argmin()]:.3e}); per |f|>110: {Ry[np.abs(f) > 110].mean():.1e}")
print(f"potenza = integrale della PSD: x {Rx.sum() * df:.4f} (varianza {x.var():.4f}); y {Ry.sum() * df:.4f} (varianza {y.var():.4f})")
print(f"potenza di y: r_y(0) = {ry[0]:.4f}; sigma^2 sum(h^2) = {sigma**2 * np.sum(h**2):.4f}; "
      f"filtro ideale sigma^2 fc/(Fs/2) = {sigma**2 * fc / (Fs / 2):.4f}")
ry_teo = sigma**2 * np.r_[np.correlate(h, h, "full")[L:], np.zeros(150 - L)]   # r_y(kT) = sigma^2 sum_m h[m] h[m+k]
ry_id = sigma**2 * 2 * fc / Fs * np.sinc(2 * fc * np.arange(151) / Fs)            # filtro ideale: 0.16 sinc(2 fc tau)
print("r_y(0..3): empirica", np.round(ry[:4], 4), "| FIR teorica", np.round(ry_teo[:4], 4), "| ideale", np.round(ry_id[:4], 4))
print(f"max |r_y empirica - FIR teorica| sui lag 0..150: {np.abs(ry - ry_teo).max():.4f}")
print(f"curtosi in eccesso di y (0 per una gaussiana): {stats.kurtosis(y):.3f}")
plt.show()

(3) Risultati e cosa si vede

• Nel tempo. $x$ è una nuvola irregolare di valori in $[-4,4]$ (dev. standard $1$), senza struttura visibile; $y$ è molto più liscia (le variazioni veloci sono state eliminate), con escursioni minori (deviazione standard $\sqrt{0{,}155}=0{,}39$).
• Correlazione dell'ingresso. $\hat r_x(0)=1{,}0037$ (teoria $\sigma^2=1$) e $\lvert\hat r_x(k)\rvert\le0{,}0076$ per $k=1,\dots,150$: un impulso in $0$ su un fondo che fluttua di $\approx\sigma^2/\sqrt N=0{,}003$.
• Correlazione dell'uscita. Non è più un impulso: è un seno cardinale di larghezza dell'ordine di $\frac1{2f_c}$, con $\hat r_y(k)$ per $k=0,\dots,3$ uguale a $0{,}1555;\ 0{,}1496;\ 0{,}1326;\ 0{,}1069$ (teoria per il FIR: $0{,}1542;\ 0{,}1482;\ 0{,}1310;\ 0{,}1052$; per l'ideale $0{,}16\ \mathrm{sinc}(0{,}16\,k)$: $0{,}1600;\ 0{,}1533;\ 0{,}1344;\ 0{,}1059$). L'errore massimo fra stima e teoria del FIR sui lag $0,\dots,150$ è $0{,}0028$.
• PSD dell'ingresso. Piatta, media $1{,}002\times10^{-3}$ (teoria $\sigma^2/F_s=1{,}000\times10^{-3}$), con una fluttuazione relativa del $9{,}1\%$ da punto a punto (il valore atteso per $126$ medie è $\approx9\%$). L'area $\sum R_x\,df=1{,}0025$ coincide con la varianza campionaria $1{,}0037$.
• PSD dell'uscita. A forma di "finestra": piatta a $1{,}015\times10^{-3}$ per $\lvert f\rvert<50$ Hz, scende a $2{,}4\times10^{-4}$ in $f=\pm f_c$ (teoria $\lvert H(f_c)\rvert^2\sigma^2/F_s=2{,}5\times10^{-4}$: $-6$ dB), e vale $\approx2\times10^{-10}$ oltre $110$ Hz, cioè $-67$ dB rispetto alla banda passante. Il rapporto $R_y/R_x$ riproduce $\lvert H\rvert^2$.
• Potenza dell'uscita. Dall'area della PSD: $0{,}1551$; dalla correlazione, $\hat r_y(0)=0{,}1555$ (uguale alla varianza campionaria $0{,}1555$); teoria per il FIR reale $\sigma^2\sum h_k^2=0{,}1542$; teoria per il passa-basso ideale $\sigma^2\frac{f_c}{F_s/2}=0{,}16$. Il valore ideale è leggermente più alto perché il FIR ha una banda di transizione e lascia passare un po' meno potenza nei dintorni di $f_c$.
• Gaussianità. $y$ è combinazione lineare di variabili gaussiane indipendenti, quindi gaussiana (Vettori gaussianiX = (X₁, ..., Xₙ) è un vettore gaussiano N(m, Σ) se ogni combinazione lineare a·X è gaussiana (equivalentemente X = m + AZ con Z gaussiane standard indipendenti); se Σ è invertibile ha densità exp(−½(x−m)ᵀΣ⁻¹(x−m)) / √((2π)ⁿ det Σ). Proprietà chiave: AX + b ~ N(Am + b, AΣAᵀ), le marginali sono gaussiane e componenti non correlate sono indipendenti.Vettori gaussiani →) con varianza $\sigma^2\sum h_k^2$: la curtosi in eccesso misurata vale $-0{,}008\approx0$.

Errori tipici.

• Scrivere la PSD del rumore bianco discreto come $\sigma^2$ invece di $\sigma^2T=\sigma^2/F_s$: la PSD è per unità di frequenza e il suo integrale su un periodo $F_s$ deve dare la potenza $\sigma^2$.
• Moltiplicare per $\lvert H\rvert$ invece che per $\lvert H\rvert^2$ nel passaggio da $R_x$ a $R_y$.
• Dimenticare che la PSD ottenuta da un solo periodogramma non è consistente (serve mediare, come fa Welch) e leggere come "righe" le fluttuazioni casuali.

Vedi anche: Esercizio - processi aleatori, media d'insieme, correlazione e potenza temporale 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 →.