El Filtro de Kalman Antonio Miguel Fumero Reverón
Temas Especiales de Análisis Numérico José Manuel Corrales Sendino
Madrid. Día 18, Febrero de 2000
________________________________________________________________ TEAN
Motivación y Objetivo La elaboración de este documento parte del interés personal del alumno por ciertos temas relacionados con las comunicaciones, y los problemas clásicos que ha planteado la comunicación entre entidades (máquinas) dotadas de cierta inteligencia (entendida en el sentido de la posesión de una capacidad de decisión determinada). Este interés surge en cierta medida tras la lectura de los artículos de Shannon acerca de la teoría matemática de la comunicación, y sobre todo, de la conocida obra de Norbert Wiener “Cibernética”, originalmente editada en 1948 por el MIT bajo el título de “Cybernetics or Control Control and Commun Communicat ication ion in the Animal Animal and the Machine” Machine” . Llama la atención la
analogía entre los problemas de la predicción estadística y el control determinístico; así como la cantidad de trabajo que ha surgido en torno a las ecuaciones normales de Wiener. El interés por el tema se ha visto incrementado tras la asistencia a un pequeño seminario acerca de los sistemas no lineales y de control adaptativo que se impartió por vez vez prim primer era a como como asig asigna natu tura ra de libr libre e elec elecci ción ón en el P94 P94 dura durant nte e el segu segund ndo o cuatrimestre del curso 1998/1999. En el seminario se presentaba el problema de Wiener y sus aplicaciones en el mundo del control adaptativo y ecualización de canales de comunicación. El objetivo principal del documento es volver a plantearse esos problemas, así como las soluciones que propuso en su día R. E. Kalman, que han dado lugar a una serie de algoritmos numéricos de amplia aplicación en diversos campos de nuestra técnica. Para ello se ha acudido a algún texto de filtrado adaptativo y ecualización, así como a los artículos originales de Kalman donde propone su solución del problema de Wiener y sus posibles aplicaciones. Se han pretendido usar las ideas y conceptos presentados en el curso de TEAN, impartido en el primer semestre del presente curso 1999/2000, para elaborar un estudio del problema y las aplicaciones del Filtro de Kalman con el objetivo de que sirva de
______________________________________________________________________
________________________________________________________________ TEAN alguna forma para la evaluación de la asignatura en términos de su mayor o menor grado de aprovechamiento. Dicho sea ya antes de nada que el autor se declara una persona más preocupada por los conceptos involucrados en el problema que se presenta y su alcance en cuanto a sus implicaciones en las técnicas de uso común en la ingeniería, que por el detalle operativo de los métodos numéricos que se emplean en la resolución de los problemas concr concret etos os.. Este Este hecho, hecho, quiz quizás ás moti motiva vado do por por la leja lejaní nía a de los los prim primero eross años años de universidad en que uno se acostumbraba a lidiar en el día a día con el utillaje y la operativa operativa matemática, matemática, podría llevar a pensar que estas líneas no constituyen constituyen otra cosa que un de ejercicio de demagogia demagog ia o disquisición metafísica sin utilidad práctica alguna.
______________________________________________________________________
________________________________________________________________ TEAN
Introducción Aunque la difusión que va a recibir este documento no va a ser en absoluto amplia, creo que no es mala idea comenzar con una pequeña introducción al problema del filtrado adaptativo, estimación lineal, y optimización que servirán al menos para establecer unos términos básicos sobre los que se desarrollará todo el documento. y(n)
Filtro Programable
x(n)
- Σ+ ŷ(n)
e(n)
h(n) Algoritmo Adaptativo Figura 1:Diagrama de Bloques de un Filtro Adaptativo Genérico (Tomada de [Mulgrew88])
Lo primero es dejar claro que cuando hablamos de filtros adaptativos estamos hablando de sistemas lineales programables; es decir de filtros lineales que operan sobre una serie temporal de entrada, pero cuyos coeficientes (que definen la operación de filtrado que realizan sobre esa entrada) se pueden variar a voluntad voluntad de un instante instante temporal al siguie siguiente nte (de una muestra muestra a la siguient siguiente) e) median mediante te un algori algoritmo tmo adaptati adaptativo vo que se basará en alguna relación recursiva: hi(n+k) = hi(n) + g{x(n),ŷ(n),e(n)} donde hi(n) es el coeficiente i-ésimo del filtro en el instante temporal n, x(n) es un vector de las entradas entradas pasadas del filtro, filtro, ŷ(n) es un vector de las salidas pasadas del filtro, y
e(n) = y(n) - ŷ(n) es un vector que representa el pasado de la señal de error. La est estruct ructur ura a bási básica ca que que se repr repres esen entta en la Figu Figura ra 1 tien tiene e tres res modo modoss de funcionamiento fundamentales:
______________________________________________________________________
________________________________________________________________ TEAN (i)
Identificación de Sistemas (Figura 2(i)) : se trata de alimentar al sistema
desconocido que se pretende identificar, y a la estructura adaptativa de antes con la misma serie temporal de entrada; de forma que tras la convergencia convergencia la salida del filtro filtro adaptativo adaptativo ŷ(n) se aproximará a la salida del sistema desconocido y(n) en un sentido óptimo (normalmente se usará el criterio de convergencia en media cuadrática). cuadrática). En esa situación situación se dirá que hemos identifica identificado do el sistema. Normalmente Normalmente tendremos tendremos ruido blanco en la entr entrad ada a (err (error ores es en las las obse observ rvac acio ione nes) s) que que nos nos alej alejar ará á del del comportamiento ideal que se pretende. (ii)
Problema Inverso de Análisis o Modelado de Sistemas(Figura 2(ii)) : es
típico de las estructuras que se usan para la ecualización de canales de comu comuni nica caci ción ón.. En este este caso caso se trat trata a de que que la sali salida da del del sist sistem ema a desconocido sirva de entrada para el filtro adaptativo. Lo que se pretende es obtener la salida del filtro adaptativo una versión retardada de la señal de entrada, con lo cual, tras la convergencia, la función de transferencia del filtro coincidirá con la inversa de la función de transferencia del sistema desconocido que lo precede (típicamente el canal de comunicaciones). (iii)
Predicción Predicción Lineal(Figura Lineal(Figura 2(iii)) : se usa mucho en los codificadores en
canales de comunicación vocal. En este caso le metemos como entrada al filtro adaptativo una versión retardada de la misma señal que deseamos tener a la salida del mismo; con lo cual se trata de predecir las muestras futuras de la entrada. Este comportamiento sólo se consigue si la señal de entrada está bastante lejos de tener un espectro plano. Con estas ideas más o menos claras y algunos conceptos de sistemas lineales lineales que se han ido destilando a lo largo de los años pasados en la Escuela, pienso que sería buena idea plantear el problema de la estimación lineal óptima y la ecuación de Wiener que es la motivación principal de toda esta es ta historia. El problema de la estimación lineal se puede expresar de una forma más o menos rigurosa con un enunciado como este: Dada una secuencia aleatoria de observaciones
______________________________________________________________________
________________________________________________________________ TEAN {x(n)}, que se corresponde con una versión distorsionada de una señal {y(n)} (otra secuencia aleatoria que transporta cierta información), encontrar un filtro lineal que opere sobre {x(n)} para obtener una estimación { ŷ(n)} de {y(n)}. La calidad de la estimación se medirá con una función f(.) de la secuencia de error {e(n)}, considerado este como la diferencia e(n) = y(n) - ŷ(n), que asignará un precio o penalización a las estimaciones incorrectas (de ahí que habitualmente se la conozca como función de pérdidas). Se trata de una función que debe ser positiva (o no negativa), y decreciente. Puesto que la secuencia de error también es una secuencia aleatoria, se suele elegir el filtro de forma que minimice una función de coste l(.) que sea la esperanza matemática de la función de pérdidas: l(e(n)) = E[f(e(n))] . La función de coste más habitual suele ser el error cuadrático medio (MSE) ξ(n) = E[e2(n)]. Sistema Desconocido
x(n)
y(n)
- Σ+
Filtro Adaptativo
Figura 2(i):Identificación de Sistemas (Tomada de [Mulgrew88])
El filtro FIR óptimo se obtendría siguiendo estas considera ciones:
ŷ(n) es la salida de un filtro causal, y como tal se puede expresar como la suma de convolución de la secuencia de entrada y la respuesta al impulso, que se puede truncar para obtener un filtro FIR de orden N-1
ŷ(n) = Σhix(n-i) ; (0≤i≤N-1) que si se expresa en forma vectorial es algo como
ŷ(n) = hT x(n) ; h = [h0 h1 ... h N-1]T ∧ x(n) = [x(n) x(n-1) ... x(n-N+1)]T La función de coste MSE la podemos poner como ξ = E[y2(n)] + hT xx
h – 2hT
xx
xy
= E[x(n)xT(n)] ∈ ℳN×N es la matriz de autocorrelación
______________________________________________________________________
________________________________________________________________ TEAN xy
= E[x(n)y(n)] ∈ ℳN×1 es el vector de correlación cruzada
Para minimizar la función de coste, igualaremos el gradiente a cero: ∇ = ∂ξ/∂h = 2
h–2
xx
xy
=0⇨
h
xx opt
=
xy
Resulta que si la densidad espectral de potencia de la secuencia de entrada {x(n)} no tiene nulos, la matriz de autocorrelación es definida positiva, y por tanto no singular. Y sabemos que bajo esas condiciones la respuesta al impulso óptima es única y viene dada por -1 xx
hopt =
xy⇨ξopt
= E[y2(n)] - hTopt
xy
Esa ecuación es la que se conoce como el filtro FIR de Wiener (o filtro de Levinson). Lo que que ocur ocurre re aquí aquí es que que neces necesititam amos os la matr matriz iz de auto autocor correl relac ació ión n y el vect vector or de correlación cruzada; es decir necesitamos estadísticos de segundo orden cuando lo que normalmente se conoce son las secuencias aleatorias de datos. La labor de determinar el filtro óptimo a partir de esos datos se le suele encomendar a un filtro adaptativo, que se puede ver en esos términos como un algoritmo que opera sobre las secuencias {x(n)} e {y(n)} para formar a partir de ellas un vector de respuesta al impulso h(k) variable con el tiempo que converge en media a la respuesta al impulso óptima. El filtro de Levinson será por tanto el objetivo que debe perseguir el filtrado adaptativo.
x(n)
Sistema Desconocido
Filtro Adaptativo
- Σ+
Figura 2(ii):Problema Inverso de Modelado o Análisis de Sistemas (Tomada de [Mulgrew88])
El cálculo de la hopt requiere la solución de N ecuaciones lineales simultáneas con N incógnitas. Para una matriz no singular en general el método más eficiente sabemos que es el de eliminación de Gauss, que requiere O(N 3) operaciones, siendo también N el número de coeficientes del vector de respuesta al impulso. De lo que se dio cuenta Levinson fue del hecho de que la matriz de autocorrelación tiene una estructura muy especial: (i)
es simétrica ⇨
[i,j] =
xx
[j,i]
xx
______________________________________________________________________
________________________________________________________________ TEAN (ii) siendo
y es de Toeplitz⇨
[i,j] =
xx
[i-j]
xx
[i,j] el elemento en la fila i, columna j de la matriz de autocorrelación, y {
xx
(m)}
xx
la secuencia de autocorrelación asociada con la secuencia estacionaria en sentido amplio {x(n)}. Levinson se dio cuenta de que con esta estructura sólo necesitaba la primera fila para definir toda la matriz; con lo cual este hombre ideó un algoritmo que sólo requería O(N2) operaciones. Hay también otra forma más intuitiva de atacar el problema que se conoce como SMI (de Sampled Matrix Inversion) que consta de dos fases: en una primera fase se hace una estimación de los coeficientes de correlación y autocorrelación a partir de las secuencias de datos disponibles; y en la segunda fase se usan esos coeficientes estimados para construir la ecuación de Wiener, que se resuelve luego con el algoritmo de Levinson. Cuando sólo tenemos conjuntos finitos de datos (x(n), y(n), n=0,...k), se sustituye el opera operado dorr espe esperan ranza za mate matemá mátitica ca de la funci función ón de cost coste e MSE MSE por por un suma sumato tori rio, o, obteniéndose lo que se llama la función de coste LS ( por Least Squares). Al final resulta que el problema de la minimización LS es dual del de la minimización MSE quedando algo como
r xx xx(k)h(k) = r xy xy(k) T r xx xx (k) = Σ(n=0,...,k) x(n)x (n) ∈ ℳN×N es la matriz de autocorrelación
r xy xy (k) = Σ(n=0,...,k) x(n)y(n) ∈ ℳN×1 es el vector de correlación cruzada Existirá solución única si existe r xx xx (k) y es no singular. Para k≤N-1, r xx xx (k) será siempre singular. Para k≥N, la singularidad de r xx xx (k) dependerá de los conjuntos de datos particulares que tengamos disponibles. El problema de la singularidad aveces se puede soslayar usando la pseudoinversa de Moore-Penrose, y hay autores que lo han propuesto. Hay varias formas de construir estimaciones LS dependiendo de las hipótesis que se hagan hagan sobre sobre los datos datos dispon disponibl ibles: es: se puede puede usar direct directamen amente te la covaria covarianza, nza, la autocorrelación, o el enventanado (pre- o post- enventanado). También existe un algoritmo recursivo RLS para los casos en que se hace necesario actualizar la estimación LS según vamos teniendo nuevos datos disponibles.
______________________________________________________________________
________________________________________________________________ TEAN Hay toda una serie de algoritmos, algoritmos, que se conocen como algoritmos algoritmos rápidos. Surgieron Surgieron a partir de la publicación del FKA (Fast Kalman Algorithm) que trataba de aprovechar la estructura estructura peculiar de la matriz r xx estimación LS; le siguieron siguieron más algoritmos algoritmos xx (k) de la estimación buscando una mejora en la eficiencia computacional: FAEST (Fast A poSteriori Error Technique), FTF (Fast Transversal Filter). x(n)
Retardo
- Σ+
Filtro Adaptativo
Figura 2(iii):Predicción Lineal (Tomada de [Mulgrew88])
También los métodos estocásticos de gradiente dan técnicas iterativas donde el número de operaciones operaciones que se deben llevar llevar a cabo para encontrar la solución solución no se conoce de antemano. El método de descenso en el sentido del gradiente es una técnica iterativa que se usa en programación lineal y problemas de optimización para encontrar una variable que minimice una función de coste objetivo. Para aplicarlo a la ecuación de Wiener: se hace una estimación de la solución hi (siendo i el paso en el proceso de iteración), que tendrá asociada un valor de la función de coste MSE, ξ(hi). Luego se hace una mejora de la estimación en dos fases: en la primera se calcula el gradiente de la misma,
i
que será un vector que irá en la dirección de la máxima pendiente de la
func funció ión n de coste coste MSE; MSE; en la segu segunda nda fase fase se cons constr truy uye e una una nueva nueva estim estimac ació ión n restándole a la primera una versión escalada del gradiente µ (hi). Queda algo así:
hi+1 = hi - µ (hi) i
= ∂ξ/∂h(hi) = 2
h –2
xx i
xy
como condición de convergencia tenemos que 0<µ<1/λmáx siendo λmáx el mayor autovalor de la matriz de autocorrelación Este método evita la necesidad de la inversión de la matriz de autocorrelación, pero sigue necesitando conocer esa matriz y el vector de correlación cruzada. Pero hay un
______________________________________________________________________
________________________________________________________________ TEAN algoritmo que no necesita el conocimiento de esos estadísticos, es el LMS (por Least Mean Squares). Este algoritmo sustituye el proceso de iteración por una recursión temporal y el gradiente por una estimación del mismo.
h(k+1) = h(k) - µ (k) (k) = ∂ξ/∂h(h(k)) = ∂/∂h E[(y(n) – hT(k)x(n))2] = -2E[ x(n)(y(n) – hT(k)x(n))] (k) = -2x(k+1)e(k+1) e(k+1) = y(k+1) – hT(k)x(k+1) Al final el algoritmo se puede expresar como
h(k+1) = h(k) + 2µx(k+1)e(k+1) que es bastante más simple que el RLS. Existe una variante, el BLMS, en donde se mantiene constante la estimación de la respuesta al impulso en un bloque de L puntos de la secuencia de datos. Da muy buenos resultados cuando se usa junto con la FFT.
El Filtro de Kalman Una vez que estamos de acuerdo sobre donde comienza toda esta historia, si nos vamo vamoss ahor ahora a a mira mirarr en cual cualqu quie ierr text texto o dedi dedica cado do al filt filtra rado do adap adapta tatitivo vo y sus sus
______________________________________________________________________
________________________________________________________________ TEAN aplicaciones, y buscamos alguna referencia al filtro de Kalman, encontraremos en muchos su utilización como Ecualizador IIR adaptativo justificada más o menos en los siguientes términos. El filt filtro ro IIR IIR de Wien Wiener er,, que que resul resulta ta ser ser bast bastan ante te mejo mejorr que el FIR FIR para para canal canales es conocidos, comienza a dar problemas cuando, ante un canal desconocido o variable en el tiempo exige la utilización de estructuras adaptativas. Si la estimación se realiza con el mismo algoritmo LMS que se usa para el caso FIR, existe cierta probabilidad de que los polos del filtro se salgan de la circunferencia unidad; pudiendo dar lugar a una situación de inestabilidad si se mantienen así durante duran te mucho tiempo. La idea dea es que que la real realiz izac ació ión n del del ecua ecuallizad izador or IIR óptim ptimo o de Wien Wiener er impli mplica ca consideraciones y conceptos relacionados con el espacio de estados, y eso es lo que nos lleva al planteamiento original que Kalman presenta en sus artículos de principios de los años 60. Lo que se hace es usar el Filtro de Kalman como ecualizador; y para darle carácter adaptativo a la estructura se usa en paralelo con un algoritmo de identificación de sistemas que nos dará los M coeficientes de la respuesta al impulso del canal (una estimación que se minimizará con un criterio MSE). También hay autores que han propuesto soluciones donde meten los coeficientes del canal en la información asociada al estado del modelo de Kalman (es lo que llaman el Filtro de Kalman Extendido, EKF). El modelo de Kalman se puede poner, en tiempo discreto y con una notación más o menos estándar como sigue.
S(k+1) = A(k+1)s(k) + w(k+1) z(k) = H(k)s(k) + u(k)
s(k) es lo que se llama el vector de estado y contiene los valores de los M parámetros que definen el estado del sistema en el instante ins tante k.
______________________________________________________________________
________________________________________________________________ TEAN
A(k) es una matriz de dimensiones M ×M, que se conoce como la matriz de transición de estado. (k) es una una matri atrizz de dime dimens nsiiones ones L×M, que se conoce com como la matriz de H(k) observaciones.
w y u son dos vectores de M y L componentes respectivamente; y representan a dos procesos blancos incorrelados y de media nula de covarianzas W y U, de forma que: E[w(k)wT(k)] = W(k) E[u(k)uT(k)] = U(k) Se supone que las matrices A, W, H, U, aunque pueden variar con el tiempo, son conocidas a priori. Las estimaciones óptimas del vector de estado s(k) se generan de forma recursiva a part partir ir de la secu secuenc encia ia de obser observac vacio iones nes ruid ruidosa osass {z(k)}. (k)}. Las ecuaci ecuaciones ones pueden pueden ponerse como sigue:
ESTIMACIÓN Ŝ(k/k) = Ŝ(k/k-1) + K(k)[z(k) – H(k) Ŝ(k/k-1)] PREDICCIÓN Ŝ(k/k-1) = A(k) Ŝ(k-1/k-1) GANANCIA K(k) = V(k/k-1)HT(k)[H(k)V(k/k-1)HT(k) + U(k)]-1 COVARIANZA DEL ERROR V(k/k-1) = A(k)V(k-1/k-1)AT(k) + W(k) V(k-1/k-1) = V(k-1/k-2) – K(k-1)H(k-1)V(k-1/k-2) V(k/k) = E[( s(k) – Ŝ(k/k))(s(k) – Ŝ(k/k))T] V(k/k-1) = E[(s(k) – Ŝ(k/k-1))(s(k) – Ŝ(k/k-1))T]
Al final, ¿Qué es lo que pasa aquí? ¿Qué es lo más molesto en el filtro IIR de Wiener? Cons Conseg egui uirr la facto actori riza zaci ción ón de fase ase míni mínim ma del del espec spectr tro o de pot potenci encia a de las las observaciones.
______________________________________________________________________
________________________________________________________________ TEAN Si todos los procesos son estacionarios, y el ruido en las observaciones es blanco, el filtro de Kalman para el estado estacionario y el filtro IIR de Wiener son idénticos. Un filtro FIR nos lleva directament directamente e a una representación representación sobre la base del espacio de estados. Pero hay que tener cuidado: un vector de estado con (M-1) estados nos llevará a una formulación formulación donde el ruido de planta y el de las observaciones observaciones están correlados, y un filtro de Kalman en esa situación es condicionalmente estable. Con M estados representando representando al canal, tendríamos sin embargo una formulación formulación en la que el ruido de planta y el de las observaciones estarían incorrelados, y el filtro de Kalman sería incondicionalmente estable. Para Para mont montar ar un ecua ecualiliza zado dorr IIR IIR con con el filt filtro ro de Kalm Kalman an tendr tendría íamo moss la sigu siguie ient nte e formulación:
s(k) = [s(k) s(k-1) ... s(k-M+1)]T es el estado del modelo FIR de canal s(k) = as(k-1) + bs(k) a es una matriz de desplazamiento de M ×M bT = [1 0 0 ... 0]
↪x(k) = hT(k)s(k) + n(k) {s(k)} y {n(k)} son secuencias de ruido blanco, incorreladas, y de media nula con varianzas σs2(k) y σn2(k). Este es un planteamiento válido para cuando sólo necesitamos una estimación del estado actual del canal. Pero si tratamos con canales que no son de fase mínima, necesitaremos montarnos otra historia:
s(k) = [s T (k) s T (k-1) ... s T (k-d)]T El vector de estado se segmenta para contener (d+1) elementos
s T (k) = [s(k) s(k-1) ... s(k-M+1) ... s(k-d)]T Los vectores a y b se construyen igual que antes, pero con (d+1) elementos La matriz de observaciones se construye rellenando con ceros la respuesta al impulso del canal:
H(k) = [h0(k) h1(k) ... hM-1(k) 0 0 ... 0]
______________________________________________________________________
________________________________________________________________ TEAN Con estas consideraciones se pueden escribir las ecuaciones del ecualizador IIR de Kalman, cuya estructura se representa en la Figura 3:
ESTIMACIÓN Ŝ(k/k) = Ŝ(k/k-1) + K(k)[x(k) – H(k) Ŝ(k/k-1)] PREDICCIÓN Ŝ(k/k-1) = aŜ(k-1/k-1) GANANCIA K(k) = V(k/k-1)HT(k)[H(k)V(k/k-1)HT(k) + σn2(k)]-1 K(k) = [k0(k) k1(k) ... kd(k)]T COVARIANZA DEL ERROR V(k/k-1) = aV(k-1/k-1)aT(k) + bbTσs2(k) V(k/k) = [I – K(k)H(k)]V(k/k-1)] X(k)
k0
k1
z-1
0
h0
h1
kM-1
z-1
Kd
z-1
hM-1
Figura 3:Estructura del Ecualizador Toma Tomada da de Mul Mul rew8 rew88 8
La idea del montaje adaptativo adaptativo (Figura 4) es que funcione de la siguiente siguiente forma. Habrá un peri period odo o inic inicia iall de entr entren enam amie ient nto o en el cual cual se tran transm smititir irá á una una secu secuen enci cia a prede predete term rmin inad ada a que servi servirá rá para para que que el algor algorititmo mo de ident identifific icaci ación ón de sist sistem emas as
______________________________________________________________________
Ŝ(k-d)
________________________________________________________________ TEAN construya una estimación ĥ de la respuesta al impulso del canal h. A continuación se le pasa pasa esa esa esti estima maci ción ón al Filt Filtro ro de Kalm Kalman an para para que que se inic inicia ialilice ce,, y comi comien enza za la transmisión de datos. El hecho de utilizar una secuencia de entrenamiento supone que la inicialización del Filtro será exacta. El vector de estado en el instante cero estará constituido por las d últimas muestras de la secuencia de entrenamiento; la matriz de covarianza de error será la matriz nula puesto que el vector de estado se conoce con probabilidad uno. El Filtro de Kalman proporciona la varianza del error de las estimaciones en la diagonal principal de la matriz de covarianza de error. Una vez que el Filtro de Kalman ha alcanzado su estado estacionario, se usará como entrada para el decisor el elemento del vector de estado que alcance cierto umbral según el criterio que se haya decidido utilizar. n(k) ECUALIZADOR DE KALMAN s(k)
X(k) FIR
ŝ(k-d)
K
m(k-d) DECISOR
— CANAL z --1 z --d ĥ
a
Entrenamiento
IDENTIFICACIÓN DE SISTEMA —
FILTRADO ADAPTATIVO LMS Entrenamiento
Durante la transmisión de datos se puede utilizar el ecualizador adaptativo en la
Figura 4:Montaje del Ecualizador Adaptativo modalidad de decisión directa (la otra Toma Tomada posición da de Mul Mulde rew8 rew88 los8 conmutadores); posición en la
cual se usará la salida del decisor como entrada para el algoritmo de identificación de
______________________________________________________________________
________________________________________________________________ TEAN sistemas. Con esta forma de funcionamiento se supone que si el decisor toma las decisiones correctas, el ecualizador será capaz de seguir las variaciones temporales del canal. En esas condiciones la salida del decisor representará la entrada del canal con un retardo de d muestras. De forma que para que el algoritmo de identificación de sistemas funcione bien hay que meter un retardo de d muestras entre la salida del canal y la entrada del algoritmo de identificación de sistemas. Estos dibujos y esas ecuaciones son las que uno se encuentra en los libros de filtrado adaptativo; y es lo que se mete en la máquina y funciona. La idea ahora es ir a los artículos de R. E. Kalman donde se presentan los conceptos en los que se basan estas estructuras y modelos...
Este hombre, R. E. Kalman, se propone en su artículo de 1960 abordar el problema del filtrado y la predicción lineal usando la representación de Bode-Shannon de procesos estocásticos y el método de transición de estados para el análisis del modelo de
______________________________________________________________________
________________________________________________________________ TEAN sistema dinámico que propone como solución al mismo. De hecho la aportación más importante importante que hace este autor autor es la formulación formulación del problema problema de Wiener desde un punto de vista diferente al que aportaba la teoría clásica de análisis de sistemas: un planteamiento en términos del espacio de estados; planteando también la dualidad entre el problema de la estimación lineal óptima y el problema de control determinístico del regulador sin ruido. Segú Según n el plan plante team amie ient nto o orig origin inal al que que plan plante tea a el auto autorr, se trat trata a de obte obtene nerr la especificación de un sistema lineal dinámico para la predicción, separación o detección de una señal aleatoria. El artículo presenta un nuevo enfoque del problema de Wiener que ayuda a superar algu alguna nass limi limita taci cione oness de los los méto método doss dispo disponi nibl bles es en la época época que que coar coarta taban ban su potencialidad y utilidad práctica en aplicaciones de ingeniería. Hay una serie de suposiciones suposiciones básicas de común acuerdo que están debajo de toda la exposición de Kalman y el desarrollo de su modelo: • Todos los cálculos que se realizan y los resultados que se presentan se apoyan
en el conocimiento a priori de los estadísticos de primer y segundo orden. aleatorias arbitrarias se representan como la salida salida de un sistema sistema • Las señales aleatorias lineal dinámico excitado por señales aleatorias incorreladas (la típica suposición que se suele adoptar adoptar en el estudio de sistemas sistemas lineales que asume la existencia existencia de ruid ruido o blanc blanco o gaus gausia iano no adit aditivo ivo a la entr entrad ada a del del sist sistem ema a provo provoca cand ndo o un determinado error en las observaciones). Es en este sentido donde el autor ya difiere del enfoque clásico para la descripción de los sistemas: Kalman propondrá un sistema de ecuaciones en diferencias (o diferenciales en el caso contínuo) de primer orden). El modelo se basa en los conceptos de ESTADO del sistema, y TRANSICIÓN de estado. • La solución del problema de Wiener nos proporcionará los coeficientes del filtro
a partir de la matriz de covarianza. • Prácticamente la totalidad de su trabajo trata con sistemas discretos: que en los
térmi términos nos en que que los los vamo vamoss a trat tratar ar se pued pueden en cons consid ider erar ar como como sist sistem emas as contínuos que han sido discretizados mediante muestreo.
______________________________________________________________________
________________________________________________________________ TEAN El autor formula sus resultados en forma de teoremas expresados en un lenguaje matemático más o menos formal; utilizando una notación, que en mi opinión no es la más afortunada que se le podía haber ocurrido. De esos resultados que aparecen en los artículos, puede que aparezca alguno a lo largo del documento, pero en la lectura que yo he hecho me ha parecido más destacable el planteamiento que propone del problema de Wiener, los conceptos que subyacen a todo el aparato estadístico que soporta sus razonamientos, y el alcance de los mismos; no me gustaría que después de leer estos artículos sólo pudiera extraer un par de algoritmos numéricos muy eficientes para meter en el MatLab, o cualquier otra herramienta al uso, y jugar con ellos para sacar un montón de gráficas más o menos vistosas donde se ponga de manifiesto, por ejemplo, la diferencia en la velocidad de convergencia de los mismos. Es por esto que nos fijaremos en las bases conceptuales que soportan el razonamiento que sigue Kalman para construir construir su modelo. modelo. En este sentido, sentido, lo primero que podemos podemos destacar son una serie de consideraciones acerca del problema de la construcción de estimaciones óptimas. El problema parte de la suposición suposición de que tenemos tenemos una señal x1(t) a la que se le suma un ruido x2(t); partiremos por tanto de una serie temporal de observaciones y(t) = x 1(t) + x2(t), i.e. supondremos conocidos conocidos los valores de y(t 0), ... , y(t). De lo que se trata es de ver qué podemos inferir sobre los valores de la señal en un instante t 1. Este planteamiento nos lleva directamente a tratar con un problema de estadística y teoría de la probabilidad, que se podría enunciar en los siguientes términos: Dado un conjunto conjunto de observaciones observaciones de y(t) (una variable variable aleatoria aleatoria V.A.), V.A.), η(t0), ... , η(t) deberíamos poder calcular la probabilidad de la ocurrencia simultánea de varios valores ξ(t) de la V.A. x 1(t1), para lo cual propone echar mano de la función de probabilidad
condicional F(ξ1) = Pr [x1(t1) ≤ ξ1|y(t0)=η(t0), ... , y(t)= η(t)]. La idea es que cualquier estimación será una función de esa distribución, es decir, una función determinística de las V.A. y(t 0), ... ,y(t) que podemos poner como X 1(t1|t).
______________________________________________________________________
________________________________________________________________ TEAN Se definirá una función de pérdidas pérdidas para penalizar penalizar las estimaciones estimaciones incorrectas; incorrectas; y será una función positiva y no decreciente del error cometido en la estimación ε = x1(t1) – X1(t1). El criterio para elegir X 1 será normalmente la minimización de las pérdidas medias, o el “riesgo” en la estimación. E{L[x1(t1) – X1(t1)]} = E[E{L[x 1(t1) – X1(t1)]|y(t0), ... ,y(t)}] = E{L[x 1(t1) – X1(t1)]|y(t0), ... ,y(t)}
Como concepto importante que soporta el razonamiento que sigue Kalman a la hora de atacar el problema de la estimación lineal, aparece la proyección ortogonal. La idea es que los resultados resultados obtenidos obtenidos mediante estimación estimación lineal se pueden mejorar mejorar con una estimación no lineal sólo cuando los procesos aleatorios son no gausianos y además considerando funciones de distribución de tercer orden. Para Para justif justifica icarse, rse, Kalman Kalman busca busca una interpr interpreta etació ción n geomét geométric rica, a, y nos plante plantea a las siguientes consideraciones. Si tenemos las V.A. y(t 0), ... ,y(t) el conjunto de todas las combinacione combinacioness lineales lineales ∑(i=t0,...,t)aiy(i) formará un espacio vectorial (variedad lineal) ℒ(t), que se puede considerar como un subespacio de dimensión finita del espacio
formado por todas las posibles observaciones. En esas condiciones, dados dos vectores u,v∈ℒ(t) diremos que son ortogonales sii Euv = 0. Usando el método de Schmidt podremos seleccionar una base ortonormal de ℒ(t) {et0, ... ,e t} ; Eeie j = δij = 1 si i = j con i,j = t 0, ... ,t Χ = ∑(i=t0,...,t)aiei, ∀ Χ∈ℒ(t)
Cualquier V.A. se podrá descomponer de forma unívoca en una parte Χ que caerá en ℒ(t), y otra parte X’ que será ortogonal a ℒ(t)
X = Χ + X’ = ∑(i=t0, ... ,t) (Exei)ei + X’ Con lo cual podremos escribir E X’ei = E(x – Χ )e )ei = Exei - E Χ ei = 0 :(i=t 0, ... ,t) Puesto que las coordenadas son únicas. Se dice entonces que Χ es la proyección ortogonal de x en ℒ(t).
______________________________________________________________________
________________________________________________________________ TEAN A Χ se le puede ver como el vector de ℒ(t) (es decir una función lineal de las V.A. y(t0), ... ,y(t)) que minimiza la función de péridas cuadrática. Eso quiere decir que si cogemos cualquier otro vector de ℒ(t), w, tendremos E(x-w)2 = E(x’ + x – w) 2 = E[(x –x) + (x – w)] 2 Como resulta que x’ es ortogonal a todos los vectores de ℒ(t), y en particular particular a (x – w), tendremos que E(x-w)2 = E(x –x)2 + (x – w) 2 ≥ E(x – x)2 Si además w minimiza minimiza la función cuadrática de pérdidas, E(x – w)2 = 0, y eso significa que x = w Lo que viene a decir este razonamiento es que si tenemos {x(t)} e {y(t)} procesos aleatorios de media nula, y un conjunto de observaciones y(t0), ... ,y(t) podemos decir que si ocurre que los procesos son gausianos, o bien que la estimación óptima se restringe a una función lineal de las V.A. observadas y la función de pérdidas L( ε) = ε2 entonces entonces x*(t1|t) la estimación óptima de x(t 1) dados y(t0), ... ,y(t) coincide con x(t 1|t) la proyección ortogonal de x(t1) en ℒ(t). Esto lo expresa Kalman de una forma más compacta compacta como x(t1|t) = x*(t1|t) = Ê[x(t1)| ℒ(t)] ; que también se puede interpretar interpretar como que en las condiciones mencionadas la mejor estimación que se puede hacer es una combinación lineal de las observaciones previas. Estamos diciendo que el estimador óptimo es un filtro al que le metemos los valores reales de las V.A. V.A. observables. La aportación importante que hace Kalman a la resolución del problema del filtrado lineal es su propuesta propuesta de un modelo para los procesos procesos estocásticos estocásticos que le servirá para introducir un nuevo enfoque en el e l Análisis de Sistemas. La aprox aproxim imac ació ión n habi habitu tual al part parte e de la consi conside dera ració ción n de los los proce procesos sos alea aleato tori rios os macroscópicos como procesos gausianos independientes. Cuando tal aproximación no es aplicable se suele explicar la dependencia estadística entre las observaciones de señales aleatorias en instantes temporales distintos con la suposición de la existencia de un sistema dinámico entre la fuente y el observador. Se puede ver por tanto una señal aleatoria como la salida de un sistema dinámico excitado por una serie de procesos estocásticos gausianos independientes.
______________________________________________________________________
________________________________________________________________ TEAN Con esto se trata de aprovechar la buena costumbre que tienen las señales aleatorias gausianas de seguir siendo gausianas al pasar por un sistema lineal. Si hemos supuesto que las fuentes son gausianas e independientes, y las señales aleatorias que forman el conjunto de observaciones son también gausianas, podremos suponer que el sistema dinámico que hay entre la fuente y el observador es LINEAL. Adem Además ás sabe sabemo moss que que dado dado un proc proces eso o esto estocá cást stic ico o del del que que se cono conoce cen n sus sus estadísticos de primer y segundo orden, podemos encontrar un proceso gausiano con las mismas propiedades. Con estas ideas iniciales, Kalman se propone buscar la forma de describir el sistema en base al concepto de ESTADO . El autor entiende por estado un conjunto de información cuantitativa que constituye la cantidad mínima de datos que necesitamos conocer acerca del comportamiento del sistema en el pasado para poder predecir su futuro. La dinámica del sistema quedará descrita en términos de las TRANSICIONES entre estados. El sistema quedará descrito por un modelo como el siguiente: dx/dt = F(t)x(t) + D(t)u(t) y(t) = M(t)x(t) donde x =[x1 ... xn]T es el estado del sistema, y además u(t)∈ℳm×1 ; m≤n F(t)∈ℳn×n ; D(t)∈ℳn×m y(t)∈ℳp×1 ; M(t)∈ℳn×p ; p≤n Es un sistema en tiempo contínuo. Si suponemos que es un sistema estacionario (es decir que todos los coeficientes de F, D y M son constantes: un sistema LTI), y que u(t) es constante en el periodo de muestreo, u(t+ τ) = u(t) ; 0 ≤τ≤1 , t=0,1, t=0,1, ... podemos podemos transformarlo en un sistema discreto. Kalman lo expresa así: X(t+1)=Φ(1)x(t) + ∆(1)u(t) ; t=0,1,... Φ(1) = exp F = ∑(i=0..∞)(Fi/i!) ; F 0=I
∆(1) = ( ∫0,1expFτdτ)D
______________________________________________________________________
________________________________________________________________ TEAN Si no es estacionario, podremos ponerlo como X(t+1)=Φ(t+1;t)x(t) + ∆(t)u(t) ; t=0,1,... y(t) = M(t)x(t) Φ(t+1;t) es la MATRIZ DE TRANSICIÓN, y es tal que: Φ(t2;t1) indica la transición desde t1 hacia t2 Φ(t;t) = I Φ(t;s) Φ(s;r) = Φ(t;r) Φ-1(t;s) = Φ(s;t)
Si el sistema es LTI, Φ(t+1;t) = Φ(t+1-t) = Φ(1) = Constante y Φ(t;τ) = exp F(t- τ) Por otro lado, {u(t)} será un proceso vectorial, independiente, y gausiano de media nula: Eu(t) = 0, ∀t Eu(t)u’(z) = 0 si t≠z Eu(t)u’(t) = Q(t) {x(t)} también será un proceso gausiano de media nula, pero no independiente tal que x(t) = ∑(r=-∞ .. t-1) Φ(t;r+1)u(r) Si t≥s : Ex(t)x’(t) = ∑(r=-∞ .. s-1) Φ(t;r+1)Q(r) Φ’(s;r+1) Partiendo del modelo de sistema hemos llegado a la matriz de covarianzas Ex(t)x’(t); sin embargo lo normal en los problemas del mundo real es que intentemos llegar a la matriz de covarianzas partiendo de las observaciones, para luego construir el modelo a partir de la matriz. Con este modelo de sistema se puede ver el problema de Wiener de forma natural. Considerando el siguiente sistema dinámico: X(t+1) = Φ(t+1;t)x(t) + u(t) Y(t) = M(t)x(t) Donde u(t) es un proceso gausiano independiente vectorial (vectores de n componentes), y de media nula; x(t) es e s un vector de n componentes; y(t) y( t) tiene p componentes (p≤n); Φ(t+1;t)∈ℳn×n ; M(t)∈ℳp×n siendo sus elementos funciones determinísticas del tiempo (señales).
______________________________________________________________________
________________________________________________________________ TEAN Dadas las observaciones y(t0), ... ,y(t) Encontrar un estimador x*(ti|t) |t) de x(t x(t1) que minimice la esperanza de la función de pérdidas. Sabemos que la solución será la proyección ortogonal de x(t 1) sobre la variedad lineal generada por las V.A. observadas. Al final, jugando con la variedad lineal generada por todas las observaciones y la proyección ortogonal, Kalman propone la siguiente expresión para el Filtro Óptimo: X*(t+1) = Φ*(t+1;t)x*(t|t-1) + ∆*(t)y(t) Φ*(t+1;t) = Φ(t+1;t) - ∆*(t)M(t)
cuyo diagrama de bloques podemos ver en la Figura 5. x*(t+s|t) (t+s;t+1) y(t)
x*(t|t-1)
(t|t-1) ŷ (t|t-1)
x*(t+1|t)
d
*(t)
y(t|t-1) M(t)
(t+1;t) x*(t+1|t-1)
-I Figura 5:Diagrama de Bloques del Filtro Óptimo de Kalman
En el artículo artículo de Kalman, Kalman, todo esto que caracterizado caracterizado en forma de Teorema, Teorema, donde se recogen las relaciones de recursión que constituyen el algoritmo de resolución del problema de Wiener:
______________________________________________________________________
________________________________________________________________ TEAN La estimación óptima x*(t+1|t) de x(t+1) dados y(t0), ,,, ,y(t) se genera con el sistema lineal dinámico x*(t+1|t) = Φ*(t+1;t)x*(t|t-1) + ∆*(t)y(t) El error de la estimación viene dado por x(t+1|t) = Φ*(t+1;t) x(t|t-1) + u(t) La matriz de covarianza del error es P*(t) = Ex(t|t-1)x(t|t-1) La pérdida en media cuadrática es ∑(i=1, ... ,n) Exi2(t|t-1) = traza {P*(t)}
Las matrices ∆*(t), Φ*(t+1;t), y P*(t) siguen las siguientes relaciones de recursión ∆*(t) = Φ(t+1;t)P*(t)M’(t)[M(t)P*(t)M’(t)]-1 Φ*(t+1;t) = Φ(t+1;t) - ∆*(t)M(t)
P*(t+1) = Φ*(t+1;t)P*(t)Φ’(t+1;t) + Q(t) Para obtener las iteraciones debemos especificar las covarianzas P*(t0) de x(t0) y Q(t) de u(t) La generalización para cualquier s≥0, si Φ es inversible, vendrá dada por x*(t+s|t) = Φ(t+s;t+1)x*(t+1|t) = Φ(t+s;t+1)Φ*(t+1;t)Φ(t;t+s-1)x*(t+s-1|t-1) + Φ(t+s;t+1)∆*(t)y(t) Si eliminamos ∆* y Φ* de las relaciones de recursión, obtendremos una ecuación en diferencias no lineal P*(t+1) = Φ(t+1;t){P*(t) – P*(t)M’(t)[M(t)P*(t)M’(t)] -1P*(t)M(t)}Φ’(t+1;t) + Q(t) ; t ≥t0
Es interesante destacar en una primera lectura del artículo el planteamiento que hace el autor del problema de control determinístico mostrando su dualidad con el de la estimación lineal óptima. Considerando el sistema dinámico x(t+1) =
(t+1;t)x(t) + M(t)u(t)
______________________________________________________________________
________________________________________________________________ TEAN donde x(t)∈ℳn×1 ; u(t)∈ℳm×1 (t)∈ℳn×n ; M(t)∈ℳn×m ; m≤n siendo sus elementos funciones no aleatorias del tiempo El problema que se plantea es el siguiente: sig uiente: Dado cualquier estado x(t) en un instante t, encontrar la secuencia u(t), ... ,u(T) de vectores de control que minimicen el índice de rendimiento V[x(t)] = ∑(τ=t .. T+1) x’(τ)Q(τ)x(τ) Dond Donde e Q(t) Q(t) es una una matri atrizz defi defini nida da posi posittiva iva cuyo cuyoss elem elemen ento toss son seña señalles determinísticas. En este problema NO HAY CONSIDERACIONES ESTADÍSTICAS. Kalman propone propone que la dualidad entre los problemas de estimación estimación y control control se puede ver en el siguiente sentido: Sea τ≥0. Si en las relaciones del problema de estimación sustituimos cualquier matriz de la forma X(t) = X(t0+τ) por X’(t) = X’(T-τ) obte obtend ndre remo moss las las rela relaci cion ones es correspondientes al problema de control. Y a la inversa, si sustituimos cada matriz X(Tlas rela relaci cion ones es de cont contro roll por por X’(t X’(t0+τ), obte obtend ndrem remos os las las relaci relacione oness para para la τ) en las estimación lineal óptima. Matemáticamente lo que significa la dualidad es que ambos problemas se reducen a la resolución de una ecuación no lineal del tipo de la de Wiener-Hopf, que es la que surgía de la eliminación en las relaciones de recursión de dos de las matrices. En la Figur Figura a 6 se puede puede ver ver una una repre represe sent ntac ació ión n esqu esquem emát átic ica a del prob proble lema ma del del regulador sin ruido tomada del artículo del propio Kalman Sistema Físico a controlar
x(t+1)
U*(t) M(t)
d
*(t)
(t+1;t)
-I ______________________________________________________________________ Figura 6:Diagrama de Bloques del Regulador Sin Ruido
________________________________________________________________ TEAN
Teniendo en cuenta que U*(t) = - *(t)x(t) Las ecuaciones para el control óptimo son x(t+1) =
*(t+1;t)x(t)
V*[x(t)] = x’(t)P*(t-1)x(t) Las relaciones de recursión quedan así
*(t) = [M’(t)P*(t) M(t)]-1 M’(t) P*(t) (t+1;t) *(t+1;t) = (t+1;t) - M(t) *(t) P*(t+1) = ’(t+1;t) P*(t) *(t+1;t) + Q(t) Para t≤T Siendo inicialmente P*(T) = Q(T+1)
Conclusiones Los dos artículos de Kalman en los que se basa este documento son unos clásicos dentro del vasto campo del filtrado y la predicción lineal; constituyeron la base para el desarrollo de un serie de algoritmos de resolución del problema de Wiener como los
______________________________________________________________________
________________________________________________________________ TEAN que se mencionaron más arriba y que son y han sido de amplia aplicación en los campos de las comunicaciones y el control. Resulta clarificador estudiar el fundamento conceptual que hay detrás de los algoritmos que metemos en las máquinas a la hora de resolver resolver los problemas problemas reales; en este caso resulta interesante interesante ver el jugo que se le ha sacado a una forma intuitiva de representación del problema en base al espacio de estados y la matriz de transición; el cómo tener un adecuado modelo de los procesos estoc estocást ástic icos os ayuda ayuda a lleg llegar ar al coraz corazón ón del del probl problem ema a de estim estimaci ación ón dándo dándole le una una interpretación geométrica a partir de la proyección ortogonal, que también da pie a plantear la Ecuación Integral de Wiener-Hopf. Sobre todo resulta estimulante pensar en la idea de la dualidad entre dos problemas que a priori parecen no tener mucho que ver: un problema, el de la ESTIMACIÓN (ya sea interpolación, filtrado o predicción), cubierto por un aparato estadístico que puede conseguir desviar nuestra atención del fondo de la cuestión; y otro, el del CONTROL dete determ rmin inís ístitico co de una una plant planta a sin sin ruid ruido o donde donde no se hace hace ningu ninguna na cons consid ider eraci ación ón probabilística. Una dualidad que matemáticamente nos lleva a la ecuación de WienerHopf, y que no encontraba significado físico bajo el enfoque de Kalman, que había estudiado ambos problemas por separado.
Bibliografía Para la elaboración de este documento se han consultado, aparte de los apuntes de la asignatura, las siguientes fuentes
______________________________________________________________________
________________________________________________________________ TEAN 1. Mulgrew Mulgrew,, Bernard Bernard and F.N. Cowan, Cowan, Colin Colin “Adapt “Adaptive ive Filters Filters and Equaliser Equalisers” s” Kluwer Kluwer Academic Publishers USA, 1988 2. Wien Wiener er,, Norb Norber ertt “CIB “CIBER ERNÉ NÉTI TICA CA o el cont contro roll y comu comuni nica caci ción ón en anim animal ales es y máquinas” 2ª edición Tusquets Editores, Barcelona 1998 3. R. E. Kalman Kalman “A New Approac Approach h to Linear Linear Filter Filtering ing and Predic Predictio tion n Proble Problems” ms” Trans. Trans. ASME, Journal of Basic Engineering pp. 35-45, (March 1960) 4. R. E. Kalman Kalman and R. S. Bucy Bucy “ New Results Results in Linear Filtering Filtering and Predicti Prediction on Theory” Theory” Trans. ASME, Journal of Basic Engineering pp. 95-108, (March 1961) 5. R. E. Lawrence Lawrence and Kaufman Kaufman “The Kalman Kalman Filter Filter for the Equaliz Equalizati ation on of a Digita Digitall Communications Channel” IEEE Transactions on Communications Technology vol. 19(6) pp. 1137-1141, 137-1141, (Dec. 1971)
______________________________________________________________________
________________________________________________________________ TEAN
Temas Especiales de Análisis Numérico José Manuel Corrales Sendino Curso 1999/2000 Primer Semestre
______________________________________________________________________
