• Author / Uploaded
• api-26068058

• Categories
• Estocástico
• Matriz (Matemáticas)
• Espacio vectorial
• Filtro de Kalman
• Linealidad

Citation preview

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 and Communication in the Animal and the 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 primera como asignatura de libre elección en el P94 durante el segundo 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 concretos. Este hecho, quizás motivado por la lejanía de los primeros años de universidad en que uno se acostumbraba a lidiar en el día a día con el utillaje y la operativa matemática, podría llevar a pensar que estas líneas no constituyen otra cosa que un de ejercicio de demagogia 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 de un instante temporal al siguiente (de una muestra a la siguiente) mediante un algoritmo adaptativo 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 pasadas del 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 estructura básica que se representa en la Figura 1 tiene tres modos 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 la salida del filtro 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). En esa situación se dirá que hemos identificado el sistema. Normalmente tendremos ruido blanco en la entrada (errores en las observaciones) que nos alejará 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 comunicación. En este caso se trata de que la salida del sistema 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 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 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 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)]. x(n)

Sistema Desconocido

Filtro Adaptativo

y(n)

- Σ+

Figura 2(i):Identificación de Sistemas (Tomada de [Mulgrew88])

El filtro FIR óptimo se obtendría siguiendo estas consideraciones: ŷ(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 ... hN-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φ xxh – 2hTφ xy φ xx = 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φ xxh – 2φ xy = 0 ⇨ φ xxhopt = φ 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 hopt = φ -1xx φ 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 ocurre aquí es que necesitamos la matriz de autocorrelación y el vector 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(N3) 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 ⇨ φ xx[i,j] = φ xx[j,i]

______________________________________________________________________

________________________________________________________________TEAN y es de Toeplitz⇨ φ xx[i,j] = φ xx[i-j]

(ii)

siendo φ xx[i,j] el elemento en la fila i, columna j de la matriz de autocorrelación, y {φ xx(m)} 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 operador esperanza matemática de la función de coste MSE por un sumatorio, 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 rxx(k)h(k) = rxy(k) rxx (k) = Σ(n=0,...,k)x(n)xT(n) ∈ ℳN×N es la matriz de autocorrelación rxy (k) = Σ(n=0,...,k)x(n)y(n) ∈ ℳN×1 es el vector de correlación cruzada Existirá solución única si existe rxx (k) y es no singular. Para k≤N-1, rxx (k) será siempre singular. Para k≥N, la singularidad de rxx (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 sobre los datos disponibles: se puede usar directamente la covarianza, 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, que se conocen como algoritmos rápidos. Surgieron a partir de la publicación del FKA (Fast Kalman Algorithm) que trataba de aprovechar la estructura peculiar de la matriz rxx (k) de la estimación LS; le siguieron más algoritmos 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 que se deben llevar a cabo para encontrar la 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 función de coste MSE; en la segunda fase se construye una nueva estimación restándole a la primera una versión escalada del gradiente µ∇(hi). Queda algo así: hi+1 = hi - µ∇(hi) ∇ i = ∂ξ/∂h(hi) = 2φ xxhi – 2φ xy como condición de convergencia tenemos que 0