Densità spettrale di potenza: lo spettro di un processo stocastico e come stimarlo

Nell'articolo sui processi stocastici abbiamo descritto un processo stazionario nel tempo, con la media e la funzione di autocovarianza. Per i segnali deterministici c'era anche una descrizione in frequenza, lo spettro. Qui facciamo lo stesso per i processi stocastici.

Dall'autocovarianza allo spettro

La DTFT non si può applicare direttamente a un processo stocastico. Una realizzazione dura per sempre e non si spegne: ha energia infinita. Ha però potenza finita, perché la sua varianza resta costante nel tempo.

L'idea è trasformare la funzione di autocovarianza $\gamma_{vv}(\tau)$. Per i processi che ci interessano decade a zero abbastanza in fretta al crescere del lag, quindi è un segnale a energia finita, e la sua DTFT esiste. Si chiama densità spettrale di potenza:

$$\Gamma_{vv}(\omega) = \sum_{\tau=-\infty}^{+\infty} \gamma_{vv}(\tau)\, e^{-j\omega\tau}$$

La corrispondente trasformata Zeta si indica con $\Phi_{vv}(z)$, e vale $\Gamma_{vv}(\omega) = \Phi_{vv}(e^{j\omega})$.

$\Gamma_{vv}(\omega)$ dice come la potenza del processo si distribuisce tra le frequenze: quanto contribuisce, in media sulle realizzazioni, ciascuna componente in frequenza alla variabilità del segnale.

Conoscere lo spettro o l'autocovarianza è equivalente, perché si torna indietro con l'antitrasformata:

$$\gamma_{vv}(\tau) = \frac{1}{2\pi} \int_{-\pi}^{\pi} \Gamma_{vv}(\omega)\, e^{j\omega\tau}\, d\omega$$

Con $\tau = 0$ si ottiene un risultato importante: la varianza del processo è l'area sotto la densità spettrale, divisa per $2\pi$.

$$\operatorname{Var}[v(t)] = \gamma_{vv}(0) = \frac{1}{2\pi} \int_{-\pi}^{\pi} \Gamma_{vv}(\omega)\, d\omega$$

La densità spettrale di potenza ha quattro proprietà:

1. è reale;
2. è non negativa: $\Gamma_{vv}(\omega) \ge 0$, come deve essere una potenza;
3. è pari: $\Gamma_{vv}(-\omega) = \Gamma_{vv}(\omega)$;
4. è periodica di periodo $2\pi$.

Dimostrazioni: le quattro proprietà

Perché è reale. Per la formula di Eulero, ogni esponenziale della somma è un coseno più un seno moltiplicato per $j$:

$$e^{-j\omega\tau} = \cos(\omega\tau) - j\,\sin(\omega\tau)$$

Il seno è la parte immaginaria. Si prendono insieme i termini in $\tau$ e in $-\tau$. Hanno lo stesso coefficiente, perché $\gamma_{vv}(-\tau) = \gamma_{vv}(\tau)$. Il coseno è pari e resta uguale, mentre il seno è dispari e cambia segno: $\sin(-\omega\tau) = -\sin(\omega\tau)$. Sommando la coppia, i due seni si cancellano:

$$\gamma_{vv}(\tau)\, e^{-j\omega\tau} + \gamma_{vv}(\tau)\, e^{j\omega\tau} = 2\,\gamma_{vv}(\tau) \cos(\omega\tau)$$

Quindi la densità spettrale è una somma di soli coseni, ed è reale:

$$\Gamma_{vv}(\omega) = \gamma_{vv}(0) + 2\sum_{\tau=1}^{\infty} \gamma_{vv}(\tau) \cos(\omega\tau)$$

Perché è pari. Dalla stessa formula, perché $\cos(-\omega\tau) = \cos(\omega\tau)$.

Perché è periodica. Ogni termine lo è: $e^{-j(\omega + 2\pi)\tau} = e^{-j\omega\tau}\, e^{-j2\pi\tau} = e^{-j\omega\tau}$, perché $\tau$ è intero.

Perché è non negativa. Non si vede dalla formula. Si ricava dal periodogramma, più avanti: $I_N(\omega) = |V(e^{j\omega})|^2 / N$ è sempre $\ge 0$, e il suo valore atteso tende a $\Gamma_{vv}(\omega)$. Il limite di quantità non negative non può essere negativo.

Per le ultime due proprietà basta guardarla tra $0$ e $\pi$. Come per la DTFT, $\omega$ è una pulsazione normalizzata, in radianti per campione: $\omega = \pi$ corrisponde alla frequenza di Nyquist, $f_s/2$.

Densità spettrale di un rumore bianco

Per il rumore bianco $e(t) \sim \text{WN}(0, \lambda^2)$ l'autocovarianza è diversa da zero solo in $\tau = 0$, dove vale $\lambda^2$. Nella somma resta un solo termine:

$$\Gamma_{ee}(\omega) = \lambda^2$$

Lo spettro è costante: tutte le frequenze contribuiscono allo stesso modo, e nessuna domina. Il nome viene da qui: come la luce bianca contiene tutti i colori in parti uguali. È lo spettro di un segnale che non ha nessuna regolarità, e che quindi non si può predire.

Processi stazionari a confronto

Il grafico mette a confronto quattro processi stazionari a media nulla, tutti ottenuti dallo stesso rumore bianco $e(t) \sim \text{WN}(0, 1)$:

• rumore bianco: $e(t)$ stesso;
• regolare: $y(t) = a\,y(t-1) + e(t)$ con $a = 0{,}9$;
• alternante: lo stesso con $a = -0{,}9$;
• risonante: $y(t) = 1{,}6\,y(t-1) - 0{,}9\,y(t-2) + e(t)$.

Per ciascuno: una realizzazione, l'autocovarianza e la densità spettrale. Ogni pannello ha la sua scala.

Realizzazione, autocovarianza e spettro raccontano la stessa cosa in tre modi.

• Regolare. Ogni campione è il 90% del precedente: il segnale cambia lentamente. L'autocovarianza decade piano, perché campioni lontani restano simili. La potenza è concentrata alle basse frequenze: $\Gamma(0) = 100$, $\Gamma(\pi) = 0{,}28$.
• Alternante. Il segno si inverte quasi a ogni passo. L'autocovarianza cambia segno a ogni lag. Lo spettro è lo specchio del precedente: dominano le alte frequenze, vicino a $\pi$.
• Risonante. Il segnale oscilla in modo irregolare, con un periodo di circa 11 campioni. Anche l'autocovarianza oscilla, smorzandosi. Lo spettro ha un picco di risonanza in $\omega \approx 2\pi/11$.

Perché ci sono i picchi. Lo spettro dice a che ritmo si muove il segnale. Se cambia piano, la potenza sta alle frequenze basse, vicino a $0$. Se salta su e giù a ogni passo, sta a quelle alte, vicino a $\pi$. Se oscilla con un suo ritmo, c'è un picco proprio a quel ritmo.

Più lo spettro è concentrato, più il processo è regolare e prevedibile. Lo spettro piatto del rumore bianco è il caso opposto.

Densità cross-spettrale. Con due processi stazionari $v$ e $x$ si trasforma allo stesso modo la cross-covarianza: $\Gamma_{vx}(\omega)$ è la DTFT di $\gamma_{vx}(\tau)$, e $\Phi_{vx}(z)$ la sua trasformata Zeta. A differenza dell'autocovarianza, la cross-covarianza non è pari, ma vale $\gamma_{vx}(\tau) = \gamma_{xv}(-\tau)$; di conseguenza $\Gamma_{vx}(\omega) = \Gamma_{xv}(-\omega)$ e $\Phi_{vx}(z) = \Phi_{xv}(z^{-1})$. Servirà per collegare ingresso e uscita di un sistema.

Stimare dai dati

Nella pratica né l'autocovarianza né lo spettro sono noti: vanno stimati. Le ipotesi sono:

• il processo è stazionario ed ergodico, così basta una sola realizzazione;
• ha media nulla. Se non è così, si stima la media con la media temporale e la si sottrae;
• sono disponibili $N$ campioni, $v(0), \dots, v(N-1)$.

Autocovarianza campionaria

Siccome $\gamma_{vv}(\tau) = \mathbb{E}[v(t)\,v(t+\tau)]$, la stima viene naturale: la media dei prodotti delle coppie di campioni a distanza $\tau$. Le coppie disponibili sono $N - |\tau|$:

$$\hat\gamma_{vv}(\tau) = \frac{1}{N - |\tau|} \sum_{t=0}^{N-|\tau|-1} v(t)\, v(t+|\tau|)$$

Si usa $|\tau|$ perché, per la simmetria, la stima è la stessa per $\tau$ e $-\tau$. Con $\tau = 0$ si stima la varianza. Questo stimatore è corretto, e consistente per ogni $\tau$ fissato.

Il problema sono i lag grandi. Con $\tau$ vicino a $N$ restano pochissime coppie da mediare, e la stima diventa molto incerta. L'alternativa è dividere sempre per $N$:

$$\hat\gamma'_{vv}(\tau) = \frac{1}{N} \sum_{t=0}^{N-|\tau|-1} v(t)\, v(t+|\tau|) = \frac{N - |\tau|}{N}\, \hat\gamma_{vv}(\tau)$$

Questo è distorto: in media vale $\frac{N-|\tau|}{N}\gamma_{vv}(\tau)$, e viene schiacciato verso zero tanto più quanto $\tau$ si avvicina a $N$. Ma è asintoticamente corretto, e ai lag grandi ha varianza molto più piccola.

Un esempio con il processo $v(t) = 0{,}8\,v(t-1) + e(t)$ dell'articolo precedente, la cui autocovarianza vera vale $2{,}78 \cdot 0{,}8^{|\tau|}$. Ripetendo la stima su 3.000 realizzazioni da $N = 200$ campioni, ecco la deviazione standard dei due stimatori:

A lag piccoli i due stimatori si equivalgono. A lag grandi lo stimatore corretto dà valori a caso: a lag 190 ha una deviazione standard di 1,7, dove il valore vero è zero, mentre quello distorto resta vicino a zero. È lo stesso compromesso tra bias e varianza visto per i modelli statici: un po' di distorsione in cambio di molta stabilità. E siccome per i processi che ci interessano $\gamma_{vv}(\tau) \to 0$, schiacciare verso zero è anche la scelta giusta.

Periodogramma

Per stimare lo spettro si prende la definizione e si mette, al posto dell'autocovarianza vera, quella stimata $\hat\gamma'_{vv}$. La somma si ferma a $|\tau| < N$, perché oltre non ci sono dati. Il risultato si chiama periodogramma:

$$I_N(\omega) = \sum_{\tau=-(N-1)}^{N-1} \hat\gamma'_{vv}(\tau)\, e^{-j\omega\tau}$$

Si dimostra che è il modulo al quadrato della DTFT dei dati, diviso per $N$:

$$I_N(\omega) = \frac{1}{N}\, \bigl|V(e^{j\omega})\bigr|^2$$

dove $V(e^{j\omega})$ è la DTFT dei campioni $v(0), \dots, v(N-1)$. In pratica lo si calcola con una FFT, sui bin della DFT $\omega_k = 2\pi k / N$.

In media il periodogramma è giusto: per $N \to \infty$ il suo valore atteso tende a $\Gamma_{vv}(\omega)$. Il problema è quanto oscilla. A ogni frequenza sbaglia, tipicamente, di quanto vale lo spettro stesso, e questo errore non scende con più dati: con più campioni si ottengono più frequenze, non stime migliori. Si dice che lo stimatore non è consistente.

Il grafico lo mostra per il processo risonante. In alto il periodogramma su 4.096 campioni (la scala verticale è logaritmica): segue lo spettro vero, ma è tutto frastagliato.

Per misurare l'errore si divide, frequenza per frequenza, il periodogramma per lo spettro vero, e si calcola la deviazione standard di questi rapporti. Con una stima perfetta varrebbe zero:

Con sedici volte più dati l'errore resta lo stesso: circa il 100% del valore vero.

Metodo di Bartlett

Il periodogramma però in media è giusto: sbaglia a caso, a volte per eccesso e a volte per difetto. Facendo la media di più periodogrammi indipendenti, gli errori si compensano.

Il metodo di Bartlett divide gli $N$ dati in $K$ blocchi da $M = N/K$ campioni, calcola il periodogramma di ogni blocco e ne fa la media:

$$\bar I_{M,K}(\omega) = \frac{1}{K} \sum_{i=1}^{K} I^{(i)}_{M}(\omega)$$

Se l'autocovarianza decade in fretta, i blocchi sono quasi indipendenti e la varianza si divide per $K$. È la curva in basso nel grafico: 16 blocchi da 256 campioni, molto più regolare.

Il prezzo è che ogni blocco è corto. Con meno campioni le frequenze stimate sono più distanti tra loro ($2\pi/M$), e i picchi vengono abbassati e allargati. Con 16.384 campioni del processo risonante, il cui picco vero vale 346:

• Pochi blocchi lunghi: il picco è risolto bene, ma la stima è ancora rumorosa (400, il 15% in più del vero).
• Tanti blocchi corti: la stima è liscia, ma il picco viene spianato (215).

Se lo spettro ha picchi stretti servono blocchi lunghi, proprio come serve una misura lunga per distinguere frequenze vicine con la DFT.

Riepilogo

• La densità spettrale di potenza $\Gamma_{vv}(\omega)$ è la DTFT dell'autocovarianza: dice come la potenza del processo si distribuisce tra le frequenze. È reale, non negativa, pari e periodica, quindi basta guardarla tra 0 e $\pi$.
• L'area sotto $\Gamma_{vv}(\omega)$, divisa per $2\pi$, è la varianza del processo.
• Il rumore bianco ha spettro piatto, $\Gamma_{ee}(\omega) = \lambda^2$. Uno spettro concentrato alle basse frequenze indica un processo regolare, alle alte un processo alternante, un picco un'oscillazione dominante.
• Dei due stimatori dell'autocovarianza, quello che divide per $N$ è distorto ma ha varianza molto minore ai lag grandi.
• Il periodogramma $|V(e^{j\omega})|^2/N$ non è consistente: la sua varianza non scende con più dati. Il metodo di Bartlett fa la media di $K$ periodogrammi e divide la varianza per $K$, al prezzo di risoluzione in frequenza.