Filtro di Wiener

Da Wikipedia, l'enciclopedia libera.

Il Filtro di Wiener è un filtro digitale per l'elaborazione numerica dei segnali proposto da Norbert Wiener negli anni 1940 e pubblicato nel 1949.

Descrizione[modifica | modifica sorgente]

L'obiettivo del filtro di Wiener è eliminare il rumore presente in un segnale. È basato su un approccio statistico.

Tipicamente i filtri sono concepiti per una ben precisa risposta in frequenza. Il filtro di Wiener approccia il problema del filtraggio da un diverso punto di vista. Si assume che si abbia conoscenza delle caratteristiche spettrali del segnale originale e del rumore, e si cerca il filtro LTI la cui uscita sia più vicina possibile al segnale originale. I filtri di Wiener sono caratterizzati dalle seguenti caratteristiche:

  1. Assunzioni: il segnale e il rumore (additivo) sono processi stocastici stazionari e lineari con caratteristiche spettrali note o autocorrelazione e mutua correlazione nota
  2. Condizioni: il filtro deve essere fisicamente realizzabile, p.e. causale (questa condizione può essere abbandonata, ottenendo così una soluzione non causale)
  3. Criteri prestazionali: minimo errore quadratico medio

Questo filtro è usato frequentemente nel processo di deconvoluzione; per questo utilizzo vedi deconvoluzione di Wiener.

Impostazione del problema[modifica | modifica sorgente]

L'ingresso di un filtro di Wiener si assume essere un segnale, s(t), alterato da un rumore additivo, n(t). L'uscita, \hat{s}(t), è calcolata per mezzo di un filtro, g(t), usando la seguente convoluzione:

\hat{s}(t) = g(t) * (s(t) + n(t))

Wienerhopf.png
dove

  • s(t) è il segnale originale (da ricavare)
  • n(t) è il rumore
  • \hat{s}(t) è il segnale stimato (che speriamo sia uguale a s(t))
  • g(t) è il filtro di Wiener

L'errore è e(t) = s(t + \alpha) - \hat{s}(t) mentre l'errore quadratico è e^2(t) = s^2(t + \alpha) - 2s(t + \alpha)\hat{s}(t) + \hat{s}^2(t) dove

  • s(t + \alpha) è l'uscita desiderata filtro
  • e(t) è l'errore

A seconda del valore di α il:

Scrivendo \hat{s}(t) come integrale di convoluzione: \hat{s}(t) = \int_{-\infty}^{\infty}{g(\tau)\left[s(t - \tau) + n(t - \tau)\right]d\tau}.

Prendendo il valore atteso dell'errore quadratico si ha

E(e^2) = R_s(0) - 2\int_{-\infty}^{\infty}{g(\tau)R_{x\,s}(\tau + \alpha)d\tau} + \int_{-\infty}^{\infty}{\int_{-\infty}^{\infty}{g(\tau)g(\theta)R_x(\tau - \theta)d\tau}d\theta}

dove

Se il signale s(t) e il rumore n(t) sono non correlati (p.e., la cross-correlazione è zero) allora si ottiene

  • R_{x\,s} = R_s
  • \,\!R_x = R_s + R_n

L'obiettivo adesso è minimizzare E(e^2) trovando la g(t) ottimale.

Soluzione stazionaria[modifica | modifica sorgente]

Il filtro di Wiener ammette soluzione per due possibili casi: il caso in cui si voglia un filtro causale, e il caso in cui un filtro non causale sia accettabile. L'ultimo è più semplice ma non si presta per applicazioni in tempo reale. L'obiettivo principale di Wiener era risolvere il caso in cui la condizione di causalità fosse valida.

Soluzione non causale[modifica | modifica sorgente]

G(s) = \frac{S_{x,s}(s)e^{\alpha s}}{S_x(s)}

Posto che g(t) sia ottimale allora l'equazione di minimo errore quadratico medio si riduce in E(e^2) = R_s(0) - \int_{-\infty}^{\infty}{g(\tau)R_{x,s}(\tau + \alpha)d\tau}

E la soluzione g(t) è l'inversa della trasformata di Laplace a due lati di G(s).

Soluzione causale[modifica | modifica sorgente]

G(s) = \frac{H(s)}{S_x^{+}(s)}

Dove

  • H(s) è la parte causale di \frac{S_{x,s}(s)e^{\alpha s}}{S_x^{-}(s)} (cioè la parte di questa frazione avente una soluzione temporale positiva per trasformazione inversa di Laplace)
  • S_x^{+}(s) è la componente causale di S_x(s) (p.e. la trasformata inversa di Laplace di S_x^{+}(s) è non nulla solo se t\ge 0)
  • S_x^{-}(s) è la componente anti causale di S_x(s) (p.e. la trasformata inversa di Laplace S_x^{-}(s) è non nulla solo per t non negativo)

Questa espressione generale è complicata e necessita di una più dettagliata spiegazione. Per ottenere la soluzione G(s) per un caso specifico, bisogna seguire questi passi (vedi LLoyd R. Welch: Wiener Hopf Theory):

1. Inizia con lo spettro S_x(s) in forma razionale e fattorizzalo nelle componenti causale e anti causale:

S_x(s) = S_x^{+}(s) S_x^{-}(s)

dove S^{+} contiene tutti gli zeri e i poli nel piano levogiro (LHP) e S^{-} contiene gli zeri e i poli nel piano destrogiro (RHP).

2. Dividi S_{x,s}(s)e^{\alpha s} per {S_x^{-}(s)} e riscrivi il risultato espandendolo in frazioni parziali.

3. Seleziona solo quei termini dell'espansione che hanno poli nel piano LHP. Chiama questi termini H(s).

4. Dividi H(s) per S_x^{+}(s). Il risultato è la funzione filtro di trasferimento desiderata G(s)

FIR filtro di Wiener per serie discrete[modifica | modifica sorgente]

Per ricavare i coefficienti del filtro di Wiener, consideriamo un signale w[n] passato a un filtro di Wiener di ordine N e con coefficienti a_i, i=0,\ldots, N. L'uscita del filtro, indicata con x[n], è data dall'espressione


x[n] = \sum_{i=0}^N a_i w[n-i]

L'errore residuo è indicato con e[n] e definito come e[n]=x[n]-s[n] (Vedi il corrispondente diagramma a blocchi). Il filtro di Wiener è concepito in maniera tale da minimizzare la deviazione standard (criterio MMSE) che può essere scritta concisamente nella maniera seguente:


a_i = arg.min ~E\{e^2[n]\}

dove E\{.\} indica l'operatore di valore atteso. In generale, i coefficienti a_i possono essere complessi e possono essere derivati anche nel caso in cui w[n] e s[n] sono complesse. Per semplicità considereremo solo il caso in cui tutte queste quantità sono reali. La deviazione standard può essere riscritta come:


\begin{array}{rcl}
E\{e^2[n]\} &=& E\{(x[n]-s[n])^2\} \\
&=& E\{x^2[n]\} + E\{s^2[n]\} - 2E\{x[n]s[n]\}\\
&=& E\{\big( \sum_{i=0}^N a_i w[n-i] \big)^2\} + E\{w^2[n]\} -2 E\{  \sum_{i=0}^N a_i w[n-i]s[n]\}
\end{array}

Per trovare il vettore [a_0,\ldots,a_N] che minimizza l'espressione sopra, calcoliamo la sua derivata rispetto ai coefficienti a_i


\begin{array}{rcl}
\frac{\partial}{\partial a_i} E\{e^2[n]\} &=& 2E\{ \big( \sum_{j=0}^N a_j w[n-j] \big) w[n-i] \} - 2E\{s[n]w[n-i]\} \quad i=0, \ldots ,N \\
&=& 2 \sum_{j=0}^N E\{w[n-j]w[n-i]\} a_j  - 2 E\{ w[n-i]s[n] \}
\end{array}

Se supponiamo che w[n] e s[n] sono stazionarie, possiamo introdurre le seguenti sequenze R_w[m] ~\textit{ and }~ R_{ws}[m] note rispettivamente come l'autocorrelazione di w[n] e la cross-correlazione tra w[n] e s[n] definite come segue


R_w[m] = E\{w[n]w[n+m]\} \quad \textit{ and } \quad R_{ws}[m] = E\{w[n]s[n+m]\}

La derivata della MSE può dunque essere riscritta come (da notare che R_{ws}[-i] = R_{sw}[i])


\frac{\partial}{\partial a_i} E\{e^2[n]\} = 2 \sum_{j=0}^{N} R_w[j-i] a_j - 2 R_{sw}[i] \quad i=0, \ldots ,N

Imponendo la derivata uguale a zero otteniamo


\sum_{j=0}^N R_w[j-i] a_j = R_{sw}[i] \quad i=0, \ldots ,N

che può essere riscritta in forma matriciale


\begin{bmatrix}
R_w[0] & R_w[1] & \cdots & R_w[N] \\
R_w[1] & R_w[0] & \cdots & R_w[N-1] \\
\vdots & \vdots & \ddots & \vdots \\
R_w[N] & R_w[N-1] & \cdots & R_w[0]
\end{bmatrix}
\begin{bmatrix}
a_0 \\ a_1 \\ \vdots \\ a_N
\end{bmatrix}
=
\begin{bmatrix}
R_{sw}[0] \\R_{sw}[1]  \\ \vdots \\ R_{sw}[N]
\end{bmatrix}

Queste equazioni sono note come equazioni di Wiener-Hopf. La matrice che appare nell'equazione è una matrice di Toeplitz simmetrica. Queste matrici sono note essere definite positive e dunque non singolari, portando ad una sola soluzione nella determinazione dei coefficienti del filtro di Wiener. Inoltre esiste un algoritmo efficiente per risolvere le equazioni di Wiener-Hopf noto come algoritmo di Levinson-Durbin.

Filtro di Wiener per il riconoscimento del Pile Up[modifica | modifica sorgente]

Il filtro di Wiener è molto utile anche nel campo del riconoscimento di segnale, in particolare del Pile Up (sovrapposizione di due segnali nella stessa finestra temporale). Vediamo come affrontare il problema e l'espressione della funzione di trasferimento del filtro.

Si supponga di avere un sistema con una risposta r(t) al segnale unitario e un segnale non corrotto u(t) in ingresso (es: \delta (t-t_0)). L'uscita s(t) del sistema sarà pertanto una convoluzione tra la funzione di risposta e il segnale di ingresso u(t):

s(t) = \int_{-\infty}^{+\infty}{ r(\tau) u(t-\tau)d\tau  }

che nel dominio di Fourier si riduce semplicemente a s(\omega) = r(\omega) u(\omega).

Nel caso reale invece, l'uscita del sistema sarà un segnale corrotto c(t) = s(t) + n(t). La u(\omega) in questo caso sarà:

u(\omega) = \frac{c(\omega)}{r(\omega)}

Lo scopo del Filtro di Wiener è quello di trovare un filtro H(\omega) tale per cui vale:

u'(\omega) = \frac{c(\omega)H(\omega)}{r(\omega)}

dove u'(\omega) è la miglior stima di u(\omega), ovvero H(\omega) minimizza la differenza quadrata tra U e U'. In formule equivale a minimizzare


L(\omega) = \int_{-\infty}^{+\infty}{ \left |u'(\omega) - u(\omega)\right |^2 d\omega  } 
          = \int_{-\infty}^{+\infty}{ \left | \frac{c(\omega)H(\omega)}{r(\omega)} - \frac{s(\omega)}{r(\omega)} \right |^2 d\omega  }

Risolvendo dL/dH = 0 si ottiene l'espressione cercata:

H(\omega) = \frac{ |s(\omega)|^2 }{ |s(\omega)|^2 + |n(\omega)|^2 }

ma se consideriamo la u(t) come una delta, l'espressione si semplifica:


H(\omega) = \frac{ r^*(\omega) }{ |r(\omega)|^2 + |n(\omega)|^2 } \qquad
u'(\omega) = \frac{ c(\omega)r^*(\omega) }{ |r(\omega)|^2 + |n(\omega)|^2 }

Questo filtro può tornare molto utile nella ricerca di segnali, come ad esempio in un algoritmo di trigger o riconoscimento di pattern, dato che la U' estrapolata rappresenta appunto una delta, quindi facilmente riconoscibile.

Voci correlate[modifica | modifica sorgente]

Bibliografia[modifica | modifica sorgente]

  • Wiener, Norbert (1949), Extrapolation, Interpolation, and Smoothing of Stationary Time Series. New York: Wiley. ISBN 0-262-73005-7
  • Brown, Robert Grover and Patrick Y.C. Hwang (1996) Introduction to Random Signals and Applied Kalman Filtering. 3 ed. New York: John Wiley & Sons. ISBN 0-471-12839-2


matematica Portale Matematica: accedi alle voci di Wikipedia che trattano di matematica