Processi stocastici e sistemi LTI: ogni processo stazionario è rumore bianco filtrato

Negli articoli precedenti abbiamo descritto un processo stocastico stazionario nel tempo e in frequenza. Ora lo facciamo passare attraverso un sistema LTI. Il risultato è il collegamento che serve per tutta l'identificazione: un processo stazionario si può vedere come rumore bianco filtrato da un sistema dinamico.

Problema

Mettiamo un processo stocastico $u(t)$ in ingresso a un sistema LTI:

l'ingresso è stazionario in senso debole, con media $m_u$ e autocovarianza $\gamma_{uu}(\tau)$ ;
il sistema è asintoticamente stabile, con funzione di trasferimento razionale $G(z)$ , risposta all'impulso $g(t)$ e guadagno statico $\mu = G(1)$ .

L'uscita si ottiene come sempre con la convoluzione:

y(t) = \sum_{i=0}^{\infty} g(i)\, u(t-i)

Se l'ingresso è casuale, anche l'uscita lo è. La domanda è: che processo esce? Ha media costante? È ancora stazionario? Che autocovarianza e che spettro ha?

Come esempio useremo un sistema con un solo polo in $0{,}8$ :

G(z) = \frac{1}{1 - 0{,}8\,z^{-1}}, \qquad y(t) = 0{,}8\,y(t-1) + u(t)

Ogni uscita è l'80% della precedente più l'ingresso nuovo. Con un impulso in ingresso l'uscita vale $1, 0{,}8, 0{,}64, \dots$ , quindi la risposta all'impulso è $g(i) = 0{,}8^i$ . Il guadagno statico si ottiene ponendo $z = 1$ : $\mu = G(1) = 1/(1 - 0{,}8) = 5$ .

Valore atteso

Si parte dalla convoluzione e si prende il valore atteso. Il valore atteso è lineare, quindi passa dentro la somma. Ogni termine diventa la media dell'ingresso, che vale sempre $m_u$ perché l'ingresso è stazionario:

\mathbb{E}[y(t)] = \sum_{i=0}^{\infty} g(i)\, \mathbb{E}[u(t-i)] = m_u \sum_{i=0}^{\infty} g(i)

Resta la somma della risposta all'impulso, che è proprio il guadagno statico. Infatti $G(z) = \sum_i g(i)\, z^{-i}$ , e in $z = 1$ diventa $G(1) = \sum_i g(i) = \mu$ . Quindi:

\mathbb{E}[y(t)] = \mu\, m_u

La media dell'uscita è la media dell'ingresso moltiplicata per il guadagno statico, e non dipende da $t$ .

Nell'esempio, con un rumore bianco di media $0{,}5$ in ingresso, l'uscita ha media $5 \cdot 0{,}5 = 2{,}5$ .

Autocovarianza e Stazionarietà

La media dell'ingresso non cambia le covarianze, quindi si può supporre $m_u = 0$ . Si moltiplica l'uscita all'istante $t + \tau$ , scritta come convoluzione, per l'ingresso all'istante $t$ , e si prende il valore atteso. Il valore atteso passa dentro la somma, e ogni termine è un'autocovarianza dell'ingresso:

\gamma_{uy}(\tau) = \mathbb{E}\bigl[u(t)\, y(t+\tau)\bigr] = \sum_{i=0}^{\infty} g(i)\, \mathbb{E}\bigl[u(t)\, u(t+\tau-i)\bigr] = \sum_{i=0}^{\infty} g(i)\, \gamma_{uu}(\tau - i)

La cross-covarianza tra ingresso e uscita dipende solo dal lag $\tau$ , perché $\gamma_{uu}$ dipende solo dal lag. Con lo stesso conto, moltiplicando per l'uscita invece che per l'ingresso:

\gamma_{yy}(\tau) = \mathbb{E}\bigl[y(t)\, y(t+\tau)\bigr] = \sum_{i=0}^{\infty} g(i)\, \mathbb{E}\bigl[y(t)\, u(t+\tau-i)\bigr] = \sum_{i=0}^{\infty} g(i)\, \gamma_{yu}(\tau - i)

Anche l'autocovarianza dell'uscita dipende solo da $\tau$ , perché $\gamma_{yu}(\tau) = \gamma_{uy}(-\tau)$ dipende solo dal lag. Media costante e autocovarianza che dipende solo dal lag: l'uscita è stazionaria. È un teorema:

Se $u(t)$ è un processo stazionario e $G(z)$ è asintoticamente stabile, anche $y(t)$ è un processo stazionario.

C'è una precisazione pratica. Il conto suppone che l'ingresso sia applicato da sempre, da $t = -\infty$ . Nella realtà parte da un istante preciso, e l'uscita diventa stazionaria solo dopo un transitorio, quando l'effetto delle condizioni iniziali si è spento. La stabilità garantisce che si spenga.

Spettro dell'uscita

Passando in frequenza, le convoluzioni diventano prodotti, e si ottiene il risultato centrale:

\Gamma_{yy}(\omega) = \bigl|G(e^{j\omega})\bigr|^2\, \Gamma_{uu}(\omega)

e in termini di trasformata Zeta

\Phi_{yy}(z) = G(z)\, G(z^{-1})\, \Phi_{uu}(z)

Dimostrazione: lo spettro dell'uscita

Le due relazioni nel tempo sono convoluzioni con $g$ , quindi in frequenza diventano prodotti con $G(e^{j\omega})$ :

\Gamma_{uy}(\omega) = G(e^{j\omega})\, \Gamma_{uu}(\omega), \qquad \Gamma_{yy}(\omega) = G(e^{j\omega})\, \Gamma_{yu}(\omega)

Per la proprietà della densità cross-spettrale, $\Gamma_{yu}(\omega) = \Gamma_{uy}(-\omega)$ . Sostituendo:

\begin{aligned} \Gamma_{yy}(\omega) &= G(e^{j\omega})\, \Gamma_{uy}(-\omega) \\ &= G(e^{j\omega})\, G(e^{-j\omega})\, \Gamma_{uu}(-\omega) \\ &= \bigl|G(e^{j\omega})\bigr|^2\, \Gamma_{uu}(\omega) \end{aligned}

Nell'ultimo passaggio, $G(e^{-j\omega})$ è il complesso coniugato di $G(e^{j\omega})$ , perché i coefficienti di $G$ sono reali, e $\Gamma_{uu}$ è pari. Per la trasformata Zeta il conto è lo stesso, con $z^{-1}$ al posto di $-\omega$ .

Il sistema modula lo spettro dell'ingresso: le frequenze dove $|G(e^{j\omega})|$ è grande vengono amplificate, le altre attenuate. È la stessa risposta in frequenza dei segnali deterministici, ma al quadrato, perché lo spettro è una potenza.

Nell'esempio, con un rumore bianco $\text{WN}(0, 1)$ in ingresso, $\Gamma_{uu}(\omega) = 1$ e

\Gamma_{yy}(\omega) = \frac{1}{\bigl|1 - 0{,}8\,e^{-j\omega}\bigr|^2} = \frac{1}{1{,}64 - 1{,}6\cos\omega}

Il denominatore si calcola con Eulero: $1 - 0{,}8\,e^{-j\omega} = (1 - 0{,}8\cos\omega) + j\,0{,}8\sin\omega$ , e il modulo quadro è $(1 - 0{,}8\cos\omega)^2 + 0{,}64\sin^2\omega = 1 - 1{,}6\cos\omega + 0{,}64\,(\cos^2\omega + \sin^2\omega) = 1{,}64 - 1{,}6\cos\omega$ .

Con $\cos\omega = 1, 0, -1$ vale $1/0{,}04 = 25$ a $\omega = 0$ , $1/1{,}64 = 0{,}61$ a $\pi/2$ e $1/3{,}24 = 0{,}31$ a $\pi$ : il polo vicino a $z = 1$ trasforma un rumore piatto in un processo dominato dalle basse frequenze. Il grafico confronta la formula con lo spettro stimato, con il metodo di Bartlett, dall'uscita di una simulazione:

La varianza dell'uscita si ottiene in due modi, che devono coincidere:

in frequenza, come area sotto lo spettro divisa per $2\pi$ . Con l'integrale noto $\frac{1}{2\pi}\int_{-\pi}^{\pi} \frac{d\omega}{a - b\cos\omega} = \frac{1}{\sqrt{a^2 - b^2}}$ si ottiene $1/\sqrt{1{,}64^2 - 1{,}6^2} = 1/0{,}36 = 2{,}78$ ;
nel tempo, perché con un rumore bianco in ingresso i campioni $u(t-i)$ sono incorrelati, e le varianze dei termini della convoluzione si sommano: $\lambda^2 \sum_i g(i)^2$ , che con $\lambda^2 = 1$ vale $\sum_i 0{,}64^i$ , una serie geometrica: $1/(1 - 0{,}64) = 2{,}78$ .

La simulazione dà 2,80.

Ogni processo stazionario è rumore bianco filtrato

Il teorema si può leggere al contrario. Se in ingresso c'è un rumore bianco $e(t) \sim \text{WN}(0, 1)$ , lo spettro dell'uscita è esattamente $|G(e^{j\omega})|^2$ . Allora un processo stazionario $y(t)$ con spettro $\Gamma_{yy}(\omega)$ si può rappresentare come l'uscita di un sistema $G(z)$ asintoticamente stabile, alimentato da rumore bianco, purché

\bigl|G(e^{j\omega})\bigr|^2 = \Gamma_{yy}(\omega)

Un $G(z)$ razionale così esiste sempre quando anche lo spettro è razionale, cioè un rapporto di polinomi in $\cos\omega$ . È il caso di tutti i processi che useremo.

Scrivendo la convoluzione, il processo è una combinazione lineare di infiniti campioni di rumore bianco, pesati dalla risposta all'impulso:

y(t) = \sum_{j=0}^{\infty} g(j)\, e(t-j) = g(0)\,e(t) + g(1)\,e(t-1) + g(2)\,e(t-2) + \dots

Questa forma si chiama $\text{MA}(\infty)$ , e sarà il punto di partenza delle famiglie di modelli del prossimo articolo.

C'è anche un uso pratico: simulare un processo. Se si conosce, o si stima, lo spettro di un fenomeno e si trova un $G(z)$ stabile e causale con quel modulo quadro, si possono generare al computer quante realizzazioni si vogliono. Basta produrre una sequenza di numeri casuali incorrelati e filtrarla con $G(z)$ .

Esempio in Python: simulare un rumore bianco filtrato

python

import numpy as np
from scipy.signal import lfilter

e = np.random.default_rng().normal(0, 1, 10_000)  # rumore bianco WN(0, 1)
y = lfilter([1], [1, -0.8], e)                    # y(t) = 0.8 y(t-1) + e(t)

Modello generale

Ora si può tornare all'impostazione dell'articolo sui processi stocastici e darle una forma precisa. L'ingresso che non si può misurare sarà proprio un rumore bianco $e(t)$ .

Serie temporali. La serie $y(t)$ è un processo stazionario, che si rappresenta come rumore bianco filtrato: $y(t) = H(z)\, e(t)$ . Il filtro del rumore si indica con $H(z)$ , per distinguerlo dal sistema fisico $G(z)$ .
Sistemi ingresso/uscita. L'uscita è la somma dell'effetto dell'ingresso noto e di un disturbo stazionario, a sua volta rumore bianco filtrato:

y(t) = G(z)\, u(t) + H(z)\, e(t)

L'obiettivo dell'identificazione sarà stimare $G(z)$ e $H(z)$ dai dati. Le due funzioni hanno però una natura diversa. $G(z)$ descrive un sistema fisico reale. $H(z)$ ed $e(t)$ invece non esistono fisicamente: sono uno strumento matematico per descrivere tutto quello che $G(z)$ non riesce a spiegare del legame tra $u(t)$ e $y(t)$ .

Riepilogo

Se un processo stazionario entra in un sistema LTI asintoticamente stabile, l'uscita è stazionaria, dopo un transitorio.
La media dell'uscita è $\mathbb{E}[y] = G(1)\, m_u$ : la media dell'ingresso per il guadagno statico.
Lo spettro dell'uscita è $\Gamma_{yy}(\omega) = |G(e^{j\omega})|^2\, \Gamma_{uu}(\omega)$ : il sistema modula lo spettro dell'ingresso con il quadrato della sua risposta in frequenza.
Con rumore bianco in ingresso, $\Gamma_{yy}(\omega) = \lambda^2\, |G(e^{j\omega})|^2$ . Viceversa, ogni processo stazionario con spettro razionale si può rappresentare come rumore bianco filtrato da un $G(z)$ stabile: $y(t) = \sum_j g(j)\, e(t-j)$ .
Il modello generale sarà $y(t) = G(z)\,u(t) + H(z)\,e(t)$ : $G$ descrive il sistema fisico, $H$ ed $e$ il disturbo.