FILTRADO LINEAL OPTIMO: FILTRO DE WIENER
TRATAMIENTO AVANZADO DE SEÑAL EN COMUNICACIONES
Grupo de Tratamiento Avanzado de Señal (GTAS) Dpto. de Ingeniería de Comunicaciones (DICOM), Universidad de Cantabria
Indice • Derivación del filtro de Wiener. • Un prim primer er ejem ejempl plo o en co comu muni nica caci cion ones es:: Ig Igua uala laci ción ón.. • Cálculo teórico y estima del fi filtro de Wiener en en Ma Matlab. • Principio de or ortogonalidad. • Superficie de error. Propiedades. • Exte Extens nsio ione nes: s: IIR IIR cau causa sall y no caus causal al,, com compl plej ejo. o. • Filtrado óptimo con restricciones: “LCMV beamformer”, ejemplos. • Filtro adaptado estocástico (“eigenfilter”). • Conclusiones.
Indice • Derivación del filtro de Wiener. • Un prim primer er ejem ejempl plo o en co comu muni nica caci cion ones es:: Ig Igua uala laci ción ón.. • Cálculo teórico y estima del fi filtro de Wiener en en Ma Matlab. • Principio de or ortogonalidad. • Superficie de error. Propiedades. • Exte Extens nsio ione nes: s: IIR IIR cau causa sall y no caus causal al,, com compl plej ejo. o. • Filtrado óptimo con restricciones: “LCMV beamformer”, ejemplos. • Filtro adaptado estocástico (“eigenfilter”). • Conclusiones.
Planteamiento del problema: CASO FIR x(n)
−1
x(n-1)
z
w0
• Entrada: • Coeficientes:
−1
x(n-2)
−1
z
x(n-M+1)
z
w1
w M − 2
+
+
w M −1 +
T
x n = [ x(n), x(n − 1),..., x (n − M + 1)] T
w = [w0 , w1,..., w −1 ] M −1
• Salida:
d(n)
y (n) =
∑
T
x( n − i ) wi = x n w = w T x n
i =0
+
y(n) e(n)
Función de coste I • Error:
T
e( n) = d (n) − y (n) = d (n) − x n w
• La esti estima ma de mínimo mínimo error error cuadrá cuadrátic tico o medi medio o (MM (MMSE SE), ), es la solución que minimiza J ( w ) = E [e 2 ( n)]
• Desa Desarr rrol olla land ndo o la la fun funci ción ón de co cost ste e se se obt obtie iene ne T J ( w ) = E [d 2 (n) + w T x n x T 2 d ( n ) − w w x n ] n
T 2 T T w w x w R w w p J ( w ) = E [d 2 (n)] + w T E [ x n x T ] 2 E [ d ( n ) ] 2 − = σ + − n n d x
Función de coste II • Siendo
R x
la matriz de autocorrelación de la entrada
E [ x(n) x(n − 1)] E [ x(n) x(n)] E [ x (n − 1) x(n)] T R x = E x n x n = … E [ x(n − M + 1) x(n)]
[
]
R x (1) R x (0) R (1) R x (0) x R x = R x ( M − 1) R x ( M − 2)
… …
E [ x(n) x(n − M + 1)]
E [ x(n − 1) x(n − M + 1)] E [ x(n − M + 1) x(n − M + 1)]
… R x ( M − 1) …
R x (0) M × M
R x ( M − 2)
MATRIZ TOEPLITZ
M es la longitud del filtro FIR
Hemos tenido en cuenta que R x (k ) = R x (−k )
Función de coste III • p es el vector de correlación cruzada entre la entrada y la salida deseada E [d (n) x(n)] E [d (n) x(n − 1)] p = E [d (n) x n ] = [ − + ] E d ( n ) x ( n M 1 )
• Por ejemplo, para un filtro con dos coeficientes, la función de coste sería R x (0) R x (1) w0 p(0) 2 − 2(w0 w1 ) J ( w) = σ d + (w0 w1 ) p(−1) R x (1) R x (0) w1 J (w) = σ d 2 + w02 + w12 R x (0) + 2w0 w1 R x (1) − 2w0 p(0) − 2w1 p(−1)
Filtro de Wiener • La función de coste es cuadrática en los coeficientes del filtro
Solución única (filtro de Wiener) ∂ J ∇ J ( w ) = ∂w0 w
∂ J ∂w1
∂(σ d + w R x w − 2w p ) ∂ J = 2 R x w-2 p = 0 = ∂w M −1 ∂w T
…
2
T
Ec’s normales: R x w = p
T
w opt = R x p -1
¿Qué pasa si Rx es singular? • Sustituyendo la solución en la función de coste se obtiene el error cuadrático medio mínimo 2 − pT wopt J ( wopt ) = σ d
Ejemplo 1 retardo
• Igualación de un canal de comunicaciones digitales
s (n) símbolos
r (n)
canal
h
d (n) = s (n − d )
ruido
⊕
igualador
x(n)
w
y (n) -
+
e(n)
⊕
• Para resolver teóricamente el filtro de Wiener son necesarios los estadísticos de segundo order del problema, en concreto hay que conocer: – 1.- Modelo estadístico de la señal de entrada: i.i.d BSPK {-1,+1}. 2
– 2.- Modelo estadístico del ruido: AWGN de media cero y varianza σ r . – 3.- El canal h.
Ejemplo 1 • Para el modelo considerado: R x (k ) = Rhh ( k ) + σ r 2δ (k )
con
Rhh (k ) = h(k ) * h(− k )
p (−k ) = E [ s (n − d ) x(n − k )] =
= E s (n − d ) ∑ h(r ) s (n − k − r ) + r (n − k ) = h(d − k ) r • Ecuaciones normales: Rhh (0) + σ r 2 Rhh (1) 2 + σ r R ( 1 ) R ( 0 ) hh hh Rhh ( M − 1) Rhh ( M − 2)
… …
Rhh ( M − 1) w0
h(d ) R x ( M − 2) w1 h(d − 1) = 2 Rhh (0) + σ r w M −1 h(d − M + 1)
Ejemplo 1 • Fichero: Wiener_eq.m Wiener filter (blue), original BPSK (red)
1.5
1
h = [0.3 − 0.7
0.6
0.2]
0.5
SNR = 20 dB M = 15
d = 7
0
-0.5
-1
-1.5 0
10
20
30
40
50
60
70
80
90
100
Ejemplo 1 • Fichero: Wiener_eq.m Wiener filter (blue), original BPSK (red) 1 0.8
h = [0.3 − 0.7 SNR = 20 dB M = 15
d = 0
0.6
0.2]
0.6 0.4 0.2 0 -0.2 -0.4 -0.6 -0.8 -1
0
10
20
30
40
50
60
70
80
90
100
¿Por qué es necesario un retardo? • Suponemos una situación sin ruido y un canal causal y estable. • Si el canal es de fase no mínima (alguno de sus ceros fuera del círculo unidad) el inverso estable es no causal. hi [n] n=0
• Si no incluimos un retardo, la aproximación del filtro de Wiener al inverso del canal será muy mala: • Si permitimos un cierto retardo:
d=0
w[n] n=0 w[n]
d≠0 n = d
Estima del filtro de Wiener I • En un caso práctico no se conoce el canal ni la varianza de ruido, tan sólo se conoce la señal deseada s(n), en un intervalo n=1,…N. (Secuencia de entrenamiento). • Por ejemplo, en GSM
0.577 mseg.
Rb = 270.833 kb / s Hay 26 símbolos conocidos por el receptor
Estima del filtro de Wiener II ˆ = R x
• La matriz de autocorrelación se estima como • La correlación cruzada
pˆ =
1
N
1
N
∑
N n =1
T
x n x n
s (n) x ∑ N
n
n =1
• Si queremos incluir un retardo d pˆ =
1
N
s (n) x ∑ N
n + d
n =1
ˆ = R x
1
N
∑ N
T
x n + d x n + d
n =1
• El igualador óptimo bajo un criterio MMSE es ˆ −1 pˆ wˆ opt = R x
“Direct Matrix Inversion” (DMI)
Ejemplo 2
h = [1 − 0.3 0.6] SNR = 25 dB
• Fichero: Wiener_eq_est.m - 0.4800 0.6000 1.4546 - 0.4800 1.4546 - 0.4800 1.4546 R x = 0.6000 - 0.4800 0.6000 - 0.4800 0 0 0 0.6000
M = 5
0.60 0 - 0.48 0.6 1.4546 - 0.4800 - 0.4800 1.4546 0
Fase mínima
0
1.4734 - 0.5413 0.5827 - 0.0655 - 0.0346 - 0.5413 1.4650 - 0.5400 0.6066 - 0.0770 ˆ = 0.5827 - 0.5400 R 1.4738 - 0.5518 0.5925 x 0.0655 0.6066 0.5518 1.4132 0.5050 - 0.0346 - 0.0770 0.5925 - 0.5050 1.4253
d = 0
N = 50
1 0 p = 0 0 0
0.9452 0.2623 wopt = - 0.4150 0.2117 0.1013
1.0065 - 0.0415 pˆ = - 0.0430 0.0188 - 0.0508
0.9380 0.2521 wˆ opt = - 0.4275 0.2074 0.1050
Ejemplo 2 teórico (azul), estimado (rojo) 1.5
1
0.5
0
-0.5
-1
-1.5
0
5
10
15
20
25
30
35
40
45
50
Principio de ortogonalidad I • El filtro de Wiener se puede también derivar a partir del principio de ortogonalidad: EL ERROR ES ORTOGONAL A LOS DATOS.
E [e( n) x n ] = E [ x n ( d ( n) − x nT w )] = 0
E [ x n d ( n)] = E [ x n x nT ]w −1
w opt = Rx p
p
R x
Principio de ortogonalidad II • Interpretación geométrica e(n)
d (n)
x −∞ x n − 2
x n x n −1
T
w x n
• Corolario: el error es ortogonal a la salida estimada E [e( n) w T x n ] = w T E [e(n) x n ] = 0
Superficie de error I • La función de coste es un paraboloide con un único mínimo J ( w ) = σ d 2 + w T R x w − 2w T p
• Si consideramos un filtro de 2 coeficientes J ( w ) =
2 σ d
1 r w0 − 2(w0 w1 ) r 1 w1
+ (w0
6
4
w opt
2
0
-2
-4
-6 -6
-4
-2
0
2
4
6
p (0) w1 ) p(−1)
Superficie de error II • El eje del paraboloide es perdendicular al plano J(w)=cte. • El mínimo está en
−
w opt = Rx 1 p
• Su orientación y forma dependen exclusivamente de R x • La función de coste puede escribirse como J = J min + ( w − w opt )T R( w − w opt ) ~ T R w ~ J = J min + w x
siendo
~ = w−w w opt
Superficie canónica I ~ • w = w − w opt define una traslación: el mínimo de la función de ~ =0 coste está ahora en w opt
~) J ( w
J ( w ) 6
8 5 8 7 . 0 = ) 1 ( t
p o
w
6
4 4
w opt
2
0 = ) 1 ( t
0
p o
-2
w
-4
~ w opt
2
0
-2
-4
-6 -6
-4
-2
0
2
4
wopt (0) = 0.8360
6
-6 -6
-4
-2
0
2
~ (0) = 0 w opt
4
6
Superficie canónica II • La matriz de autocorrelación puede escribirse en términos de sus autovectores y autovalores como
R x = QΛQ = [q1 T
…
λ 1 0 q M ] 0
0
…
λ 2
…
0
0
T q 1 M = ∑ λ i qi qiT 0 T i =1 q M λ M
~ T Q Λ Q T w ~ J = J min + w
• Definimos ahora los vectores transformados (linealmente) T ~ v =Q w
Superficie canónica III • Se obtiene así la forma o superficie canónica J = J min + v Λ v = T
M
∑
λ k
v k 2
i =1
~ w 1
~ v =Q w T
~ w 0
v1
v0
• Las direcciones en el espacio transformado (canónico) están alineadas con las direcciones principales del paraboloide. • La excentricidad depende de los autovalores.
Un ejemplo • Matriz de autocorrelación
2 1 Rx = 1 2
• Autovalores y autovectores 1 (1 − 1) λ 1 = 1 ⇒ v1 = 2 R x − λ I = 0 ⇒ 1 λ 2 = 3 ⇒ v 2 = (1 1) 2
Superficie de error: ~ ) = J + (w ~ J ( w min 0
Variables acopladas
1 1
1 1
0 1 − 1
2 − 1 1 0 3 1
1
Superficie canónica:
~ w 2 1 0 ~ ~ = w1 ) 1 2 w1
~ 2 + 2w ~w ~ + 2w ~2 = 2w 0 0 1 1
R x = QΛQ = T
J (v ) = J min + (v0
1 0 v0 2 = v0 + 3v12 0 3 v1
v1 )
v1
~ w 1 ~ w 0
Variables desacopladas
v0
Caso general • Matriz de autocorrelación
siendo r =
R x (1) R x (0)
σ 2 R x = 2 r σ
r σ 2
(1 + r )σ 2 0 ⇒ Λ = 2 2 (1 − r )σ σ 0
el coeficiente de correlación.
r=0
r=0.7
r=0.3
v1
v1
v1
v0
• Dispersión de autovalores:
v0
λ max λ min
=
v0
1 + r 1 − r
• Un valor alto indica mal condicionamiento de la matriz de autocorrelación.
Estimación lineal vs. no lineal • Hemos visto hasta ahora que el filtro de Wiener es el estimador LINEAL que minimiza el MSE
ˆ ( n) = w T x d opt n
con
−1
w opt = Rx p
Es necesario conocer los estadísticos de segundo orden. • El estimador que minimiza el MSE es, en general, no lineal y se obtiene como ˆ (n) = E [ d ( n) | x ] d n
Es necesario conocer:
f d | x n (d ( n) | x n )
Estimación lineal vs. no lineal: ejemplo Est. lineal wopt =
equiprobables s(n) ∈{+1,−1} h0
AWGN
h0
sˆ(n) = wopt x(n)
ho2 + σ r 2
r (n)
⊕
x(n) x(n)h0 Est. no lineal sˆ(n) = tanh
h0 / σ r 2
tanh (⋅)
2
σ r
Algunas extensiones • Hasta ahora nos hemos concentrado en el caso del filtrado óptimo FIR con secuencias reales. • A continuación extendemos las ideas expuestas sobre filtrado óptimo a: 1.- Filtros IIR (no causales) 2.- Filtro IIR (causales) 3.- Filtros FIR complejos
Filtro de Wiener IIR d(n) y(n) - +
wn
x(n)
wn ≠ 0
e(n)
−∞ < n < ∞ ∞
• La salida del filtro IIR es
y (n) =
∑ w x(n − k ) k
k = −∞
• El principio de ortogonalidad sigue siendo válido min E e 2 ( n) ⇒ E [e( n) x ( n-l )] = 0
-∞ < l < ∞
∞
Ec’s Normales:
∑ w R
k xx
k = −∞
(l − k ) = p (l )
− ∞ < l < ∞
Filtro de Wiener IIR • Es posible escribir la solución en el dominio de la frecuencia:
W (ω ) =
S xd (ω ) S xx (ω )
• Para el problema de igualación/deconvolución visto anteriormente r (n) s(n)
h
⊕
• Si no hay ruido:
x(n)
w
W (ω ) =
W (ω ) =
y(n)
S xd (ω ) S xx (ω )
=
H (ω )* H (ω )
2
=
S xd (ω ) S xx (ω )
=
H (ω )* 2
H (ω ) + σ r 2
1 H (ω )
El filtro de Wiener es el filtro inverso (en general, IIR no causal)
Filtro de Wiener IIR (causal) • Es posible también obtener el filtro óptimo IIR causal ∞
∑ w R (l − k ) = p(l ) k x
l ≥ 0
k = 0
• La solución de este problema es algo más complicada S xx ( z ) = σ x2Q ( z )Q ( z −1 )
Q ( z )
es fase mínima
• Entonces: S xd ( z ) 1 W ( z ) = −1 Q ( z ) + σ x2Q( z )
Parte causal
Filtro de Wiener complejo • En el caso de secuencias complejas (filtro FIR) y ( n) =
M −1
∑ x(n − i)w
*
i =0
i
= w H x n
H denota hermítico (conjugado y transpuesto) • La solución del filtro óptimo sigue siendo −1
w opt = Rx p
pero ahora R x = E [ x n x nH ] p = E [ x n d ( n)] *
Rˆ x = pˆ =
1
N
∑ N 1
H
x n x n
n =1 N
∑ N
*
x n d ( n)
n =1
Filtrado lineal óptimo con restricciones • En ciertos problemas de comunicaciones, por ejemplo en conformación de haz (“beamforming”), es necesario minimizar el MSE sujeto a un restricción lineal.
Señal deseada
d
x1 ( n )
*
x2 ( n )
*
y(n) i
Señal interferente
x M (n )
*
M
minimizar
E [| y (n) |2 ] = w H R x w
s.t.
w s (θ 0 ) = g H
Modelo de propagación ( e n s d
d
s1 (t ) = s0 (t − τ (θ ))
El retado entre dos elementos del array es
s0 (t ) τ (θ )
=
d sen(θ ) c
Hipótesis: onda plana, banda estrecha, array equiespaciado: s0 (t ) = u (t )e jω 0t ⇒ s1(t ) = u (t − τ (θ ))e jω 0t e − jω 0τ (θ ) ≈ u(t)e jω 0t e − jω 0τ (θ ) f d − j 2π 0 d sen (θ ) − jπ sen (θ ) − jω 0τ (θ ) c λ / 2 = s0 (t )e = s0 (t )e s1 (t ) = s0 (t )e
“Linearly constrained minimum variance” (LCMV) • PROBLEMA:
minimizar E [| y (n) |2 ] = w H R x w siendo s (θ 0 ) = (1 e − j
θ 0
…
s.t.
w s (θ 0 ) = g H
)
j ( M −1)θ 0 T
e−
• Aplicamos el método de los multiplicadores de Lagrange J ( w ) = w H R x w + λ R Re w H s (θ 0 ) − g R + λ I Im w H s (θ 0 ) − g I J ( w ) = w H R x w + Re λ * w H s (θ 0 ) − g
“Linearly constrained minimum variance” (LCMV) • Tomando derivadas e igualando a cero ∇ J = 2 R x w +
*
λ
2
w opt = −
s (θ 0 ) = 0
*
λ
2
−1
Rx s (θ 0 )
• Para obtener el multiplicador de Lagrange w s (θ 0 ) = g = − H opt
λ
2
H
−1
λ =
s (θ 0 ) Rx s (θ 0 )
−
2 g s (θ 0 ) R x− s (θ 0 ) H
1
• El LCMV o “Minimum variance distorsionless response” (MVDR) beamformer es wopt =
g * R x−1 s (θ 0 ) H
−1
s (θ 0 ) R x s (θ 0 )
MVDR Spectrum • Tomando g=1, la potencia mínima a la salida del conformador es J min =
que θ 0
1 s H (θ 0 ) R x−1 s (θ 0 )
mantiene la respuesta minimizando el resto.
(potencia)
en
la
dirección
• Considerando la expresión anterior como función de la frecuencia espacial θ , o temporal ω S (ω ) =
1 H
−1
s (ω ) R x s (ω )
“MVDR spectrum” o “Capon spectrum”
Ejemplo LCMV • Programa: Lcmv.m
-Señal QPSK deseada en
θ 0
= −20
-Interferencia 1 en
θ 1
= 20 con SIR=-10 dB
-Interferencia 2 en
θ 2
= −60 con SIR=-15 dB
-Ruido gaussiano y espacialmente blanco con SNR=0 dB - Array Array lineal con M=10 elementos, separación entre antenas d=λ /2.
Ejemplo LCMV N=25 “snapshots” Una realización
10 realizaciones ind. 10
10
0
-10 0 -20 -10 -30 -20 -40
-50
-30
-60 -40
-80
-60
I2
-40
-20
D
0
20
I1
40
60
80
-80
-60
I2
-40
-20
D
0
20
I1
40
60
80
Ejemplo LCMV N=250 “snapshots” Una realización
10 realizaciones ind.
10
10
0
0
-10
-10
-20
-20
-30
-30
-40 -40 -50 -50 -60 -80
-60
I2
-40
-20
D
0
20
I1
40
60
80
-80
-60
I2
-40
-20
D
0
20
I1
40
60
80
Filtro adaptado estocástico I • Otro problema filtrado óptimo que aparece frecuentemente en comunicaciones es el siguiente:
v(n) s (n)
⊕
x (n)
w
y (n)
¿Cúal es el filtro que maximiza la SNR a la salida? H
• Potencia de señal a la salida del filtro:
w R s w
• Potencia de ruido a la salida del filtro:
w Rv w
SNR máxima
H
H
maximizar
SNR =
w R s w H
w Rv w
Filtro adaptado estocástico II • Si suponemos que el ruido es blanco SNR =
maximizar
Rv = σ v I 2
1 w H R s w 2
σ v
H
w w
• El problema es equivalente a
maximizar
H
w R s w H
w w
s.t.
w w =1 H
• La solución es el autovector correspondiente al máximo autovalor de : R s w opt = qmax
siendo
R s qmax = λ max qmax
y la
SNRmax =
λ max 2
σ v
• Puede verse como un filtro adaptado estocástico, también se denomina en la literatura “eigenfilter”.
