TESIS DOCTORAL NUEVAS TÉCNICAS DE DECONVOLUCIÓN DISPERSA ROBUSTA. Autor: CARLOS EUGENIO MARTÍNEZ CRUZ. Director: DR. JOSÉ LUIS ROJO ÁLVAREZ


Save this PDF as:
 WORD  PNG  TXT  JPG

Tamaño: px
Comenzar la demostración a partir de la página:

Download "TESIS DOCTORAL NUEVAS TÉCNICAS DE DECONVOLUCIÓN DISPERSA ROBUSTA. Autor: CARLOS EUGENIO MARTÍNEZ CRUZ. Director: DR. JOSÉ LUIS ROJO ÁLVAREZ"

Transcripción

1 TESIS DOCTORAL NUEVAS TÉCNICAS DE DECONVOLUCIÓN DISPERSA ROBUSTA Autor: CARLOS EUGENIO MARTÍNEZ CRUZ Director: DR. JOSÉ LUIS ROJO ÁLVAREZ Codirector: DR. MANEL MARTÍNEZ RAMÓN Departamento de TEORÍA DE LA SEÑAL Y COMUNICACIONES UNIVERSIDAD CARLOS III DE MADRID Leganés, Junio 27

2

3 Tesis Doctoral: NUEVAS TÉCNICAS DE DECONVOLUCIÓN DISPERSA ROBUSTA Autor: Director: Codirector: CARLOS EUGENIO MARTÍNEZ CRUZ DR. JOSÉ LUIS ROJO ÁLVAREZ DR. MANEL MARTÍNEZ RAMÓN El tribunal nombrado para juzgar la Tesis Doctoral arriba citada, compuesto por los doctores: Presidente: Vocales: Secretario: Acuerda otorgarle la calificación de Leganés, a

4

5 Agradecimientos Al ver en retrospectiva mi vida en estos últimos tres años, y verme en la situación de tener que escribir unos agradecimientos, solo puede surgir un nombre: Eva Salvador, gracias por haber compartido estos años conmigo. Quiero agradecer principalmente a las personas que durante el tiempo que me llevó realizar la presente Tesis estuvieron conmigo en todo momento y en todas las circunstancias: Diego Gómez, Geydy Gutiérrez, María Venegas, Rosa Barrio y Abdel, muchas gracias por estar ahí. También, a la gente que he conocido en el laboratorio 4.2.A.3 muy especialmente a Jesús Requena; y a todos los demás, M a Luz, Felipe, Javier, Oscar, Raúl, y tantos otros, muchas gracias. Así mismo agradezco la amistad de la gente con la que en estos últimos años he compartido muchas vivencias, Carlos Guerrero, Leonardo, Daniela, Hugo, Orlando, Alejandro, Omar, y en general al grupo hispanoamericano de estudiantes de Doctorado de la Carlos III. Quiero extender mi agradecimiento a un antiguo profesor, y más tarde compañero de trabajo en la Escuela de Ingeniería Eléctrica de la Universidad de El Salvador: gracias Salvador German por motivarme y apoyarme a continuar con mis estudios. Finalmente, esta Tesis ha sido finalizada gracias a la voluntad y el empeño de sus directores. Gracias José Luis y gracias Manel. Institucionalmente, el inicio, el desarrollo y la finalización de la presente Tesis ha sido posible gracias al apoyo del Departamento de Teoría de la Señal y Comunicaciones de la Universidad Carlos III de Madrid y del Programa Alβan, Programa de becas de alto nivel de la Unión Europea para América Latina.

6

7 Resumen La deconvolución dispersa consiste en estimar una señal de reducido número de muestras no nulas a partir de un conjunto de observaciones ruidosas, que se suponen originadas por la convolución entre una señal de aspecto impulsivo y la respuesta al impulso de un sistema lineal e invariante. En la introducción de esta investigación se comentan las aplicaciones más relevantes, que van desde la propuesta original en sísmica de reflexión a extensiones en procesamiento de imágenes. La primera de las dos líneas de investigación presentadas en esta Tesis introduce nuevos métodos de deconvolución dispersa no ciega basados en una técnica emergente de procesado de señal a partir de algoritmos de Máquinas de Vectores Soporte (SVM). Estos algoritmos están regularizados, por lo que son adecuados al problema de deconvolución dispersa, el cual está mal condicionado, y pueden proporcionar soluciones dispersas. En esta Tesis se introducen dos algoritmos de deconvolución no ciega utilizando SVM. El primer algoritmo se formula en el funcional SVM primal, y proporciona una solución regularizada no dispersa. El segundo se formula en el dual, y produce una solución dispersa utilizando como función núcleo la función de ambigüedad del sistema lineal invariante. Se utiliza una función de coste robusta y su grado de dispersión se ajusta a través de uno de sus parámetros libres. La segunda línea parte de la modificación del algoritmo de deconvolución homomórfica (DH), a través de la introducción de un método constructivo de transformación de la señal que permite analizar, mediante representaciones cepstrales, señales limitadas en banda. Las prestaciones de los métodos de deconvolución dispersa propuestos se analizan mediante simulaciones sobre datos sintéticos y reales. El algoritmo basado en Máquinas de Vectores soporte con función núcleo de función de ambigüedad se prueba en datos reales obtenidos en un experimento en sísmica de reflexión. Así mismo, el algoritmo de deconvolución homomórfica propuesto se prueba en un problema práctico de hemodinámica.

8

9 Abstract Sparse deconvolution is the estimation of a low nonzero samples signal from a set of noisy observations. The observations are supposed to be generated as the convolution between a sparse signal and a linear invariant system s impulse response. In this research s introduction we comment some relevant applications, from reflection seismology original applications to image processing. The first of the two research lines presented in this Thesis introduces new nonblind sparse deconvolution methods. They are based on an emerging signal processing techniques called Support Vector Machine (SVM). SVM algorithms are regularized and produce sparse solutions fitting sparse deconvolution requirements. In this Thesis we introduce two nonblind sparse deconvolution algorithms based on SVM. The first algorithm is formulated in the SVM primal functional and produces a regularized non sparse solution. The second algorithm is formulated in the dual signal model and produces a sparse solution when the kernel function is the linear time invariant system s ambiguity function. A robust cost function is utilized which allows to tune the sparsity by one of its free parameter. The second research line uses a homomorphic deconvolution (DH) modification through the introduction of a constructive method, which transforms band limited signal into full band signals. Performance of the sparse deconvolution methods are analyzed through simulations on synthetic and real data. Support Vector Machine algorithm is tested with reflection seismic real data. Also, homomorphic deconvolution algorithm is tested in a hemodynamic practical problem.

10

11 Nomenclatura y(t): señal de salida analógica x(t): señal de entrada analógica h(t): respuesta al impulso e(t): señal de error y n : señal de salida de tiempo discreto x n : señal de entrada de tiempo discreto h n : respuesta al impulso de tiempo discreto e n : señal de error muestreada P: número de muestras de la señal de entrada Q: número de muestras de la respuesta al impulso N: número de observaciones H: matriz de convolución L: lagrangiano λ: parámetro libre en función de coste L ε: parámetro de insensibilidad en función de coste ε-huber γ: parámetro de coste cuadrático en función de coste ε-huber C: parámetro de coste lineal en función de coste ε-huber ξ ( ) k : variables de holgura β ( ) k : multiplicadores de Lagrange asociados a variables de holgura α ( ) k : multiplicadores de Lagrange asociados a vectores soporte K: matriz de núcleos w: vector de pesos φ(x) : X F transformación no lineal : operador de convolución de tiempo discreto R: conjunto de los números reales R n : espacio n-dimensional de los reales X : espacio de entrada H: espacio de características () T : operador de trasnposición : norma 2 : norma 2, : producto interno

12

13 Acrónimos ARMA: BW: BG: CMPG: CSPG: db: DFT: DD: DH: DP: FFT: FIR: Hz: IIR: KKT: LS: MAP: MED: MG: MDS: ML: MPS: MSKA: MSE: PSF: PL: RKHS: SLIT: SNR: SVM: SVR: TF: TIF: WT: WTVA: Modelo Autoregresivo de Medias Móviles (Auto Regressive Moving Average) Ancho de Banda (Bandwidth) Bernoulli-Gausiano Conjunto de Trazas con Punto Medio Común (Common Midpoint Gather) Conjunto de Trazas de Disparo Común (Common Shotpoint Gather) Decibelio Transformada Discreta de Fourier Deconvolución Dispersa Deconvolución Homomórfica Deconvolución Predictiva Transformada Rápida de Fourier (Fast Fourier Transform) Respuesta al Impulso Finita (Finite Impulse Response) Hercios Respuesta al Impulso Infinita (Infinite Impulse Response) Karush-Khun-Tucker Mínimos Cuadrados (Least Squares) Máximo a posteriori Deconvolución de mínima entropía Mezcla de Gausianas Modelo Dual de Señal Máxima Verosimilitud Modelo Primal de Señal Modelo de Señal con Kernel de Ambigüedad o Autocorrelación Error Cuadrático Medio (Mean Square Error) Función Punto-Dispersión (Point Spread Function) Problema Lineal Espacio de Hilbert con Núcleo Reproductor (Reproducive Kernel Hilbert Space) Sistema Lineal Invariante en el Tiempo Relación Señal a Ruido (Signal Noise Ratio) Máquina de Vectores Soporte (Support Vector Machine) Máquina de Vectores Soporte en Regresión (Support Vector Regression) Transformada de Fourier Transformada Inversa de Fourier Trazo de Oscilación (Wavelet Trace) Trazo de Oscilación Temporal de Área Variable (Wavelet Trace Variable Area)

14

15 Índice general. Introducción.. Interés y aplicaciones de las técnicas de deconvolución La deconvolución en problemas de señal Motivación e interés de la investigación Planteamiento y objetivos de la Tesis Deconvolución de señales dispersas DD no ciega Solución usando la transformada de Fourier Solución de mínimo error cuadrático Solución mediante filtro de Wiener Deconvolución predictiva Deconvolución como problema de optimización Deconvolución norma L Deconvolución con coste Huber Deconvolución de máxima verosimilitud Deconvolución mediante mezcla de Gausianas Método de la ventana iterativa DD ciega Técnicas de DH Algoritmo de módulos constantes Métodos basados en aprendizaje máquina Consideraciones sobre la señal dispersa Consideraciones sobre el ruido a la salida Consideraciones sobre la fase Conclusiones Deconvolución de señales dispersas mediante SVM SVM en clasificación y regresión El clasificador binario con clases linealmente separables El problema de regresión lineal Núcleos de Mercer Ejemplos de funciones núcleo Núcleos invariantes al desplazamiento Núcleos de Mercer en clasificación SVM Núcleos de Mercer en regresión SVM xiii

16 XIV ÍNDICE GENERAL 3.2. Marco lineal SVM para señales de tiempo discreto Análisis espectral SVM Identificación de sistemas SVM-ARMA SVM para interpolación Algoritmo MPS para interpolación Interpolación no uniforme usando espacios invariantes Algoritmo MDS para interpolación SVM para deconvolución no ciega Formulación de la DD en el MPS Formulación de la DD en el MDS Formulación de la DD en el MDS modificado Simulaciones y Experimentos Experimento : parámetros libres de la SVM Experimento 2: señal de entrada determinista Experimento 3: señal de entrada aleatoria Experimento 4: efecto del ruido impulsivo Experimento 5: efecto de la fase no mínima Experimento 6: efecto del ruido coloreado Experimento 7: aplicación práctica Conclusiones Deconvolución homomórfica Generalidades del cepstrum complejo Definición del CC Generalización del principio de superposición Sistema característico para la convolución Deconvolución homomórfica de señales DH de señales limitadas en banda DH como identificador de sistemas DH en el Dominio Suplemental Simulaciones y Experimentos MD de DH Señales limitadas en banda Señales periódicas limitadas en banda Señales periódicas limitadas en banda con reflexiones Conclusiones Caracterización de la Impedancia de Entrada Aórtica mediante DH El Sistema Cardiovascular Ciclo cardíaco Hemodinámica Presión y flujo sanguíneo y teoría de propagación de ondas en el sistema arterial Modelos de propagación de ondas en el sistema arterial Modelo Windkessel

17 ÍNDICE GENERAL XV Modelo elástico tubular Modelos en el dominio de la frecuencia Método basado en series de Fourier Método basado en modelado ARMA Método híbrido basado en modelos históricos Métodos temporales Limitaciones de los modelos Simulaciones y Experimentos Estimación de la IEA y la AEA usando análisis cepstral Estimación de la IEA y la AEA en el DS Conclusiones Conclusiones y líneas futuras Aportaciones Líneas Futuras Apéndices 75 A. Condiciones KKT 77 B. DH utilizando raíces polinómicas 79 C. Análisis cepstral de un tren de impulsos general 83

18 XVI ÍNDICE GENERAL

19 Lista de Figuras 2.. El problema de sísmica de reflexión: (a) el subsuelo terrestre modelado como un conjunto de estructuras geológicas estimulado por una onda sísmica; (b) su representación como un SLIT Señales de prueba: (a) respuesta al impulso; (b) señal dispersa; (c) amplitud del espectro de la respuesta al impulso; (d) amplitud del espectro de la señal dispersa; (e) amplitud del espectro de la señal de salida con ruido; (f) señal dispersa estimada Filtrado de Wiener en DP: (a) diseño del filtro; (b) aplicación a un conjunto de observaciones Funciones de coste utilizadas en la deconvolución como problema de optimización: (a) función coste valor absoluto (norma L ); (b) función coste de Huber Modelo convolucional utilizado en la formulación del problema de deconvolución de ML. La entrada esta forma de la suma de dos señales aleatorias, una de estadística BG y la otra de estadística Gausiana Esquema que ilustra la deconvolución con regularización sobre el ruido. El filtro de Fourier atenúa el ruido amplificado durante la operación de inversión y el filtro Wavelet atenúa el ruido residual El problema de clasificación binario de clases linealmente separables: (a) hiperplano separador (w,b); (b) margen de una muestra al hiperplano; (c) margen máximo del hiperplano Función de regresión lineal unidimensional Efecto de la traslación φ(x, x 2 ) = (x 2, 2x x 2, x 2 2) de un espacio de entrada X a un espacio de características H Hiperplano separador calculado mediante SVM para: (a) clases linealmente separables; (b) clases linealmente separables con muestras en el margen; (c) clases no separables linealmente.. 45 xvii

20 XVIII LISTA DE FIGURAS 3.5. Esquema de regresión (figura adaptada de [8]) SVR. Las muestras en el espacio original se trasladan a un espacio de características donde es posible aplicar la SVR lineal. Las muestras fuera del intervalo definido por ε son penalizadas y se convierten en vectores soporte. La penalización se lleva a cabo a través de una función de coste: (a) función ε-insensible; (b) función ε-huber El problema de muestreo, la localización del instante de muestreo se ha marcado con una X y el valor de la muestra con un círculo: (a) una función x(t) R ha sido muestreada uniformemente y no uniformemente, paneles superior e inferior, respectivamente; (b) dos tipos de muestreo bidimensional (cartesiano y polar) con muestreo uniforme y no uniforme, paneles superior e inferior, respectivamente Representación esquemática en tres pasos del teorema de muestreo utilizando espacios invariantes al desplazamiento. El intervalo de muestreo es uniforme con T =. () La señal de entrada x(t) se ha filtrado con un filtro antisolapamiento. (2) El proceso de muestreo conduce a la obtención de los coeficientes w k. (3) La reconstrucción se obtiene a través del filtrado de los coeficientes con ϕ(t) Esquema para el MPS. Los multiplicadores de Lagrange y las observaciones son la entrada y la salida, respectivamente, de un SLIT Esquema para el MDS, donde se representa la relación entre las observaciones y los multiplicadores de Lagrange Modelo modificado dual de señal: (a) modificación de las observaciones; (b) relación entre las observaciones y los multiplicadores de Lagrange Señales de prueba: (a) señal dispersa; (b) respuesta al impulso; (c) señal de salida sin ruido; (d) amplitud de la respuesta en frecuencia del filtro; (e) fase de la respuesta en frecuencia del filtro Variación del MSE para los métodos de deconvolución basados en los modelos MSKA y MPS en función de los parámetros libres C y γ, para: SNR = 4 db, (a)-(b) y SNR = 6 db (c)-(d) MSE en función de la SNR para los métodos: MSKA, MPS, MG, L (trazo continuo, discontinuo, punteado y punteado discontinuo, respectivamente) en presencia de ruido: (a) Gausiano; (b) Laplaciano; (c) uniforme F ssd en función del parámetro libre para SNR = 4, 8, 2 y 2 db (trazo continuo, discontinuo, punteado y punteado discontinuo, respectivamente ) en: (a) MPS; (b) MSKA; (c) MG; (d) L

21 LISTA DE FIGURAS XIX 3.5. Figura de mérito F pick (trazo continuo) y F null (trazo discontinuo) para SNR = 6 db en: (a) MPS; (b) MSKA; (c) MG; (d) L Promedio de la probabilidad de detección para la señal determinista en: (a) MPS; (b) MSKA; (c) MG; (d) L Sensibilidad (trazo continuo) y especificidad (trazo discontinuo) para una SNR = 4dB en: (a) MPS; (b) MSKA; (c) MG; (d) L MSE en función de SNR para los métodos: MSKA; MPS; MG; y L (trazo continuo, discontinuo, punteado y punteado discontinuo, respectivamente) en presencia de ruido Gausiano y entrada aleatoria con distribución BG ( σ a =, p =.5) F ssd para una entrada aleatoria, con distribución BG, en función del parámetro libre para SNR = 4, 8, 2 y 2 db (trazo continuo, discontinuo, punteado y punteado discontinuo, respectivamente) en: (a) MPS; (b) MSKA; (c) MG; (d) L Promedio de la probabilidad de detección para la señal aleatoria con distribución BG obtenido con: (a) MPS; (b) MSKA; (c) MG; (d) L Ejemplo que ilustra las capacidades de los métodos de deconvolución en presencia de ruido impulsivo (los círculos indican el valor verdadero de la señal dispersa): (a) se introduce en la muestra 65 de las observaciones una espícula de magnitud 2; (b) filtrado anticausal de las observaciones como primer paso del MSKA (el trazo continuo representa la señal filtrada y el discontinuo la señal sin filtrar) ; (c) resultado de aplicar el algoritmo MSKA; (d) resultado del método MG; (e) resultado del método L MSE en función de SNR para todos los métodos: MSKA, MPS, MG, L (trazo continuo, discontinuo, punteado y punteado discontinuo, respectivamente) en presencia de ruido impulsivo para diferentes valores de probabilidad Bernoulli: (a) p=.; (b) p=.5; (c) p= F ssd para una entrada determinista con ruido impulsivo (p=.), con distribución BG, en función del parámetro libre para SNR = 4, 8, 2 y 2 db (trazo continuo, discontinuo, punteado y punteado discontinuo, respectivamente) en: (a) MPS; (b) MS- KA; (c) MG; (d) L Promedio de la probabilidad de detección para la señal determinista con ruido impulsivo, y probabilidad p =., obtenido a través de: (a) MPS; (b) MSKA; (c) MG; (d) L

22 XX LISTA DE FIGURAS F ssd del algoritmo MSKA calculado para sistemas de fase mínima, c=.6, y no mínima, c=2, (trazo continuo y discontinuo, respectivamente): (a) distribución geométrica de los polos y ceros de la respuesta al impulso de un sistema de fase mínima al que se le desplaza el cero hasta sacarle fuera de la circunferencia de radio unidad, haciéndole de fase no mínima; (b) F ssd para una SNR = 4dB; (c) F ssd para una SNR = 8dB Resultados obtenidos al aplicar el algoritmo de deconvolución MSKA a dos señales con respuesta al impulso de fase no mínima: (a) señal determinista de fase no mínima ruidosa, c =.4; (b) señal determinista de fase no mínima ruidosa, c = 2.; (c) señal dispersa estimada, c =.4; (d) señal dispersa estimada, c = Esquema que ilustra las modificaciones previas realizadas sobre las muestras de ruido, filtradas con el inverso de la respuesta al impulso Procesamiento de señales sísmicas: (a) esquema que ilustra la obtención del CSPG en seis puntos diferentes, así como un esquema para el punto medio común m 6 ; (b) el concepto de multiplicidad; (c) esquema para un único trazo fuente-receptor del punto m 6 e ilustración del concepto de apilamiento; (d) datos reales obtenidos de un conjunto CSPG [7]; (e) trazas más cercana y más lejana del conjunto CSPG y dos formas estándar de representación: trazo tipo WT y WTVA Esquema general que muestra las tres técnicas más importantes en el procesamiento de señales sísmicas Estimación de la onda sísmica: (a) obtenida traza a traza utilizando DP; (b) onda sísmica obtenida promediando las ondas sísmicas estimadas de cada una de las trazas del CSPG; (c) valor normalizado para la señal dispersa estimada; (d) comparación del valor verdadero de la traza con su valor calculado a partir de la convolución de la onda y señal dispersa estimadas Resultados de la deconvolución aplicada sobre un conjunto CSPG utilizando ε = 3 : (a) conjunto de señales dispersas estimadas; (b) conjunto de señales reconstruidas; (c) conjunto de señales originales Resultados de la deconvolución aplicada sobre un conjunto CSPG: (a) señal dispersa obtenida mediante DP; (b) señal dispersa obtenida mediante deconvolución Wiener; (c) señal dispersa obtenida con método de frecuencia; (d) señal dispersa obtenida mediante MG ciego

23 LISTA DE FIGURAS XXI 4.. Principio generalizado de superposición: (a) el concepto de transformación lineal ligado al principio de superposición se ha generalizado a través de una transformación no lineal, H, definida como un sistema homomórfico; (b) un sistema homomórfico representado a través de tres subsistemas Operador de convolución como sistema característico: (a) representación canónica; (b) aproximación usando la DFT Caso ideal en DH: la componente de señal de la izquierda, de variación lenta, y la de la derecha, de variación rápida, están claramente separadas Filtros en el dominio cepstral: (a) paso corto; (b) paso largo SLIT con entrada y salida periódicas de banda limitada Señales sintéticas: (a) respuesta al impulso; (b) señal dispersa; (c) convolución de ambas; (d) magnitud de Y (e jω ); (e) valor principal del ángulo de Y (e jω ); (f) valor principal del ángulo sin componente de fase lineal de Y (e jω ); (g) fase desenrrollada de Y (e jω ) CC de: (a) la respuesta al impulso, ȟ[n]; (b) la señal dispersa, ˇx[n] Señal sintética de prueba: (a) CC; (b) respuesta al impulso verdadera y estimada; (c) señal dispersa de entrada verdadera y estimada Señal sintética de banda limitada: (a) onda Ricker; (b) señal dispersa formada por tres impulsos; (c) señal sísmica de salida; (d) amplitud del espectro; (e) fase del espectro CC para las señales: (a) onda Ricker (MD); (b) onda Ricker (MR y DS); (c) onda sísmica (MD); (d) onda sísmica (MR y DS);(e) onda Ricker original y estimada; (f) señal dispersa original y estimada Modelo de parábola invertida: (a) representación de un ciclo; (b) tren de dos impulsos; (c) señal formada por dos ciclos; (d) amplitud del espectro; (e) fase del espectro CC del modelo de parábola invertida: (a) MD (trazo continuo) y basado en las raíces (trazo discontinuo); (b) CC de dos ciclos del modelo de parábolas calculado con el MD; (c) utilizando rellenado con ceros Señales originales y resultados de la DH en DS: (a) modelo parabólico de señal (trazo continuo) y su estimación (trazo discontinuo); (b) señal dispersa original (círculos) y su estimación (trazo discontinuo) Modelo de parábola invertida: (a) representación de un ciclo; (b) tren de seis impulsos con reflexiones; (c) señal formada por seis ciclos; (d) amplitud del espectro; (e) fase del espectro. 35

24 XXII LISTA DE FIGURAS 4.5. Señales originales y resultados de la DH en DS : (a) modelo parabólico de señal (trazo continuo) y su estimación (trazo discontinuo); (b) señal dispersa original (círculos) y su estimación (trazo discontinuo) Curvas de presión y flujo que ilustran la secuencia de hechos mecánicos que se producen durante un único ciclo cardíaco (figura adaptada de []) Ondas de presión y velocidad: (a) la onda de presión es causa del ensanchamiento de la pared arterial; (b) ondas de presión y velocidad en diferentes puntos de la aorta (Figuras adaptadas de []) Esquema de un circuito eléctrico correspondiente a un modelo Windkessel: (a) dos elementos; (b) tres elementos Esquemas para determinar el efecto de las reflexiones (la línea punteada representa los promedios): (a) las medidas experimentales de presión y flujo están relacionadas por la impedancia de entrada, Z in ; (b) en un modelo idealizado sin reflexiones la relación entre presión y flujo es directamente proporcional, el valor de la constante es la impedancia característica Problema hemodinámico: (a) Problema hemodinámico directo, dados los valores de compliancia, área de sección transversal y longitud de las arterias se estima el valor de la impedancia de entrada; (b) problema hemodinámico inverso, dada la impedancia de entrada (o valores de presión y flujo) se estiman los valores de compliancia, área y longitud Señales de presión y flujo en la entrada de la aorta: (a) un solo ciclo que muestra los periodos diastólico y sitólico; (b) en el cuarto ciclo ya no hay una recuperación de la señal de presión Presión y flujo para el cerdo : (a) onda de presión; (b) onda de flujo; (c) espectro de la onda de presión; (d) espectro de la onda de flujo; (e) raíces de la onda de presión; (f) raíces de la onda de flujo; (g) distribución de las raíces de la onda de presión; (h) distribución de las raíces de la onda de flujo Presión y flujo para el cerdo 2 : (a) onda de presión; (b) onda de flujo; (c) espectro de la onda de presión; (d) espectro de la onda de flujo; (e) raíces de la onda de presión; (f) raíces de la onda de flujo; (g) distribución de las raíces de la onda de presión; (h) distribución de las raíces de la onda de flujo CC en registros del cerdo : (a) CC de la presión (MD); (b) CC del flujo (MD); (c) CC de la presión (MR); (d) CC del flujo (MR); CC en registros del cerdo 2 : (a) CC de la presión (MD); (b) CC del flujo (MD); (c) CC de la presión (MR); (d) CC del flujo (MR);

25 LISTA DE FIGURAS XXIII 5.. IEA y AEA calculada sobre los registros del cerdo : (a) IEA (MR); (b) IEA (MD); (c) IEA (MR) con registros filtrados; (d) IEA (MD) con registros filtrados; (e) IEA (MR, -5Hz); (f) IEA (MD, -5Hz); (g) AEA (MR, -5Hz); (h) AEA (MD, -5Hz) IEA y AEA calculada sobre los registros del cerdo 2: (a) IEA (MR); (b) IEA (MD); (c) IEA (MR) con registros filtrados; (d) IEA (MD) con registros filtrados; (e) IEA (MR) (-5Hz); (f) IEA (MD, -5Hz); (g) AEA (MR, -5Hz); (h) AEA (MD, -5Hz) Respuesta al impulso (cerdo ) calculada a partir de: (a) IEA (MR); (b) IEA (MD); (c) AEA (MR); (d) AEA (MD) Respuesta al impulso (cerdo 2) calculada a partir de: (a) IEA (MR); (b) IEA (MD); (c) AEA (MR); (d) AEA (MD) Señal suplemental: (a) magnitud del espectro; (b) distribución de polos y ceros IEA, AEA y respuesta al impulso del cerdo : (a) IEA calculada mediante transformación al DS; (b) AEA calculada mediante transformación al DS; (c) respuesta al impulso de la IEA; (d) respuesta al impulso de la AEA IEA y AEA del cerdo2 : (a) IEA calculada mediante transformación al DS; (b) AEA calculada mediante transformación al DS; (c) respuesta al impulso de la IEA; (d) respuesta al impulso de la AEA IEA y RI promedios calculados mediante segmentación en 4 segmentos de los registros del cerdo y 2 : (a) IEA (cerdo ); (b) RI calculada a partir de (a); (c) AEA (cerdo ); (d) RI calculada a partir de (c); (e) IEA (cerdo 2); (f) RI calculada a partir de (e); (g) AEA (cerdo 2); (h) RI calculada a partir de (g) B.. Raíces de dos señales reales: (a) distribución de las raíces de la señal de presión aórtica en el plano complejo; (b) histograma del radio de las raíces de la señal de presión; (c) distribución de las raíces de la señal de flujo; (d) histograma del radio de las raíces de la señal de flujo B.2. Rejilla de búsqueda del método de cálculo de raíces basado en el algoritmo de la transformada rápida de Fourier

26 XXIV LISTA DE FIGURAS

27 Lista de Tablas 3.. Ejemplos de funciones núcleo de Mercer Problemas primal y dual en clasificación SVM con margen suave lineal y no lineal Problemas primal y dual para SVR con función de coste ε- insensible en los casos lineal y no lineal, y con función de coste ε-huber en el caso no lineal Problemas SVM primal y dual para el marco lineal con función de coste ε-huber Algoritmo propuesto para la estimación espectral con SVM dentro del marco lineal general Algoritmo propuesto para el modelado ARMA con SVM utilizando el marco lineal general Problemas primal y dual para interpolación mediante SVM, con función de coste ε-huber, en el modelo primal de señal Definición de las variables V p, F p, F n y V n MSE±STD (x) para señales determinista con ruido (en negrita se ha resaltado el algoritmo que obtiene las mejores prestaciones y en cursivas el segundo mejor): (a) Gausiano; (b) Laplaciano; (c) uniforme MSE±STD(x) para señales aleatoria con ruido Gausiano (en negrita se ha resaltado el algoritmo que obtiene las mejores prestaciones y en cursivas el segundo mejor) Parámetros libres calculados a partir de diferentes figuras de mérito: (a) F ssd ; (b) F pick y F null ; (c) Se y Es MSE±STD(x) para señal determinista con ruido impulsivo (en negrita se ha resaltado el algoritmo que obtiene las mejores prestaciones y en cursivas el segundo mejor): (a) p=.; (b) p=.5; (c) p= Parámetros libres calculados para el ruido impulsivo a partir de diferentes figuras de mérito y probabilidad p =.: (a) F ssd ; (b) F pick y F null ; (c) Se y Es MSE obtenido a partir de la DD de una señal determinista con el método MSKA, donde se compara el efecto de modificar previamente las muestras de ruido Procesamiento de datos sísmicos xxv

28 XXVI LISTA DE TABLAS 4.. Algoritmo de DH usando el método de la señal suplemental Resultados del experimento con la señal sísmica Determinación de la señal dispersa mediante MD y DS a partir de varios ciclos del modelo de señal de parábola invertida Determinación de la señal dispersa con reflexiones mediante MD y DS a partir de varios ciclos del modelo de señal de parábola invertida Contribución a la presión y el flujo medido en cualquier punto del sistema circulatorio, debido a cualquiera de los cuatro posibles frentes de onda. El símbolo significa aumento y disminución Relación cualitativa entre los parámetros del sistema arterial y la impedancia de entrada en tres rangos diferentes de frecuencias Algunos parámetros para el CC de las señales de presión y flujo, utilizando el MD y el MR Constantes del método de DH en el DS B.. Algoritmo para calcular la factorización de polinomios

29 Capítulo Introducción La convolución de señales es uno de los pilares fundamentales del análisis y procesado de señal, y se basa en los principios de linealidad e invarianza temporal [87]. El problema opuesto de la convolución, la deconvolución, es igualmente importante debido a la gran cantidad de aplicaciones que requieren la separación de señales que pueden modelarse mediante una mezcla convolucional [76]. La convolución de señales puede interpretarse como un problema directo, es decir, se busca la señal de salida en un sistema lineal e invariante en el tiempo (SLIT) a partir de los valores conocidos de la señal de entrada y de la respuesta al impulso que caracteriza a dicho sistema. De la misma forma, la deconvolución puede interpretarse como un problema inverso, esto es, conocidas la salida y en principio, también la respuesta al impulso, se busca estimar la señal de entrada [76]. Existen muchas aplicaciones que se pueden modelar mediante la convolución entre una señal de entrada, formada por unas pocas espículas, denominada señal dispersa, y la respuesta al impulso de un SLIT. Al problema de encontrar las espículas a partir de la señal de salida se le ha llamado deconvolución dispersa (DD). Si se conoce a priori la respuesta al impulso, se habla de DD no ciega, y en caso contrario se hablará de DD ciega. La presente Tesis introduce nuevos métodos de DD tanto en los casos donde se conoce la respuesta al impulso, como en aquellos donde se desconoce. De una parte, en la sección 3.4 se proponen algoritmos para DD no ciega basados en la utilización de una de las técnicas emergentes en procesado de señal, esto es, las máquinas de vectores soporte (SVM, Support Vector Machine) [42]. Las SVM se amoldan muy bien a los requerimientos en la solución del problema de deconvolución dispersa no ciega, ya que conllevan una regularización implícita y, además, producen soluciones dispersas [4]. De otra parte, en la sección 4.4 se propone una modificación al algoritmo de deconvolución homomórfica (DH) propuesto en [7], motivada por el estudio de un problema práctico. El problema en cuestión es la estimación de la impedancia de entrada aórtica a partir de registros cardíacos de presión y flujo medidos en la arteria aorta. La modificación de la DH se basa en una

30 2 CAPÍTULO. INTRODUCCIÓN transformación de la señal a deconvolucionar, trasladándola a un espacio de señal donde se garantiza una separación adecuada entre las señales que forman el conjunto de observaciones. En el presente capítulo se hará una breve introducción al problema de la DD, analizando, a través del gran número de aplicaciones, su importancia y su trascendencia en problemas de señal. También se describirán las motivaciones principales de la presente investigación. Por último, se definirán el plantamiento y los objetivos de la Tesis... Interés y aplicaciones de las técnicas de deconvolución Fue la necesidad de encontrar yacimientos de combustibles fósiles la que dio origen a los primeros métodos de DD [2]. Las primeras aplicaciones se centraron en la caracterización del subsuelo terrestre, a través de la sísmica de reflexión [3]. Los registros obtenidos de estas aplicaciones se modelaban como la convolución entre una onda sísmica y los reflectores correspondientes a las interfases entre estructuras geológicas. Dado que la onda sísmica se propaga a través del subsuelo, originando nuevas ondas como producto de las reflexiones, este modelo resulta ser una buena aproximación [75]. Los primeros métodos de DD que surgieron requerían algún conocimiento de la onda sísmica, por ejemplo, su autocorrelación [9, 3]. Sin embargo, en esas primeras aplicaciones no era posible conocer la onda sísmica y, por tanto, fue necesario utilizar algunas suposiciones sobre las características de la señal y sobre la naturaleza de la onda que se propaga. La suposición de que la distribución de los reflectores puede modelarse como un proceso aleatorio de ruido blanco, condujo al resultado de que la autocorrelación de la onda sísmica es igual, excepto por un factor de escala, a la autocorrelación de las observaciones [3]. Esta suposición permitió utilizar la DD sin la necesidad de conocer a priori la onda sísmica, dando origen a lo que se conoce como DD ciega. Sin embargo, recientemente se ha cuestionado la suposición que se hizo sobre la distribución de ruido blanco de los reflectores [5, 58]. En [57, 58] se propuso la medición directa de la onda sísmica, y si bien tal propuesta generó mucho rechazo [89], no obstante abrió el camino a la exploración de métodos de DD no ciega. Recientemente, la DD se ha empezado a utilizar en aplicaciones más complejas, las cuales demandan métodos más robustos de estimación. La búsqueda de yacimientos geotérmicos, donde la estructura del subsuelo es mucho más caótica que la de aquellos que contienen yacimientos petrolíferos, debido a la naturaleza volcánica del subsuelo, requiere métodos de DD mucho más robustos [94]. A través de métodos de DD se puede estimar los tiempos de llegada de la onda sísmica y, a partir de éstos, se construyen mapas de velocidades que permiten identificar la presencia de yacimientos geotérmicos [93].

31 .2. LA DECONVOLUCIÓN EN PROBLEMAS DE SEÑAL 3 Esta técnica, junto con los recientes avances en algoritmos de interpolación y visualización, ha permitido la construcción de representaciones del subsuelo terrestre en dos y tres dimensiones, las cuales ofrecen mejores herramientas de exploración [24]. En el campo de la ingeniería de comunicaciones, la deconvolución se encuentra en aplicaciones de igualación de canal. Por un lado, los métodos no ciegos se desarrollan a partir de la transmisión de una secuencia conocida de datos a través del canal, la cual sirve como señal de entrada [92]. Por otro lado, los algoritmos ciegos utilizan como única información disponible un conjunto de observaciones a la salida del canal [47, 48]. Los primeros algoritmos para deconvolución ciega o igualación surgen como resultado del intento de resolver las interferencias entre símbolos en comunicaciones vía modem, originadas al aumentar la velocidad de comunicación [68, 69]. Otro hito importante en el desarrollo de algoritmos de deconvolución e igualación se establece, en la década de 98, con la introducción del algoritmo de módulos constantes [4, 35], el cual tiene una amplia aplicación en modulación digital de señales [32]. De muy especial interés resultan las aplicaciones en formación de imágenes radar de apertura sintética [23]. A partir del modelo de señal propuesto en [77], donde se modela la energía recibida por el radar a través de un modelo convolucional, ha sido posible introducir técnicas de deconvolución en el mejoramiento de la calidad de las imágenes. Otra aplicación más reciente y muy importante de los métodos de deconvolución no ciegos se encuentra en el área de procesamiento digital de señales biomédicas, ya sea en una o en varias dimensiones. En el caso particular de imágenes médicas (caso bidimensional) se pueden considerar tres áreas donde la deconvolución ha tenido un gran impacto: ultrasonido, resonancia magnética y tomografía computarizada. Las técnicas de restauración de imágenes se basan en modelar la imagen degradada como la convolución entre la imagen real y la función punto-dispersión (PSF, Point Spread Function). Esta aplicación es un ejemplo de deconvolución no ciega, ya que la restauración de la imagen requiere del conocimiento a priori de la PSF [29]. Todas estas aplicaciones plantean nuevos retos que requieren la renovación de viejos paradigmas y la introducción de nuevas ideas en el desarrollo de técnicas robustas de deconvolución..2. La deconvolución en problemas de señal La deconvolución ha sido un campo activo de investigación durante más de cincuenta años. En [] y [3] se propusieron métodos que, sin hacer uso del término deconvolución, permitían determinar, con cierta precisión,

32 4 CAPÍTULO. INTRODUCCIÓN las reflexiones en un problema de DD. Sin embargo, el fundamento teórico que impulsó el desarrollo de técnicas de deconvolución se encuentra, unos años más atrás, en el trabajo de [5] y en la implementación algorítmica desarrollada en [65]. Norbert Wiener proponía en [5]: (...) unir la teoría y práctica de dos campos de trabajo (...) las series temporales en estadística y la ingeniería de comunicaciones. La implementación condujo al diseño de filtros de mínimos cuadrados, aunque el trabajo de [65] es ampliamente reconocido por la solución eficiente de sistemas normales de ecuaciones asociados al diseño de filtros de mínimos cuadrados. La relación explícita entre series temporales y deconvolución fue establecida en [3], dando origen al método de deconvolución de mínimos cuadrados más importante: la deconvolución predictiva (DP) [89]. Los métodos de solución por mínimos cuadrados y mínimos absolutos han estado presentes durante más de dos siglos [28], y su aplicación a problemas de DD surgieron con el nacimiento del procesamiento de señales en tiempo discreto. Los métodos de deconvolución desarrollados durante las décadas de 95 y 96 fueron principalmente mínimos cuadrados [48]. No fue sino hasta la década 97 que se introdujo la solución por mínimos absolutos o norma L [25]. En [3] se utilizó un factor de penalización que permitía la regularización del problema, y además se formuló, de una forma muy precisa la DD como problema de programación lineal (PL). En [8] se exploraron las ventajas del método Simplex como técnica de solución del problema de PL en DD. Así mismo, los trabajos en norma L estimularon investigaciones en normas más generales [7, 54]. A finales de la década de 97 y durante toda la década de 98, se inició y desarrolló uno de los métodos de deconvolución más importantes de los últimos 25 años, la deconvolución de máxima verosimilitud (ML, Maximum Likelihood). El método de deconvolución ML se inició utilizando modelos de variables de estado [6], pero su amplia utilización en problemas de procesado de señal condujo a su formulación convolucional [76]. ML incorporó conocimiento de áreas tales como la identificación de sistemas y la teoría de detección y estimación. Así mismo, se incorporó un nuevo modelo estadístico para la señal dispersa [75, 76]. En la década de 99 empezaron a publicarse trabajos que demostraban que una de las hipótesis más importantes en los métodos de DD clásicos, el modelado de la señal dispersa como una secuencia de ruido blanco aleatorio, carecía de fundamento real [2, 4, 57]. Desde entonces, se han propuesto nuevos modelos para la señal dispersa, y en [3] se analizan los más importantes. Ricker, en [], acuñó el término compresor de ondas, ya que interpretaba las observaciones como la superposición de ondas originadas por el efecto de reflexión a través de discontinuidades. Dichas reflexiones podían calcularse aplicando un filtro capaz de contraer cada una de las ondas hasta hacer de ellas impulsos que correspondiesen con los coeficientes de reflexión. Por otro lado, Robinson, en [3], utilizó el término descomposición predictiva, que más tarde sustituiría por deconvolución predictiva.

33 .3. MOTIVACIÓN E INTERÉS DE LA INVESTIGACIÓN 5 A modo de reflexión cabe decir que las modificaciones actuales sobre la naturaleza de la señal dispersa no han conducido a la desestimación de los métodos clásicos, basados en el filtrado de mínimos cuadrados, como por ejemplo, la solución mediante el filtro de Wiener. De hecho en [3, 4] se resalta que la DD mediante el filtrado de Wiener es uno de los métodos más robustos frente a una gran variedad de condiciones, como por ejemplo, la estadística del ruido o la señal de entrada..3. Motivación e interés de la investigación Los algoritmos propuestos en la presente Tesis se basan en dos técnicas muy diferentes, que se corresponden con el caso donde se conoce la respuesta al impulso y en el que se desconoce completamente. El caso donde se conoce la respuesta al impulso tiene como problemas principales: el mal condicionamiento, la sensibilidad al tipo de ruido y la limitante impuesta de poder encontrar soluciones únicamente en los casos donde la fase de la respuesta al impulso es mínima. El problema del mal condicionamiento se da cuando pequeñas variaciones en las observaciones conducen a grandes variaciones en las estimaciones. La sensibilidad al tipo de ruido se da cuando variaciones en la estadística del ruido o la presencia de valores atípicos conducen a malas estimaciones de la señal dispersa. Por último, el problema de la fase mínima se da cuando no se puede encontrar una solución directa a través de un filtro inverso estable. El caso donde no se conoce la respuesta al impulso adolece de los mismos problemas del caso anterior, pero la situación es aún más compleja, y es necesaria la introducción a priori de información sobre la señal, como por ejemplo, algún conocimiento sobre la estadística de las observaciones o sobre el ancho de banda de las señales implicadas. Por tanto, aquí se exploran en dos direcciones las posibilidades para nuevos algoritmos de DD: para el caso no ciego, se propone extender las propiedades de robustez y dispersión de las SVM a la solución del problema; y para el caso ciego, se proponen modificaciones del algoritmo de la DH de acuerdo con la naturaleza conocida de la aplicación práctica a resolver en la presente Tesis, es decir, la estimación de la impedancia de entrada aórtica. Si bien las SVM han sido propuestas inicial y fundamentalmente como procedimiento de clasificación en problemas biclásicos [42], recientemente se han desarrollado algoritmos SVM para problemas de procesado digital de señal, que requieren de técnicas robustas de procesado, tales como el modelado ARMA (Auto Regressive Moving Average), el análisis espectral clásico, el filtrado gamma o la interpolación [7,,, 8]. Dadas las buenas prestaciones observadas en la extensión de las SVM a problemas clásicos de procesado digital de señal, se extenderán las propiedades de robustez de la SVM a la solución del problema de DD. La DH fue desarrollada de forma simultánea y por fuentes separadas en

34 6 CAPÍTULO. INTRODUCCIÓN [, 83]. Su aplicación como técnica ciega ha tenido un amplia aceptación en muchos campos de la ingeniería [23, 82]. Las modificaciones recientes del método [7], así como los desarrollos de métodos numéricos para el cálculo de raíces [22], y a partir de éstas los coeficientes cepstrales, han estimulado el uso de la DH como método robusto de deconvolución ciega. Su incorporación en la presente Tesis ha sido motivada, además de sus buenas prestaciones, por la necesidad de tener una técnica robusta de deconvolución que permita resolver una aplicación particular, esto es, la caracterización de la impedancia de entrada aórtica..4. Planteamiento y objetivos de la Tesis En la presente Tesis se desarrollarán métodos de DD basados en dos técnicas diferentes, como son la deconvolución mediante SVM y la DH. La primera sirve como alternativa a propuestas no ciegas y se basa en dos conceptos claves: la minimización de funciones de coste robustas y las técnicas de optimización. La segunda es una alternativa a métodos ciegos propuestos hasta la fecha y se basa en modificaciones del algoritmo de DH. La extensión de las SVM como método de DD requiere el desarrollo de nuevos modelos algorítmicos. El desarrollo de un modelo adecuado para la solución del problema ha conducido a explorar dos propuestas SVM dentro del marco de procesado de señal: En la primera propuesta, se estudiarán los trabajos en SVM realizados en aplicaciones de procesado de señal [7,,, 8] y generalizados dentro de un marco lineal general en [2, 7]. Esta generalización permitirá incluir a la DD como un caso particular dentro del marco SVM lineal general. Este método de deconvolución será llamado deconvolución SVM basada en el modelo primal de señal (MPS). Del análisis teórico y de la simulación sobre el MPS se obtendrán importantes conclusiones sobre sus ventajas y limitaciones, que conducirán a la búsqueda de otros algoritmos que superen sus limitaciones. La segunda propuesta se basará en el modelo dual de señal (MDS), introducido en [8] como modelo alternativo en problemas de interpolación, que utiliza el algoritmo SVM desarrollado en regresión (SVR, Support Vector Regression). Este modelo formulará el problema de DD como una regresión no lineal e incorporará la utilización de los espacios de características. El análisis teórico del MDS mostrará ventajas sobre el MPS, pero a la vez dificultades teóricas en su implementación, por ejemplo en sistemas causales. La necesidad de un algoritmo para aplicaciones de DD ciega ha llevado en esta Tesis a la utilización y modificación del algoritmo para la DH propuesto en [86]. En este sentido, se explorarán aquí las siguientes propuestas:

35 .4. PLANTEAMIENTO Y OBJETIVOS DE LA TESIS 7 Se desarrollará, a través del estudio de los algoritmos precedentes de DH, una nueva versión que se adaptará a los modelos de señal más comúnmente encontrados en problemas de DD ciega. Esta modificación conducirá a la transformación previa de la señal a un dominio que se ha llamado dominio suplemental (DS), que mejorará los resultados de la DD ciega. Se aplicará el algoritmo modificado de DH a una aplicación de bioingeniería, la caracterización de la impedancia de entrada aórtica. Esta aplicación requerirá de algoritmos robustos que permitirán hacer inferencias sobre las características del sistema cardiovascular a partir de mediciones realizadas en el flujo sanguíneo. La presente investigación se estructurará en seis capítulos, los cuales se detallan a continuación: El presente capítulo ha servido de introducción elemental al problema de la DD. El Capítulo 2 realizará una revisión bibliográfica de los métodos clásicos y las aportaciones más recientes en deconvolución ciega y no ciega. En el Capítulo 3 se introducirán las SVM para los problemas clásicos de clasificación y regresión, y se estudiarán las aplicaciones más recientes de las SVM a problemas de procesado de señal. Además, se presentará la primera aportación de la presente Tesis: el desarrollo de dos algoritmos de DD no ciega basados en SVM. Por último, se presentarán los resultados, a través de varios experimentos, que ilustrarán la robustez del método para diferentes tipos de ruido, de señales de entrada y de respuestas al impulso. El Capítulo 4 describirá brevemente la DH y presentará modificaciones realizadas sobre la señal que permitirán mejorar las prestaciones del algoritmo de DH. Una de estas modificaciones es la segunda aportación de la presente Tesis, esto es, la DH en el Dominio Suplemental. También al final del capítulo, se realizarán pruebas sobre señales sintéticas. El Capítulo 5 desarrollará una aplicación práctica de los métodos de DH en bioingeniería, esto es, la caracterización de la impedancia de entrada aórtica. Finalmente, el Capítulo 6 se resumirán las conclusiones de la investigación realizada en la presente Tesis, y se propondrán las líneas futuras de investigación.

36 8 CAPÍTULO. INTRODUCCIÓN

37 Capítulo 2 Deconvolución de señales dispersas En muchas situaciones prácticas se desea caracterizar la estructura interna de un sistema físico a partir de medidas del comportamiento del mismo. También, es posible que a partir de estas medidas se tenga interés no en la estructura interna, sino en la estimulación que dio origen a dichas medidas. Estas dos situaciones se pueden formular como un único problema [44] que se conoce con el nombre de problema inverso, y que matemáticamente se formula a través de la ecuación integral siguiente: y(t) = t K(t, τ)x(τ)dτ + e(t) (2.) donde las señales y(t), K(t, τ) y x(t) representan la salida, la estructura interna y la entrada del sistema físico en cuestión, respectivamente, y e(t) es un proceso aleatorio que modela la componente de ruido. La Ecuación (2.) se conoce con el nombre de ecuación integral de Fredholm de primera clase [9, 36], y la función K(t, τ) en la misma es la función núcleo o kernel la cual debe cumplir con las condiciones de la sección Un caso particularmente importante se da en los problemas donde el núcleo es una función del tipo convolucional, es decir, el núcleo depende solamente de la diferencia entre dos variables independientes. Este tipo particular de problema inverso se conoce con el nombre de deconvolución [45] y matemáticamente se escribe como: y(t) = t h(t τ)x(τ)dτ + e(t) (2.2) donde h(t τ) es la respuesta al impulso del sistema, la cual se desplaza en función de la variable τ. Así, el problema de deconvolución se refiere a la solución de una ecuación integral lineal del tipo Fredholm de primera clase con núcleo invariante [55]. Cabe hacer notar que, en la Teoría de los Sistemas Lineales, los límites de la integral (2.2) están definidos en el intervalo 9

38 CAPÍTULO 2. DECONVOLUCIÓN DE SEÑALES DISPERSAS (, ), con lo que se consigue incluir aquellas funciones definidas en todo el intervalo de tiempo. Desde el punto de vista de la Teoría de Operadores, la integral (2.2) se interpreta como una transformación lineal sobre la señal de entrada x(t). Desde el punto de vista de los sistemas lineales, la Ecuación (2.2) se conoce con el nombre de integral de convolución, donde la señal de salida se origina a partir de la convolución de la señal de entrada con la respuesta al impulso. Así, la deconvolución se interpreta como una operación que cancela los efectos de una convolución efectuada sobre una señal [46]. El caso particular de la DD se da cuando el sistema tiene una repuesta al impulso formada por una serie de impulsos distribuidos de forma dispersa y aleatoria. A modo de ejemplo, y para facilitar la descripción del problema, se utiliza a continuación la aplicación que más ha estimulado el desarrollo de algoritmos de deconvolución: la sísmica de reflexión [75, 3]. Como puede verse en la Figura 2.(a), en esta aplicación la señal de salida corresponde generalmente a una señal sísmica generada artificialmente. La onda sísmica se origina a través del uso de explosivos o vibradores mecánicos que hacen propagar ondas a través del subsuelo. La onda sísmica se propaga hacia abajo y, cuando encuentra discontinuidades geológicas, parte de ésta se refleja hacia la superficie, donde es registrada por los equipos de medida. El modelo geológico de capas del subsuelo terrestre, tal y como se observa en la Figura 2.(a), corresponde a un modelo de parámetros distribuidos [75], y cualquier señal medida corresponderá a un conjunto de réplicas retrasadas y escaladas de la onda sísmica de entrada. Matemáticamente, la señal sísmica puede modelarse a través de la siguiente ecuación: y(t) = h j x(t t j ) + e(t) (2.3) j= donde x(t) e y(t) representan a la onda que se hace propagar a través del subsuelo y a la onda sísmica total medida a la salida, respectivamente; h j y t j representan la amplitud y el instante de tiempo correspondiente a la reflexión en el instante i-ésimo, respectivamente. Es posible reescribir la respuesta al impulso como un sumatorio de impulsos: h(t) = h j δ(t t j ) (2.4) j= donde δ(t) representa la función delta de Dirac. Debido a que la convolución es un operador conmutativo, es posible intercambiar el orden de las señales a convolucionar, y suponer que el tren de impulsos de (2.4) sirve de entrada a un sistema de respuesta al impulso igual a la onda sísmica. Este cambio se realiza debido a dos cuestiones importantes: en primer lugar se debe a que la serie de impulsos no posee una longitud finita definida y, en segundo lugar, la serie de impulsos se modela en muchas aplicaciones como una secuencia aleatoria [75].

39 Fuente Medidor Medidor x n y n Onda Transmitida x n h n en + y n (a) (b) Figura 2.: El problema de sísmica de reflexión: (a) el subsuelo terrestre modelado como un conjunto de estructuras geológicas estimulado por una onda sísmica; (b) su representación como un SLIT. La Ecuación (2.3) representa un modelo de señal de tiempo continuo de gran importancia en el modelado de aplicaciones que involucran señales dispersas, pero su utilidad es solamente teórica. Las aplicaciones que requieren la estimación de la señal dispersa a partir de un número limitado de observaciones se basan en un modelo equivalente de tiempo discreto donde la señal de entrada y la respuesta al impulso tienen longitud finita, el cual se escribe como sigue: P y[n] = x[j]h[n j] + e[n] = h[n] x[n] + e[n] (2.5) j= donde n =,,..., N. La señal dispersa de entrada x[n] y la respuesta al impulso h[n] tienen longitud finita igual a P muestras y Q=N-P + muestras, respectivamente. También, es posible simplificar la notación de sumatorio a través de la introducción del operador, que se conoce con el nombre de operador de convolución de tiempo discreto. Resulta de mucha utilidad escribir la suma de convolución (2.5) de forma matricial, y para ello se utiliza una representación algorítmica, dada por: y n = P x j h n j+ + e n (2.6) j= donde n =,..., N. Obsérvese que, a diferencia del modelo de tiempo discreto (2.5), los índices empiezan en la posición n =, es decir x n = x[n ]. Esta notación se mantendrá a lo largo de la presente Tesis cuando se utilice la representación algorítmica. Para tener en cuenta (2.6) en todos los valores del índice n, se puede

40 2 CAPÍTULO 2. DECONVOLUCIÓN DE SEÑALES DISPERSAS utilizar la formulación matricial siguiente: h... h 2 h... y y 2. y N.. h 2... =. h Q.... h Q..... x x 2. x P + e e 2. e N (2.7)... h Q Reescribiendo de forma simplificada, se tiene la expresión: y = Hx + e (2.8) donde la matriz H tiene dimensión N P, y corresponde a la matriz de convolución generada a partir de la respuesta al impulso h n [7]. El problema de la DD se describe desde dos puntos de vista claramente diferenciados, según el conocimiento que se tenga sobre la respuesta al impulso [5]. El primer enfoque supone conocida la respuesta al impulso h n, además de la señal de salida y n, y busca calcular la señal dispersa x n. Este caso es conocido como DD no ciega. El segundo enfoque supone que tanto la señal dispersa x n como el filtro h n son desconocidos, y se denomina DD ciega. En el presente capítulo se realizará una revisión bibliográfica de los diferentes métodos de DD. Se incluirán métodos básicos, como el de la transformada discreta de Fourier, de mucha utilidad en la introducción de algunos problemas que suelen encontrarse en este tipo aplicaciones. También se estudiarán métodos clásicos basados en mínimos cuadrados y mínimos absolutos, así como técnicas de estimación máximo a posteriori, entre las cuales destaca la deconvolución de máxima verosimilitud. En cuanto a las técnicas no ciegas, se revisará el método de la DH y el método de módulos constantes. Por último, se hará un análisis de varios aspectos relevantes encontrados al aplicar las técnicas de DD, entre los cuales están el modelo de señal dispersa y la fase de la respuesta al impulso. 2.. DD no ciega La DD es un problema clásico dentro del procesamiento digital de señales. De hecho, fue la deconvolución de señales dispersas, tal como se introdujo en [3], el primer método de procesamiento digital de señales sin una contraparte analógica [2]. Los métodos analógicos, anteriores a la introducción de la deconvolución, buscaban la estimación de señales dispersas a través de la aplicación directa de filtros de ranura sobre la la señal de salida. Por

41 2.. DD NO CIEGA 3 otro lado, la deconvolución se introdujo como una técnica de tiempo discreto [2]. A continuación, en esta sección se hará una revisión de las técnicas no ciegas más comunes Solución usando la transformada de Fourier En [28, 27, 5] se introdujo el algoritmo de la transformada rápida de Fourier (FFT, Fast Fourier Transform) en el cálculo de la transformada de Fourier de tiempo discreto (DFT, Discrete Fourier Transform). El problema de la deconvolución tiene una solución teórica directa en el dominio transformado a través del uso de la DFT, siempre que las señales sean de longitud finita [2]. La señal de entrada en principio puede calcularse directamente a través de: ( Y [k] ) x[n] = DFT H[k] ( E[k] ) + DFT H[k] (2.9) donde DFT representa la DFT inversa, y H[k], Y [k] y E[k] representan las DFT de h[n], y[n] y e[n], respectivamente. La estimación de la señal de entrada obtenida a través de (2.9) es muy sensible a errores en los datos, introducidos por la componente de ruido. En particular, cuando se divide por un valor de H[k] con un valor de módulo pequeño, se amplifica el efecto del ruido, obteniéndose una mala estimación de la señal de entrada [98]. Al segundo término del miembro de la derecha de (2.9) se le conoce como ruido de alta frecuencia, debido a qué es el responsable de la distorsión de alta frecuencia [32]. El ruido de alta frecuencia se produce debido al carácter paso bajo de la mayoría de sistemas físicos, esto es, la señal de entrada se ve severamente atenuada en sus componentes de frecuencias altas y, al realizar la operación de deconvolución, éstas se ven amplificadas junto con las componentes de ruido. Para ilustrar la solución al problema de deconvolución mediante la DFT, se introduce a continuación un ejemplo en el cual se recupera la señal de entrada mediante (2.9). Las ecuaciones de este ejemplo se presentan en la sección de experimentos del capítulo 4 y corresponde a la convolución de una onda Ricker con una secuencia de tres impulsos. La Figura 2.2(a) y 2.2(b) muestran la respuesta al impulso de un filtro paso banda y una señal dispersa formada por tres impulsos, respectivamente. En las Figuras 2.2(c) y 2.2(d) se muestran las amplitudes de los espectros de la respuesta al impulso y la señal dispersa, respectivamente. La amplitud de la DFT de la señal de salida con ruido aditivo Gausiano se muestra en la Figura 2.2(e). La amplitud de la DFT de la señal de entrada obtenida a partir de (2.9) se muestra en la Figura 2.2(f), donde se observa como el método es muy sensible al ruido aditivo sobre la señal de salida, no solamente en alta frecuencia, sino también, en este caso, en baja frecuencia, dado el carácter paso banda de la respuesta al impulso utilizada.

42 4 CAPÍTULO 2. DECONVOLUCIÓN DE SEÑALES DISPERSAS Amplitud n (a) Amplitud n (b) 2.2 Amplitud rad/seg (c) Amplitud rad/seg (d) Amplitud rad/seg (e) Amplitud rad/seg (f) Figura 2.2: Señales de prueba: (a) respuesta al impulso; (b) señal dispersa; (c) amplitud del espectro de la respuesta al impulso; (d) amplitud del espectro de la señal dispersa; (e) amplitud del espectro de la señal de salida con ruido; (f) señal dispersa estimada. Este enfoque basado en la DFT fue analizado en [98], donde se propuso la regularización para mejorar el problema de la alta sensibilidad al ruido aditivo. Dicha regularización condujo a la deconvolución mediante filtro de Wiener [3, 99], de la cual se hablará más adelante.

43 2.. DD NO CIEGA Solución de mínimo error cuadrático El método de mínimos cuadrados (LS, Least Squares) ha cumplido ya su bicentenario, y su incorporación en la solución del problema de DD surgió con el nacimiento mismo de la deconvolución. Fueron los trabajos, primero el teórico en [5] y luego la implementación en [65], los que condujeron al desarrollo de filtros de deconvolución de LS [89]. En esta sección solo se considerarán los métodos LS que se basan en la minimización directa del error cuadrático, con un enfoque puramente matricial tal como fue planteado en []. Solución de mínimo error cuadrático no regularizado Utilizando notación matricial, la solución LS consiste en encontrar un vector x que minimice la norma cuadrática del error, es decir se busca minimizar el funcional: J LS = y Hx 2 2 (2.) La solución viene dada por: ˆx = (H T H) H T y = H + y (2.) La expresión H + se conoce con el nombre de matriz generalizada de H, o pseudoinversa de Moore-Penrose [46]. Un resultado similar se obtiene con el estimador de máxima verosimilitud [49], en presencia de ruido Gausiano, donde se busca maximizar la probabilidad condicional p(y x). Obsérvese que H T H es una matriz cuadrada QxQ de tipo Toeplitz. Si se denota la autocorrelación de la respuesta al impulso por: γ k = Q h i h i+k, k =, 2,..., Q (2.2) i= se tiene que: γ γ 2... γ Q H T γ 2 γ... γ Q H =... γ Q γ Q... γ (2.3) es la matriz de autocorrelación de la respuesta al impulso h n, donde γ i es el elemento que forma i-ésima diagonal. La solución (2.) tuvo inicialmente alguna aceptación en la comunidad científica [5, 26, 37, ] fundamentalmente por su sencillez, pero al igual que en el método basado en la DFT, la solución de (2.) es muy sensible a las variaciones debido al ruido en la señal de salida. Esta alta sensibilidad El método de mínimos cuadrados fue utilizado por Gauss en 795, sin embargo fue publicado por Legendre en 85 [58].

44 6 CAPÍTULO 2. DECONVOLUCIÓN DE SEÑALES DISPERSAS significa que pequeñas perturbaciones en la señal de salida producen grandes variaciones en la señal dispersa estimada, a esto se le conoce como problema mal condicionado [44]. Formalmente, un problema está bien condicionado si cumple con tres condiciones [42, 44]: () tiene solución, (2) la solución es única y (3) la solución depende continuamente de datos y parámetros, es decir, pequeñas perturbaciones en los datos no deben causar grandes perturbaciones en los resultados. El problema está mal condicionado cuando no se cumple alguna de estas condiciones. A continuación se ilustrará con un ejemplo simple el problema de mal condicionamiento. Supónganse los valores de señal siguientes: [ ]. y =.25, A =.7., x =, e = donde y, A, x y e corresponden con el vector de observaciones, la matriz mal condicionada, el vector de entrada y el vector de ruido, respectivamente. La estimación de x obtenida a través de (2.) es: [ ] 7. ˆx = -8.4 El resultado obtenido utilizando la solución LS de (2.) conduce a una mala estimación del valor verdadero de la señal, x = [, ] T, a pesar de tener un valor pequeño para la norma del error, Ax y 2 =.22. La mala estimación del valor de x se debe al mal condicionamiento de la matriz A. En [54] se estudió el problema del mal condicionamiento de la matriz de convolución en problemas de deconvolución, el cual se atribuyó a la dependencia cuasi-lineal de sus filas. A través de la expansión en serie de Taylor, se demostró que cualquier fila de la matriz H se aproxima como la combinación lineal de dos filas adyacentes. La dependencia cuasi-lineal de las filas de H se debe a la suavidad de la respuesta al impulso. Así, un aumento de la frecuencia de muestreo empeora el condicionamiento de la matriz, ya que se incrementa la dependencia entre filas [54]. Una caracterización más explícita del problema de mal condicionamiento en problemas de deconvolución se realizó en [32], a través de la utilización de un modelo convolucional diferente de (2.8) y dado por: y + δy = H(x + δx) (2.4) donde δy y δx, representa las distorsiones en la señal de salida y entrada respectivamente. Utilizando la representación basada en la descomposición en valores singulares [32], se obtiene una solución para la señal de entrada dada por: x + δx = N i= y T u i λ i v i + N i= δy T u i λ i v i (2.5)

45 2.. DD NO CIEGA 7 donde {λ i, u i } y {λ i, v i } son los autovalores/autovectores de las matrices HH T y H T H, respectivamente. La solución dada por (2.5) caracteriza el problema de mal condicionamiento en un problema de deconvolución, y permite identificar el origen de las distorsiones en la estimación de la señal [32]. Del segundo término del miembro de la derecha de (2.5), se observa que, incluso en las situaciones donde la componente de distorsión es pequeña, el valor de δx puede ser grande, si el autovalor λ i correspondiente es pequeño, en comparación con el autovalor de valor máximo λ max. El problema de mal condicionamiento hace necesaria la búsqueda de alternativas que permitan mejorar las prestaciones del método de deconvolución LS. A continuación se presenta una solución a dicho problema. Solución de mínimo error cuadrático regularizado El mal condicionamiento puede corregirse a través de un proceso llamado regularización, el cual reemplaza la matriz pseudoinversa H + de (2.) por una matriz C γ, donde el parámetro γ se denomina parámetro de regularización [8]. La forma de regularización más utilizada se lleva acabo modificando el número de condición de la matriz H T H, dado por el cociente entre sus autovalores mayor y menor, agregando una valor constante a la diagonal principal. Esta modificación conduce a la solución LS regularizada siguiente: ˆx = (H T H + γi) H T y (2.6) donde γ > es el parámetro de regularización y la matriz C γ viene dada por: ( ) H C γ = H T T H + γi (2.7) La ecuación (2.6) es la misma solución obtenida por Tikonov en [3] a partir de la optimización del funcional: J LSR = y Hx γ x 2 2 (2.8) donde γ x 2 2 es el término de regularización y γ es el parámetro de regularización, interpretado en este caso como una constante de compromiso entre la minimización del error y la norma de la señal dispersa. A través de un enfoque puramente probabilístico, en [5, 28] se muestra que es posible obtener la solución regularizada de (2.6) a partir de la maximización de la probabilidad p(x y). Este principio se formula como el estimador máximo a posteriori siguiente: { ( ˆx MAP = max exp y ) ( Hx 2 2 exp 2σ 2 )} x 2 2 2σx 2 (2.9)

46 8 CAPÍTULO 2. DECONVOLUCIÓN DE SEÑALES DISPERSAS donde σ y σ 2 x son las varianzas del ruido y de la señal dispersa, respectivamente. Maximizar (2.9) equivale a minimizar el funcional: J MAP = y Hx σ2 x 2 σx 2 2 (2.2) que conduce a la expresión de LS regularizada (2.6), con γ = σ 2 /σ 2 x. El enfoque probabilístico muestra que la deconvolución de LS regularizada adolece de dos problemas fundamentales: la norma L 2 del error no favorece la obtención de señales dispersas, debido a la suposición de estadística Gausiana, y es necesario conocer a priori un valor adecuado para el parámetro de regularización γ. Otras técnicas de deconvolución basadas en LS se han propuesto en la literatura, a saber: mínimos cuadrados con restricciones [38, 78] y mínimos cuadrados total [74, 4]. Estas técnicas, si bien es cierto que realizan algún tipo de regularización que resuelve en cierta medida el problema del mal condicionamiento, no ofrecen resultados que se adapten a los requerimientos en la solución de problemas de DD, debido fundamentalmente a que presuponen una estadística Gausiana para la señal dispersa Solución mediante filtro de Wiener La solución del problema de deconvolución usando filtrado de Wiener busca la obtención de un filtro lineal óptimo, mediante un enfoque de series temporales. La idea fundamental del filtrado de Wiener consiste en que a partir de una versión ruidosa de la señal deseada, es posible obtener, mediante filtrado lineal óptimo, una estimación de ésta. La estimación óptima se obtiene a través de la minimización del error cuadrático medio, dado por: E{ x x opt 2 2} (2.2) donde x opt es la estimación óptima de la señal dispersa. La estimación lineal de x se obtiene multiplicando las observaciones y por una matriz de estimación [49], esto es: x opt = M opt y (2.22) donde M opt es la matriz de estimación óptima. La matriz M opt no se puede obtener a partir de la minimización de (2.2), debido a que el valor de x es desconocido, y en su lugar se busca un aproximación ˆx = x opt de x mediante el principio de ortogonalidad para filtrado lineal óptimo [46], y que consiste en minimizar el funcional: J W = E{(x x opt )y T } (2.23)

47 2.. DD NO CIEGA 9 Sustituyendo (2.8) en (2.23), se obtiene la expresión siguiente: E{xx T H T + xe T } E{(M opt Hx + M opt e)(x T H T + e T )} = (2.24) donde E{xx T } = C x y E{ee T } = C e son las covarianzas de la señal dispersa y del ruido, respectivamente, y E{xe T } =. Sustituyendo estos valores en (2.24) y despejando para la matriz óptima, se obtiene: M opt = (H T H + C e C x ) H T (2.25) Sustituyendo (2.25) en (2.22) se obtiene la expresión cerrada para la señal estimada: x opt = (H T H + C e C x ) H T y (2.26) La solución mediante filtro de Wiener requiere el conocimiento a priori de la covarianza de la señal dispersa, y debido a esto no es posible obtener una solución directa. En [9] se demuestra que es posible utilizar la covarianza de las observaciones y en [5] se utilizan métodos iterativos para encontrar mejores soluciones. La solución del problema de DD a través del filtro de Wiener obtenida en (2.25) coincide con: el estimador MAP (2.2) cuando x y e son realizaciones de procesos aleatorios Gausianos de media nula [49, 5], la solución LS regularizada (2.6), y la solución regularizada en el sentido Tikhonov [3]. El filtro de Wiener fue una de las primeras técnicas utilizadas en DD [3], y su aplicación dio origen a la DP, de la cual se hablará en la siguiente sección Deconvolución predictiva Si bien es cierto la DP se considera un método de deconvolución ciego, su enfoque basado en el filtro de Wiener permite que se le incluya con los métodos de deconvolución dispersa no ciegos. El filtro de predicción es un operador que actúa sobre un registro de entrada en un tiempo n y estima la amplitud del trazo en un tiempo futuro, n+k, donde n es la variable independiente que representa el tiempo discreto. Es razonable definir la salida deseada como una versión adelantada en el tiempo de la entrada, con lo cual el funcional a minimizar está dado por: J DP = E{ y n y n+k 2 2} (2.27)

48 2 CAPÍTULO 2. DECONVOLUCIÓN DE SEÑALES DISPERSAS Entrada Salida deseada y d n n Cálculo del operador de predicción a n Operador de predicción (a) Entrada y n Operador predictor de error f=[,,...,, - a ] x ^n Salida (b) Figura 2.3: Filtrado de Wiener en DP: (a) diseño del filtro; (b) aplicación a un conjunto de observaciones. En [9] se demostró que el sistema de ecuaciones para la deconvolución, basada en el filtro de Wiener, y el del filtro predictor de longitud unitaria, son idénticos, excepto por un factor de escala. Así, el uso de distancias de predicción mayores que uno conduce a una solución más general que el filtro óptimo de Wiener, ya que permite controlar la longitud de la señal de salida, y por lo tanto especifica el grado de resolución deseado de la señal dispersa (el caso ideal es el impulso unitario). El modelo general basado en el filtro de Wiener involucra tres señales, como se muestra en la Figura 2.3(a). La señal de entrada y, la señal deseada d y la señal de salida, que corresponde a los coeficientes del filtro predictor. La minimización del funcional (2.27) conduce a la solución mediante filtro de Wiener [46, 3], dado por el sistema normal de ecuaciones siguiente: r r... r Q r r 2... r Q r Q r Q 2... r a a. a Q = r k r k+. r k+q la cual se puede reescribir de forma simplificada como a continuación: (2.28) Ra = p (2.29) donde a es el vector con los coeficientes del predictor, R es la matriz de autocorrelación de la respuesta al impulso, y que para la suposición de señal dispersa como secuencia aleatoria blanca es igual a la autocorrelación de las observaciones [9], y p es la correlación cruzada entre la señal de entrada y la señal deseada. Los elementos de la matriz de autorcorrelación y de la correlación cruzada se calculan a partir de: r k = Q y i y i+k, k =, 2,..., Q (2.3) i=

49 2.. DD NO CIEGA 2 La solución de (2.29) permite obtener el operador de predicción de la Figura 2.3(a) el cual se utiliza como filtro de predicción de error, Figura 2.3(b). La salida de esta operación de filtrado es la señal de error, que bajo este modelo corresponde a la estimación de la señal dispersa. A pesar de las ventajas que introduce el parámetro k del predictor, fundamentalmente en el control sobre la resolución de la salida, el método adolece de los mismos problemas que los métodos LS y filtrado de Wiener, ya que DP es en sí un método LS. En las secciones siguientes se introducen nuevos enfoques en la solución del problema de deconvolución Deconvolución como problema de optimización El problema de DD no ciega se puede formular como un problema de optimización [8] que consiste en minimizar el funcional: J P O = e p p + λ q ˆx q q (2.3) donde e = y Hˆx es el vector de residuos, p y q son las normas p y q, respectivamente, y λ q es un parámetro de regularización, que se denomina peso relativo o factor de amortiguamiento. El parámetro λ q se selecciona como compromiso entre dos objetivos antagónicos [8], la minimización del error e y la conservación del carácter disperso de la señal de entrada ˆx. El caso particular p = 2 y λ p = corresponde a la solución de LS, mientras que la combinación p = 2 y λ 2 corresponde a LS regularizado, descrito anteriormente. Otras variantes del funcional (2.3) se obtienen con la incorporación de restricciones a la solución [78]. Otros casos particulares se han presentado en la literatura [25, 8, 3, 54], destacando el caso que corresponde con la minimización del valor absoluto de los residuos, también conocido como deconvolución norma L, y del cual se hablará en la próxima sección Deconvolución norma L La DD norma L fue introducida en [25] como alternativa a la deconvolución de LS. Al igual que LS, el método de mínimos absolutos, o norma L, fue introducido hace más de dos siglos y su incorporación en la solución de problemas de DD nace junto con el desarrollo de otras disciplinas, como por ejemplo las investigaciones operativas y las técnicas de optimización. En [25] se propuso la solución de (2.3), con λ = y p =, utilizando programación lineal, es decir la minimización del valor absoluto del error sin regularización. En [3] se incorporó un factor de regularización, λ, quedando el funcional a minimizar como: J L = e + λ ˆx (2.32)

50 22 CAPÍTULO 2. DECONVOLUCIÓN DE SEÑALES DISPERSAS sujeto a e = y Hˆx (2.33) donde ˆx es la estimación de la señal dispersa y λ es una constante de regularización. El método utilizaba programación paramétrica para calcular el valor óptimo. En [3] se introdujo el método Simplex como método de solución del problema de deconvolución, pero sin la utilización de un parámetro de regularización, y sin seguir directamente la formulación establecida por el funcional (2.32). Sin embargo, en [8] se utilizó el funcional regularizado (2.32) junto con las restricciones (2.33) y se escribió como problema de PL, donde se buscaba minimizar: e + + e + λ ˆx + + ˆx (2.34) sujeto a [ ] [ ] T H H I I x + x e + e = y (2.35) donde x +, x, e + y e son los valores absolutos de la parte positiva y negativa de los vectores ˆx = x + x y e = e + e, reescritos de esta forma para poder ser compatibles con la formulación estándar del método Simplex. Además de aplicar el método Simplex en la solución de (2.34)-(2.35), se estableció un criterio para la determinación del valor más adecuado de la constante de penalización λ, en función de la relación señal a ruido y de un valor máximo de este parámetro [8]. La deconvolución norma L es uno de los métodos de deconvolución dispersa de señales más robustos, debido a que es menos sensible a la presencia de muestras atípicas. La insensibilidad a los valores atípicos tiene una interpretación estadística simple: métodos robustos están asociados a funciones de densidad con distribución extendida, mientras que la función Gausiana tiene una distribución concentrada [28]. Una generalización de la deconvolución norma L para funcionales de norma L p fue presentada en [7, 54], donde el problema se presenta como la solución de un sistema lineal de ecuaciones a través de métodos iterativos a saber: mínimos cuadrados recursivos (IRLS, Iterative Recursive Least Squares) y métodos basados en gradiente conjugado Deconvolución con coste Huber En la Figura 2.4(a) se observa que la función valor absoluto es singular en el origen, es por eso que su minimización se hace difícil debido a este punto de singularidad. Las soluciones propuestas en [3, 8, 3] se basan generalmente en métodos de programación lineal que son computacionalmente muy costosos. Las limitaciones de los métodos LS y norma L han conducido a la utilización de métodos híbridos que combinan las normas L y L 2 en el tratamiento de residuos grandes y pequeños, respectivamente [].

51 2.. DD NO CIEGA 23 La función de coste Huber [52] se define matemáticamente como : e 2 n L H (e n ) = 2ε, e n ε e n ε 2, e n ε (2.36) donde ε es el parámetro que controla la frontera entre la zona cuadrática y la zona lineal. De esta forma, la zona cuadrática resuelve el problema de la singularidad y penaliza cuadráticamente los residuos pequeños, mientras que la zona lineal penaliza los residuos atípicos. En [43], y apoyándose en la función de Huber [52], se formalizó matemáticamente el problema de deconvolución que combinaba el uso de las normas p = y p = 2. El funcional a minimizar viene dado por: J H = e 2 n + 2ε n I n I 2 ( e n ε ) 2 (2.37) donde I e I 2 corresponden a las zonas cuadrática y lineal, respectivamente, determinadas por el valor de los residuos respecto a la constante ε. Debido a que (2.37) es una función no lineal, requiere métodos no lineales para su optimización. En [43] se utilizó un método recursivo no lineal basado en el algoritmo de Newton, donde el valor actual de la señal dispersa viene dado por: x k+ = x k µ k F k J ε(e k ) (2.38) donde es el operador gradiente de derivadas parciales, F k es el Hessiano del funcional (2.37), y µ k es una constante de aprendizaje. La deconvolución con función de coste Huber es un método robusto que incorpora las ventajas de los métodos individuales basados en norma L y L 2. Este método mejora sus prestaciones en presencia de ruido y valores atípicos, manteniendo la suavidad en los pequeños residuos que son penalizados con norma cuadrática. Las prestaciones de la función de coste de Huber servirán de base en los algoritmos que se desarrollarán en el capítulo siguiente Deconvolución de máxima verosimilitud Uno de los métodos de deconvolución más importantes desarrollados en los últimos 25 años ha sido el método de máxima verosimilitud (ML, Maximum Likelihood). La deconvolución de ML se introdujo usando modelos de ecuaciones de espacio de estado [6], y posteriormente se reformuló usando el modelo de convolución [76]. En la Figura 2.5 se muestra el modelo convolucional utilizado en la deconvolución de ML. Pueden identificarse tres componentes principales: la entrada x, la respuesta al impulso h y el ruido aditivo a la salida y. La entrada está formada el producto de las secuencias r y q de distribución Bernoulli y

52 24 CAPÍTULO 2. DECONVOLUCIÓN DE SEÑALES DISPERSAS L(e n) H L (e n) (a) Punto de discontinuidad e - (b) e Figura 2.4: Funciones de coste utilizadas en la deconvolución como problema de optimización: (a) función coste valor absoluto (norma L ); (b) función coste de Huber. Gausiana, respectivamente, más una componente de ruido aditivo Gausiano e B, que modela el efecto de las refracciones [76]. La respuesta al impulso se modela como un proceso ARMA. El ruido a la salida se supone blanco, aditivo, Gausiano y de media cero. La deconvolución de ML modifica el modelo para la señal de entrada, lo que lleva a reescrir (2.8) como sigue: y = HQr + He B + e (2.39) donde Q es la matriz diagonal Q = diag[q, q 2,..., q N ] formada por los elementos de la secuencia Bernoulli q, el vector columna r está formado por los elementos de la secuencia Gausiana, e B es un vector columna con las muestras de ruido aditivo Gausiano a la entrada y e es el vector columna con las muestras de ruido aditivo Gausiano a la salida. La relación directamente proporcional que existe entre la función de verosimilitud y la función de probabilidad [76] se representa por: L { h, s, q, r, e B y } = p(y, q, r, e B h, s) (2.4) donde h, q, r, e B son los parámetros del modelo descritos en la Figura 2.5, y s = {σ q, σ B, σ e, λ} es un vector con los parámetros estadísticos de las diferentes procesos aleatorios. A partir de (2.4), y utilizando el teorema de Bayes, se formula el funcional de la deconvolución de ML [76]. Para el caso en el que la respuesta al impulso es conocida, el funcional a minimizar queda como: J MLD = 2σ r r 2 2 2σ B e B 2 2 2σ n y HQr He B cte (2.4) donde el término cte depende de la distribución de la secuencia Bernoulli [76]. La optimización del funcional (2.4) se lleva a cabo a través de algoritmos de búsqueda exhaustiva. Durante la primera mitad de la década de 98 se

53 2.. DD NO CIEGA 25 (r. q) n + e B n x n h n en + y n Figura 2.5: Modelo convolucional utilizado en la formulación del problema de deconvolución de ML. La entrada esta forma de la suma de dos señales aleatorias, una de estadística BG y la otra de estadística Gausiana. desarrollaron muchos algoritmos de deconvolución de ML, entre los cuales destaca el algoritmo basado en la detección y reemplazo de la espícula más probable (SMLR, Single Most-Likely Replacement Detector). La optimización basada en el algoritmo SMLR ha sido una de las claves más importantes en la optimización del funcional (2.4), pero su implementación es computacionalmente muy costosa, y además presenta mínimos locales [57] Deconvolución mediante mezcla de Gausianas La deconvolución basada en la mezcla de Gausianas (MG) es un método basado en la minimización de la norma L 2 del error y un término de penalización [5]. El término de penalización se obtiene a partir de la modelización de la señal dispersa como una mezcla de dos Gausianas ponderadas de media nula y con distintas varianzas. La distribución Gausiana de mayor varianza, σ b, modela las muestras correspondientes a las espículas de la señal dispersa y la distribución Gausiana de menor varianza, σ n, modela las muestras nulas de la misma. La distribución de probabilidad de una señal con distribución MG se representa a través de: p(x) = π n 2πσn e x 2 2σ 2 n + π b 2 e x 2σ b 2 (2.42) 2πσb donde p(x) es la función densidad de probabilidad de la variable x, σ n y σ b son las desviaciones estándar de las Gausianas estrecha y ancha, respectivamente; y π n y π b son las probabilidades a priori de las Gausianas, cuya suma es igual a uno. A partir del modelo (2.42), y utilizando estimación MAP, se deriva el funcional para la deconvolución de MG como sigue: J MG (x) = z Hx 2 2 α MG N i= log 2 π j p j (x i ) (2.43) donde α MG es un parámetro de regularización que se determina con métodos heurísticos, por ejemplo, utilizando el modelo de señal, el grado de dispersividad o la varianza del ruido [5]. j=

54 26 CAPÍTULO 2. DECONVOLUCIÓN DE SEÑALES DISPERSAS El método de MG puede interpretarse como un modelo general, que tiene como casos particulares la distribución Gausiana del término de penalización (π n = π b y σ n = σ b ) y la Bernoulli-Gausiana (σ n = ) de la deconvolución ML propuesto en [6, 76]. En [5] se propusieron dos métodos para la minimización de (2.43). El primero actualiza la señal muestra a muestra, mientras que el segundo actualiza toda la señal de forma simultánea. La deconvolución MG presenta problemas en la convergencia debido a la alta sensibilidad de sus algoritmos de optimización, y sobre esto se hablará en la sección de experimentos del siguiente capítulo Método de la ventana iterativa La deconvolución mediante la utilización iterativa de ventanas (IWM, Iterative Window Method) surgió como una modificación al algoritmo SLMR, y fue introducido por [56], basándose en la maximización de la densidad de probabilidad a posteriori, esto es: p(q r, y) p(y q, r) p(q r) p(r) (2.44) donde q y r son los vectores de amplitud y posición, respectivamente, de la señal dispersa. El método de deconvolución se basa en la maximización de (2.44), o estimación MAP, a través de la selección adecuada de los valores de q y r. Para cada valor de r se calculan los valores óptimos de ˆq, mediante la solución de un sistema lineal de ecuaciones. Sin embargo, p(ˆq, r y) es una función no lineal, y su maximización con respecto a r se lleva a cabo iterativamente. La búsqueda empieza comparando un valor de referencia r con un número de vecinos, dentro de una ventana. Si alguno de los vecinos incrementa el valor de p(ˆq, r y) se adopta este nuevo valor de r como valor de referencia, la búsqueda se detiene cuando ningún vecino mejora el valor de probabilidad a posteriori. A continuación, se selecciona otra ventana y se repite el procedimiento. Las iteraciones se detienen cuando ya no hay ningún cambio ni en la amplitud ni en la localización de las espículas. Si el número de observaciones es grande, la estimación de las amplitudes se vuelve computacionalmente muy costosa. El método IWM reduce el coste computacional calculando las estimaciones sobre un subconjunto de componentes en cada iteración. La maximización se restringe a una ventana donde se busca localmente la posición de las espículas que maximizan (2.44). Así la deconvolución mediante IWM se formula partiendo r en dos partes, una para los componentes dentro de la ventana (r w ) y otra para los que están fuera (r w ). Lo mismo se hace con la señal dispersa y la matriz de convolución del filtro; H w representa las columnas de H que corresponden a las espículas dentro de la ventana y H w representa las columnas correspondientes a las espículas fuera de la ventana.

55 2.2. DD CIEGA 27 El valor del funcional a minimizar viene dado por: J IW (r w ) = (v w ) T ( (H w ) T H w + γi ) v w + g(r w ) (2.45) donde v w = (H w ) T y (H w ) T H w q w y el término g(r w ) se calculan a partir de los estadísticos de la señal dispersa [56]. g(r) = 2σ 2 eln ( p(r) ) σ 2 eln ( 2πσ q ) M (2.46) donde σ e, representa la varianza del ruido de salida e, σ q representan la varianza de la magnitud de la señal dispersa de entrada q, respectivamente, y M representa el número total de espículas [56, 57] A pesar de las mejoras introducidas por la deconvolución IWM, su coste computacional sigue siendo alto, debido principalmente a la búsqueda iterativa del modelo que máximiza la probabilidad a posteriori (2.44). También cabe decir que IWM ha sido comparado con otros algoritmos utilizando modelos de señal sintéticos, pero muy pocas pruebas se han realizado con señales reales DD ciega Los métodos de deconvolución ciegos parten del desconocimiento de la respuesta al impulso. Para poder realizar la deconvolución es necesario conocer, además del conjunto de observaciones, información a priori sobre la señal a procesar. Los métodos LS, como por ejemplo el filtro de Wiener y la DP, suelen considerarse a veces como algoritmos ciegos, aunque en su formulación aparece la respuesta al impulso como conocida 2. En esta sección se revisarán brevemente tres de las técnicas ciegas más relevantes que han surgido en la literatura Técnicas de DH La técnica de DH surgió casi simultáneamente de dos fuentes diferentes en la década de 96 [, 83], pero fue en [86, 83] donde se sentaron las bases teóricas y en [86, 7, 34] se desarrollaron aplicaciones de lo que se llamó superposición generalizada, y donde la DH es un caso particular. El principio de la DH consiste en convertir las operaciones de convolución en sumas, para lo cual se recurre a un operador que combina varias operaciones, y se utiliza un modelo de observaciones que no incluye el ruido blanco aditivo a la salida 2 Tradicionalmente, la mayoría de métodos LS resuelven el desconocimiento de la respuesta al impulso asumiendo que la señal dispersa se puede modelar como ruido aleatorio blanco, ya que con esta suposición se obtiene que la autocorrelación de las observaciones es igual a la de la respuesta al impulso, excepto por un factor de escala. Se ha demostrado que tal suposición sobre la señal dispersa carece de validez [2], y la necesidad de conocer la respuesta al impulso surge de forma natural [57].

56 28 CAPÍTULO 2. DECONVOLUCIÓN DE SEÑALES DISPERSAS del modelo convolucional (2.8). Sobre este modelo se aplica la secuencia de operaciones siguiente: la transformada Z (TZ), el logaritmo complejo y la TZ inversa. Este conjunto de operaciones traslada las observaciones a un dominio diferente, llamado dominio cepstral, quedando la expresión para las observaciones como: ˇy = ˇx + ȟ (2.47) donde ˇx y ȟ representan la transformada Z inversa de log[x(z)] y log[h(z)], respectivamente. Teóricamente, el método de DH permite separar señales que se han originado a partir de una convolución pero en la práctica el método presenta algunos problemas. Quizá el problema que más se ha tratado en la literatura es el de desenrrollado de la fase [26, 85], y del que siguen publicándose métodos mejores y más eficientes [6, 6]. Otro problema de la DH es el cálculo del logaritmo de valores cercanos o iguales a cero, muy común en señales limitadas en banda. Algunas soluciones proponen ampliar el ancho de banda de la señal [34] o construir un espacio de señal de banda completa, donde no existan valores en frecuencia de la señal iguales o cercanos a cero [7, 72, 4]. Sobre estas soluciones y la aplicación de las mismas se hablará en los Capítulos 4 y 5 de la presente Tesis Algoritmo de módulos constantes El algoritmo de deconvolución basado en el criterio de módulos constantes (MC) fue introducido en [4], y en [47] como parte de un conjunto más general de algoritmos conocidos como técnicas Bussgang [47]. En [4] se propuso la minimización de la función de coste de dispersión: J p = 2p E { ( y n p γ MC ) 2 } donde la constante γ MC se define como: } } γ MC = E { x n 2p /E { x n p (2.48) (2.49) de tal forma que se asegura la existencia de un mínimo local para el funcional (2.48). El caso particular de p = 2 se conoce como criterio de módulos constantes, y la sustitución de (2.49) en (2.48) conduce al funcional: { J MC = ( 4 E y n 2 E{ x n 4 } ) } (2.5) σx 2 El funcional (2.5) muestra como el criterio de módulos constantes penaliza el cuadrado del módulo de salida y n 2 que se alejan de la constante γ MC = E{ x n 4 }/σ 2 x.

57 2.3. CONSIDERACIONES SOBRE LA SEÑAL DISPERSA Métodos basados en aprendizaje máquina En la presente revisión del estado del arte no se encontró ningún método SVM en DD, siendo el método más parecido el de igualación ciega tipo bussgang propuesto en [6]. En [6] se considera el problema de una secuencia de símbolos, pertenecientes a un alfabeto, por ejemplo {s i ±} para el caso binario, de módulos constantes { s i = }, como señal de entrada de un SLIT, y a cuya salida se agrega ruido blanco aditivo Gausiano. Formalmente se puede enunciar este problema como sigue: supónganse que, a partir del conjunto de N observaciones en la salida del canal (y,..., y N ), donde la observación y i = (y i, y i,..., y i M+ ) T es de dimensión M. Se desea recuperar la propiedad de módulos constantes a través de la aplicación del filtro: (w T y i ) 2 = z 2 i = (2.5) donde i =,..., N, y w es el vector de coeficientes del filtro de longitud M. El problema se puede formular como la minimización del funcional: J(w) = 2 w 2 + C N (w T y i ) 2 ε (2.52) i= donde C > es una constante de penalización y (w T y i ) 2 ε es la función de coste ε-insensible [42] definida en el capítulo siguiente. La función de coste establece un compromiso entre la complejidad del igualador y el nivel de desviación por encima de la salida módulo constante unitaria. Dentro de la formulación SVM, el funcional (2.52) se reescribe como un problema de optimización con restricciones, a través de la introducción de variables positivas de holgura, y se aplica la técnica de los multiplicadores de Lagrange para encontrar su solución [6] Consideraciones sobre la señal dispersa La mayoría de métodos de DD desarrollados hasta la fecha suponen que la señal dispersa puede modelarse como un señal aleatoria de ruido blanco. Esta suposición permite estimar la función de autocorrelación de la respuesta al impulso a partir de las observaciones. En trabajos más recientes se demostró que los modelos de señal dispersa asociados a ciertos fenómenos físicos no tienen un comportamiento aleatorio de ruido blanco [2, 43]. De hecho, en [57] se alentó a abandonar completamente la suposición de ruido blanco y orientar los esfuerzos a la medición precisa de la onda. Esta posición generó mucho debate en su momento [89] y estimuló la investigación de nuevos modelos de señal dispersa. En [4] se formuló el problema de deconvolución considerando la señal dispersa como la convolución de dos señales: y[n] = h[n](x m [n] x a [n]) (2.53)

58 3 CAPÍTULO 2. DECONVOLUCIÓN DE SEÑALES DISPERSAS y n H + ~ x n Filtro Fourier Filtro Wavelet ^x n Figura 2.6: Esquema que ilustra la deconvolución con regularización sobre el ruido. El filtro de Fourier atenúa el ruido amplificado durante la operación de inversión y el filtro Wavelet atenúa el ruido residual. donde x m [n] es la componente de fase mínima y ruido no blanco y x a [n] es la componente paso todo y de ruido blanco. El filtro de deconvolución estimado a partir de (2.53) viene dado por: f[n] = (h[n] x m [n]) (2.54) y al aplicarlo a las observaciones se obtiene como resultado solamente la componente paso todo de ruido blanco de la señal dispersa, x a [n]. Al usar (2.54) como filtro de deconvolución en problemas de DD se obtiene una estimación de la señal dispersa con características de ruido blanco. Esto no significa que la señal estimada posea realmente esa característica, sino que (2.54) extrae únicamente la componente blanca de la señal dispersa. La suposición de ruido blanco condujo a resultados equivocados [2, 4, 57]. En [3, 4] se propuso modelar la señal dispersa como ruido aleatorio fraccionadamente integrado (FIN, Fractionally Integrated Noise). A partir de éste nuevo modelo, y de considerar que la señal dispersa se divide en una componente coloreada de fase mínima y otra blanca de fase no mínima, se generalizó la deconvolución mediante filtro de Wiener. La idea de introducir el ruido aleatorio FIN ha sido una contribución importante en la mejora de los métodos de deconvolución. Sin embargo, el modelo tiene como principal limitación el no incluir el ruido aditivo a la salida del modelo de convolución Consideraciones sobre el ruido a la salida En general, muchos métodos de deconvolución no consideran el efecto que el ruido coloreado, originado a partir de la aplicación del filtro de deconvolución a las observaciones, tiene sobre la estimación de la señal dispersa. En [79] se introdujo un método híbrido de deconvolución no ciego, que regulariza el ruido a la salida del modelo convolucional a través de la combinación de un doble filtrado. Los filtros se obtienen utilizando el criterio de Wiener en los dominios de frecuencia y Wavelets. Si se aplica la inversa o pseudoinversa de la matrix de convolución a (2.8) se obtiene una mala estimación de la señal dispersa debido al término de

59 2.4. CONSIDERACIONES SOBRE EL RUIDO A LA SALIDA 3 ruido coloreado, como se muestra a continuación: ˆx = H + Hx + H + e ˆx = x + e (2.55) donde e = H + e es la componente de ruido coloreado y presenta una varianza muy grande cuando la matriz H esta mal condicionada, originando un error cuadrático medio entre x y ˆx muy alto. La estimación de la señal dispersa x a partir de ˆx se lleva acabo a través de la representación adecuada de esta última en un dominio que permite minimizar el efecto del ruido coloreado. Así, dado un conjunto de funciones base ortogonales {z k } N k= la estimación preliminar de la señal dispersa puede representarse como: ˆx = N k= ( ) x, z k + H e, z k z k (2.56) Una mejor estimación de la señal dispersa se consigue si se modifica cada una de las componentes de ˆx, a través de la multiplicación por un escalar λ f k. La nueva ecuación queda como: ˆx = N k= ( ) x, z k + H e, z k λ f k z k (2.57) donde λ f k se interpreta como una regularización del problema de deconvolución. El valor de la constante de regularización λ f k se escoge como compromiso entre la distorsión de la señal x y la atenuación del ruido coloreado e. Si el conjunto de funciones base esta dado por funciones exponenciales complejas, la estimación obtenida a través de (2.57) se interpreta como una operación de un filtro en el dominio de la frecuencia. Los coeficientes del filtro están dados por los diferentes valores de las constantes de regularización λ f k. El método es equivalente a la deconvolución basada en la DFT y por tanto la señal dispersa estimada se ve seriamente distorsionada debido al mal condicionamiento de la matriz de convolución. Como solución a este problema, en [79] se propone una nueva regularización a través de la representación de (2.57) mediante transformada Wavelet. La representación de (2.57) en el dominio Wavelet y su regularización a través de la multiplicación de cada una de sus componentes por un escalar λ w k, de forma similar a como se hizo en el dominio de la frecuencia, conduce a mejores estimaciones de la señal dispersa [79]. Sin embargo, es muy difícil escoger los parámetros adecuados de regularización de los filtros de posprocesado. Esta limitación es muy común en los métodos de deconvolución que utilizan regularización.

60 32 CAPÍTULO 2. DECONVOLUCIÓN DE SEÑALES DISPERSAS 2.5. Consideraciones sobre la fase El concepto de fase mínima surge dentro del análisis de circuitos a partir del desarrollo de amplificadores realimentados, y posteriormente se introduce en la teoría de predicción de series temporales [9]. La idea parte del concepto de lazo de realimentación en una red eléctrica donde la salida de ésta depende de la entrada y de la señal de realimentación. El proceso de realimentación tiene que darse a cierta velocidad, de tal forma que la salida del sistema sea estable. La velocidad a la cual se da la realimentación es controlada por la función de fase, y ésta tiene que tener un retraso de fase mínimo para que el sistema realimentado sea estable [2]. Desde el punto de vista del dominio de la frecuencia, un sistema de fase mínima, dentro de un conjunto de sistemas con igual amplitud para la respuesta en frecuencia, es aquél para el cual la componente de fase posee el menor desplazamiento de fase en función de la frecuencia [33]. Por otro lado, desde el punto de vista temporal, un sistema de fase mínima es aquel para el cual la energía de la respuesta al impulso está muy concentrada al inicio, y en el caso de señales discretas se traduce en valores grandes para los primeros coeficientes. Así, una señal de fase mínima tiene la forma decreciente exponencialmente con el tiempo, donde la parte mayor está al inicio [2]. Desde un punto de vista geométrico, la condición de fase mínima se interpreta a partir de funciones de transferencia racionales. Para que un SLIT sea de fase mínima sus polos y ceros deben estar dentro de la circunferencia de radio unidad. En deconvolución, el concepto de fase mínima se introdujo en [3] como una forma de garantizar que el filtro de predicción de la DP fuese estable. A partir de ese trabajo, la mayoría de métodos de deconvolución consideran la respuesta al impulso del SLIT como de fase mínima [3]. En [37] se utilizó el método de la DH en la deconvolución de señales con respuesta al impulso de fase no mínima, demostrándose que, bajo ciertas condiciones, era posible deconvolucionar este tipo de señales. La limitación viene impuesta por la separación entre el primer reflector y la componente de fase mínima de la respuesta al impulso; si estas se solapan no es posible separarles, obteniéndose una estimación distorsionada de la señal dispersa. Hasta la publicación del trabajo en [62], no se sabía si era posible la estimación de la respuesta al impulso de una onda de fase no mínima a partir de un conjunto de observaciones [75], pero bajo la suposición de que la señal dispersa corresponde a una señal aleatoria de ruido blanco con distribución BG, es posible su estimación. Otro método que permite calcular un filtro de deconvolución con independencia de la fase de la respuesta al impulso es la deconvolución de mínima entropía (MED, Minimum Entropy Deconvolution) [52]. La MED calcula un filtro que maximiza una función objetivo que mide el carácter no Gausiano

61 2.6. CONCLUSIONES 33 de la señal, a través de la maximización de la kurtosis de la señal dispersa. El método no hace ninguna suposición sobre la fase de la respuesta al impulso pero algunas dificultades en cuanto a su pobre robustez lo hicieron rápidamente poco atractivo como técnica de deconvolución [52] Conclusiones En la revisión bibliográfica presentada en este capítulo se ha tratado de caracterizar el estado del arte y de la práctica de los métodos de deconvolución de señales dispersas. No cabe duda de que la mayor parte de técnicas se basan mayoritariamente en los métodos de LS. Las propuestas más recientes incorporan variaciones en el modelo de señal dispersa, la fase de la respuesta al impulso y el tipo de ruido a la salida del modelo convolucional. Para los algoritmos no ciegos propuestos en la presente Tesis, resultan de especial interés tres contribuciones realizadas en el campo de la deconvolución. La primera es la formulación introducida en [25], donde el problema de deconvolución se escribe como un problema de optimización. La segunda se refiere a los trabajos que incluyen las restricciones al funcional de optimización, como una forma de regularización [8, 3]. Por último, resulta de valiosa importancia la introducción de funciones de coste robustas realizadas en [43]. Estas tres contribuciones servirán de base en el desarrollo de las nuevas técnicas de DD no ciega propuestas en esta Tesis. Tal como se ha presentado en la presente revisión bibliográfica, las técnicas ciegas de DD, en particular la DH, siguen líneas muy diferentes a las técnicas no ciegas. La interpretación que se le puede dar a algunos de estos algoritmos es que se basan en trasladar la señal de salida a un espacio de funciones donde es posible separar la señal dispersa de la respuesta al impulso. Esta será una de las claves importantes en cuanto a las técnicas ciegas propuestas en la presente Tesis.

62 34 CAPÍTULO 2. DECONVOLUCIÓN DE SEÑALES DISPERSAS

63 Capítulo 3 Deconvolución de señales dispersas mediante SVM En el presente capítulo se desarrollan algoritmos SVM para deconvolución de señales dispersas. Con el objeto de facilitar la introducción de estos algoritmos, se estudian conceptos importantes del aprendizaje máquina. Entre los conceptos introductorios a tratar se pueden destacar: la clasificación y la regresión mediante algoritmos SVM, en particular haciendo énfasis en conceptos y definiciones básicas fundamentales en el área de las máquinas de aprendizaje; los núcleos de Mercer, estudiando en particular los núcleos invariantes al desplazamiento; las soluciones no lineales de los problemas de clasificación y regresión mediante algoritmos SVM, permitiendo la introducción de las formulaciones del problema primal y dual; el marco lineal general de los algoritmos SVM para procesado digital de señal, introduciendo aplicaciones como el análisis espectral y el modelado ARMA; y el problema de interpolación no uniforme de señales ruidosas, revisando la aplicación del algoritmo basado en el MDS. El estudio de estos temas tiene un doble objetivo. El primero es introducir los conceptos que son más relevantes para esta Tesis, relacionados con el aprendizaje mediante SVM. El segundo se refiere a la introducción de modelos de señal que permiten ampliar las capacidades de los algoritmos SVM en aplicaciones de procesado de señal. En este capítulo se proponen dos algoritmos de DD desarrollados a partir de los modelos de señal MPS, y MDS. El MPS fue introducido en [7] y se 35

64 36 CAPÍTULO 3. DECONVOLUCIÓN DE SEÑALES DISPERSAS MEDIANTE SVM ha aplicado a la solución de diferentes problemas de procesamiento de señales [2, 7,, ]. Por otro lado, el MDS fue introducido un poco más tarde [8] en la solución de problemas de interpolación. El capítulo se divide en seis secciones. La primera sección introduce los conceptos más importantes en SVM mediante la solución de los problemas de clasificación y regresión, para los casos lineal y no lineal. El caso no lineal se trata mediante el estudio previo de los núcleos de Mercer. En la segunda sección se estudia el MPS y sus aplicaciones en la solución de problemas de procesado de señal. La tercera sección presenta el problema de interpolación mediante SVM, y que permite desarrollar un nuevo algoritmo basado en el MDS. La cuarta sección presenta los resultados más novedosos de la presente Tesis, esto es, los algoritmos SVM para DD. La quinta sección muestra los resultados obtenidos en DD a partir de la comparación de los nuevos métodos SVM propuestos con algoritmos de DD ampliamente utilizados en la literatura. Los experimentos se realizan para diferentes tipos de ruido, con señales disperas deterministas y aleatorias y con respuesta al impulso de fase mínima y no mínima. También se ha incluido una aplicación práctica que resuelve el problema de DD en un problema de geofísica. 3.. SVM en clasificación y regresión En esta sección se describen dos problemas del aprendizaje máquina, la clasificación y la regresión, y su resolución mediante algoritmos SVM. Ambos problemas han motivado el desarrollo teórico y algorítmico de técnicas estadísticas y matemáticas que permiten su solución de una forma eficiente. La solución de ambos problemas no es ni mucho menos un tema cerrado, y nuevos resultados se siguen presentando continuamente en la literatura El clasificador binario con clases linealmente separables El problema de clasificación más simple es aquél en el que hay que decidir entre dos clases claramente diferenciables en el sentido de que pueden separarse mediante un hiperplano. A este problema se le llama clasificación binaria con conjuntos linealmente separables, y formalmente puede enunciarse como sigue. Sea el conjunto de muestras de entrenamiento: S = {(x, y ),..., (x N, y N )} (X Y) N (3.) donde N es el número de muestras, y X e Y denotan los espacios de entrada y salida, respectivamente. Además, los vectores x i X R n del espacio de entrada se denominan ejemplos de entrenamiento, mientras que los valores y i Y = {, } del espacio de salida se conocen como etiquetas.

65 3.. SVM EN CLASIFICACIÓN Y REGRESIÓN 37 Sea una función real f() : X R n R, la entrada x = [x,..., x n ] T pertenece a la clase positiva si f(x), y a la clase negativa en caso contrario. La función f(x) es una función lineal con respecto a x X, y corresponde a la ecuación de un hiperplano que se escribe como: f(x) = w, x + b (3.2) donde el vector de pesos w R n y el sesgo b R son los parámetros que controlan la función f(x), y deben de calcularse a partir de los datos, y, representa el operador producto escalar. La regla de decisión está dada por la función sgn(f(x)). Una interpretación geométrica del problema de clasificación binaria es la que considera al espacio de entrada X partido en dos por un hiperplano definido por la ecuación w, x + b = [29]. En la Figura 3.(a) se muestra cómo el hiperplano definido por una recta divide el espacio de entrada R 2 en dos regiones, separando así el espacio de entrada en dos clases distintas. Bajo esta interpretación geométrica, ya que el vector w define un vector perpendicular al hiperplano, y la constante b mueve al hiperplano en una dirección paralela, es necesario calcular n + parámetros. Existen diferentes criterios para seleccionar el hiperplano separador. En esta sección se presenta el criterio de máximo margen, pues es el criterio en el que se basan los algoritmos SVM para clasificación. Se define el margen como la distancia entre el hiperplano y la muestra más cercana x i, y se requiere que la ecuación del hiperplano sea w T x i + b =. Puede demostrarse que la distancia entre la muestra y el hiperplano está dada por [29]: d = w (3.3) donde d es la distancia entre la muestra x i y el hiperplano, y el problema de maximizar el margen es equivalente a minimizar la norma w. Las Figuras 3.(a), 3.(b) y 3.(c) muestran dos hiperplanos separadores para un mismo conjunto de muestras. En general para un conjunto de entrenamiento con muestras linealmente separables existen infinito número de hiperplanos separadores. De este número infinito de hiperplanos el clasificador de margen máximo corresponde al hiperplano con margen geométrico máximo, tal como se observa en la Figura 3.(c), sujeto a las restricciones: y i [ w, xi + b ] (3.4) donde i =, 2,..., N. Dado que el margen geométrico es igual al inverso de la norma, ζ = / w, el problema puede escribirse como un problema clásico de programación cuadrática con restricciones de desigualdad, donde se busca minimizar la norma cuadrática del vector lineal del hiperplano separador sujeta a la restricción (3.4). Éstos problemas de optimización se resuelven

66 38 CAPÍTULO 3. DECONVOLUCIÓN DE SEÑALES DISPERSAS MEDIANTE SVM b w i _ w (a) (b) (c) Figura 3.: El problema de clasificación binario de clases linealmente separables: (a) hiperplano separador (w,b); (b) margen de una muestra al hiperplano; (c) margen máximo del hiperplano. usando la técnica de los multiplicadores de Lagrange [29], de donde se obtiene el funcional Lagrangiano siguiente: 2 wt w N i= [ α i {y i w, xi + b ] } (3.5) donde los α i son los multiplicadores de Lagrange correspondientes a (3.4). La solución de (3.5) conduce al hiperplano de máximo margen con margen geométrico ζ = / w. Como se verá más adelante, el clasificador de máximo margen es la base en la que se apoyan los algoritmos SVM para la construcción de clasificadores más complejos. Además, como se explicará a continuación, el clasificador de máximo margen tiene un gran parecido con el problema de ridge regression, y ha sido este parecido el que ha permitido una de las primeras ampliaciones de las SVM a la solución de problemas de regresión [29] El problema de regresión lineal El problema de regresión lineal consiste en encontrar una función lineal de la forma: f(x) = w, x + b (3.6) que se ajuste lo mejor posible a los datos de entrenamiento S = {(x, y ),..., (x N, y N )}, donde a diferencia del conjunto (3.), aquí el espacio de salida Y R. El problema de regresión se corresponde con el ajuste de un hiperplano al conjunto dado de entrenamiento. En la Figura 3.2 se muestra como ejemplo una función de regresión lineal unidimensional, donde la distancia e i es el error para una muestra particular de entrenamiento. La solución por mínimos cuadrados selecciona los parámetros (w, b) que

67 3.. SVM EN CLASIFICACIÓN Y REGRESIÓN 39 minimizan el siguiente funcional: e 2 i = N ( yi w, x i b ) 2 i= (3.7) que se conoce como función cuadrática de pérdidas, y que se reescribe de forma matricial como a continuación: donde y, X y w se definen como: y y 2 y =., X = y N e 2 = ( y X w ) T ( y X w ) x T, x T 2,. x T N, (3.8), w = (wt, b) T (3.9) Obteniendo el gradiente e igualando a cero (3.8) se obtiene la solución por mínimos cuadrados siguiente: w = ( XT X) XT y (3.) A menudo la matriz X T X, de tamaño (N+) (N+), es una matriz mal condicionada y el cálculo de su inversa conduce a problemas de estabilidad [29]. Una forma de enfrentar dicho problema es a través de la introducción de un factor de regularización en el funcional (3.7), como se escribe a continuación: N ( ) 2 λ w, w + w, x i + b y i (3.) i= donde el parámetro λ R + controla el compromiso entre mantener un nivel bajo de pérdidas y una norma de la solución suave. La solución de (3.) se conoce como ridge regression y viene dada por: w = ( XT X + λin+ ) XT y (3.2) donde I N+ es la matriz identidad con el elemento (N+) (N+) igual a cero y cuyos otros elementos de la diagonal son iguales a la constante λ. Obsérvese el parecido entre (3.) y el funcional para la clasificación lineal separable de máximo margen (3.5), y que ambos presentan un funcional formado por dos términos, uno que controla la complejidad de la hipótesis y el otro su exactitud sobre los datos de entrenamiento. Hasta este punto se han introducido los problemas de clasificación y regresión más simples. En la mayoría de problemas prácticos se encuentra que, en el caso de la clasificación, las clases no son linealmente separables y, en el

68 4 CAPÍTULO 3. DECONVOLUCIÓN DE SEÑALES DISPERSAS MEDIANTE SVM y b e i x Figura 3.2: Función de regresión lineal unidimensional. caso de la regresión, no existen hiperplanos lineales que se ajusten de forma apropiada a las observaciones. Las SVM surgieron como algoritmos robustos propuestos para resolver problemas de clasificación, y más tarde se extendieron a la solución de problemas de regresión [42]. La capacidad de las SVM de resolver problemas de clasificación y regresión, que incluyen datos ruidosos, posiblemente de estadística desconocida, y con presencia de valores atípicos, se basa en buscar su solución en un espacio posiblemente diferente del espacio de entrada. Estos espacios se conocen con el nombre de espacios de características, y su incorporación conduce a la utilización de funciones núcleo o kernel. Las funciones núcleo han dotado a las SVM de la capacidad de extender las soluciones a problemas no lineales, preservando las ventajas de las soluciones a problemas lineales Núcleos de Mercer Frecuentemente, las aplicaciones reales de clasificación y regresión requieren hipótesis más complejas que las que se limitan a funciones lineales. Una de estas hipótesis se basa en la idea de que los datos del espacio de entrada X se pueden representar en un espacio de mayor dimensión (espacio de características) H mediante una transformación no lineal φ(): R N H sobre los datos, es decir: x = (x,..., x N ) φ(x) = (φ (x),..., φ NH (x)) (3.3) donde N H es la dimensión del espacio de características H y ésta puede ser de dimensión infinita. Como ejemplo, considérese el problema de separar dos conjuntos de muestras linealmente no separables, tal como se muestra en la Figura 3.3. Claramente, en el espacio de entrada el problema no es linealmente separable,

69 3.. SVM EN CLASIFICACIÓN Y REGRESIÓN 4 ya que no existe ningún plano en R 2 que pueda separar ambas clases. Si se transforman los datos del espacio de entrada utilizando la transformación no lineal siguiente: φ(x) = (x 2, 2x x 2, x 2 2) = (z, z 2, z 3 ) (3.4) entonces el espacio de entrada en R 2 se traslada a R 3. Si se supone que la frontera de separación en R 2 de las clases corresponde a una elipse, ésta se transformará en R 3 en la ecuación de un plano, como se muestra a continuación: ( x ) 2 ( x2 ) 2 z φ : + = a b a + z 3 2 b = (3.5) 2 que corresponde a la ecuación de un plano en R 3. Así, el propósito de la transformación φ(x) es linealizar, en el espacio de características, las estructuras no lineales. Sin embargo, los espacios de características son más complicados y de mayor dimensión que el espacio de entrada y, además, tienen la dificultad añadida de requerir el desarrollo de la operación producto interno en ese espacio [29]. Es posible calcular este producto interno de una forma implícita, solamente con los datos del espacio de entrada. Para ilustrar como se puede realizar este cálculo se desarrolla el producto interno usando la transformación no lineal (3.4), así: φ(x i ), φ(x j ) = (x i, 2x i x i2, x 2 i), (x j, 2x j x j2, x 2 j) = (x i x j + 2x i x i2 x j x j2 + x 2 ix 2 j) = x i, x j 2 = K(x i, x j ) (3.6) donde K(x i, x j ) es la función núcleo. Obsérvese que, para calcular el producto escalar φ(x i ), φ(x j ) en el espacio de características, no es necesario calcular la transformación no lineal φ(x) = (x 2, 2x x 2, x 2 2), sino que se calcula de una forma indirecta a través del producto interno cuadrático de las muestras del espacio de entrada, x i, x j 2. Esta forma indirecta de calcular el producto interno en el espacio de características es muy bien conocida en el mundo de las máquinas de aprendizaje, y se ha denominado el truco del núcleo [8]. Es decir, la función núcleo permitirá calcular el producto escalar K(x i, x j ) = φ(x i ), φ(x j ) = φ(x i ) T φ(x j ) (3.7) directamente a partir de las muestras x i y x j, sin tener que calcular las transformaciones. El uso de funciones núcleo sobre conjuntos finitos de datos conduce a la obtención de matrices de tipo Gram como a continuación: K(x, x ) K(x, x 2 )... K(x, x N ) K(x 2, x ) K(x 2, x 2 )... K(x 2, x N ) G = K(x i, x j ) = (3.8).... K(x N, x ) K(x N, x 2 )... K(x N, x N )

70 42 CAPÍTULO 3. DECONVOLUCIÓN DE SEÑALES DISPERSAS MEDIANTE SVM x 2 b a x z 2 b 2 z 3 a 2 z Figura 3.3: Efecto de la traslación φ(x, x 2 ) = (x 2, 2x x 2, x 2 2) de un espacio de entrada X a un espacio de características H. que son simétricas y definidas positivas, es decir, todos sus autovalores son mayores que cero [29]. La validez de la función núcleo está determinada por las condiciones de Mercer [29]. Formalmente, las condiciones de Mercer establecen que dada una aplicación bivariada K(x i, x j ) continua y simétrica, se puede asegurar la existencia de la aplicación φ : X X H tal que: K(x i, x j ) = φ(x i ), φ(x j ) = λ n φ n (x i )φ n (x j ) (3.9) donde λ n R + y la función de transformación no lineal φ(x) = (φ (x),..., φ NH (x)) permite la generalización del producto interno en un espacio de Hilbert. La expansión es posible si y solo si para toda función g(x) de energía finita, tal que: g 2 (x)dx < (3.2) n= se cumple la condición: K(x i, x j )g(x i )g(x j )dx i dx j (3.2) Las condiciones de Mercer llevan de forma implícita dos propiedades necesarias que deben cumplir las funciones núcleo [29], las cuales se escriben a continuación: La función núcleo de Mercer es siempre simétrica de tal forma que: K(x i, x j ) = φ(x i ), φ(x j ) = φ(x j ), φ(x i ) = K(x j, x i ) (3.22) La función núcleo de Mercer debe cumplir con la desigualdad triangular, es decir debe satisfacer la desigualdad de Cauchy-Schwartz definida por:

71 3.. SVM EN CLASIFICACIÓN Y REGRESIÓN 43 K(x i, x j ) 2 = φ(x i ), φ(x j ) 2 φ(x i ) 2 φ(x j ) 2 = φ(x i ), φ(x i ) φ(x j ), φ(x j ) = K(x i, x i )K(x j, x j ) (3.23) En la siguiente sección se darán algunos ejemplos de las funciones núcleo más utilizadas Ejemplos de funciones núcleo La Tabla 3. presenta un resumen de las funciones núcleo más habitualmente utilizadas, de las cuales la función núcleo lineal K(x i, x j ) = x i, x j es la más simple. También diferentes transformaciones no lineales pueden conducir a la misma función núcleo. Por ejemplo, las transformaciones no lineales siguientes: ( x = (x, x 2 ) R 2 φ(x) = 2 (x 2 x 2 2), ) 2x x 2, (x 2 + x 2 2) R 3 2 x = (x, x 2 ) R 2 φ(x) = (x 2, x x 2, x x 2, x 2 2) R 4 conducen a la función núcleo (3.6) del ejemplo de la Figura 3.3. Otra función núcleo comúnmente utilizada es la función de núcleo polinómico de grado d=2, definida en la Tabla 3.. La función de transformación para este tipo de núcleo es una función no lineal que transforma observaciones en R 2 a un espacio de dimensión 6, en R 5 a través de: φ(x) = (, 2x, 2x 2, 2x x 2, x 2, x 2 2) (3.24) y desarrollando el producto escalar, se tiene: φ(x i ), φ(x j ) = + 2x i x j + 2x i2 x j2 + 2x i x i2 x j x j2 + x 2 ix 2 j + x 2 i2x 2 j2 = ( x i, x j + ) 2 (3.25) En general, la dimensión del espacio de características de la función núcleo polinomial es ( ) N+d d, donde N es la dimensión del espacio de origen [2]. Otra función núcleo ampliamente utilizada es la función de núcleo Gausiano K(x i, x j ) = e x i x j 2 2σ 2, donde σ representa la anchura de la función núcleo. La transformación no lineal φ() para esta función núcleo es desconocida, y la dimensión del espacio de características es infinita [73] Núcleos invariantes al desplazamiento En [56] se propuso el uso de una función núcleo basada en la transformada Wavelet, demostrando que las funciones invariantes al desplazamiento

72 44 CAPÍTULO 3. DECONVOLUCIÓN DE SEÑALES DISPERSAS MEDIANTE SVM Función núcleo K(x i, x j ) = x i, x j K(x i, x j ) = ( x i, x j + ) d K(x i, x j ) = e x i x j 2σ 2 K(x i, x j ) = e x i x j σ K(x i, x j ) = sinc ( (x i x j )σ ) K(x i, x j ) = tanh(γ x i, x j + µ) K(x i, x j ) = x i x j 2 c 2 2 Tipo Lineal Polinómico Gausiano Laplaciano Sinc Sigmoidal Multicuadrática Tabla 3.: Ejemplos de funciones núcleo de Mercer. cumplen con las condiciones de Mercer para funciones núcleo. Las funciones invariantes a desplazamientos forman espacios vectoriales denominados espacios invariantes al desplazamiento. Las funciones invariantes al desplazamiento se encuentran en aplicaciones muy conocidas en análisis de señal como, por ejemplo, la interpolación [38]. La combinación de prefiltrado y muestreo puede escribirse en términos de productos internos [4, 38] como sigue: (h f)(t) t=n = f(t)h(n t)dt = f, ϕ n (3.26) donde ϕ n = h(n t), y el periodo de muestreo es T=. Obsérvese como sustituyendo ϕ(t) por la función sinc se tiene la ecuación de muestreo uniforme. El conjunto de funciones {ϕ n = ϕ(t n) n Z } expande un espacio de Hilbert. En [56] se demostró que las funciones núcleo invariantes al desplazamiento son núcleos de Mercer si cumplen con la condición, necesaria y suficiente, de tener una transformada de Fourier no negativa, es decir: (2π) /2 + t= ϕ(t)e j2πft dt f R (3.27) Un tipo de función núcleo invariante al desplazamiento de mucha importancia, en este capítulo, en el desarrollo de algoritmos para deconvolución en SVM, es la función de ambigüedad de señales finitas. La función de ambigüedad R h n de una señal h n de longitud finita se define como: R h n = h n h n (3.28) La transformada de Fourier de (3.28) es no negativa y por tanto puede utilizarse como una función núcleo de Mercer del tipo invariante al desplazamiento. La función originada como la convolución entre la respuesta al impulso y su inverso temporal es la función de ambigüedad [6], y que no debe confundirse con la función de autocorrelación. Sin embargo, en el presente capítulo se utiliza el término función de autocorrelación como sinónimo de función de ambigüedad.

73 3.. SVM EN CLASIFICACIÓN Y REGRESIÓN 45 wx + b = wx + b = wx + b = - (a) wx + b = wx + b = wx + b = - j wx + b = wx + b = wx + b = - i (b) (c) Figura 3.4: Hiperplano separador calculado mediante SVM para: (a) clases linealmente separables; (b) clases linealmente separables con muestras en el margen; (c) clases no separables linealmente Núcleos de Mercer en clasificación SVM En la Sección 3. se introdujo el problema de clasificación binaria de clases linealmente separables, y se definió el criterio de margen máximo como criterio de optimización. Además, la solución se planteó como un problema de optimización con restricciones que se reescribe a continuación: minimizar 2 w 2 (3.29) ( sujeto a y i w, xi + b ), i =,..., N (3.3) donde el cuadrado de la norma se utiliza con el objetivo de formular un problema de optimización cuadrática. Este problema se conoce como problema primal, y su solución conduce al hiperplano de máximo margen con margen geométrico γ = / w. El Lagrangiano del problema primal se escribe como sigue: 2 wt w N i= [ α i {y i w, xi + b ] } (3.3) Los multiplicadores de Lagrange α i, tales que las restricciones (3.3) se verifican con igualdad, se denominan vectores soporte y pueden interpretarse como los puntos más próximos a la superficie de decisión. La Figura 3.4(a) representa la solución para un problema linealmente separable de dimensión 2, obsérvese la existencia de un hiperplano que maximiza el margen entre ambos conjuntos de clases.

74 46 CAPÍTULO 3. DECONVOLUCIÓN DE SEÑALES DISPERSAS MEDIANTE SVM Derivando parcialmente con respecto a las variables primales w k y b, y aplicando las condiciones Karush-Kuhn-Tucker (KKT), se obtiene el problema dual del problema primal que consiste en: donde α i. maximizar sujeto a N α i 2 i= N i= N y i y j α i α j x i, x j (3.32) j= N y i α i = (3.33) i= El Apéndice A presenta, dentro del marco de un problema de optimización general, un estudio más amplio sobre las condiciones KKT. Las condiciones KKT permiten hacer algunas conclusiones sobre la formulación de los problemas primal y dual. Éstas establecen que en la solución óptima (α i, w, b ), el producto entre variables duales y restricciones es igual a cero, α i [ yi ( w, x i + b ) ] = (3.34) lo que implica que solamente las entradas x i para las que el margen es igual a uno tienen asociados valores de α i diferentes de cero, todos los otros valores de α i son iguales a cero. Para los diferentes algoritmos SVM que se presentan en este capítulo, las formulaciones de los problemas primal y dual se realizarán a través de tablas de doble columna, donde los problemas primal y dual estarán representados en las columnas de la izquierda y derecha, respectivamente. Clases no separables linealmente. En muchos problemas reales, los datos no son linealmente separables, por ejemplo, con datos ruidosos, y no es posible encontrar el hiperplano separador tal que se verifique la restricción (3.3). La solución a este tipo de problemas motivó el uso de medidas más robustas de la distribución del margen. Esta nueva distribución considera posiciones de otros puntos de entrenamiento, además de los cercanos a la frontera de decisión, es decir, se suaviza el margen. Para suavizar el margen se introducen variables de pérdidas ξ i. Dicha modificación conduce a las formulaciones primal y dual mostradas en la parte superior de la Tabla 3.2, donde C es un parámetro de compromiso entre maximizar el margen y minimizar el error de entrenamiento. La Figura 3.4(c) muestra el problema de clasificación no separable con dos muestras i y j mal clasificadas. En general, si x i es clasificada correctamente por el hiperplano y está fuera del margen, ξ i = ; si está correctamente clasificada y en el margen, < ξ i < ; y si está incorrectamente clasificada, ξ i >. Clasificación no lineal. Cuando se introdujo el concepto de función núcleo se dijo que era posible linealizar la estructura de datos a través de transformaciones no lineales a espacios de mayor dimensión. En este espacio es posible

75 3.. SVM EN CLASIFICACIÓN Y REGRESIÓN 47 Problema primal minimizar 2 w 2 + C N i= ξ i sujeto a Problema lineal Problema dual maximizar N i= α i N 2 i,j= α iα j y i y j x i, x j sujeto a y i ( w, x i + b) ξ i ξ i α i C N i= α iy i = Problema primal minimizar 2 w 2 + C N i= ξ i sujeto a Problema no lineal Problema dual maximizar N i= α i N 2 i,j= α iα j y i y j K(x i, x j ) sujeto a y i ( w, φ(x i ) + b) ξ i ξ i α i C N i= α iy i = Tabla 3.2: Problemas primal y dual en clasificación SVM con margen suave lineal y no lineal. aplicar la teoría de clasificación separable a problemas que no lo eran en el espacio de entrada. Así, el problema consiste en encontrar un hiperplano en el espacio de características tal que ( φ(x i ), w + b) > para cada muestra x i con clase y i = y ( φ(x i ), w + b) < para cada muestra x i con clase y i =. La formulación de los problemas primal y dual se muestran en la parte inferior de la Tabla Núcleos de Mercer en regresión SVM A diferencia de los métodos LS, donde se busca el hiperplano f(x) = w, x + b que mejor se ajusta al conjunto de entrenamiento {(x i, y i ), i =,..., N donde x i R N e y i R}, minimizando el error cuadrático medio, el criterio de error del algoritmo SVR busca el hiperplano que mejor se ajusta a los datos, con una tolerancia igual a ε. A la diferencia entre la función estimada y las observaciones se le llama residuo de la salida, y se utiliza como indicador de la exactitud en la aproximación. Por consiguiente es necesario medir la importancia de esta exactitud

76 48 CAPÍTULO 3. DECONVOLUCIÓN DE SEÑALES DISPERSAS MEDIANTE SVM a través de la función de coste. El tipo de función de coste utilizado incide directamente en los resultados de la regresión. La regresión SVM fue propuesta en [42] y utiliza la función de coste ε-insensible, que permite encontrar hiperplanos con valores óptimos para w ignorando al mismo tiempo valores de residuos menores que el parámetro ε. La función de coste ε-insensible se define como: { e i < ε L(e i ) = (3.35) e i ε e i > ε El problema de regresión se escribe como un problema de optimización donde se busca minimizar la norma del hiperplano sujeto a unas restricciones, como se muestra a continuación: minimizar 2 w 2 (3.36) sujeto a y i w, x i b ε (3.37) donde i =, 2,..., N. w, x i + b y i ε (3.38) La formulación del problema (3.36)-(3.38) supone que la optimización es factible, es decir la función f(x i ) aproxima todos los pares (x i, y i ) con tolerancia ε. Esta situación no es muy frecuente en la práctica, y es necesario suavizar el margen a través de la introducción de variables de holgura, quedando el problema de optimización como sigue: minimizar 2 w 2 + C N (ξ i + ξi ) (3.39) sujeto a y i w, x i b ε + ξ i (3.4) i= w, x i + b y i ε + ξ i (3.4) donde C es una constante que determina el compromiso entre el valor óptimo de la norma de w y la función de pérdidas, ξ i y ξ i son las variables de holgura o pérdidas por exceso o por defecto, respectivamente, sobre la solución respecto al valor de ε. Aplicando el teorema de Lagrange y las condiciones KKT se obtiene el problema dual. La formulación del problema dual, aparte de relacionar los multiplicadores de Lagrange con los datos de entrenamiento y los parámetros libres de la SVM, facilita la extensión del algoritmo SVR lineal a su versión no lineal. En la parte superior izquierda de la Tabla 3.3 se ha vuelto a escribir el problema primal (3.39)-(3.4), y a su derecha se encuentra su respectivo problema dual. La expresión de la solución del hiperplano o vector de pesos estimada es: w = N (α i αi )x i (3.42) i=

77 3.. SVM EN CLASIFICACIÓN Y REGRESIÓN 49 (a) H n L (e ) x i j ( x ) j + - (b) * i j - e L H (e n) x j ( ) x i * i * i j - ec e Figura 3.5: Esquema de regresión (figura adaptada de [8]) SVR. Las muestras en el espacio original se trasladan a un espacio de características donde es posible aplicar la SVR lineal. Las muestras fuera del intervalo definido por ε son penalizadas y se convierten en vectores soporte. La penalización se lleva a cabo a través de una función de coste: (a) función ε-insensible; (b) función ε-huber. donde α i y αi son los multiplicadores de Lagrange que corresponden a las restricciones (3.36) y (3.37), respectivamente. Sustituyendo en la ecuación del hiperplano se tiene la expresión conocida como expansión vectores soporte, dada por: N ŷ = f(x) = (α i αi ) x i, x + b (3.43) i= que describe la solución como una combinación lineal de las muestras de entrenamiento x i. Nótese como (3.43) es función solamente de los multiplicadores de Lagrange y el producto interno entre los datos. Es decir, no es necesario calcular explícitamente el vector de pesos w. Tal como se hizo con el problema de clasificación, el algoritmo SVR se puede generalizar a problemas no lineales a través de una transformación no lineal. Esto se lleva a cabo trasladando los datos de entrada a un espacio de características a través de la transformación no lineal φ : X H, y luego resolviendo para el modelo lineal siguiente: ŷ = f(x) = φ(x), w + b (3.44) donde ŷ i son los valores estimados de y i, w es el vector de pesos en el espacio de características y b es el sesgo. Similarmente al caso anterior, el algoritmo SVR solo depende de los productos sobre los datos transformados φ(x). Por lo que solo es necesario conocer K(x i, x).

78 5 CAPÍTULO 3. DECONVOLUCIÓN DE SEÑALES DISPERSAS MEDIANTE SVM Problema primal Problema lineal, coste ε-insensible Problema dual minimizar maximizar 2 w 2 + C N i= (ξ i + ξi ) N N 2 i= j= (α i αi ) x i, x j (α j αj)+ + N i= ((α i αi )y i (α i + αi )ε) sujeto a sujeto a y i w, x i b ε + ξ i w, x i + b y i ε + ξi ξ i, ξi (α i αi ) C N i= (α i αi ) = Problema no lineal, coste ε-insensible Problema primal Problema dual minimizar maximizar 2 w 2 + C N i= (ξ i + ξi ) N N 2 i= j= (α i αi )K(x i, x j )(α j αj)+ + N i= ((α i αi )y i (α i + αi )ε) sujeto a sujeto a y i w, φ(x i ) b ε + ξ i w, φ(x i ) + b y i ε + ξi ξ i, ξi (α i αi ) C N i= (α i αi ) = Problema primal Problema no lineal, coste ε-huber Problema dual minimizar maximizar 2 w 2 + N 2γ i= (ξ2 i + ξi 2 )+ 2 i,j (α i αi ) ( ) K(x i, x j ) + γδ ij (αj αj)+ +C i I 2 (ξ i + ξi ) + N i= ((α i αi )y i (α i + αi )ε) i I 2 γ C2 2 sujeto a y i w, φ(x i ) b ε + ξ i w, φ(x i ) + b y i ε + ξi ξ i, ξi sujeto a (α i αi ) C N i= (α i αi ) = Tabla 3.3: Problemas primal y dual para SVR con función de coste ε- insensible en los casos lineal y no lineal, y con función de coste ε-huber en el caso no lineal.

79 3.2. MARCO LINEAL SVM PARA SEÑALES DE TIEMPO DISCRETO 5 En la parte de en medio de la Tabla 3.3 se muestra la formulación para los problemas primal y dual para el caso no lineal, y de donde la solución SVR viene dada por: ŷ = f(x) = N (α j αj) φ(x j ), φ(x) + b j= A pesar de las buenas prestaciones de la formulación SVR, existen limitaciones debido principalmente a su función de coste que de una forma implícita asume una estadística específica para los residuos del modelo [8]. En muchas aplicaciones resulta más real asumir que las observaciones han sido contaminadas con diferentes fuentes de ruido de estadística diferente. Por ello, es necesaria la inclusión de funciones de coste más robustas. La función ε-huber introducida en [] ha demostrado ser una función de coste robusta, con muy buenas prestaciones en problemas que involucran diferentes tipos de ruido. La función ε-huber se define como:, e i ε L εh (e i ) = 2γ ( e i ε) 2, ε e i e C (3.45) C( e i ε) 2 γc2, e i e C donde ε, C y γ son los parámetros libres de la función de coste y tienen que ser calculados a priori. La función de coste ε-huber combina tres regiones: una región insensible, una cuadrática y una lineal. Las tres regiones permiten tratar diferentes tipos de ruido. Así, la zona insensible ignora errores menores que el valor absoluto ε; la zona de coste cuadrática utiliza la norma L 2 de los errores, apropiada para perturbaciones con estadística Gausiana; y la zona de coste lineal limita el efecto de valores atípicos en la solución [8, 9, ]. Los problemas SVR primal y dual para la función de coste ε-huber se muestran en la parte inferior de la Tabla 3.3. De donde se observa que la formulación del problema dual es muy similar a la obtenida con función de coste ε-insensible excepto por el factor γδ ij, donde δ ij es la función delta de Kronecker. El factor γ añadido a la diagonal principal de la matriz núcleo funciona como parámetro de regularización [], por ello si se compara el funcional del problema dual ε-huber con el obtenido para la función ε- insensible, se observa que el primero es una versión regularizada del último [] Marco lineal SVM para señales de tiempo discreto Durante los últimos años se han planteado muchos algoritmos para problemas de procesado digital de señal mediante el principio de las SVM. Muchas

80 52 CAPÍTULO 3. DECONVOLUCIÓN DE SEÑALES DISPERSAS MEDIANTE SVM de esas aplicaciones se basan en la idea de buscar una relación directa entre los datos y las observaciones, y así poder utilizar el algoritmo SVR. Además, en la mayoría de los casos se aprovecha el truco del núcleo para desarrollar versiones no lineales de los algoritmos propuestos. En esta sección se estudia el marco general para la solución de problemas supervisados de señales de tiempo discreto desarrollados a partir del algoritmo SVR y que se basa en el concepto de MPS introducido en [2]. Sea Y el espacio de Hilbert N dimensional, donde cada uno de sus elementos es una secuencia discreta {y n } de N elementos y que puede representarse de forma vectorial como y = [y, y 2,..., y N ] T. Además, sea el conjunto de vectores {z,..., z P } una base que expande el espacio Y, y donde cada vector de la base se define como z i = (z i,..., z Ni ) T. Cada vector de señal observado y podrá representarse como una combinación lineal de elementos de esta base, más un término de error e = [e,..., e N ] T mediante: y y 2. y N = w z z 2. z N + + w P z P z 2P. z NP + e e 2. e N (3.46) Para un instante de tiempo dado n, el modelo lineal de series temporales puede escribirse como: P y n = w p z np + e n = w, v n + e n (3.47) p= donde w = [w,..., w P ] es el vector de pesos del modelo, que debe estimarse, y v n = [z n,..., z np ] representa el espacio de entrada en el instante n. Los coeficientes del modelo se estiman minimizando una función de coste residual más un término de regularización, es decir, se minimiza: N 2 w 2 + L p (e n ) (3.48) n= Para establecer una conexión entre los residuos y los parámetros, se introducen restricciones adicionales basadas en el modelo (3.47). Siguiendo la misma metodología de desarrollo algorítmico SVM, se presentan en la Tabla 3.4 los problemas primal y dual. A partir del problema primal, derivando parcialmente su Lagrangiano con respecto a las variables primales y aplicando las condiciones KKT, se obtiene una expresión para el vector de pesos en función de los multiplicadores de Lagrange, dada por: N w p = (α n αn)z pn (3.49) n=

81 3.2. MARCO LINEAL SVM PARA SEÑALES DE TIEMPO DISCRETO 53 Problema primal Marco lineal SVM Problema dual minimizar 2 w 2 + 2γ maximizar N n= (ξ2 n + ξn 2 )+ (α 2 α ) T (R v + γi)(α α )+ +C n I 2 (ξ n + ξ n) n I 2 γ C2 2 +(α α ) T y ε T (α + α ) sujeto a y n w T v n ε + ξ n y n + w T v n ε + ξn ξ n, ξn sujeto a α ( ) C Tabla 3.4: Problemas SVM primal y dual para el marco lineal con función de coste ε-huber. Además, del desarrollo del problema dual, y tal como se muestra en la parte derecha de la Tabla 3.4, es posible identificar la matriz de correlaciones del espacio de entrada: R v = v s, v t (3.5) Tras encontrar los multiplicadores de Lagrange, que se representan en forma vectorial mediante α ( ) = [α ( ),..., α ( ) N ]T, donde α ( ) i se utiliza para referirse a cualquiera de los dos multiplicadores α i ó αi, el modelo de series temporales se reescribe como: y n = f(v n ) = N (α l αl ) v l, v n (3.5) l= En la aplicación de este marco lineal general es importante identificar el conjunto de funciones base que expande el espacio de observaciones Y. A través de algunos ejemplos, en la siguiente sección se muestra la solución de algunos problemas clásicos de procesado de señal como casos particulares dentro del marco propuesto en esta sección Análisis espectral SVM El algoritmo SVM para análisis espectral fue propuesto en []. Formalmente se puede enunciar como sigue. Sea {y tn } una serie temporal obtenida a través de muestreo en los correspondientes instantes de tiempo {t,, t N } de una función continua y(t), se puede representar como un modelo de señal basada en sinusoides mediante: P y tn = A i cos(ω i t n + φ i ) + e tn (3.52) i=

82 54 CAPÍTULO 3. DECONVOLUCIÓN DE SEÑALES DISPERSAS MEDIANTE SVM donde los parámetros desconocidos son las amplitudes A i, las fases φ i y las frecuencias ω i para un número P de componentes sinusoidales, y e tn es el error del modelo. Aplicando identidades trigonométricas para la función coseno y definiendo los coeficientes c i = A i cos(φ i ) y d i = A i sin(φ i ), es posible reescribir (3.52) como: y tn = P ( ci cos(ω i t n ) + d i sin(ω i t n ) ) + e tn (3.53) i= Si se hace v n = [cos(ω t n ),..., cos(ω P t n ), sin(ω t n ),..., sin(ω P t n )] T y los coeficientes c T = [c,..., c P ] y d T = [d,..., d P ] se concatenan en w T = [c T, d T ], es posible reescribir (3.53) como: y tn = w, v n + e n (3.54) El problema de estimación espectral se ha reformulado dentro del marco SVM lineal de manera similar a (3.47), con lo cual se puede proceder a la formulación de los problemas primal y dual de forma automática. En la Tabla 3.5 se muestra un algoritmo, en tres pasos, que resuelve el problema de estimación espectral mediante SVM. Como se observa a partir de la reformulación de (3.52) en (3.54), es posible resolver el problema de la estimación espectral dentro del marco lineal general de las SVM. La robustez de este método de estimación espectral ha sido demostrada a través de ejemplos sintéticos y en la solución de problemas de estimación de canal en problemas de comunicaciones digitales [] Identificación de sistemas SVM-ARMA En [] se introdujo el uso de las SVM como método de identificación de sistemas basado en modelado ARMA. Formalmente, el problema puede escribirse como sigue. Sean {x n } y {y n } las secuencias de entrada y salida, respectivamente, de un sistema invariante en el tiempo caracterizado por la ecuación en diferencias siguiente: y n = M Q a i y n i + b j x n j+ + e n (3.55) i= j= donde {a i } y {b j } son los coeficientes AR y MA del sistema, respectivamente, y e n es la componente de ruido del sistema. Se definen los vectores yn T = [y n, y n 2,..., y n M ] y x T n = [x n, x n,..., x n +Q ], y su concatenación como v n = [yn, T x T n] T. Además, se definen los vectores de coeficientes a T = [a, a 2,..., a M ] y b T = [b, b 2,..., b Q ], y su concatenación como w T = [a T, b T ]. Se puede reescribir (3.55) como: y n = w, v n + e n (3.56)

83 3.2. MARCO LINEAL SVM PARA SEÑALES DE TIEMPO DISCRETO 55 Análisis espectral mediante SVM Se construye la matriz de correlaciones del espacio de entrada, la cual está formada por la suma de las matrices R cos y R sin que se definen como: R cos (m, k) = P i= cos(w it m ) cos(w i t k ) R sin (m, k) = P i= sin(w it m ) sin(w i t k ) Utilizando la Tabla 3.4, se formula y resuelve el problema dual, que consiste en maximizar: 2 (α α ) T [R cos + R sin + γi] (α α ) + (α α ) T y ε T (α + α ) sujeto a α ( ) C. Los coeficientes del modelo (3.53) se obtienen a partir de los multiplicadores de Lagrange calculados en el paso anterior, y están dados por: c i = N k= (α k αk ) cos(ω it k ) d i = N k= (α k αk ) sin(ω it k ) Tabla 3.5: Algoritmo propuesto para la estimación espectral con SVM dentro del marco lineal general. De (3.56) se observa que el conjunto de vectores base está formado por réplicas desplazadas de las secuencias de entrada y salida, y será necesario conocer las condiciones iniciales n = k o,..., N, donde k o = máx (M +, Q). La estimación de los coeficientes del modelo se realiza siguiendo los mismos pasos que en el marco lineal general. La formulación de los problemas primal y dual se realiza de manera similar al método general. En la Tabla 3.6 se muestra un algoritmo, en tres pasos, que permite estimar los coeficientes del modelo ARMA de la Ecuación (3.55). Identificación de sistemas SVM-ARMA no lineal. Las capacidades de las SVM se han extendido a la formulación de nuevos algoritmos que incluyen modelado no lineal. Esto es posible a través de la transformación de los datos de entrada a un espacio de características que cumple con las propiedades de los espacios de Hilbert con núcleo reproductor (RKHS, Reproducing Kernel Hilbert Space) [73]. En [73] se hace referencia a dos técnicas de modelado no lineal. Una primera solución, implícita al problema, se obtiene si se aplica una transformación no lineal, φ( ), al espacio de entrada, conduciendo a un modelo más general

84 56 CAPÍTULO 3. DECONVOLUCIÓN DE SEÑALES DISPERSAS MEDIANTE SVM Modelado ARMA mediante SVM Se construyen las matrices de correlación del espacio de entrada y salida, R Q x y R M y, respectivamente, y que se definen como: R Q x (m, k) = Q j= x m j+x k j+ R M y (m, k) = M i= y m iy k i Se formula y resuelve el problema dual que consiste en maximizar: 2 (α α ) T [ R Q x + R M y ] (α α ) + (α α ) T y ε T (α + α ) sujeto a α ( ) C. Los coeficientes del modelo (3.55) se obtiene a partir de los multiplicadores de Lagrange calculados en el paso anterior y están dados por: a i = N n=k o (α n α n)y n i b j = N n=k o (α n α n)x n j+ Tabla 3.6: Algoritmo propuesto para el modelado ARMA con SVM utilizando el marco lineal general. dado por: y n = w, φ(v n ) + e n (3.57) El otro método se basa en un solución explícita a través de la construcción del modelo ARMA en un espacio de características del tipo RKHS que permite un análisis más profundo de la estructura de la serie temporal. El modelo de señal viene dado por: y n = a T φ(y n ) + b T φ(x n ) + e n (3.58) La formulación de los problemas primal y dual para los modelos de señal (3.57) y (3.58), y ejemplos que miden sus prestaciones, se tratan con detalle en [73] SVM para interpolación En [8] se introdujo el uso de las SVM en la solución de problemas de interpolación sobre muestras no uniformes y ruidosas. Además de utilizar el marco lineal desarrollado en la sección anterior, en [9, 8] se introdujo un

85 3.3. SVM PARA INTERPOLACIÓN 57 nuevo algoritmo de interpolación basado en el MDS, que plantea la solución del problema de interpolación como una regresión no lineal sobre los instantes de tiempo de las observaciones, utilizando funciones núcleo adecuadas. En esta sección se presenta el problema de interpolación desde dos enfoques, uno basado en el MPS y otro en el MDS. Desde el punto de vista de procesado de señal, el problema de interpolación trata sobre la reconstrucción de señales a partir de sus muestras. La solución teórica a dicho problema fue formulada de manera independiente por tres fuentes diferentes, y por eso el teorema que da la solución se conoce con el nombre de Whitaker-Shannon- Kotel nikow (WSK), en honor a sus autores [39]. El teorema de muestreo establece que una señal limitada en banda puede reconstruirse a partir de sus muestras si éstas han sido obtenidas a través de muestreo uniforme y con una frecuencia igual a dos veces la frecuencia máxima de la señal [9]. Matemáticamente, la fórmula de reconstrucción se escribe como: x(t) = n= ( n ) x sinc ( W t n ) (3.59) W donde πw representa el ancho de banda de la señal, y las muestras equidistantes de x(t) se interpretan como los coeficientes de un conjunto de funciones base formadas a través de réplicas desplazadas de la función sinc(t) = sin(πt). Es πt decir, desde el marco lineal SVM se tiene que las funciones base que expanden el espacio de salida (espacio de reconstrucción) están dadas por réplicas de la función sinc centradas en los instantes de muestreo, {z n } = {sinc(w (t n))}, para n =,..., N. En algunas aplicaciones es posible suponer que el conjunto de muestras se ha obtenido a través de muestreo uniforme, es decir, los instantes de muestreo forman una rejilla uniforme en una o varias dimensiones. Sin embargo, en muchas situaciones prácticas las observaciones forman un conjunto de datos que corresponden a muestras obtenidas en instantes de tiempo no uniformemente espaciados. En [55] no solo se introdujo la prueba de la unicidad de la recuperación de señales limitadas en banda no uniformemente muestreadas, sino que además se desarrollaron ecuaciones cerradas para su cálculo. Así, solamente bajo la suposición de que la señal es limitada en banda con energía mínima es posible reconstruir la señal original a partir de un conjunto de muestras de tamaño 2W T, no uniformemente muestreadas a través de: { N N } x(t) = x(t n ) b mn sinc(w (t t m )) (3.6) n= m= donde b mn es el (m, n)-ésimo elemento del inverso de la matriz cuyos elementos están dados por sinc(w (t n t m )), con n, m =,..., N. Si se agrupan de manera diferente los términos del doble sumatorio como sigue: { N N } x(t) = b mn x(t n ) sinc(w (t t m )) (3.6) m= n=

86 58 CAPÍTULO 3. DECONVOLUCIÓN DE SEÑALES DISPERSAS MEDIANTE SVM x x x x x x x x x x x x x x xx xx x x xx xx x x x (a) (b) 5 Figura 3.6: El problema de muestreo, la localización del instante de muestreo se ha marcado con una X y el valor de la muestra con un círculo: (a) una función x(t) R ha sido muestreada uniformemente y no uniformemente, paneles superior e inferior, respectivamente; (b) dos tipos de muestreo bidimensional (cartesiano y polar) con muestreo uniforme y no uniforme, paneles superior e inferior, respectivamente. el término entre llaves en (3.6) corresponde a los coeficientes del modelo, a m = N b mn x(t n ) (3.62) n= Reescribiendo (3.62) utilizando notación matricial se tiene: a = S x (3.63) donde a = [a, a 2,..., a N ] T es el vector de coeficientes del modelo, x = [x(t ),..., x(t N )] T es el vector de observaciones, correspondientes al muestreo no uniforme, y la matriz S se define como a continuación: S = sinc[w (t t )]... sinc[w (t t N )] sinc[w (t 2 t )]... sinc[w (t 2 t N )]... sinc[w (t N t )]... sinc[w (t N t N )] (3.64) La solución (3.63) es óptima en el sentido LS, pero requiere la inversión de la matriz S, que normalmente adolece de problemas de mal condicionamiento. Otro problema comúnmente encontrado en la práctica se debe a que las observaciones raramente representan muestras exactas de la señal x(t) debido, entre otros problemas, a la presencia de ruido. Así, para ser más precisos se

87 3.3. SVM PARA INTERPOLACIÓN 59 asume un modelo más general que incluye el efecto de ruido aditivo con una estadística normalmente conocida a priori, esto es: x (t) = x(t) + e(t) = N a m sinc(w (t t m )) + e(t) (3.65) m= de donde el problema de interpolación se puede escribir dentro del marco lineal SVM de forma directa, tal como se detallará a continuación Algoritmo MPS para interpolación La formulación a partir de un conjunto {x (t n ) = x(t n )+e n, n =,..., N} de observaciones de (3.65) dentro del marco SVM lineal fue introducida en la solución de problemas de interpolación en [8]. La ecuación (3.65) se reescribe como: x (t n ) = N a m sinc(w (t n t m )) + e(t n ) = a, v n + e n (3.66) m= donde v n = [sinc(w (t n t )),..., sinc(w (t n t N ))] y a es el vector de coeficientes del modelo, tal como se definió en (3.63). Al igual que en las aplicaciones en identificación de sistemas, a través del modelado ARMA, y la estimación espectral de la sección anterior, se sigue la misma metodología en el desarrollo de los problemas primal y dual. Es decir, se busca minimizar una función de coste robusta, la ε-huber en este caso, sujeta a unas restricciones lineales formuladas a partir del MPS. La formulación de los problemas primal y dual se muestra en la parte superior de la Tabla 3.7. De la formulación del Lagrangiano para el problema primal y de las condiciones KKT se obtiene la matriz del espacio de entrada [8], que está dada por: T(k, m) = N sinc(w (t n t k ))sinc(w (t n t m )) (3.67) n= El problema dual del MPS se resuelve, de nuevo, utilizando técnicas de programación cuadrática, y a partir de su solución se obtienen los multiplicadores α ( ) = [α ( )... α ( ) N ]T con los que se calculan los coeficientes {a j } del modelo N a j = (α i αi )sinc(w (t i t j )) (3.68) i= Obsérvese cómo los coeficientes del modelo son proporcionales a la correlación cruzada entre los multiplicadores de Lagrange y las funciones base.

88 6 CAPÍTULO 3. DECONVOLUCIÓN DE SEÑALES DISPERSAS MEDIANTE SVM Problema primal Modelo Primal de Señal x n = N i= a isinc(w (t n t i )) + e n Problema dual minimizar 2 a 2 + 2γ maximizar N n= (ξ2 n + ξn 2 )+ (α 2 α ) T (T + γi)(α α )+ +C n I 2 (ξ n + ξ n) n I 2 γ C2 2 +(α α ) T x ε T (α α ) sujeto a sujeto a x n N i= a isinc(σ (t n t i )) ε + ξ n α ( ) C N i= a isinc(σ (t n t i )) x n ε + ξ n ξ n, ξ n Problema primal Modelo Dual de Señal x n = w, φ(t n ) + e n Problema dual minimizar 2 w 2 + 2γ maximizar N n= (ξ2 n + ξn 2 )+ (α 2 α ) T (G + γi)(α α )+ +C n I 2 (ξ n + ξ n) n I 2 γ C2 2 +(α α ) T x ε T (α α ) sujeto a x n w, φ(t n ) ε + ξ n w, φ(t n ) + x n ε + ξn ξ n, ξn sujeto a α ( ) C Tabla 3.7: Problemas primal y dual para interpolación mediante SVM, con función de coste ε-huber, en el modelo primal de señal Interpolación no uniforme usando espacios invariantes La generalización del Teorema de Muestreo (3.59) a espacios de reconstrucción más generales fue posible a través de la introducción de los espacios invariantes [4, 39]. Su ampliación a la interpolación no uniforme ha sido desarrollada recientemente [2, 2, 66]. Las primeras consideraciones de no uniformidad estaban enfocadas a pequeñas perturbaciones a la hora de obtener las muestras [66]. En trabajos más recientes se ha utilizado la teoría de espacios invariantes, desarrollada en principio para muestreo uniforme, obteniendo una teoría única que incluye interpolación uniforme y no uniforme con ruido en las observaciones [3].

89 3.3. SVM PARA INTERPOLACIÓN 6 Si se hace ϕ k (t) = sinc(t k), x n = w k y T= el teorema de muestreo (3.59) puede reescribirse como: x(t) = k Z w k ϕ k (t) (3.69) donde w k son los coeficientes del modelo y ϕ k (t) representa un conjunto de funciones base que expande un espacio V de aproximación, { V (ϕ(t)) = x(t) = } w k φ k (t) : w k l 2 (3.7) k Z En general, un espacio invariante al desplazamiento es un espacio de funciones en R d [3] y se define por: { N H } V (ϕ,..., ϕ NH ) = wkϕ i i (t k) i= k Z d (3.7) donde las funciones {ϕ i (t k)} pueden formar una base, ya no estrictamente ortogonal, de funciones limitadas o no en banda [39, 38]. Los coeficientes w k óptimos de (3.69) en el sentido de LS y se obtienen a partir de los operadores de proyección ortogonal y proyección oblicua para bases ortogonales y no ortogonales, respectivamente. P V (ϕ) x(t) = k Z x(t), ϕ k (t) ϕ k (t) (3.72) donde ϕ k (t) es el conjunto de funciones base dual de ϕ k (t) que cumplen con la condición de biortogonalidad ϕ k (t), ϕ l (t) = δ k l (t). El producto interno w k = x(t), ϕ k (t) representa las contribuciones de la señal a lo largo de la dirección especificada por ϕ k (t). Éste producto interno, como se muestra en la Figura 3.7, es equivalente a, primero, filtrar la señal de entrada con un filtro paso bajo ideal (filtro antisolapamiento), y luego muestrear la señal filtrada, así: (x h)(t) t=k = x(t)h(k t)dt = x(t), ϕ k (t) (3.73) donde ϕ(t) = h T (t k) y para el caso formulado en (3.59) el prefiltrado y el postfiltrado son ambos filtros paso bajo con respuesta al impulso h(t) = ϕ(t) = ϕ( t) y el producto interno resulta en las muestras uniformemente espaciadas x(k) = x(t), ϕ k (t). La representación de señal utilizando los operadores de proyección es de mucha utilidad cuando se introduce el concepto de RKHS. Si reescribimos (3.72) como:

90 62 CAPÍTULO 3. DECONVOLUCIÓN DE SEÑALES DISPERSAS MEDIANTE SVM x(t) Prefiltrado Muestreo Postfiltrado h(t) (-t) + (t) x(t) (t-n) Figura 3.7: Representación esquemática en tres pasos del teorema de muestreo utilizando espacios invariantes al desplazamiento. El intervalo de muestreo es uniforme con T =. () La señal de entrada x(t) se ha filtrado con un filtro antisolapamiento. (2) El proceso de muestreo conduce a la obtención de los coeficientes w k. (3) La reconstrucción se obtiene a través del filtrado de los coeficientes con ϕ(t). P V (ϕ) x(t) = k Z x( ), ϕ( k) ϕ(t k) (3.74) = x( ), k Z ϕ( k)ϕ(t k) (3.75) = x( ), K(t, ) = x(t) (3.76) de donde para el espacio invariante V (ϕ) especificado por (3.7) el núcleo reproductor es igual a: K(t, s) = k Z ϕ(t k) ϕ(s k) (3.77) El resultado anterior dice que el espacio de reconstrucción V (ϕ) es un RKHS cuyo núcleo reproductor es (3.77). La discusión anterior ha servido para establecer la validez de un nuevo modelo de señal presentado en [8] y que se discutirá a continuación Algoritmo MDS para interpolación Un algoritmo SVM diferente puede obtenerse para la interpolación si se usa el principio de la SVR. El problema se plantea como a continuación. Dado un conjunto de observaciones {x n } en los instantes de tiempo {t n } se traslada los instantes de tiempo a un espacio de características H a través de una transformación no lineal t R φ(t) H. En el espacio de características, se tiene que las transformaciones sobre los instantes temporales aproximan linealmente las observaciones mediante: x n = w, φ(t n ) + e n (3.78)

91 3.4. SVM PARA DECONVOLUCIÓN NO CIEGA 63 donde x n es el conjunto de observaciones de una señal limitada en banda con ruido aditivo Gausiano. La formulación de los problemas primal y dual se muestra en la parte inferior de la Tabla 3.7. De la formulación para el problema dual se identifica la siguiente matriz: G(k, m) = φ(t k ), φ(t m ) = K(t k, t m ) (3.79) donde K(t k, t m ) es un núcleo de Mercer y permite obviar el conocimiento explícito de la transformación no lineal φ( ). El vector de pesos en H está dado por: w = N (α j αj)φ(t n ) (3.8) n= La solución final se expresa como: x n = N η j K(t j, t n ) (3.8) j= donde η j = (α j αj) y el núcleo se define como K(t k, t n ) = sinc(σ (t k t n )). Es posible demostrar que la función núcleo anterior cumple con las propiedades de un núcleo de Mercer [56] y se conoce como núcleo sinc. De (3.8) se observa que es posible utilizar otros tipos ) de núcleo de Mercer. Por ejemplo, si se define K(t k, t n ) = exp ( t k t n 2, tendríamos un interpolador 2σ 2 Gausiano. También, (3.8) puede interpretarse como la convolución entre un filtro cuya respuesta al impulso es la función y la señal de entrada es la secuencia de multiplicadores de Lagrange correspondientes a cada instante de tiempo. Así, si se asume que {η n } son las observaciones de un proceso de tiempo discreto η[n] y que {K(t n )} es la versión de tiempo discreto de una función de autocorrelación dada por K[n] [9, 8], la solución x[n] puede escribirse como: x[n] = η[n] K[n] (3.82) Obsérvese que (3.82) tiene validez en el caso de la interpolación con núcleo sinc, ya que la función núcleo sinc cumple con la condición de tener una transformada de Fourier no negativa [56] SVM para deconvolución no ciega En esta sección se desarrolla la principal aportación de la presente Tesis en DD no ciega. Se proponen dos algoritmos SVM desarrollados a partir de los modelos MPS Y MDS analizados en las secciones 3.2 y 3.3. Tal como se

92 64 CAPÍTULO 3. DECONVOLUCIÓN DE SEÑALES DISPERSAS MEDIANTE SVM describe en esta sección, la DD no ciega mediante SVM presenta ventajas sobre los métodos hasta ahora conocidos, debido fundamentalmente a las siguientes cuestiones: la solución se expresa en función de un subconjunto de observaciones, conduciendo a soluciones dispersas; la solución conlleva una regularización implícita ofreciendo buenas prestaciones en problemas mal condicionados; y la incorporación de funciones de coste robustas mejora las capacidades de los algoritmos SVM frente a ruido con estadística no Gausiana y de carácter impulsivo. Estas ventajas se amoldan a los requerimientos en la solución de problemas de DD no ciega. Debido a la estrecha relación existente entre los problemas de interpolación y deconvolución [39, 38, 34], la formulación de algoritmos SVM utilizando los modelos de señal MPS y MDS surgen de forma natural. La solución SVM del problema de DD no ciego basada en MPS se obtiene aplicando el marco lineal de la sección 3.2, y de una forma simple se referirá a éste como algoritmo MPS. La formulación de la solución mediante el algoritmo MPS no conduce a una solución dispersa, pero descubre uno de los conceptos fundamentales en DD no ciega mediante SVM: la relación implícita entre la respuesta al impulso y los núcleos de Mercer. Por otra parte, la solución SVM basada en MDS conduce a soluciones dispersas, pero la respuesta al impulso debe cumplir con un requerimiento fundamental; ésta debe ser un núcleo de Mercer. A partir de la combinación de las ventajas de los algoritmos MPS y MDS, se propone y se desarrolla un nuevo algoritmo. Debido a que el algoritmo se basa en una modificación previa de las observaciones, que consiste en filtrado con el inverso temporal de la respuesta al impulso, se ha denominado con el nombre de algoritmo basado en un Modelo de Señal con Kernel de Ambigüedad o Autocorrelación (MSKA). La metodología a seguir en la presentación de cada uno de los algoritmos es como sigue: se introduce el modelo de señal correspondiente; se introduce el funcional a optimizar; se formula el problema primal; se formula el funcional Lagrangiano; se formula el problema dual; y

93 3.4. SVM PARA DECONVOLUCIÓN NO CIEGA 65 se obtiene una expresión cerrada para la señal dispersa en función de los multiplicadores de Lagrange. Cada uno de los algoritmos propuestos parten de un modelo de señal diferente, en cambio, el funcional a optimizar puede escribirse de una forma genérica, adaptando los cambios a cada uno de los casos particulares. En el Capítulo 2 se estableció que el problema de DD no ciega puede formularse como un problema de optimización a través del funcional (2.3). Esta formulación permite amoldar el problema de deconvolución dentro del marco SVM, y con ello la introducción del funcional general, con funciones de coste más robustas y con la incorporación de términos de penalización más adecuados. El funcional (2.3) se puede reescribir como sigue: J SV M = N L εh (e n ) + τ 2 2 (3.83) n= donde L εh (e n ) es la función de coste ε-huber de los residuos y τ 2 2 es un término genérico de penalización, definido de forma diferente para cada uno de los algoritmos propuestos. A continuación se desarrollan cada uno de los algoritmos de DD no ciega mediante SVM, utilizando el funcional general (3.83) y sus respectivos modelos de señal Formulación de la DD en el MPS El algoritmo MPS se obtiene reformulando el problema de DD, tal como se hizo en los ejemplos de análisis espectral e identificación de sistemas de las secciones 3.2. y 3.2.2, respectivamente. Siguiendo esa misma línea se reescribe el modelo de señal utilizando la representación algorítmica de la convolución: P y n = x j h n j+ + e n (3.84) j= donde, de nuevo, x j es la señal dispersa de entrada que ha de ser estimada. El conjunto de funciones base se forma a partir de desplazamientos de la respuesta al impulso, la cual es una función invariante en el tiempo. De la ecuación (2.7) se observa que la respuesta al impulso está dada por el vector columna siguiente: h = [ h, h 2,..., h Q,,..., ] T (3.85) de tal manera que el conjunto de funciones base está formado por P funciones correspondientes a cada una de las columnas de la matriz de convolución H, formada por réplicas desplazadas de la respuesta al impulso.

94 66 CAPÍTULO 3. DECONVOLUCIÓN DE SEÑALES DISPERSAS MEDIANTE SVM El espacio de entrada se obtiene de forma similar a como se obtuvo en el ejemplo de identificación de sistemas de la Sección 3.2.2, y está formado por cada una de las filas de la matriz de convolución: v n = [ ] h n, h n, h n 2..., h n (P ) (3.86) donde h n = en los valores de índice P <n<. También, si se considera la señal dispersa como el conjunto de coeficientes de un problema de identificación de sistemas, la ecuación (3.84) puede escribirse como: y n = x, v n + e n (3.87) La estimación de la señal dispersa puede realizarse siguiendo los pasos del marco lineal general de la Sección 3.2 pero, a diferencia de los ejemplos de estimación espectral e identificación de sistemas, en esta sección se desarrolla paso a paso la obtención de los problemas primal y dual, permitiendo evaluar con mayor precisión las capacidades del algoritmo. El algoritmo MPS se obtiene de la optimización del funcional (3.83) sujeto a restricciones lineales determinadas a partir del modelo original, en este caso, la suma de convolución. La función de coste residual es la ε-huber, mientras que el parámetro de regularización está dado por la norma de la estimación de la señal dispersa, y se incorpora en (3.83), haciendo τ = ˆx. Sustituyendo en (3.83), se tiene que el funcional a optimizar está dado por: J MP S = N L εh (e n ) + 2 ˆx 2 2 (3.88) n= Como normalmente se hace en el desarrollo SVM, se sustituye (3.45) en (3.88), conduciendo a la minimización del funcional siguiente: 2 P k= ˆx 2 k + 2γ (ξn 2 + ξn 2 ) + C (ξ n + ξn) γc 2 2 n I n I 2 n I 2 (3.89) con respecto a ˆx k, ξ n y ξ n, donde ξ n y ξ n son variables de holgura, sujeto a las restricciones lineales siguientes: y n y n + P ˆx j h n j+ ε + ξ n (3.9) j= P ˆx j h n j+ ε + ξn (3.9) j= ξ n (3.92) ξ n (3.93) donde n =,,..., N e I e I 2 son las regiones para los cuales los residuos tienen un coste cuadrático y lineal, respectivamente. Obsérvese que las restricciones (3.9) y (3.9) se han escrito utilizando el modelo de señal (3.84)

95 3.4. SVM PARA DECONVOLUCIÓN NO CIEGA 67 pero los mismos resultados se obtienen si se utiliza el modelo de señal (3.87). Utilizando los multiplicadores de Lagrange se obtiene el funcional Lagrangiano siguiente: 2 P k= ˆx 2 k + 2γ N α n ( y n + n= N αn(y n n= (ξn 2 + ξn 2 ) + C (ξ n + ξn) n I n I 2 P j= ˆx j h n j+ + ε + ξ n ) P ˆx j h n j+ + ε + ξn) j= N (β n ξ n + βnξ n) γc 2 n I 2 n= 2 (3.94) donde α n, α n, β n y β n son los multiplicadores de Lagrange correspondientes a (3.9), (3.9), (3.92) y (3.93), respectivamente. Derivando parcialmente el funcional Lagrangiano con respecto a las variables primales ˆx k y aplicando las condiciones KKT se obtiene la siguiente expresión genérica para la señal dispersa: ˆx n = N (α i αi )h i n+ = i= N η i h i n+ (3.95) donde de nuevo η i = α i α i. La Ecuación (3.95) se puede reescribir como sigue: ˆx n = η n h n δ n+q (3.96) donde δ n representa la secuencia impulso unitario de tiempo discreto definida como δ n = para n y δ n = para n = [33]. De (3.96) se observa que la señal dispersa estimada es igual a la convolución entre los multiplicadores de Lagrange y la inversión temporal desplazada de la respuesta al impulso. Sobre este hecho se harán algunos comentarios relevantes al final de esta sección. El problema dual se obtiene derivando parcialmente el funcional Lagrangiano (3.94) con respecto a las variables primales y aplicando las condiciones KKT. Como paso previo, se obtiene la matriz de entrada sustituyendo (3.95) en (3.94) como a continuación: i= R = R(i, l) = P h i k+ h l k+ (3.97) k= El problema dual consiste en maximizar: 2 (α α ) T (R + γi)(α α ) + (α α ) T y ε T (α + α ) (3.98)

96 68 CAPÍTULO 3. DECONVOLUCIÓN DE SEÑALES DISPERSAS MEDIANTE SVM sujeto a α ( ) C (3.99) donde, al igual que en los algoritmos SVM anteriores, los vectores α ( ) = [α ( ), α ( ) 2,, α ( ) N ]T, y = [y, y 2,, y N ] T, y = [,,..., ] T representan el conjunto de multiplicadores de Lagrange, las observaciones y un vector columna de unos de longitud N, respectivamente. Resulta relevante en el análisis del MPS la relación entre los residuos y los multiplicadores de Lagrange, que fue demostrada en [], y está dada por: sign(e n )C, e n e C η n = sign(e n ) γ ( e n ε), ε e n e C (3.), e n < ε De (3.) se observa que se puede controlar convenientemente la dispersión de los multiplicadores de Lagrange a través del parámetro libre ε y, además, el efecto de los valores atípicos está limitado por el valor de la constante C. Utilizando (3.96) y (3.) se puede construir un esquema para el MPS donde los multiplicadores de Lagrange y la estimación de las observaciones son, respectivamente, la entrada y la salida de un sistema lineal e invariante con el tiempo. Esta relación se ilustra en la Figura 3.8, donde el rectángulo a trazo discontinuo muestra la existencia implícita de una función núcleo de ambigüedad, la cual está dado por: K n = R h n = h n h n (3.) En el esquema de la Figura 3.8 se observa que la señal dispersa ˆx n se encuentra entre los SLTI correspondientes a la respuesta al impulso y la repuesta al impulso invertida temporalmente. También el mismo esquema muestra un sistema de retraso de tiempo discreto δ n+q que hace compatible el índice de muestras con el índice de señal, compensando el desplazamiento temporal introducido por la convolución con la respuesta al impulso con inversión temporal. El esquema del MPS de la Figura 3.8 permite hacer tres observaciones importantes. Primero, es posible controlar la dispersión de los multiplicadores de Lagrange a través del parámetro libre ε. Segundo, a pesar de que la secuencia η n es una secuencia dispersa, de acuerdo a la naturaleza de los algoritmos SVM, el valor estimado de la señal de entrada no lo es, debido a que la solución (3.96) incluye el filtrado en sentido anticausal, perdiéndose el carácter disperso de la secuencia η n. Finalmente, esta situación conduce a explorar otros modelos de señal que aprovechen mejor el carácter disperso de los multiplicadores de Lagrange, por ello en la siguiente sección se desarrolla un algoritmo SVM para DD basado en el MDS.

97 3.4. SVM PARA DECONVOLUCIÓN NO CIEGA 69 n n + Q C R h n = Kn e h- n - C x n > y n e n + y n > h n Figura 3.8: Esquema para el MPS. Los multiplicadores de Lagrange y las observaciones son la entrada y la salida, respectivamente, de un SLIT Formulación de la DD en el MDS La formulación del problema de DD no ciega utilizando el MDS se formula siguiendo la misma línea que en la Sección 3.3.3, introducida en [8]. Sea una secuencia {y n } discreta en un espacio de Hilbert de la que se disponen de N observaciones consecutivas, agrupadas en un vector y = [y, y 2,..., y N ] T, en los instantes de tiempo {n}, se aplica una transformación no lineal φ: Z H sobre estos últimos. Dicha transformación establece una correspondencia entre el espacio de entrada, dado por los instantes de tiempo, y un RKHS o espacio característico de dimensión N H. Cada vector de señal observado y puede representarse como una combinación lineal de un conjunto de funciones base {φ(n)} y un término de error e = [e,..., e N ] T. Así, para un instante de tiempo dado n el modelo lineal puede escribirse como: y n = N H j= w j φ j (n) + e n = w, φ(n) + e n (3.2) donde n =,, N, φ j es la componente j ésima de la función de transformación no lineal φ(), w j es el j-ésimo elemento del vector de pesos del modelo w y N H es la dimensión del espacio de características. En general, la dimensión del espacio de características es mucho mayor que la del espacio de entrada. Sin embargo, las SVM permiten trabajar con espacios de dimensiones muy altas. Esto es posible si se incorporan en los algoritmos SVM, como se vio en la Sección 3..3, el uso de núcleos de Mercer [29], lo que facilita la formulación matemática y el cálculo computacional.

98 7 CAPÍTULO 3. DECONVOLUCIÓN DE SEÑALES DISPERSAS MEDIANTE SVM A partir del modelo de señal (3.2) y del funcional generalizado (3.83), es posible construir el funcional para el MDS. El parámetro de regularización en el MDS está dado por la norma al cuadrado del vector de pesos, y se incorpora en (3.83) haciendo τ = w. Sustituyendo en (3.83) se tiene: J MDS = Así, se busca minimizar: N L εh (e n ) + 2 w 2 2 (3.3) n= 2 sujeto a N H k= w 2 k + 2γ (ξn 2 + ξn 2 ) + C (ξ n + ξn) γc 2 2 n I n I 2 n I 2 (3.4) y n y n + N H j= N H j= w j φ j (n) ε + ξ n (3.5) w j φ j (n) ε + ξ n (3.6) ξ n (3.7) ξ n (3.8) Siguiendo la metodología SVM se obtiene, de manera similar a como se hizo con el algoritmo MPS, el funcional Lagrangiano siguiente: 2 N H k= w 2 k + 2γ N ( α n y n + n= N n= α n ( y n N β n ξ n n= n I (ξ 2 n + ξ 2 n ) + C n I 2 (ξ n + ξ n) n I 2 γc 2 2 N H j= N H j= w j φ j (n) + ε + ξ n ) ) w j φ j (n) + ε + ξn N β n ξn (3.9) n= donde α n, α n, β n y β n son los multiplicadores de Lagrange correspondientes a (3.5), (3.6), (3.7) y (3.8), respectivamente. Derivando parcialmente el funcional Lagrangiano con respecto a las variables primales w k se obtiene el vector de pesos w k como función de los multiplicadores de Lagrange, como se escribe a continuación: w k = N (α n αn)φ k (n) (3.) n=

99 3.4. SVM PARA DECONVOLUCIÓN NO CIEGA 7 Sustituyendo (3.) en (3.9) y simplificando, se obtiene la siguiente matriz de Gram: G = G(i, l) = φ(i), φ(l) = K(i, l) (3.) donde la matriz G(i, l) se obtiene evaluando la función núcleo para los valores i, l =, 2,..., N. Sustituyendo (3.), (3.) en el funcional Lagrangiano y utilizando las condiciones KKT se formula el problema dual de (3.4)-(3.8), que consiste en maximizar: [ ] 2 (α α ) T G + γi (α α ) + (α α ) T y ε T (α + α ) (3.2) sujeto a α ( ) C (3.3) Sustituyendo (3.) en el primer término del miembro derecho de (3.2), se obtiene una expresión para la estimación de las observaciones en función de los multiplicadores de Lagrange: ŷ n = = = N H j= i= ( N ) (α i αi )φ j (i) φ j (n) = i= N ( N H ) (α i αi ) φ j (i)φ j (n) = j= N η i K(i, n) (3.4) i= donde, de nuevo, α i, αi son los multiplicadores de Lagrange, e ŷ n es la estimación de la observación n-ésima. De (3.4) y tal como se muestra en la Figura 3.9 se observa que se pueden relacionar los multiplicadores de Lagrange, la función núcleo y las observaciones mediante un SLIT como sigue: ŷ n = η n K n = η n h n = ˆx n h n (3.5) donde la señal dispersa de entrada se define a partir de los multiplicadores de Lagrange, la respuesta al impulso corresponde a la función núcleo y la señal de salida a la estimación de las observaciones. Además, al igual que en el MPS, y tal como se mostró en (3.), el grado de dispersión se controla a través del parámetro ε. Este enfoque requiere que la respuesta al impulso sea compatible con los núcleos de Mercer. En [56] se demostró que los núcleos invariantes K(i, l) = K(i l) son núcleos de Mercer si y solo si su transformada de Fourier es no negativa (véase Sección 3..5). Por ejemplo, el núcleo sinc, K(t k, t n ) =

100 72 CAPÍTULO 3. DECONVOLUCIÓN DE SEÑALES DISPERSAS MEDIANTE SVM n = x n > x n > C e K n - C y n + e n y n > Figura 3.9: Esquema para el MDS, donde se representa la relación entre las observaciones y los multiplicadores de Lagrange. sinc(σ (t k t n )), utilizado en la interpolación con muestreo no uniforme de señales de banda limitada [8], tiene transformada de Fourier real y simétrica y, por tanto, cumple con el criterio de tener una transformada de Fourier no negativa. La respuesta al impulso encontradas en las aplicaciones prácticas tienen una naturaleza causal, y por tanto no son simétricas. Por consiguiente, no cumplen con una de las condiciones necesarias de los núcleos de Mercer. Observando el MPS representado en la Figura 3.8 se tiene que el SLIT que relaciona los multiplicadores de Lagrange y las observaciones es la función de autocorrelación, y como se ha demostrado esto es un núcleo de Mercer. De esta observación surge la posibilidad de aprovechar las ventajas de los modelos MPS y MDS. A continuación se presenta un modelo que aprovecha las ventajas de los métodos precedentes y supera sus limitaciones, introduciendo previamente un preprocesado sobre el conjunto de observaciones Formulación de la DD en el MDS modificado En la sección anterior se determinó que la principal desventaja del algoritmo MDS es el requerimiento de tener una respuesta al impulso que sea un núcleo de Mercer. En esta sección se introduce un nuevo algoritmo que permite resolver la limitación antes citada [9]. A través del preprocesado

101 3.4. SVM PARA DECONVOLUCIÓN NO CIEGA 73 de las observaciones es posible modificar la respuesta al impulso del sistema SLIT, convirtiéndose ésta en la función de autocorrelación, la cual cumple con las condiciones de los núcleos de Mercer [56]. Debido a la incorporación de la función de autocorrelación de la respuesta al impulso en la formulación del algoritmo SVM, se ha denominado a éste algoritmo basado en el MSKA. A partir de los resultados obtenidos en el MPS, se modifica el modelo convolucional filtrando la salida y n con la inversión temporal desplazada de la respuesta al impulso h Q n. La Figura 3.(a) muestra cómo se modifican las observaciones a través de filtrado y cuyo resultado puede verse como la convolución de la señal dispersa con una respuesta al impulso simétrica, es decir: z n = y n h n δ n+q = = ˆx n [ h n h n ] δn+q + e n = = ˆx n R h n δ n+q + e n = = ˆx n K h n δ n+q + e n (3.6) donde e n = e n h n δ n+q son los residuos modificados y K h n es un núcleo invariante. Las observaciones modificadas z n pueden interpretarse como la convolución entre la señal dispersa y el núcleo invariante dado por la autocorrelación de la respuesta al impulso. Este nuevo modelo permite aprovechar el carácter disperso intrínseco en el MDS sin los inconvenientes causados por la falta de simetría de la respuesta al impulso. En ese sentido, las observaciones modificadas z n se utilizan, tal como se hizo en MDS, para construir el modelo de señal modificado: z n = N H i= v i ψ i (n) + e n (3.7) donde n =,, N, ψ i es la componente i ésima de la función de transformación no lineal ψ(), v i es el i-ésimo elemento del vector de pesos del modelo v y N H es la dimensión del espacio de características. A partir del modelo de señal (3.7) y del funcional generalizado (3.83) es posible construir el funcional a optimizar en el algoritmo MSKA, donde el parámetro de regularización está dado por la norma cuadrática del vector de pesos, τ = v, escribiéndose como sigue: J MSKA = N L εh (e n ) + v 2 2 (3.8) n= El problema primal consiste en la minimización de: 2 N H k= v 2 k + 2γ (ξn 2 + ξn 2 ) + C (ξ n + ξn) γc 2 2 n I n I 2 n I 2 (3.9)

102 74 CAPÍTULO 3. DECONVOLUCIÓN DE SEÑALES DISPERSAS MEDIANTE SVM x n y n z n x n h n e n + y n h Q -n z n (a) n = x n > n + Q C - C e K n = R n h h + n h-n z n + e n z > n (b) Figura 3.: Modelo modificado dual de señal: (a) modificación de las observaciones; (b) relación entre las observaciones y los multiplicadores de Lagrange. sujeto a z n z n + N H j= N H j= v j ψ j (n) ε + ξ n (3.2) v j ψ j (n) ε + ξ n (3.2) ξ n (3.22) ξn (3.23) El problema de optimización se resuelve aplicando los multiplicadores de

103 3.4. SVM PARA DECONVOLUCIÓN NO CIEGA 75 Lagrange, de donde se obtiene el correspondiente Lagrangiano: 2 N H k= v 2 k + 2γ N α n ( z n + n= N αn(z n n= n I (ξ 2 n + ξ 2 n ) + C n I 2 (ξ n + ξ n) n I 2 γ 2 C2 N H j= N H j= v j ψ j (n) + ε + ξ n ) v j ψ j (n) + ε + ξ n) N β n ξ n n= N βnξ n (3.24) donde α n, αn, β n y βn son los multiplicadores de Lagrange correspondientes a (3.2), (3.2), (3.22) y (3.23), respectivamente. Derivando parcialmente (3.24) con respecto a las variables primales v k se obtiene la siguiente expresión para los pesos del modelo: v k = n= N (α n αn)ψ k (n) (3.25) n= Sustituyendo (3.25) en (3.24) y aplicando las condiciones KKT se obtiene la matriz N N siguiente, A = A(i, l) = N H k= ψ k (i)ψ k (l) = K(i, l) (3.26) definida por un producto escalar en el RKHS, y como consecuencia, puede sustituirse por un núcleo de Mercer. Al igual que en los casos MPS y MDS, se utilizan las condiciones KKT para formular el problema dual, que consiste en maximizar: 2 (α α ) T (A + γi)(α α ) + (α α ) T z ε T (α + α) (3.27) sujeto a α ( ) C (3.28) Sustituyendo (3.25) en el primer término del miembro derecho de (3.7) se obtiene una expresión para las observaciones modificadas, dada por: ẑ n = = = N H j= i= ( N ) (α i αi )ψ j (i) ψ j (n) = i= N ( N H ) (α i αi ) ψ j (i)ψ j (n) = N η i K(i, n) = i= j= N η i Ri n h (3.29) i=

104 76 CAPÍTULO 3. DECONVOLUCIÓN DE SEÑALES DISPERSAS MEDIANTE SVM donde ẑ n es la estimación de las observaciones obtenidas como función de los multiplicadores de Lagrange de manera similar a como se hizo en el caso MDS, se puede reescribir como sigue: ẑ n = η n K n = η n R h n = ˆx n R h n (3.3) En la Figura 3.(b) se representa, a través de un esquema, las relaciones establecidas en (3.3), de donde se observa que el grado de dispersión de la señal estimada ˆx n = η n se ajusta a través del parámetro ε. De esta forma, el algoritmo MSKA incorpora las ventajas de los métodos MPS y MDS, ya que la señal dispersa estimada se obtiene directamente a partir de los multiplicadores de Lagrange. La solución del algoritmo MSKA es dispersa fundamentalmente como consecuencia de la formulación SVM. Así mismo es robusta debido a la regularización implícita en el funcional (3.8), que incorpora la función de coste ε-huber. En la siguiente sección se desarrollan experimentos que permiten evaluar las prestaciones de los algoritmos SVM propuestos, estos son: el MPS y el MSKA. Las prestaciones se miden utilizando dos conjuntos de figuras de mérito, uno que mide las prestaciones en la detección de espículas y otro que miden las prestaciones en la estimación de las mismas. También, se incluye, como una forma de evaluar las capacidades de los algoritmos SVM en aplicaciones prácticas, la deconvolución de señales dispersas en problemas de sísmica de reflexión Simulaciones y Experimentos En esta sección se realizan experimentos que, en general, permiten evaluar las prestaciones de los algoritmos SVM propuestos en la presente Tesis. Los experimentos se orientan a medir prestaciones en lo referente a diferentes estadísticas de ruido aditivo; se consideran, además, situaciones donde la respuesta al impulso es de fase mínima y no mínima; y se valora, también, el efecto del ruido coloreado. Principalmente, los experimentos irán orientados a medir los siguientes aspectos: Experimento. Busca determinar el valor de los parámetros libres óptimos γ, C y ε a utilizar en los algoritmos SVM. Experimento 2. Mide las prestaciones de los algoritmos cuando la señal dispersa de entrada es determinista. Experimento 3. Mide las prestaciones de los algoritmos cuando la señal dispersa de entrada es aleatoria con distribución BG. Experimento 4. Mide las prestaciones de los algoritmos cuando el ruido aditivo a la salida del modelo convolucional es de tipo impulsivo.

105 3.5. SIMULACIONES Y EXPERIMENTOS 77 Experimento 5. Mide las prestaciones de los algoritmos cuando la respuesta al impulso es de fase no mínima. Experimento 6. Mide el efecto del ruido coloreado propio del algoritmo MSKA. Experimento 7. Se aplica el algoritmo MSKA en la estimación de señales dispersas en un problema de sísmica de reflexión. Además, en los experimentos 2, 3, 4 y 5 se aplican a las señales sintéticas otros métodos de deconvolución, lo que permite comparar las prestaciones obtenidas con métodos ampliamente aceptados y utilizados en DD no ciega. Los métodos con los cuales se realizan las comparaciones son: MG y deconvolución norma L, véase las Secciones 2..6 y Estos algoritmos han sido seleccionados debido a que ambos métodos son muy robustos frente a señales dispersas de entrada y ruido aditivo a la salida de estadística no Gausiana [76, 5] Experimento : parámetros libres de la SVM Los algoritmos SVM utilizados en la solución de problemas de DD no ciegos requieren del conocimiento previo de los parámetros libres C, γ y ε. En [22, 64] se han introducido métodos para estimar de una forma analítica los parámetros ε y C basados en el conocimiento a priori de la estadística de la señal y el ruido. También, y a pesar de su elevado coste computacional, se han utilizado los métodos de remuestreo en la estimación de los parámetros libres [6]. Búsqueda de los parámetros C y γ. Debido a que el parámetro ε, tal como se estableció en (3.), controla el grado de dispersión de la señal este experimento se enfoca a analizar el efecto de los parámetros restantes C y γ en la solución del problema de DD no ciega. Dejando el valor de ε fijo se realiza búsqueda exhaustiva sobre los parámetros C y γ. El valor óptimo de dichos parámetros viene dado por aquellos valores que minimizan, en el sentido LS, el valor del error cuadrático medio (MSE, Mean Square Error) como figura de mérito, el cual se define como: MSE = N N (ˆx n x n ) 2 (3.3) n= Descripción de los datos. La búsqueda de parámetros se realiza utilizando la señal determinista propuesta en [35, 36], donde la señal dispersa está formada por 28 muestras con cinco picos en las posiciones x 2 = 8., x 25 = 6.845, x 47 =-5.4, x 7 =4. y x 95 =-3.6. Así mismo, se utiliza el filtro de respuesta al impulso de longitud infinita (IIR, Infinite Impulse Response) cuya

106 78 CAPÍTULO 3. DECONVOLUCIÓN DE SEÑALES DISPERSAS MEDIANTE SVM Amplitude n (a) Amplitude n (b) Amplitud n (c) Amplitud Rad/seg (d) Fase Rad/seg (e) Figura 3.: Señales de prueba: (a) señal dispersa; (b) respuesta al impulso; (c) señal de salida sin ruido; (d) amplitud de la respuesta en frecuencia del filtro; (e) fase de la respuesta en frecuencia del filtro. función de transferencia se escribe a continuación: H(z) = z.6 z 2.44z +.64 (3.32) donde el filtro tiene un cero en z =.6 y dos polos en z =.8 exp(±j5π/2). Como se muestra en la Figura 3.(c) el filtro tiene una banda de paso muy estrecha que elimina muchas de las componentes de la señal dispersa, haciendo el problema de DD no ciega más difícil. Consideraciones sobre la simulación. Utilizando un valor fijo de ε= se realizaron simulaciones para diferentes valores de parámetros C y γ con

107 3.5. SIMULACIONES Y EXPERIMENTOS 79 diferentes valores de relación señal a ruido (SNR, Signal to Noise Ratio). Para cada valor de SNR y dejando fijo dos de los tres parámetros ({ε, C} o {ε, γ}) se obtuvo el promedio MSE (3.3) de realizaciones para cada valor de parámetro libre. Los valores óptimos obtenidos para C y γ fueron muy similares para el rango de SNR de 2 db. Utilizando los valores óptimos estimados para C y γ, como parámetros fijos, se realizó una búsqueda de los valores óptimos, para cada SNR, del parámetro ε. Resultados. Las simulaciones se realizaron utilizando los algoritmos SVM desarrollados para MPS y MSKA. Para cada valor óptimo de ε se calculó el MSE utilizando los rangos de valores para C=[ 4, 4 ] y γ=[ 6, 2 ]. En la Figura 3.2 se observa que tanto para MPS como para MSKA existe un rango de valores de C y γ para los que el MSE tiene un valor mínimo. Comparando los resultados obtenidos para diferentes valores de SNR, como por ejemplo los obtenidos en las Figuras 3.2(a) y 3.2(c), correspondientes a 4dB y 6dB, respectivamente, se observa que los rangos de valores para C y γ son independientes del nivel de ruido. Así, de los resultados obtenidos en ambos casos, MPS Y MSKA, es posible escoger los valores C= 2 y γ= 2 como parámetros óptimos. Estos resultados facilitan la aplicación de los algoritmos SVM, ya que en general la determinación de los parámetros óptimos es un problema difícil y aún no resuelto. Siguiendo el trabajo desarrollado en [64], donde se estudia la dependencia lineal entre ε y las estimaciones del ruido, la estimación de parámetros se centra en la determinación del parámetro ε óptimo Experimento 2: señal de entrada determinista Utilizando la señal sintética determinista del experimento anterior se determinaron las prestaciones de los algoritmos de DD basados en SVM. Las prestaciones de los métodos se midieron a través de dos tipos de figuras de mérito, a saber: figuras de mérito de estimación y figuras de mérito de detección. Las prestaciones de los algoritmos SVM se compararon con la de los métodos MG y L en todos los casos. Figuras de mérito para estimación de señal. Además de utilizar el valor MSE de (3.3) para medir las prestaciones de los algoritmos de DD en la estimación de la señal dispersa, se utilizaron las figuras de mérito propuestas en [8], que se definen como a continuación: F ssd = n (ˆx n x n ) 2 (3.33) F pick = x n (ˆx n x n ) 2 (3.34) F null = (ˆx n x n ) 2 = ˆx 2 n (3.35) x n = x n =

108 8 CAPÍTULO 3. DECONVOLUCIÓN DE SEÑALES DISPERSAS MEDIANTE SVM MSKA MPS γ C (a) γ.5 C (b) MSKA MPS γ C (c) C (d) γ Figura 3.2: Variación del MSE para los métodos de deconvolución basados en los modelos MSKA y MPS en función de los parámetros libres C y γ, para: SNR = 4 db, (a)-(b) y SNR = 6 db (c)-(d). donde F ssd es la suma de las desviaciones al cuadrado de cada instante de señal de su valor verdadero, F pick compara los resultados solamente en localizaciones donde haya picos verdaderos, y F null estima las desviaciones al cuadrado en los instantes nulos de señal. Valores bajos para cada una de las figuras de mérito anteriores indican buenas prestaciones [8]. Figuras de mérito para detección de señal. Las figuras de mérito (3.33)- (3.35) miden las prestaciones en la estimación, si por otro lado se consideran solo las capacidades de detección, se puede usar la sensibilidad y la especificidad como figuras de mérito, las cuales se definen como: V p Se = V p + F n (3.36) Es = V n V n + F p (3.37) donde Se y Es son la sensibilidad y la especificidad, respectivamente. Las variables V p, F p, F n y V n se definen en la Tabla 3.8. Consideraciones sobre la simulación. Para cada método y para cada valor de SNR en el intervalo 4 a 2 db se generaron realizaciones para

109 3.5. SIMULACIONES Y EXPERIMENTOS 8 Señal dispersa de prueba Resultado de la det. Pico Nulo detectado Verdadero positivo(v p ) Falso positivo (F p ) no detectado Falso negativo (F n ) Verdadero negativo (V n ) Tabla 3.8: Definición de las variables V p, F p, F n y V n. MSE SNR [db] (a) MSE MSE SNR [db] (b) SNR [db] (c) Figura 3.3: MSE en función de la SNR para los métodos: MSKA, MPS, MG, L (trazo continuo, discontinuo, punteado y punteado discontinuo, respectivamente) en presencia de ruido: (a) Gausiano; (b) Laplaciano; (c) uniforme. cada parámetro libre. Además, las simulaciones se repitieron para tres tipos diferentes de ruido: Gausiano, Laplaciano y uniforme. MSE y SNR. A partir de la selección previa de los parámetros óptimos se calculó el valor MSE promedio para diferentes valores de SNR. El parámetro óptimo libre de cada uno de los métodos se escogió de acuerdo con la figura de mérito (3.33), es decir, el parámetro que conduce a valores promedios mínimos de F ssd. Una vez seleccionado el parámetro libre para cada uno de los métodos de deconvolución, se aplicó cada uno de los algoritmos obteniéndose el MSE promedio de realizaciones. Las Figuras 3.3(a), 3.3(b) y 3.3(c) muestran los resultados obtenidos para todos los métodos de deconvolución en presencia de ruido Gausiano,

110 82 CAPÍTULO 3. DECONVOLUCIÓN DE SEÑALES DISPERSAS MEDIANTE SVM Laplaciano y uniforme, respectivamente. Es evidente que el algoritmo MSKA tiene mejores prestaciones que los otros métodos. Solo para relaciones señal a ruido altas, superiores a los db, los métodos MG y L ofrecen prestaciones comparables a los resultados obtenidos con el algoritmo MSKA. Además, para obtener una mayor precisión en el análisis de los resultados, se presentan en la Tabla 3.9 los promedios y la desviación estándar del MSE para todos los métodos y tipos de ruido. En general los resultados obtenidos de MSE y su respectiva desviación estándar son similares en cada uno de los tres tipos de ruido analizados, pero en cada uno de éstos se observa que el algoritmo MSKA supera a los otros algoritmos al obtener valores MSE inferiores, así como también valores de desviación estándar relativos a la media más pequeños. Así mismo se observa que todos los algoritmos presentan mejores prestaciones de estimación en presencia de ruido uniforme, comparadas con las obtenidas en los casos de ruido Gausiano y Laplacio, siendo las prestaciones en éstos últimos bastante similares. El algoritmo L ocupa el segundo lugar en prestaciones, pero con una desventaja muy evidente: su alto coste computacional. De hecho las prestaciones de L mejoran en presencia de ruido Laplaciano, comparadas con las obtenidas con ruido Gausiano, pero éstas no superan los resultados obtenidos con el algoritmo MSKA. Cabe notar que el método MG obtiene resultados pobres con SNR bajas, y que además, el valor de su desviación estándar supera el valor promedio de MSE debido a la alta variabilidad de éste método. En resumen, para los tres tipos de ruido y para los diferentes valores de SNR, el algoritmo MSKA presenta las mejores prestaciones en cuanto a la obtención de niveles de MSE menores que los obtenidos con los otros algoritmos. Desviación cuadrática del error. La Figura 3.4 muestra los resultados obtenidos con la figura de mérito F ssd para cuatro valores de SNR. En todos los casos puede observarse que el valor del parámetro libre óptimo que minimiza su valor depende del nivel de ruido. Tanto en los métodos SVM como en L se puede identificar que valores grandes del parámetro libre conducen a la solución cero. La solución cero para la señal estimada puede identificarse como un valor constante, igual a la norma cuadrática de la señal verdadera de prueba, así: F ssd = F pick = N ( x n ) 2 = n= N x 2 n (3.38) El valor de parámetro libre que conduce a la solución cero puede interpretarse como el valor de una cota máxima. Por una parte, de la Figura 3.4(d) se observa que la deconvolución L tiene un valor de cota máximo muy bien definido y que no depende del valor de SNR. Por otra parte, y como pue- n=

111 3.5. SIMULACIONES Y EXPERIMENTOS 83 Método 4 db 8 db 2 db 6 db 2 db L 23.2 ± ± ±.9.6 ± ±.36 MG 85.2±2..± ± ± ±.8 MPS 43.4± ± ±3. 7.3± ±.85 MSKA.8± ±.8.9±.9.74±.3.33±.5 (a) Método 4 db 8 db 2 db 6 db 2 db L 7.7 ± ± ±.7.9 ±.6.45 ±.25 MG 59.9± ±8. 3.7± ± ±.4 MPS 5.7± ± ± ± ±. MSKA 3.±6. 5.2±2.2 2.±.9.85±.4.35±.5 (b) Método 4 db 8 db 2 db 6 db 2 db L 3.4 ±.2.3 ±.7.5 ±.2.2 ±.7.9 ±.3 MG 2.3± ±.7 3.±. 2.56± ±.32 MPS 7.±2.2. ±.5 5.3±.8 2.8±.32.66±.7 MSKA.4±.6.6 ±.2.2±..±.4.4±. (c) Tabla 3.9: MSE±STD (x) para señales determinista con ruido (en negrita se ha resaltado el algoritmo que obtiene las mejores prestaciones y en cursivas el segundo mejor): (a) Gausiano; (b) Laplaciano; (c) uniforme. de verse en las Figuras 3.4(a) y (3.4)(b), en los algoritmos SVM valores grandes de ε conducen de forma suave a la solución cero. También, en la Figura 3.4 se observa que valores excesivamente pequeños del parámetro óptimo conducen a estimaciones muy malas de la señal dispersa. En el caso de deconvolución L, esto se debe a que valores pequeños de λ conducen a una mala regularización, mientras que, para el caso de los algoritmos SVM, valores pequeños de ε conducen a soluciones poco dispersas. De lo anterior se concluye que el valor óptimo del parámetro libre se encuentra en un punto intermedio entre la cota máxima y valores no muy pequeños del parámetro libre. Las Figuras 3.4(a) y 3.4(b) muestran que es posible seleccionar un valor de ε que produzca un valor mínimo de F ssd en los métodos MPS y MSKA. En el caso del método MG no es posible identificar un mínimo para la figura de mérito F ssd. Analizando a través de simulaciones el comportamiento diferente del método MG, se extendieron los rangos de búsqueda del parámetro libre α mg, determinándose que no es posible alcanzar la solución cero y que siempre existe algún valor de reflector detectado. Otras figuras de mérito de estimación. Las figuras de mérito F pick y F null

112 84 CAPÍTULO 3. DECONVOLUCIÓN DE SEÑALES DISPERSAS MEDIANTE SVM F ssd 5 F ssd ε (a) ε (b) F ssd F ssd α mg (c) λ (d) Figura 3.4: F ssd en función del parámetro libre para SNR = 4, 8, 2 y 2 db (trazo continuo, discontinuo, punteado y punteado discontinuo, respectivamente ) en: (a) MPS; (b) MSKA; (c) MG; (d) L. pueden utilizarse también para determinar los parámetros libres óptimos. La función F pick determina la desviación cuadrática de los picos reales y su valor se incrementa con el valor del parámetro libre hasta alcanzar un valor constante, determinado por la solución cero, mientras que F null decrece con el valor del parámetro libre hasta que, en la solución cero, alcanzar el valor de cero. La intersección de ambas figuras de mérito conduce a un punto óptimo, ya que corresponde con un valor de compromiso entre la estimación de los reflectores y la estimación de los valores de señal igual a cero. En la Figura 3.5 se presentan las figuras de mérito F pick y F null para una SNR igual a 4dB. Al igual que en el caso F ssd, se pueden identificar valores de parámetros libres que conducen a la solución cero, excepto en MG. Una buena elección del parámetro óptimo en los métodos SVM es el punto de intersección de ambas figuras de mérito, y se observa que éstos son muy parecidos a los obtenidos con la figura de mérito F ssd. Probabilidad de detección. Se calculó la probabilidad de detección en función del parámetro libre en cada uno de los métodos de deconvolución. Las pruebas de detección se hicieron con una SNR igual a 4dB. Las Figuras 3.6(a)-(d) muestran los promedios de detección obtenidos

113 3.5. SIMULACIONES Y EXPERIMENTOS F pick y F null 2 5 F pick y F null ε (a) ε (b) F pick y F null F pick y F null α mg (c) λ (d) Figura 3.5: Figura de mérito F pick (trazo continuo) y F null (trazo discontinuo) para SNR = 6 db en: (a) MPS; (b) MSKA; (c) MG; (d) L. para la señal determinista. Nuevamente se puede comprobar, en todos los métodos, salvo MG, la existencia de un valor de parámetro libre que sirve de cota máxima, y superado dicho valor se obtiene la solución cero. En promedio, todos los métodos poseen algún valor de parámetro libre que permite alcanzar valores de detección cercanos al %. Sensibilidad y especificidad. Las figuras de mérito Se y Es se calculan en función de los parámetros libres. A diferencia de las figuras de mérito (3.33)- (3.36), la sensibilidad y especificidad miden la capacidad de detección de los métodos y no la de estimación. El punto de intersección entre ambas figuras de mérito conduce a un valor de compromiso entre la sensibilidad y la especificidad. En la Figura 3.7, para una SNR igual a 4dB, se ha calculado Se y Es para todos los métodos. Se observa que el método MSKA supera a los demás métodos, al conducir a valores más altos de sensibilidad y especificidad en el punto de cruce.

114 86 CAPÍTULO 3. DECONVOLUCIÓN DE SEÑALES DISPERSAS MEDIANTE SVM Probabilidad de detección Probabilidad de detección ε (a) ε (b) Probabilidad de detección Probabilidad de detección α mg (c) λ (d) Figura 3.6: Promedio de la probabilidad de detección para la señal determinista en: (a) MPS; (b) MSKA; (c) MG; (d) L Experimento 3: señal de entrada aleatoria El propósito de este experimento es el de mostrar las capacidades de los diferentes algoritmos SVM en la DD no ciega de señales dispersas aleatorias. Se utilizó para este experimento un tipo de señal aleatoria con distribución Bernoulli-Gausiana (BG) [76], donde las muestras se generan de acuerdo con el siguiente modelo: x k = r k q k (3.39) donde r k es ruido blanco Gausiano, de media cero y varianza σ a, y q k es una secuencia Bernoulli, para la cual: P r (q k ) = { p, qk = p, q k = (3.4) Se seleccionaron 28 muestras de una señal dispersa aleatoria con distribución BG. La localización de los valores de señal diferente de cero tienen una distribución Bernoulli con una probabilidad p =.5, mientras que sus magnitudes tienen una distribución Gaussiana de media cero y varianza σ a =. La respuesta al impulso es la misma que se usó para el caso de la señal determinista.

115 3.5. SIMULACIONES Y EXPERIMENTOS ε (a) ε (b) α mg λ (c) (d) Figura 3.7: Sensibilidad (trazo continuo) y especificidad (trazo discontinuo) para una SNR = 4dB en: (a) MPS; (b) MSKA; (c) MG; (d) L. Método 4 db 8 db 2 db 6 db 2 db L.76 ±.8.37 ±.48.6 ±.2.7±.5.3±.3 MG.64±.65.39±.59.32±.47.29±.44.22±.29 MPS.27±.27.65±.49.34±.37.2±.7.3±.8 MSKA.4±.39.9±.9.2±.36.8 ±.6.7 ±. Tabla 3.: MSE±STD(x) para señales aleatoria con ruido Gausiano (en negrita se ha resaltado el algoritmo que obtiene las mejores prestaciones y en cursivas el segundo mejor). MSE y SNR. En la Figura 3.8 se observa que el método MSKA consigue los niveles de MSE más bajos que los otros métodos, siendo superado únicamente por la deconvolución L para SNR altas. Siguiendo el mismo procedimiento que en el caso determinista, se presentan en la Tabla 3. los promedios y la desviación estándar para todos los métodos de deconvolución en presencia de ruido Gausiano. Los valores de desviación estándar muestran un problema común en todos los métodos: la presencia de valores de desviación estándar tan grandes como el valor promedio MSE. A través de simulaciones se ha comprobado que el MSE calculado para realizaciones

116 88 CAPÍTULO 3. DECONVOLUCIÓN DE SEÑALES DISPERSAS MEDIANTE SVM x MSE SNR [db] Figura 3.8: MSE en función de SNR para los métodos: MSKA; MPS; MG; y L (trazo continuo, discontinuo, punteado y punteado discontinuo, respectivamente) en presencia de ruido Gausiano y entrada aleatoria con distribución BG ( σ a =, p =.5). tiene una estadística no Gausiana debido a la alta variabilidad de las soluciones, y este es un problema bastante común en deconvolución de señales dispersas [76]. Figuras de mérito de estimación. Las figuras de mérito de estimación son las establecidas en las Ecuaciones (3.33)-(3.35). Los resultados obtenidos son muy similares a los que se obtuvieron en el caso de una señal de entrada determinista. En la Figura 3.9 se muestra la variación de la desviación cuadrática del error, representado a través de la figura de mérito F ssd, donde se observa que el valor del parámetro óptimo que minimiza el error en los algoritmos L y SVM se encuentra en un valor intermedio entre el parámetro que determina la solución cero y valores cercanos a cero. Como puede observarse en las Figuras 3.9(a)-3.9(d) los valores obtenidos para el parámetro ε son más pequeños que los obtenidos en la DD con señal determinista. Cabe hacer una observación con respecto a los resultados obtenidos con los métodos MG y L. El primero, al igual que el resultado obtenido cuando la señal de entrada era determinista, obtuvo promedios constantes de F ssd para valores elevados del parámetro libre, mientras que el segundo condujo a la solución cero para valores de parámetro libre similares al del experimento con entrada determinista. Además, y como puede observarse de la Figura 3.9(d), el método L no presenta un trazo suave en su curva de F ssd, debido a que, para algunas realizaciones, el algoritmo Simplex produce valores anormalmente elevados de MSE. Las figuras de mérito F pick y F null condujeron a valores de parámetros óptimos muy parecidos a los obtenidos con la desviación cuadrática del error. En la Tabla 3. (a)-(b) se muestran los resultados de calcular el parámetro óptimo: a través del mínimo de F ssd y a través de la intersección entre F pick y F null. Es evidente que el valor de parámetro libre óptimo en ambos casos

117 3.5. SIMULACIONES Y EXPERIMENTOS F ssd F ssd ε (a) ε (b) F ssd 4 3 F ssd α mg (c) λ (d) Figura 3.9: F ssd para una entrada aleatoria, con distribución BG, en función del parámetro libre para SNR = 4, 8, 2 y 2 db (trazo continuo, discontinuo, punteado y punteado discontinuo, respectivamente) en: (a) MPS; (b) MSKA; (c) MG; (d) L. es muy parecido, dando un voto de confianza al valor de parámetro óptimo calculado. Figuras de mérito de detección. La probabilidad de detección se calculó en función del parámetro libre, donde para cada valor del parámetro se promedió el resultado de realizaciones. Las simulaciones se realizaron para diferentes valores de SNR. Las Figuras 3.2(a)-3.2(d) muestran los promedios de detección obtenidos para cada valor de parámetro libre para SNR = 4dB. Se observa que, para este tipo de entrada, y con el nivel de ruido alto producido por una SNR = 4dB, la probabilidad de detección es baja para todos los métodos. Los valores del parámetro ε óptimo obtenidos en las Figuras 3.2(a)-3.2(d), es decir, aquellos para los que se tienen probabilidades de detección más altas, se corresponden con los valores óptimos obtenidos con las figuras de mérito de estimación (3.33)-(3.35). Las otras figuras de mérito relacionadas a la detección vienen dadas por (3.36)-(3.37), y los resultados de éstas conducen a la obtención de un valor óptimo de parámetros libres, de manera similar a como se obtuvo con las figuras de mérito F pick y F null. La Tabla 3.(c) muestra los valores de

118 9 CAPÍTULO 3. DECONVOLUCIÓN DE SEÑALES DISPERSAS MEDIANTE SVM Met(par) 4 db 6 db 8 db db 2 db 4 db 6 db 2 db L (λ) MG(α) MPS (ε) MSKA(ε) (a) Met(par) 4 db 6 db 8 db db 2 db 4 db 6 db 2 db L (λ) MG(α) MPS (ε) MSKA(ε) (b) Met(par) 4 db 6 db 8 db db 2 db 4 db 6 db 2 db L (λ) MG(α) MPS (ε) MSKA(ε) (c) Tabla 3.: Parámetros libres calculados a partir de diferentes figuras de mérito: (a) F ssd ; (b) F pick y F null ; (c) Se y Es. parámetros óptimos obtenidos a partir de las intersecciones de las curvas de sensibilidad (3.36) y especificidad (3.37) que, coinciden con los valores obtenidos con las otras figuras de mérito Experimento 4: efecto del ruido impulsivo Los algoritmos de DD se ven seriamente comprometidos cuando las observaciones se ven distorsionadas por ruido impulsivo. A través de simulaciones, en el presente experimento se midieron las prestaciones de los diferentes algoritmos frente al ruido impulsivo, el cual se generó como una señal aleatoria blanca con distribución BG, de manera similar a como se hizo con la señal de entrada aleatoria con distribución BG del experimento anterior. La Figura 3.2 muestra, a través de un ejemplo, el comportamiento de los diferentes métodos de DD cuando las observaciones han sido distorsionadas con ruido impulsivo. A través de una distorsión controlada, utilizando un único impulso en diferentes posiciones y para diferentes magnitudes, se puede inferir las capacidades y las limitaciones de los diferentes métodos. El ejemplo de la Figura

119 3.5. SIMULACIONES Y EXPERIMENTOS 9 Probabilidad de detección Probabilidad de detección ε (a) ε (b) Probabilidad de detección Probabilidad de detección α (c) λ (d) Figura 3.2: Promedio de la probabilidad de detección para la señal aleatoria con distribución BG obtenido con: (a) MPS; (b) MSKA; (c) MG; (d) L. 3.2 ilustra la distorsión introducida por un único impulso de magnitud 2 en la posición y 65. Este sencillo experimento condujo a una importante conclusión sobre las limitaciones del algoritmo L : la presencia de ruido impulsivo ralentiza la convergencia del algoritmo Simplex, y en la mayoría de los casos diverge. Este problema se resuelve agregando una componente de ruido Gausiano, produciendo una SNR alta igual a 2 db. El método MG conduce a estimaciones razonables pero interpreta valores grandes de ruido impulsivo como espículas de la señal dispersa. De los resultados de la Figura 3.22(c) se observa que el método MSKA es el que mejor resuelve el problema de ruido impulsivo, a pesar de ser inferior a MG en cuanto a estimación de la amplitud. A continuación se presentan experimentos más exhaustivos, que utilizan la señal determinista del Experimento 2, los que permiten valorar más de cerca las capacidades de los diferentes métodos. MSE y SNR. Se calculó el MSE en función de la SNR para todos los métodos de DD, de manera similar a como se hizo en los experimentos anteriores. Los parámetros libres óptimos fueron previamente seleccionados de acuerdo a los resultados obtenidos con las diferentes figuras de mérito. La Figura 3.22 muestra los resultados que se obtuvieron para tres valores diferentes de ruido impulsivo, generados a partir de la variación del paráme-

120 92 CAPÍTULO 3. DECONVOLUCIÓN DE SEÑALES DISPERSAS MEDIANTE SVM Amplitud 8 6 Ruido n (a) 8 Amplitud n (b) (d) (c) l(4db),ε =.9 x xest (e) Figura 3.2: Ejemplo que ilustra las capacidades de los métodos de deconvolución en presencia de ruido impulsivo (los círculos indican el valor verdadero de la señal dispersa): (a) se introduce en la muestra 65 de las observaciones una espícula de magnitud 2; (b) filtrado anticausal de las observaciones como primer paso del MSKA (el trazo continuo representa la señal filtrada y el discontinuo la señal sin filtrar) ; (c) resultado de aplicar el algoritmo MSKA; (d) resultado del método MG; (e) resultado del método L. tro de probabilidad Bernoulli (p). Se observa que los métodos L y MSKA condujeron a los mejores resultados, obteniendo el primero mejores estimaciones para SNR bajas. Sin embargo, en ambos casos el nivel de MSE es muy similar en la mayor parte del rango de SNR, en general a partir de valores

121 3.5. SIMULACIONES Y EXPERIMENTOS MSE SNR [db] (a) MSE MSE SNR [db] (b) SNR [db] (c) Figura 3.22: MSE en función de SNR para todos los métodos: MSKA, MPS, MG, L (trazo continuo, discontinuo, punteado y punteado discontinuo, respectivamente) en presencia de ruido impulsivo para diferentes valores de probabilidad Bernoulli: (a) p=.; (b) p=.5; (c) p=.. de SNR superiores a los 8 db. Para analizar con más detenimiento los resultados, la Tabla 3.2 presenta los valores promedios de MSE así como sus respectivas desviaciones estándar. Se observa que los algoritmos L y MS- KA, a diferencia de los experimentos anteriores, donde el valor de desviación estándar es grande en relación con el valor de la media, obtienen valores de desviación estándar más pequeños. Figuras de mérito de estimación. La Figura 3.23 muestra la variación de la desviación cuadrática del error, F ssd, en función del parámetro libre para los diferentes métodos de DD, para un valor de probabilidad p =.. En general, todos los algoritmos mostraron respuestas similares a las obtenidas en los experimentos 2 y 3. De esta manera, los métodos MSKA, MPS y L presentaron valores de parámetros libres para los cuales el valor de F ssd es mínimo, así como, también, valores a partir de los cuales la respuesta F ssd alcanza un valor constante e igual a la solución cero. Además, en todos los algoritmos la estimación de la señal dispersa resultó ser menos dependiente del nivel de ruido, los valores de F ssd presentan promedios muy similares para diferentes valores de SNR, como se observa en la Figura Más importan-

122 94 CAPÍTULO 3. DECONVOLUCIÓN DE SEÑALES DISPERSAS MEDIANTE SVM Método db 4 db 8 db 2 db 6 db 2 db L.7±.3.7 ±.4.6 ±.3.7 ±.3.7 ±.4.6 ±.3 MG 3.4± ±. 2.3±.6 2.4±.2 2.5±.7 2.4±.2 MPS 7.5±7 5.3±3 4.2±.8 3.7± 3.5± ±.8 MSKA.3 ±2..6±.6.4±.2.3±.2.3±..3±. (a) Método db 4 db 8 db 2 db 6 db 2 db L.9±2..4±.8.3±.4.2 ±..2 ±.. ±.6 MG 7.5±9.2 3.± ±.3 2.5±.4 2.4±.6 2.3±.4 MPS 7.3±.8 8.5± ± ±.3 2.2±.6.9±.3 MSKA 3. ±2..3 ±..6 ±.5.2±.6.±.5.7±.3 (b) Método db 4 db 8 db 2 db 6 db 2 db L 2.2 ±2.8.5 ±.4.±..9 ±.4.8 ±.5.7 ±.3 MG 3.8±9. 3.9± ± ±.9 2.3±. 2.4±.3 MPS 3.4± ± ± ±.6 4.2±. 3.7±.8 MSKA 4.7 ± ±.4. ±.7.6±.3.4±.3.7±. (c) Tabla 3.2: MSE±STD(x) para señal determinista con ruido impulsivo (en negrita se ha resaltado el algoritmo que obtiene las mejores prestaciones y en cursivas el segundo mejor): (a) p=.; (b) p=.5; (c) p=.. te aún, fue el hecho de que en los algoritmos MPS, MSKA y L los valores de parámetro libre obtenidos fueron bastante similares a los obtenidos en el experimento 2, pudiendo usarse éstos en cualquiera de los dos tipos de ruido. Los resultados que se obtuvieron usando las figuras de mérito F pick y F null no se muestran en el presente experimento, pero al igual que con F ssd condujeron a prestaciones muy similares a las obtenidas en el experimento 2. En la Tabla 3.3 se muestran los valores de parámetros libres óptimos obtenidos con F pick y F null, los cuales son muy similares a los conseguidos con F ssd. Figuras de mérito de detección. Se calcularon los promedios de la probabilidad de detección como función del parámetro libre. Para cada valor de parámetro libre se aplicaron los algoritmos de DD sobre realizaciones. Los experimentos se repitieron para diferentes valores de SNR. La Figura 3.24 muestra los promedios de detección obtenidos para cada valor de parámetro libre para una SNR = 4dB. El valor alto de probabilidad de detección debe interpretarse con cuidado, ya que el ruido impulsivo conduce a falsas detecciones, por lo que se presta especial atención a las figuras de mérito que miden la sensibilidad y la especificidad dadas por (3.36)-(3.37).

123 3.5. SIMULACIONES Y EXPERIMENTOS F ssd F ssd ε (a) ε (b) F ssd 8 F ssd α mg (c) λ (d) Figura 3.23: F ssd para una entrada determinista con ruido impulsivo (p=.), con distribución BG, en función del parámetro libre para SNR = 4, 8, 2 y 2 db (trazo continuo, discontinuo, punteado y punteado discontinuo, respectivamente) en: (a) MPS; (b) MSKA; (c) MG; (d) L. Teniendo en cuenta las capacidades de detección, el parámetro libre óptimo se calculó a través de las figuras de mérito Es y Se. En la Tabla 3.3 se comparan los valores de parámetros libres óptimos obtenidos con diferentes figuras de mérito. Se observa que los valores obtenidos a partir de Es y Se son mayores que los obtenidos a partir de las otras figuras de mérito Experimento 5: efecto de la fase no mínima La mayoría de métodos de DD no ciegos se basan fundamentalmente en calcular el inverso de la respuesta al impulso, suponiendo para ello que ésta es de fase mínima. El siguiente experimento tiene como objetivo ilustrar las capacidades del algoritmo MSKA cuando la respuesta al impulso es de fase no mínima. Consideraciones sobre la simulación. Se realizó un conjunto de simulaciones con el sistema (3.32), y cuya función de transferencia se reescribe a

124 96 CAPÍTULO 3. DECONVOLUCIÓN DE SEÑALES DISPERSAS MEDIANTE SVM Met(par) -4 db db 4 db 8 db 2 db 6 db 2 db L (λ) MG(α) MPS (ε) MSKA(ε) (a) Met(par) -4 db db 4 db 8 db 2 db 6 db 2 db L (λ) MG(α) MPS (ε) MSKA(ε) (b) Met(par) -4 db db 4 db 8 db 2 db 6 db 2 db L (λ) MG(α) MPS (ε) MSKA(ε) (c) Tabla 3.3: Parámetros libres calculados para el ruido impulsivo a partir de diferentes figuras de mérito y probabilidad p =.: (a) F ssd ; (b) F pick y F null ; (c) Se y Es. continuación: H(z) = z c 5π j (.8 e 2 z )(.8 e j 5π 2 z ) (3.4) donde c es el valor del único cero de la respuesta al impulso (3.4). Como señal dispersa de entrada se utilizó la señal determinista de 5 impulsos del experimento 2. El ruido aditivo utilizado fue de tipo Gausiano. El experimento consistió fundamentalmente en evaluar las prestaciones del algoritmo MSKA para algunos valores de c que están fuera de la circunferencia de radio unidad, es decir aquellos para los que (3.4) es de fase no mínima. La influencia de la posición geométrica del cero, (z c), de (3.4) en la estimación del parámetro óptimo se muestra en la Figura Las Figuras 3.25(b) y 3.25(c) muestran la F ssd para SNR = 4 db y SNR = 8 db, respectivamente, en ambos casos se muestran los resultados para los valores de c =.6 y c = 2.. En ambos casos la figura de mérito F ssd alcanzó el mismo nivel de error cuadrático medio, pero a diferencia del caso de fase mínima, el sistema con respuesta al impulso de fase no mínima requirió valores de parámetro libre más grandes.

125 3.5. SIMULACIONES Y EXPERIMENTOS 97 Probabilidad de detección Probabilidad de detección ε (a) ε (b) Probabilidad de detección Probabilidad de detección α (c) λ (d) Figura 3.24: Promedio de la probabilidad de detección para la señal determinista con ruido impulsivo, y probabilidad p =., obtenido a través de: (a) MPS; (b) MSKA; (c) MG; (d) L. DD de señales deterministas. En las Figuras 3.26(c) y 3.26(d) se muestran señales dispersas estimadas al aplicar el algoritmo MSKA a las señales deterministas de las Figuras 3.26(a) y 3.26(b), respectivamente. En ambos casos se seleccionaron respuestas al impulso de fase no mínima, con ceros en c =.4 y c = 2., y valor del parámetro libre ε =. y SNR = 4dB. En el primero caso, para c =.4, el valor de ε =. obtuvo buenas prestaciones en detección y estimación. Por otro lado, cuando el cero se alejó más, hacia la posición c = 2., la estimación de la señal dispersa introdujo algunos valores de espículas espurias. En este punto cabe decir que el ejemplo de la Figura 3.26(c) es bastante exagerado si se considera que, en general, las señales correspondientes a fenómenos físicos tienen sus ceros situados muy cerca de la circunferencia de radio unidad [53, 6, 22] Experimento 6: efecto del ruido coloreado En el presente experimento se estudia uno de los problemas intrínsecos en el desarrollo del algoritmo MSKA. En la Sección se formuló el problema

126 98 CAPÍTULO 3. DECONVOLUCIÓN DE SEÑALES DISPERSAS MEDIANTE SVM.8 Parte Imaginaria Fase no mínima Parte Real (a) F ssd 5 F ssd ε (b) ε (c) Figura 3.25: F ssd del algoritmo MSKA calculado para sistemas de fase mínima, c=.6, y no mínima, c=2, (trazo continuo y discontinuo, respectivamente): (a) distribución geométrica de los polos y ceros de la respuesta al impulso de un sistema de fase mínima al que se le desplaza el cero hasta sacarle fuera de la circunferencia de radio unidad, haciéndole de fase no mínima; (b) F ssd para una SNR = 4dB; (c) F ssd para una SNR = 8dB. modificando previamente las observaciones como sigue: z n = ˆx n h n δ n+q + e n (3.42) donde donde e n = h Q n e n es la componente de ruido aditivo modificado, y se le denomina ruido coloreado. El ruido coloreado de (3.42) tiene un ancho de banda definido por la respuesta al impulso y sus muestras ya no están incorrelacionadas. En la literatura este problema suele ignorarse, recurriendo a un modelo más simplificado que no incluye el ruido aditivo a la salida [4], ignorándose el efecto que éste pueda tener en los resultados. Trabajos recientes han propuesto minimizar el efecto del ruido coloreado a través de un doble filtrado sobre los resultados de la deconvolución [59, 79]. En este experimento, se demuestra la robustez de la función de coste ε- Huber que, a través de su regularización implícita, permite minimizar el efecto del ruido coloreado. Para mostrar su robustez se realizó una comparación del

127 3.5. SIMULACIONES Y EXPERIMENTOS Amplitud 2 4 Amplitud n (a) n (b) Parte Imaginaria Parte Imaginaria Parte Real Parte Real (c) (d) Figura 3.26: Resultados obtenidos al aplicar el algoritmo de deconvolución MSKA a dos señales con respuesta al impulso de fase no mínima: (a) señal determinista de fase no mínima ruidosa, c =.4; (b) señal determinista de fase no mínima ruidosa, c = 2.; (c) señal dispersa estimada, c =.4; (d) señal dispersa estimada, c = 2. desempeño del método MSKA frente a dos tipos de ruido: el coloreado, obtenido a partir de (3.42), y una versión modificada que recupera el carácter blanco de las observaciones modificadas. Las características espectrales en magnitud del ruido blanco se recuperaron a través de la modificación previa de las muestras de error, tal como se ilustra en la Figura Esto se logró filtrando las muestras de ruido aditivo con el inverso de la respuesta al impulso. e mod n = e n h n (3.43) donde h n es el inverso de la respuesta al impulso. Las muestras de ruido modificadas se sumaron a la salida del modelo convolucional. Al aplicar el algoritmo MSKA, el conjunto de observaciones y n se modificó filtrando con el inverso temporal h Q n obteniéndose las observaciones: z n = { } x n h n + e mod n h Q n = x n K h n + e blanco n (3.44)

128 CAPÍTULO 3. DECONVOLUCIÓN DE SEÑALES DISPERSAS MEDIANTE SVM e + n - h n x n h n + y n h Q -n z n Figura 3.27: Esquema que ilustra las modificaciones previas realizadas sobre las muestras de ruido, filtradas con el inverso de la respuesta al impulso. donde e blanco n = e n h n h Q n (3.45) es la componente con características espectrales en amplitud de ruido blanco y K h n es el kernel de autocorrelación. A las observaciones modificadas según (3.42) y (3.44), se aplicó el algoritmo MSKA sobre realizaciones. Los valores de parámetros óptimos utilizados fueron los que se determinaron en el experimento de búsqueda de parámetros libres, es decir, C = y γ = 6. Para las simulaciones se utilizó la señal dispersa determinista y la respuesta al impulso del Experimento 2, con lo cual el parámetro ε utilizado, calculado para cada SNR, es el que minimiza la figura de mérito F ssd. En la Tabla 3.4 se presentan los resultados obtenidos para el valor promedio de MSE calculado para los tipos de modificaciones sobre las observaciones. Los valores promedios MSE son bastante parecidos en ambas situaciones, mientras que la desviación estándar del MSE es más pequeña para el caso de ruido coloreado. Estos resultados permiten determinar que el algoritmos MSKA conduce a buenas estimaciones de la señal dispersa incluso en situaciones donde se modifica las características estadísticas del ruido. Ruido Gausiano Ruido Gausiano modificado 4 db 8 db 2 db 2 db 4 db 8 db 2 db 2 db.7± ±2. 2.±.9.3±..8± ±3.4.9±..3±.2 Tabla 3.4: MSE obtenido a partir de la DD de una señal determinista con el método MSKA, donde se compara el efecto de modificar previamente las muestras de ruido Experimento 7: aplicación práctica Uno de los métodos más comúnmente utilizados en estudios de sísmica de reflexión consiste en muestrear varias veces un mismo punto del subsuelo, teniendo como referencia un punto medio común. En la Figura 3.28(a) se

129 3.5. SIMULACIONES Y EXPERIMENTOS Intervalo de grupo (IG) Array de medidores R R R R R5 R R R R R5 R6 m m 2 m 3 m 4 m 5 m 6 R8R7R6R5 R3 R2 R S2 S3 S5 S6 S7S8 4 S S4 R S S6 S 8 Puntos fuente Puntos medios (CMP s).5 IG R R R R R5 Puntos medios R6 S S6 S8 Punto medio m 6 Desplazamiento [mts] r r 2 r 3 r 4 r 5 r 6 r 7 r 8 Traza Apilada Tiempo [seg] t Multiplicidad (a) (b) S Superficie Plano Reflector m M 6 R vto/2 vt/2 h t o/2 S eje t/2 M hipérbola t/2 h R eje h Reflector (c) m 6.2 Segundos Tiempo [segs].4.6 Trazo tipo WTVA Metros Distancia [mts] (d) (e) Trazo tipo WT Figura 3.28: Procesamiento de señales sísmicas: (a) esquema que ilustra la obtención del CSPG en seis puntos diferentes, así como un esquema para el punto medio común m 6 ; (b) el concepto de multiplicidad; (c) esquema para un único trazo fuente-receptor del punto m 6 e ilustración del concepto de apilamiento; (d) datos reales obtenidos de un conjunto CSPG [7]; (e) trazas más cercana y más lejana del conjunto CSPG y dos formas estándar de representación: trazo tipo WT y WTVA. observa cómo, bajo la suposición de un único estrato plano, la adquisición se realiza moviendo a lo largo de una línea imaginaria una fuente sísmica y un array de medidores. La principal ventaja del método radica en que a partir de la obtención de varios registros sobre un mismo punto se mejora la SNR [4]. Al número de veces que se muestrea un mismo punto se le conoce como multiplicidad, y depende del número de medidores utilizados, de la separación

130 2 CAPÍTULO 3. DECONVOLUCIÓN DE SEÑALES DISPERSAS MEDIANTE SVM x n Deconvol. NMO Apilamiento Migración y n Figura 3.29: Esquema general que muestra las tres técnicas más importantes en el procesamiento de señales sísmicas. entre éstos y de la separación de los puntos fuente [53, 4]. Para simplificar el análisis, se asume que la onda se propaga a través de medios homogéneos, isótropos y sin perdidas, es decir la velocidad en cada una de las capas que forman el modelo del subsuelo es constante. La Figura 3.28(a) muestra un array de 6 medidores separados uniformemente por una distancia llamada intervalo de grupo (IG). El subsuelo se estimula en 8 puntos fuente S n, con una separación entre sí igual a la mitad del IG. Los puntos de reflexión se conocen con el nombre de puntos medios y tienen una separación igual a la de los puntos fuente. Cada estimulación genera un grupo de trazas con fuente común (Common Source Point Gathering, CSPG) medidas por el array de sensores. El grupo de trazas comunes a un punto medio forman el grupo de trazas de punto medio común (Common Midpoint Gather, CMPG). En la Figura 3.28(b) se muestra que para la configuración del experimento se tiene que a partir del punto medio común m 6, la multiplicidad es igual a 6. En general los experimentos se realizan con configuraciones de multiplicidad alta con las que se consiguen SNR altas. En la Figura 3.28(c) se muestra, a través de un esquema simple formado por un único par fuente-receptor, la distribución hiperbólica que adquieren las trazas agrupadas entorno a un punto medio común [4]. La distancia h, del punto medio M al punto R, se conoce como distancia de desplazamiento. El valor de tiempo t o, que emplea la onda recorrer la distancia M-m 6 -M, se conoce con el nombre de tiempo de la distancia de desplazamiento cero. A partir de la Figura 3.28(c), y utilizando trigonometría básica, se obtiene la expresión siguiente: ( vt ) 2 ( vto ) 2 = + h 2 (3.46) 2 2 donde v representa la velocidad de la onda, t el tiempo que emplea ésta en cubrir la distancia S-m 6 -R, y h es la distancia de desplazamiento u offset. Si se reescribe (3.46) como a continuación: ( t 2) 2 h 2 v 2 = ( t 2 o 2 ) (3.47) y ya que se ha supuesto que tanto v como t o son constantes, la Ecuación (3.47) corresponde a la ecuación de una hipérbole con coordenadas (h, t/2). El conjunto de datos obtenidos a partir de un experimento de sísmica de reflexión, tal como se describe en la Figura 3.28, involucra tres variables

131 3.5. SIMULACIONES Y EXPERIMENTOS 3 Coordenada Tiempo (t) Desplazamiento (h) Punto medio común (m) Técnica de procesamiento de señal Deconvolución Corrección CMP (NMO) y apilamiento Migración Tabla 3.5: Procesamiento de datos sísmicos. importantes: el tiempo t, la distancia h y el punto medio común m. Cada una de estas variables se corresponde con tres métodos de procesamiento digital de señal, a saber: deconvolución, apilamiento y análisis de velocidad, y migración, respectimamente. La Tabla 3.5 muestra por una parte la correspondencia entre las coordenadas {t, h, m} y sus respectivas técnicas de procesamiento, mientras que, por otra parte, la Figura 3.29 muestra la secuencia de ejecución de estas técnicas. Previo a la aplicación de cualquiera de éstas técnicas de procesado de señal, se realizan sobre los registros varios tipos de preprocesado, entre los cuales se pueden mencionar: el demultiplexado, la edición de trazas, la aplicación de ganancias o la corrección estática. El interés de la presente aplicación se centra en la deconvolución de registros sísmicos del tipo CSPG. En la Figura 3.28(d) se muestra un conjunto CSPG de señales sísmicas medidas a lo largo de una línea de sensores de un kilómetro de longitud [7]. Para mejorar la presentación visual del conjunto de trazas, éstas se representan como una imagen, donde la amplitud de la traza se corresponde con una escala de grises previamente definida [7]. Debido a que las trazas que forman el CSPG presentan rangos dinámicos muy altos, se hace necesario, para mejorar la representación visual, modificar el registro a través de la amplificación y truncamiento de cada una de las trazas. Así, la Figura 3.29(d) muestra el conjunto CSPG con una amplificación y truncamiento de 4 db, es decir, cada una de las trazas se amplifica 4 db truncando aquellos valores que superan ese valor. La Figura 3.28(e) muestra dos trazas correspondientes a las mediciones realizadas por los sensores situados a y metros, respectivamente. Además, para la traza más lejana se han dibujado los dos tipos de representación más comúnmente encontrados en la literatura, a saber: el trazo de oscilación (Wiggle Trace, WT) y el trazo de oscilación temporal de área variable (Wiggle Trace Variable Area, WTVA). Resultados de la deconvolución. La deconvolución sísmica require la estimación previa de la onda sísmica. Los métodos de deconvolución tradicionales tales como la DH y la DP suelen conducir a muy buenas estimaciones de la onda sísmica, no así de la señal dispersa. Las Figuras 3.3(a) y 3.3(b) muestran las ondas sísmicas de cada una de las trazas y el promedio de éstas, respectivamente, obtenidas a partir de las trazas de la Figura 3.28(d) utilizando DP. Utilizando la onda sísmica estimada,y aplicando el algoritmo SVM- MSKA a cada traza del CSPG, se estimó la señal dispersa para cada una

132 4 CAPÍTULO 3. DECONVOLUCIÓN DE SEÑALES DISPERSAS MEDIANTE SVM amplitud 5 5 (a) Tiempo [seg] (b).5. Amplitud.5 Amplitud Tiempo [seg] (c) Tiempo [seg] (d) Figura 3.3: Estimación de la onda sísmica: (a) obtenida traza a traza utilizando DP; (b) onda sísmica obtenida promediando las ondas sísmicas estimadas de cada una de las trazas del CSPG; (c) valor normalizado para la señal dispersa estimada; (d) comparación del valor verdadero de la traza con su valor calculado a partir de la convolución de la onda y señal dispersa estimadas. de las trazas. Los parámetros libres utilizados se ajustaron de acuerdo a los resultados obtenidos en los experimentos con señales sintéticas, γ= 6 y C=. El valor de ε empleado fue de ε=x 3. La Figura 3.3(c) muestra los valores normalizados de una única traza del CSPG junto con el valor estimado de su señal dispersa. Como una forma de comprobar si el resultado para la señal dispersa es adecuado la Figura 3.3(d) compara una traza estimada, calculada como la convolución de la onda sísmica y la señal dispersa estimada, y su valor verdadero. La similitud entre ambas señales indica la congruencia de la solución obtenida. Las Figuras 3.3(a), 3.3(b) y 3.3(c) muestran la señal dispersa, la señal verdadera y la señal estimada para varias trazas, respectivamente. Con el objetivo de mejorar la visualización, estas Figuras solamente presentan un segmento, entre los.3 y los.9 segundos. Además, se utiliza la forma de representación estándar WTVA que rellena de color la parte positiva de la onda.

133 3.5. SIMULACIONES Y EXPERIMENTOS 5.3 Longitud [mts] Tiempo [seg] (a) Longitud [mts] Longitud [mts] Tiempo [seg].5.6 Tiempo [seg] (b) (c) Figura 3.3: Resultados de la deconvolución aplicada sobre un conjunto CSPG utilizando ε = 3 : (a) conjunto de señales dispersas estimadas; (b) conjunto de señales reconstruidas; (c) conjunto de señales originales. Como se dijo en la introducción de la presente aplicación, la deconvolución sólo representa una parte importante en las aplicaciones de caracterización del subsuelo a través de sísmica de reflexión. Los resultados obtenidos a través de ésta permiten magnificar zonas de los registros distorsionadas por la onda sísmica y por el ruido en la medida de los registros. En la Figura 3.32 se muestran los resultados obtenidos con otros métodos de deconvolución. Los métodos de DD mediante filtro predictor, filtro de Wiener y deconvolución en el dominio de la frecuencia de las Figuras 3.3(a), 3.3(b) y 3.3(c), respectivamente, presentan resultados muy similares. Esto tiene una explicación bastante clara debido a que, como se dijo en el capítulo anterior, estos métodos conducen, bajo la suposición de estadística Gausiana, a la solución LS. Estas soluciones no tienen un carácter disperso y, como se observa, apenas se nota el efecto de compresión de onda sísmica [] propio de la DD de señales sísmicas. La Figura 3.32(d) presenta los resultados obtenidos a través de una variante del algoritmo MG para DD ciega desarrollado en [5]. La solución mediante MG ciega solo permite realizar estimaciones de aquellas zonas de la traza con energía alta, aquellas zonas donde la energía de la señal es muy baja no permiten la obtención de valores de señal dispersa.

134 6 CAPÍTULO 3. DECONVOLUCIÓN DE SEÑALES DISPERSAS MEDIANTE SVM Longitud [mts] Longitud [mts] Tiempo [seg].5.6 Tiempo [seg] (a) Longitud [mts] (b) Longitud [mts] Tiempo [seg].5.6 Tiempo [seg] (c) (d) Figura 3.32: Resultados de la deconvolución aplicada sobre un conjunto CSPG: (a) señal dispersa obtenida mediante DP; (b) señal dispersa obtenida mediante deconvolución Wiener; (c) señal dispersa obtenida con método de frecuencia; (d) señal dispersa obtenida mediante MG ciego. Como nota final puede decirse que la solución mediante SVM produce resultados dispersos de señal. Por otra parte, las otras soluciones conducen a estimaciones muy malas de la señal dispersa. Es decir, si se reinterpreta la deconvolución como un método de compresión de ondas [], los métodos de DD mediante filtro predictor, filtro de Wiener y deconvolución en la frecuencia no son capaces de efectuar una compresión adecuada de la onda sísmica. Cabe decir también, que la principal desventaja del algoritmo SVM-MSKA es su coste computacional alto Conclusiones La principal aportación de este capítulo ha sido la introducción de dos algoritmos de DD no ciega. El algoritmo MPS presenta buenas prestaciones en cuanto a regularización, pero su aplicación directa se ve limitada por que la solución no tiene la dispersión buscada en este tipo de aplicaciones. El algoritmo MSKA surge como una modificación del algoritmo basado en el MDS. Las principales características del MSKA son su solución dispersa,

135 3.6. CONCLUSIONES 7 intrínseca en la naturaleza de los algoritmos SVM, y su robustez frente a ruido de estadística no Gausiana. A través de varios experimentos se ha mostrado las capacidades de los algoritmos frente a diferentes situaciones: tipos de ruido (Gausiano, Laplaciano, uniforme e impulsivo), tipo de señal de entrada (aleatoria y determinista) y características en la fase (mínima y no mínima). También, se han probado las capacidades del algoritmo MSKA en datos reales correspondientes a un experimento de geofísica: la reflexión sísmica. Los resultados obtenidos a través de dicho experimento resaltan el carácter disperso de la solución del algoritmo SVM para deconvolución.

136 8 CAPÍTULO 3. DECONVOLUCIÓN DE SEÑALES DISPERSAS MEDIANTE SVM

137 Capítulo 4 Deconvolución homomórfica En este capítulo se tiene interés en un tipo de señal dispersa de características peculiares, que posee una transformada de Fourier igual al cociente de dos señales paso bajo. La necesidad de disponer de un algoritmo que se adapte a este caso particular ha conducido en la presente Tesis a estudiar, aplicar, y ampliar algoritmos de DH en la estimación de señales dispersas. La DH está ligada a dos conceptos muy importantes: los sistemas homomórficos y el análisis cepstral. El concepto de sistema homomórfico fue introducido en [84], y permitió la generalización del principio de superposición de los sistemas lineales. Casi al mismo tiempo, y de forma separada, fue introducido el concepto de cepstrum [], el cual se definió como la transformada inversa de Fourier (TIF) del logaritmo del espectro de potencia de una señal. El cepstrum corresponde a un caso particular de sistema homomórfico [86], aquél para el que la operación entre señales se define por el operador de convolución. Este capítulo está formado por cinco secciones principales. La primera sección hace una introducción teórica al cepstrum complejo, y, particularmente, al algoritmo de DH, estudiando los conceptos más importantes, así como los problemas en su implementación. La segunda sección analiza el problema de la DH de señales limitadas en banda, revisando los métodos que proponen la traslación de la señal a un espacio de señal de banda completa. La tercera sección formula un problema particular de DD, donde la transformada de Fourier de la señal dispersa está formada por el cociente de dos señales limitadas en banda. La cuarta sección presenta la segunda contribución de la presente Tésis, esto es, se utiliza una modificación al algoritmo de DH propuesto en [7] para su aplicación en problemas donde la señal dispersa en el dominio de la frecuencia está formada por el cociente de dos señales paso bajo; dicha DD se realiza trasladando las señales limitadas en banda a un dominio de señales de banda completa. Finalmente se realiza un conjunto de simulaciones y experimentos sobre señales sintéticas que permiten evaluar las prestaciones de los algoritmos. 9

138 CAPÍTULO 4. DECONVOLUCIÓN HOMOMÓRFICA 4.. Generalidades del cepstrum complejo El cepstrum fue empleado en [] utilizando únicamente la magnitud de la transformada de Fourier de la señal. La incorporación de la fase condujo al concepto más general denominado cepstrum complejo (CC) [86]. El CC no debe interpretarse como que éste es una función compleja, de hecho el CC de una señal real es una señal real [87]. La presente sección está formada por cuatro subsecciones. En la primera se define formalmente el CC, introduciendo la notación a utilizar en el resto del capítulo. La segunda explica el concepto de CC dentro del marco de los sistemas homomórficos. En la tercera, y a partir de la representación canónica, se muestra la forma de calcular numéricamente el valor del CC. Finalmente, se introduce la DH desde el enfoque de la Teoría de Sistemas Definición del CC Considérese una señal de tiempo discreto x[n] cuya transformada Z (TZ) se expresa como a continuación: X(z) = + n= x[n]z n = X(z) e j X(z) (4.) donde X(z) y X(z) representan la magnitud y el ángulo de la TZ compleja X(z), respectivamente. El logaritmo complejo de X(z) se define [33] como: ˇX(z) = log[x(z)] = = log [ X(z) e j X(z)] = = log X(z) + j X(z) (4.2) El CC ˇx[n] de la señal de tiempo discreto x[n], se calcula aplicando la TZ inversa: ˇx[n] = [ ] log X(z) z n dz (4.3) 2πj donde el CC ˇx[n] debe cumplir con la condición necesaria de ser una señal estable [87], para lo cual la región de convergencia, donde se aplica la integración de contorno, debe incluir la circunferencia de radio unidad. Esta condición define la TZ del CC ˇX(z) como una señal de tiempo discreto estable, es decir, ésta posee una representación convergente en series de potencias [87]. La condición de estabilidad del CC permite representarlo utilizando la TIF, como sigue: ˇx[n] = π ˇX(e jω )e jωn dω (4.4) 2π π

139 4.. GENERALIDADES DEL CEPSTRUM COMPLEJO x[n]. y[n] H x[n] + + D L D -. x[n] ^ + + y[n] ^. y[n] (a) (b) Figura 4.: Principio generalizado de superposición: (a) el concepto de transformación lineal ligado al principio de superposición se ha generalizado a través de una transformación no lineal, H, definida como un sistema homomórfico; (b) un sistema homomórfico representado a través de tres subsistemas. donde ˇX(e jω ) es la transformada de Fourier de ˇx[n], que corresponde al caso particular z = e jω de (4.2) y que se escribe como: ˇX(e jω ) = log X(e jω ) + j X(e jω ) (4.5) donde el ángulo X(e jω ) representa la respuesta en frecuencia de la fase, la cual tiene un carácter ambiguo, es decir, se puede sumar cualquier múltiplo entero de 2π en cualquier valor de frecuencia sin que cambie ˇX(e jω ) [33]. La existencia de ˇX(e jω ) requiere que ésta sea periódica con periodo igual a 2π, y por tanto, la función de fase X(e jω ) debe ser una función de ω, continua y de periodo 2π. Además, el CC de señales reales es, también, una señal real, y por tanto ˇX(e jω ) es una función compleja simétrica conjugada. Es decir, la función de fase es una función continua de ω, periódica con periodo 2π y de simetría impar. El cálculo numérico de la función de fase es un tema abierto de investigación [6, 82] y los resultados obtenidos en DH dependen de las buenas prestaciones que posean los algoritmos al enfrentar este problema. En la siguiente sección se introduce el análisis cepstral dentro del marco de los sistemas homomórficos. Dicho enfoque permite estudiar la DH desde el punto de vista de análisis de sistemas Generalización del principio de superposición Algebraicamente, el principio de superposición se interpreta como una transformación lineal en el mismo espacio vectorial al que pertenece la señal de entrada x[n] [24]. Matemáticamente, sea T la transformación que el sistema hace sobre la señal de entrada, sean dos señales de entrada cualesquiera x [n] y x 2 [n], y sea un escalar c C, se tiene que, por el principio de superposición [87] se cumplen las dos propiedades siguientes: T { x [n] + x 2 [n] } = T { x [n] } + T { x 2 [n] } (4.6) T { cx [n] } = ct { x [n] } (4.7)

140 2 CAPÍTULO 4. DECONVOLUCIÓN HOMOMÓRFICA La extensión del principio de superposición (4.6)-(4.7) a casos no lineales fue introducida en [84, 86] y se denominó principio de superposición generalizado. Matemáticamente, si se denota por H a la operación de transformación generalizada, al operador que combina las señales de entrada, al operador que combina señales de entrada con escalares, al operador que combina las señales de salida y al operador que combina señales de salida con escalares, es posible extender las propiedades que definen el principio de superposición (4.6)-(4.7), como sigue: H { x [n] x 2 [n] } = H { x [n] } H { x 2 [n] } (4.8) H { c x[n] } = c H { x[n] } (4.9) Los sistemas que cumplen con (4.8)-(4.9) se denominan sistemas lineales generalizados o sistemas homomórficos, debido a la analogía que existe con el concepto de homomorfismo de la Teoría de Grupos, interpretándose como transformaciones algebraicas en espacios de señal [87]. Tal como se ilustra en la Figura 4., cualquier sistema homomórfico puede descomponerse en tres subsistemas homomórficos conectados en cascada, y esta configuración se conoce como representación canónica [87]. El primer sistema depende del operador de entrada, el segundo es un sistema lineal L, y el tercero es un sistema que depende solo del operador de salida. El sistema asociado al operador de entrada D [ ] se conoce como sistema característico de entrada, y al sistema asociado al operador de salida D [ ] como sistema característico de salida. En la siguiente subsección se estudia un tipo de sistema característico utilizado en la deconvolución de señales Sistema característico para la convolución El CC de una señal de tiempo discreto puede formularse, dentro del marco de los sistemas lineales generalizados, a través de la utilización de un sistema homomórfico que transforma señales convolucionadas a la entrada en la suma de éstas a la salida. De esta manera, el CC de una señal de tiempo discreto x[n] = x [n] + x 2 [n] define como operador de entrada al operador de convolución, al operador de salida como el operador suma +, y a la transformación no lineal como el sistema característico D [ ]. Aplicando las propiedades de superposición generalizada (4.8)-(4.9), se tiene: ] ] D [x[n] = D [x [n] x 2 [n] = ] [ ] = D [x [n] + D x 2 [n] = = ˇx [n] + ˇx 2 [n] = ˇx[n] (4.) donde las señales ˇx [n], ˇx 2 [n], ˇx[n] corresponden al CC de x [n], x 2 [n] y x[n], respectivamente.

141 4.. GENERALIDADES DEL CEPSTRUM COMPLEJO 3 XO x[n] OX D OX x[n] ^ L (a) y[n] ^ - D OX XO y[n] x[n] DFT log DFT - x[n] + ^ D OX x[n] ^ + DFT exp DFT - OX x[n] - D OX (b) Figura 4.2: Operador de convolución como sistema característico: (a) representación canónica; (b) aproximación usando la DFT. Utilizando el sistema característico D [ ] y su sistema inverso D [ ], la Figura 4.2(a) muestra la representación canónica de la clase de sistemas homomórficos cuyas señales de entrada y salida están combinadas a través del operador de convolución. Para una señal de entrada x[n] = x [n] x 2 [n], la salida está dada por: y[n] = y [n] y 2 [n] (4.) donde y [n] e y 2 [n] son las respuestas del sistema a las entradas x [n] y x 2 [n], respectivamente. La implementación numérica de los sistemas característicos D [ ] y D [ ] se muestra en la Figura 4.2(b), donde las DFT se calculan utilizando el algoritmo de la FFT. Además, el sistema L puede ser cualquier sistema lineal, y su papel es el de eliminar o modificar cualquiera de las componentes del CC de entrada. Este tipo de aplicaciones se introdujo en el Capítulo 2 como DH, y en las secciones siguientes se presentan más detalles sobre su implementación Deconvolución homomórfica de señales Existen algunas consideraciones a tener en cuenta antes de usar el sistema lineal L de la Figura 4.2(a) como método de DH. Primero, debe observarse que la entrada a dicho sistema no es una señal de tiempo discreto, sino el CC de ésta. En segundo lugar, L no es invariante en el tiempo, sino en la frecuencia [87], lo que modifica el concepto de filtrado de frecuencias. La

142 4 CAPÍTULO 4. DECONVOLUCIÓN HOMOMÓRFICA Magnitud Normalizada x 2 [n] ^ - ^ x [n] n Figura 4.3: Caso ideal en DH: la componente de señal de la izquierda, de variación lenta, y la de la derecha, de variación rápida, están claramente separadas. relación entre la señal de entrada y la señal de salida en este tipo de sistemas viene dada por la ecuación: ˇy[n] = l[n]ˇx[n] (4.2) donde l[n] es el sistema lineal invariante en la frecuencia, y desempeña la función de filtro homomórfico. Los sistemas invariantes en la frecuencia son de mucha utilidad práctica cuando las componentes del CC no se solapan. Por ejemplo, en el caso ideal que se muestra en la Figura 4.3, la señal de entrada x[n] = x [n] x 2 [n] tiene un CC ˇx[n] = ˇx [n] + ˇx 2 [n] cuyas componentes están suficientemente separadas. Así, se puede decir que el CC está formado por dos partes: una componente que cambia suavemente, que corresponde a ˇx [n], y otra que experimenta cambios muy rápidos, formada por ˇx 2 [n]. Cada una de las componentes del CC puede separarse únicamente si éstas no se solapan. Es decir, existe un punto n tal que: ˇx [n] =, n n ˇx 2 [n] =, n n Para el caso separable, como puede verse, es necesario encontrar el punto n que hace el papel de frecuencia de corte en los filtros homomórficos. Por ejemplo, si se quiere eliminar la componente de oscilación rápida se puede recurrir a un filtro cuya respuesta en el dominio temporal sea como la que se muestra en la Figura 4.4(a). Si, en cambio, se desea eliminar la componente de oscilación lenta, el filtro homomórfico a aplicar debe tener las características presentadas en la Figura 4.4(b). La selección de los parámetros libres N y N 2 se obtienen a partir del conocimiento a priori de las características de la señal a deconvolucionar, y algunos criterios para su selección se presentan en [2, 3, 37]. Además, debido al carácter periódico del CC, el filtro invariante

143 4.2. DH DE SEÑALES LIMITADAS EN BANDA 5 l [ n] l [ n] - N N 2 N N N N N N 2 N-N 2 N 2 N (a) (b) Figura 4.4: Filtros en el dominio cepstral: (a) paso corto; (b) paso largo. en frecuencia debe de poseer periodicidad, como puede verse en las Figuras 4.4(a) y 4.4(b) DH de señales limitadas en banda En muchas de las aplicaciones que se encuentran en la práctica de procesado de señal se trabaja con señales limitadas en banda. Esto representa un problema para el cálculo directo del CC, ya que la función logaritmo complejo calculada sobre los intervalos de respuesta nula están indeterminados. A lo largo del desarrollo y refinamiento de los métodos de DH, se han propuestos métodos que buscan solucionar este problema, a través de la conversión de las señales limitadas en banda a señales de banda completa. En [85], la transformación se realiza interpolado, diezmando y filtrado en el dominio de la frecuencia. El problema que presenta este método está relacionado con los valores a los que hay que ajustar los interpoladores, diezmadores y las frecuencias de corte de los filtros. Posteriormente, se propuso en [7] construir una señal suplemental de banda completa que se añade a la señal original [7]. Este método se describe en la presente sección. Considérese una señal paso banda x[n] = r[n] w[n] donde r[n] es un tren de impulsos y w[n] es una onda que satisface las propiedades del CC. Se definen las señales X rec (z) e Y (z) como sigue: X rec (z) = X(z) + ɛ (4.3) Y (z) = X rec (z) + S(z) (4.4) donde ɛ es una constante que se agrega a X(z) evitando valores de frecuencia cercanos o iguales a cero, y la señal S(z) es la señal suplemental que se define como: S(z) = X rec (z) + S i (z) (4.5) donde la señal S i (z) presenta las siguientes características: es una señal IIR de segundo orden, de fase mínima y de espectro suave;

144 6 CAPÍTULO 4. DECONVOLUCIÓN HOMOMÓRFICA sus polos están localizados en: π + ω 2, si ω π ω 2 ω c = 2 ω c, otro caso; 2 su ancho de banda se calcula a partir del radio del polo según la relación BW 2( r)/ r. El procedimiento constructivo se empieza reescribiendo (4.4) como sigue: Si se define Φ(z) = Xrec(z) S(z) (4.6), se tiene: [ Xrec (z) ] Y (z) = S(z) + S(z) (4.6) + y se aplica la función logaritmo complejo a ˇY (z) = log [ S(z) ] + log [ Φ(z) ] = = Š(z) + ˇΦ(z) (4.7) Considerando la condición Xrec(z) < es posible escribir ˇΦ(z) en serie S(z) de potencias como sigue: [ ] n ( ) ˇΦ(z) n X rec (z) = = n S(z) n= [ ]{ [ ] [ ] X rec (z) = X rec (z) + 2 } X rec (z)... = S(z) 2 S(z) 3 S(z) = [ X rec (z) S(z) ] Γ(z) (4.8) donde Γ(z) se define como: Γ(z) = { 2 [ ] [ ] X rec (z) + 2 } X rec (z)... S(z) 3 S(z) (4.9) Si se define la función logaritmo complejo de ˇΦ(z) como Ψ(z) = log [ˇΦ(z) ] y se aplica en (4.7) y (4.8) se tiene: Ψ(z) = log [ˇΦ(z) ] = (4.2) = log [ ˇY (z) Š(z) ] = (4.2) = log [ X rec (z) ] log [ S(z) ] + log [ Γ(z) ] = (4.22) = ˇX rec (z) Š(z) + ˇΓ(z) (4.23)

145 4.3. DH COMO IDENTIFICADOR DE SISTEMAS 7 de donde se obtiene una expresión para ˇX rec (z) en función de la TZ de las observaciones, la TZ de la señal suplemental y la serie Γ(z), como sigue: ˇX rec (z) = log [ ˇY (z) Š(z) ] + Š(z) ˇΓ(z) (4.24) La expresión para ˇX rec (z) obtenida a través del procedimiento anterior corresponde a una transformación de una señal paso banda a una de banda completa [7]. Este procedimiento resuelve el problema de la función logaritmo complejo cuando se aplica a señales de banda limitada, ya que se ha incrementado aquellos valores de potencia que producen resultados indeterminados. También, X rec (z) está bien definida ya que siempre es posible construir una señal suplemental S(z) y, además, las señales ˇY (z), Š(z) y ˇΓ(z) están bien definidas [7]. En la Tabla 4. se resume el algoritmo para la DH usando señal suplemental, tal como fue propuesto en [7, 4]. En la siguiente sección se formula un problema particular de señal donde las señales de entrada y salida del SLIT son señales periódicas y de banda limitada DH como identificador de sistemas Como se definió en el Capítulo 2, la DD no ciega parte del conocimiento de la señal de salida y la respuesta al impulso o la señal de entrada del SLIT. En el Capítulo 3 se desarrollaron ejemplos y simulaciones donde se conocía la entrada y la respuesta al impulso. Por otro lado, en esta sección se estudia un modelo de señal donde se conoce la señal de salida y la señal de entrada, con la peculiaridad añadida de que éstas son periódicas y de banda limitada. La Figura 4.5 muestra un esquema donde se ilustra el modelo de señal a utilizar, problema que también suele conocerse como identificación de sistemas. En esta aplicación se considera que el sistema a identificar corresponde a una señal dispersa de longitud finita. La función de transferencia del sistema de la Figura (4.5) está dada por: I(z) = P (z) Q(z) (4.25) donde Q(z), P (z) e I(z) corresponde a la TZ de la señal de entrada, la señal de salida y la respuesta al impulso, respectivamente. Tal como se describió en el Capítulo 2, la solución directa mediante DFT es muy sensible a errores en los datos, introducidos por la componente de ruido. Además Q(z) es de banda limitada y por tanto I(z) se vuelve indeterminada. Este problema se puede resolver si, en lugar de realizar el cociente directo de (4.25), se calcula el logaritmo complejo de ésta, como sigue: Ǐ(z) = ˇP (z) ˇQ(z) (4.26)

146 8 CAPÍTULO 4. DECONVOLUCIÓN HOMOMÓRFICA Paso. Calcular el espectro X(e jw ) de la secuencia x[n] y determinar de manera aproximada el ancho de banda de la señal suplemental (BW) que hace que ésta sea de banda completa. Paso 2. Calcular las frecuencias de los polos w c de la señal s i [n] a partir de las frecuencias de corte w y w 2. Paso 3. Calcular el radio del polo r a partir del ancho de banda. Paso 4. Con los valores de radio y frecuencia de localización de los polos construir la señal S i como la respuesta de un sistema IIR de segundo orden. Paso 5. Calcular X rec (e jw ) y S(e jw ) = X rec (e jw ) + S i (e jw ). Paso 6. Calcular las constantes G, β y β 2 y escalar las señales del paso anterior tal como se explica en [7]. Paso 7. Crear Y (z) con las señales escaladas del paso anterior y obtener ˇX rec (z). Paso 8. A partir de ˇX rec (z), obtener ˇx rec [n] mediante filtrado homomórfico. Tabla 4.: Algoritmo de DH usando el método de la señal suplemental. Aplicando la TZ inversa se obtiene la expresión del CC de la respuesta al impulso: ǐ[n] = ˇp[n] ˇq[n] (4.27) donde ˇq[n], ˇp[n] e ǐ[n] corresponde al CC de la entrada, la salida y la respuesta al impulso, respectivamente. Puede observarse que en (4.27) se ha resuelto el problema de la división por valores de frecuencias iguales o cercanos a cero de Q(z), ya que el cociente de (4.25) se ha transformado en una diferencia de señales. Pero la estimación de la respuesta al impulso no está del todo resuelta, ya que las funciones ˇQ(z) y ˇP (z) se vuelven indeterminadas para aquellos valores de Q(z) y P (z), respectivamente, con valores cercanos o iguales a cero. La solución a este problema se encuentra en la transformación previa de las señales de banda limitada a señales de banda completa, tal como se explicó en la sección anterior. Este nuevo dominio de señal se conoce con el nombre de Dominio Suplemental (DS) y su aplicación en la estimación de la respuesta al impulso (4.25) se estudia en la siguiente sección.

147 4.4. DH EN EL DOMINIO SUPLEMENTAL 9 Q P I Figura 4.5: SLIT con entrada y salida periódicas de banda limitada DH en el Dominio Suplemental El procedimiento para estimar la respuesta al impulso de (4.25) sigue los pasos establecidos en la Sección 4.2. Así, se definen las señales P rec (z) y Q rec (z) como sigue: P rec (z) = P (z) + ɛ (4.28) Q rec (z) = Q(z) + ɛ (4.29) A partir de las cuales se construye sus respectivas señales suplementales como: S P (z) = P rec (z) + S i (z) (4.3) S Q (z) = Q rec (z) + S i (z) (4.3) donde la señal S i (z) es la misma para Q(z) y P (z) y se calcula siguiendo el algoritmo descrito en la Tabla 4.. Utilizando (4.3) y (4.3) se construyen las señales de banda completa siguientes: Y P (z) = P rec (z) + S P (z) (4.32) Y Q (z) = Q rec (z) + S Q (z) (4.33) donde Y Q (z) e Y P (z) son señales de banda completa. El cociente entre (4.32) y (4.33) está dado por: donde Φ P complejo a (4.34) se obtiene: = P rec(z) S P (z) Y I (z) = Y P (z) Y Q (z) = P rec(z) + S P (z) Q rec (z) + S Q (z) = [ ] S P (z) Prec(z) + S P (z) = [ ] = S Q (z) Qrec(z) S Q + (z) = S P (z)φ P (z) S Q (z)φ Q (z) + y Φ Q = Q rec(z) S Q (z) ˇY I (z) = log[y I (z)] = ˇY P (z) ˇY Q (z) = (4.34) +. Aplicando la función logaritmo = ŠP (z) ŠQ(z) + ˇΦ P (z) ˇΦ Q (z) (4.35)

148 2 CAPÍTULO 4. DECONVOLUCIÓN HOMOMÓRFICA Desarrollando ˇΦ P (z) y ˇΦ Q (z) en series de potencias se tiene: ˇΦ P (z) = P rec(z) S P (z) Γ P (z) (4.36) ˇΦ Q (z) = Q rec(z) S Q (z) Γ Q(z) (4.37) donde las variables Γ P (z) y Γ Q (z) se definen como: Γ P (z) = [ Prec (z) ] + [ Prec (z) ] (4.38) 2 S P (z) 3 S P (z) Γ Q (z) = [ Qrec (z) ] + [ Qrec (z) ] (4.39) 2 S Q (z) 3 S Q (z) donde al igual que en (4.9), los cocientes en (4.38) y (4.39) están sujetos a P rec(z) < y Q rec(z) <. Si se define la función Ψ(z) S P (z) S Q como el logaritmo (z) complejo de Ψ(z) = log[ˇφ P (z) ˇΦ Q (z)] y se sustituye en (4.35) se tiene: Ψ = log[ˇφp (z) ˇΦ ] Q (z) = [ P rec (z) = log S P (z) Γ P (z) Q rec(z) = log { P rec (z) Q rec (z) S Q (z) Γ Q(z) = [ Γ P (z)q rec (z) Q2 rec(z)γ Q (z) S P (z) S Q (z)p rec (z) ] ]} = ˇP rec (z) ˇQ rec (z) + log[res(z)] (4.4) donde RES(z) es una variable introducida para simplificar los términos en series de la ecuación anterior, reescribiendo para ˇP rec (z) ˇQ rec (z) se tiene: Ǐ rec (z) = ˇP rec (z) ˇQ rec (z) = Ψ(z) log[res(z)] (4.4) donde Ǐrec(z) representa la función de transferencia del sistema de la Figura 4.5 en el DS. La Ecuación (4.4) corresponde a la transformación de banda completa de (4.27), es decir, la resta de las señales de presión y flujo se realiza en el DS. Este nuevo enfoque introduce dos cuestiones fundamentales: en primer lugar, el cociente espectral de señales se transforma en una resta espectral de señales y, en segundo lugar, la transformación de las señales de presión y flujo con ancho de banda paso bajo a ancho de banda completa. En una aplicación de este tipo, con señales limitadas en banda, el cambio de cociente espectral a resta espectral de señales y la transformación al DS son cuestiones fundamentales. En la siguiente sección se desarrollan varias simulaciones que permiten ilustrar muchas de las propiedades de la DH tal como han sido presentadas en este capítulo. =

149 4.5. SIMULACIONES Y EXPERIMENTOS Simulaciones y Experimentos A continuación se presentan cuatro experimentos realizados sobres señales sintéticas. Los experimentos llevan un orden incremental de dificultad y permiten evaluar las capacidades y limitaciones de los métodos de DH. El primer experimento compara resultados teóricos, a través de la obtención de una expresión matemática para el CC de una señal, con resultados obtenidos utilizando el método directo (MD). Además, éste experimento permite mostrar la mayor parte de los problemas encontrados en el cálculo numérico del CC. El segundo experimento añade un grado más de complejidad al hacer que la señal de variación suave sea de banda limitada, para ello se utiliza el modelo matemático de una onda sísmica [7]. El tercer experimento es una aproximación sintética al problema que ha motivado este capítulo, y consiste de un sistema lineal disperso cuya entrada y salida están formadas por señales periódicas de banda limitada. El último experimento consiste en una modificación del tercer experimento que incluye el efecto de reflexiones, éste modelo es una aproximación más realista del modelo que caracteriza la impedancia a la entrada de la aorta y que se estudiará en el siguiente capítulo MD de DH Este primer experimento se ha tomado de [87] y tiene por objetivo analizar el CC de señales con diferentes características de fase, mínima o mixta, y características temporales, variación suave o variación impulsiva dispersa. También el ejemplo permite mostrar los problemas que surgen cuando se calcula numéricamente el CC. Las señales a convolucionar se escogen de tal forma que es posible calcular expresiones matemáticas cerradas, del tipo polinómico racional, para sus transformadas Z. Considérese una señal de tiempo discreto que se origina a partir de la convolución de dos señales como se muestra a continuación: y[n] = x[n] h[n] (4.42) donde x[n] y h[n] representan la señal dispersa de entrada y la respuesta al impulso, respectivamente. La señal dispersa está forma de tres impulsos separados N muestras entre sí, con un valor de amplitud decreciente exponencialmente. Matemáticamente se escribe como: x[n] = δ[n] + βδ[n N ] + β 2 δ[n 2N ] (4.43) La respuesta al impulso utilizada es una señal sinusoidal amortiguada decreciente tipo FIR, dada por: h[n] = b w[n] + b w[n ] (4.44)

150 22 CAPÍTULO 4. DECONVOLUCIÓN HOMOMÓRFICA donde w[n] es igual a: w[n] = { [ ]} rn cos(θn) cos θ(n + 2) u[n], θ, π (4.45) 2sin 2 θ y que corresponde a una señal IIR con un polo complejo conjugado situado en re ±jθ. Las Figuras 4.6(a), 4.6(b), 4.6(c) muestran las señales h[n], x[n] e y[n], respectivamente, utilizando los valores de parámetros específicos siguientes: b =.98, b =, β = r =.9, θ = π/6, N = 5. Tal como se observa en las Figuras 4.6(d) y 4.6(e), la respuesta en frecuencia de h[n] es paso bajo con fase no lineal. Esta señal es una versión simplificada de modelos más complejos encontrados en aplicaciones prácticas de análisis y procesado de señal. La simplicidad del modelo puede observarse a partir de la TZ de h[n], que viene dada por: H(z) = b + b z ( re jθ z )( re jθ z ) (4.46) Esta ecuación solo presenta un pico de resonancia. Normalmente, en las aplicaciones de tratamiento digital de voz y sísmica de reflexión se tienen modelos con, al menos, diez picos de resonancia [87]. A pesar de su simplicidad, el modelo presenta todas las propiedades del CC de una señal digital cuya TZ es racional. Cálculo analítico del CC. Primero se reescribe la expresión para la TZ de h[n] de la siguiente manera: ] b z [ + b b z H(z) = (4.47) ( re jθ z )( re jθ z ) El factor z en el numerador de (4.47) conduce a dificultades teóricas, ya que introduce una componente lineal en el ángulo del logaritmo de H(z), produciendo discontinuidades en las frecuencias de fase ω = ±π de Ȟ(z), con lo cual Ȟ(z) no será analítica en el círculo unitario. Para soslayar dicho problema, se multiplica H(z) por z, cuyo efecto en el dominio temporal es un desplazamiento hacia la izquierda de la señal digital, es decir, h[n] se adelanta a h[n + ]. El adelanto también afecta a y[n], y esto deberá tenerse en cuenta al interpretar los resultados finales. Se reescribe la expresión para H(z), sin el factor z, como sigue: ] b [ + b b z H(z) = ( re jθ z )( re jθ z ) (4.48) Aplicando la función logaritmo a (4.48) y luego calculando su TZ inversa,

151 4.5. SIMULACIONES Y EXPERIMENTOS 23 4 Amplitud 3 2 Amplitud n (a) n (b) Amplitud n (c) Amplitud [db] Frecuencia [rads] (d) Radianes Frecuencia [rads] (e) Radianes Frecuencia [rads] (f) Radianes Frecuencia [rads] (g) Figura 4.6: Señales sintéticas: (a) respuesta al impulso; (b) señal dispersa; (c) convolución de ambas; (d) magnitud de Y (e jω ); (e) valor principal del ángulo de Y (e jω ); (f) valor principal del ángulo sin componente de fase lineal de Y (e jω ); (g) fase desenrrollada de Y (e jω ).

152 24 CAPÍTULO 4. DECONVOLUCIÓN HOMOMÓRFICA se tiene que el cepstrum de la respuesta al impulso ȟ[n] viene dado por: log b, n =, ] ȟ[n] = [(re jθ ) n + (re jθ ) n, n >, n (4.49) ( b ) n, n <, n b Para determinar ˇx[n], primero se calcula la TZ de x[n], la cual viene dada por: X(z) = + βz N + β 2 z 2N = β3 z 3N (4.5) βz N Aplicando la función logaritmo complejo a X(z), se obtiene: ( ) ( ) ˇX(z) = log β 3 z 3N log βz N (4.5) y ya que β <, se utiliza la expansión en series de potencias como sigue: ˇX(z) = k= β 3k k z 3N k + k= β k k z N k (4.52) Aplicando la TZ inversa se obtiene el CC para la señal dispersa de entrada, el cual está dado por: ˇx[n] = k= β 3k k δ[n 3N k] + k= β k k δ[n N k] (4.53) El CC de la señal (4.42) se obtiene mediante la suma del CC de la respuesta al impulso (4.49) y de la señal dispersa de entrada (4.53): ˇy[n] = ȟ[n] + ˇx[n] (4.54) Las Figuras 4.7(a) y 4.7(b) muestran el CC de la respuesta al impulso y de la señal dispersa de entrada, respectivamente. Es importante hacer notar las siguientes observaciones: la señal h[n] es de fase mixta y su cepstrum tiempos positivos y negativos; ȟ[n] está definido para la señal h[n] es de variación suave y, por tanto, su cepstrum se concentra entorno al origen (como puede observarse en (4.49), su valor decae exponencialmente); la señal x[n] es de fase mínima y su cepstrum ˇx[n] solo está definido para valores positivos de tiempo;

153 4.5. SIMULACIONES Y EXPERIMENTOS 25 Amplitud n (a) Amplitud n (b) Figura 4.7: CC de: (a) la respuesta al impulso, ˇx[n]. ȟ[n]; (b) la señal dispersa, la señal x[n] es de variación rápida, está formada por un tren de deltas y, por tanto, su cepstrum, como puede observarse en (4.53) corresponde a un secuencia de impulsos localizados en posiciones múltiplos de N ; la componente de fase mínima de ȟ[n] se atenúa antes de alcanzar la posición temporal del primer impulso de ˇx[n]. De los resultados teóricos obtenidos a partir de este ejemplo, se puede concluir que solo es posible separar las señales que componen a (4.42) si la señal de variación suave no se solapa con la señal de variación rápida. Para este caso particular, la separación es posible si la componente de fase mínima de ȟ[n] se atenúa lo suficiente, hacia valores cercanos a cero, antes de alcanzar la posición de la primera muestra de ˇx[n]. Cálculo numérico del CC. La señal (4.42) está formada por 248 muestras, y para visualizar mejor su representación solo se muestra en la Figura 4.6(c) el segmento donde la señal tiene valores diferentes de cero. La Figura 4.6(e) muestra la fase o valor principal del ángulo de (4.46), donde se observan discontinuidades en los valores de fase ±π. La discontinuidad se debe al carácter multivaluado de la parte imaginaria de la función logaritmo complejo [87], y representa uno de los mayores problemas en el cálculo numérico del CC. El paso fundamental consiste en calcular la función de fase verdadera. Este paso se conoce como desenrollado de fase, y ha sido ampliamente estudiado en la literatura [6, 85]. El desenrrollado de la fase requiere primero que se elimine la componente de fase lineal, en este ejemplo multiplicando por z. Las Figuras 4.6(f) y 4.6(g) muestran la fase de la respuesta al impulso sin la componente lineal y el valor verdadero de la fase, respectivamente. En la Figura 4.8(a) se muestra el CC de la señal (4.42), donde se pueden identificar el CC individual de cada componente de la señal. Las Figuras 4.8(b) y 4.8(c) muestran los resultados obtenidos después de deconvolucionar mediante DH la respuesta al impulso y la señal dispersa, respectivamente.

154 26 CAPÍTULO 4. DECONVOLUCIÓN HOMOMÓRFICA.5 Amplitud n (a) 3 h est [n] x est [n] 2 h[n].8 x[n] Amplitud Amplitud n (b) n (c) Figura 4.8: Señal sintética de prueba: (a) CC; (b) respuesta al impulso verdadera y estimada; (c) señal dispersa de entrada verdadera y estimada Obsérvese que en las Figuras 4.8(b) y 4.8(c) ya no se utiliza la representación de tiempo discreto de señales. La representación de trazo continuo facilita la comparación de las señales estimadas con las originales y se utilizará en los siguientes experimentos. En general las estimaciones tanto de la respuesta al impulso como de la señal dispersa se realizan adecuadamente, es decir, el MD de DH permite separar ambas componentes. El MD funciona bien en este ejemplo debido a las siguientes cuestiones: el cepstrum debido a la componente de fase mínima de la respuesta al impulso no se solapa con el cepstrum de la señal dispersa; es posible llevar a cabo el desenrrollado de la fase; la respuesta al impulso tiene un ancho de banda amplio; la ausencia de ruido aditivo en el modelo de señal (4.42); no existe periodicidad en el modelo de señal.

155 4.5. SIMULACIONES Y EXPERIMENTOS Amplitud.4.2 Amplitud n (a) n (b).6.4 Amplitud n (c) Y(e jω ) Frecuencia [rad] (d) Y(e jω ) Frecuencia [rad] (e) Figura 4.9: Señal sintética de banda limitada: (a) onda Ricker; (b) señal dispersa formada por tres impulsos; (c) señal sísmica de salida; (d) amplitud del espectro; (e) fase del espectro Señales limitadas en banda En la Figura 4.6(d) se observa que la señal del modelo sintético de (4.42) tiene un ancho de banda que se extiende por casi todo el espectro de frecuencias. En muchas situaciones prácticas se tiene que la señal tiene un ancho de banda muy estrecho, como por ejemplo, las señales sísmicas, y que además tienen una respuesta en frecuencia paso banda. Las señales limitadas en banda presentan dificultades añadidas en la apli-

156 28 CAPÍTULO 4. DECONVOLUCIÓN HOMOMÓRFICA cación del algoritmo de la DH. Como se estudió en las secciones precedentes, la principal fuente de dificultades proviene de la necesidad de calcular el logaritmo complejo de señales cuya TZ tiene valores cercanos o iguales a cero. Considere la señal de variación lenta definida por la onda Ricker, que es igual a la segunda derivada de la función Gausiana [7, 4], ampliamente utilizada en la generación sintética de ondas sísmicas, la cual esta dada por: h[n] = π [ (.6n 3.464) 2 ] e 2 [(.6n 3.464) 2 ] (4.55) La señal dispersa de entrada o señal de variación rápida está formada por tres impulsos, como a continuación: x[n] = δ[n] +.75δ[n 2] δ[n 24] (4.56) Las Figuras 4.9(a), 4.9(b) y 4.9(c) muestran la onda Ricker, la señal dispersa y la onda sísmica, respectivamente. La onda sísmica está formada por la convolución entre la señal dispersa y la onda Ricker. También, las Figuras 4.9(d) y 4.9(e) muestran la amplitud y la fase del espectro, respectivamente. La amplitud del espectro muestra una señal de banda limitada que dificulta o, en el peor de los casos, impide el calculo del CC. Además se observa como la fase del espectro presenta discontinuidades en aquellos valores que superan el valor de ±π, y esto hace necesaria la introducción de técnicas de desenrrollado de fase. A diferencia del experimento anterior, en este experimento no es posible obtener una expresión cerrada para el CC de la onda Ricker. Las Figuras 4.(a) y 4.(c) muestran el CC de la onda Ricker y la señal sísmica, respectivamente, obtenidos de forma directa utilizando el esquema de la sección Por otro lado, la señal dispersa (4.56) tiene una expresión cerrada dada por (4.52). En principio el cepstrum de (4.56) debería de ser identificable, ya que está formado por una secuencia de impulsos múltiplos del factor de separación 2, pero de la Figura 4.(c) se observa que no es posible determinarlo con claridad. Las Figuras 4.(b) y 4.(d) muestran el CC de la onda Ricker y de la señal sísmica, respectivamente, calculado mediante dos métodos de DH: el primero se basa en la transformación de la señal sísmica al DS y el segundo en el cálculo de las raíces. El método de DH basado en el cálculo de las raíces (MR) se explica en el Apéndice B, y en esta sección sirve para validar los resultados obtenidos con el el método de transformación al DS. Su utilización como método de validación se debe a que, para señales no ruidosas de longitud pequeña, que no superan las cuatro mil muestras, el valor del cepstrum es muy preciso, debido principalmente a que no requiere el cálculo del logaritmo complejo, ni del desenrrollado de la fase. Como se comentó en el Experimento se utiliza representación de trazo continuo de señales, se hace de esta forma para mantener coherencia con

157 4.5. SIMULACIONES Y EXPERIMENTOS MR DS Amplitud Amplitud n (a) n (b) 6 4 MR DS Amplitud Amplitud n (c) n (d).8.6 Original DS.8 Original DS Amplitud.4.2 Amplitud n (e) n (f) Figura 4.: CC para las señales: (a) onda Ricker (MD); (b) onda Ricker (MR y DS); (c) onda sísmica (MD); (d) onda sísmica (MR y DS);(e) onda Ricker original y estimada; (f) señal dispersa original y estimada. los resultados comúnmente encontrados en la literatura [85, 33, 37]. En la Figura 4.(d) se identifica con bastante claridad la componente de variación rápida ˇx[n] de la componente de variación lenta ȟ[n]. Las Figuras 4.(e) y 4.(f) muestran la onda Ricker y la señal dispersa estimadas, respectivamente. Estas señales se han estimado utilizando el algoritmo de DH con transformación DS como se describe a continuación: las señales tienen una frecuencia de muestreo de 25 Hz y una longitud

158 3 CAPÍTULO 4. DECONVOLUCIÓN HOMOMÓRFICA Retraso Señal dispersa de prueba n Valor impulso MR DS Tabla 4.2: Resultados del experimento con la señal sísmica. de 248 muestras; las frecuencias que determinan el ancho de banda son ω = y ω 2 =.6786, de donde se obtiene el ancho de banda BW = 2 Hz, la frecuencia de resonancia de polos ω c =.9 radianes y la longitud de polo r=.7783; se aplica el algoritmo de la Tabla 4. obteniéndose las señales X(e jω ) S(e jω ) e Y (e jω ); se obtiene el cepstrum de la onda sísmica; se obtiene las componentes de variación suave ȟ[n] y de variación rápida ˇx[n]. Ya que este experimento se ha llevado a cabo con un número pequeño de muestras, es posible validar los resultados utilizando el MR. La Tabla 4.2 muestra los resultados obtenidos en la estimación de la señal dispersa (4.56), utilizando ambos métodos. Los resultados son muy parecidos, pero el algoritmo basado en la transformación al DS obtiene aproximaciones ligeramente mejores y más próximas a los resultados reales Señales periódicas limitadas en banda La DH de señales periódicas presentan un grado añadido de dificultad, debido a que, desde el punto de vista de los sistemas lineales, éstas se interpretan como la convolución entre un ciclo aislado y una secuencia de impulsos unitarios separados uniformemente y de longitud igual a. Esta última secuencia tiene ceros en la circunferencia de radio unidad, imposibilitando teóricamente el cálculo del CC [6, 87]. En el presente experimento se ha escogido una señal, representada en la Figura 4.(a), modelada por un ciclo aislado de señal mediante una parábola invertida, tomando sus valores positivos, de duración No seguida de una 3 deflexion negativa, de duración igual a No. Matemáticamente, se define como 5 sigue: 36 n n, n =,..., ( h[n] = No 2 N o )N o (4.57), otro caso

159 4.5. SIMULACIONES Y EXPERIMENTOS 3 Amplitud t[segundos] (a) Amplitud t[segundos] (b).8 Amplitud t[segundos] (c) H(e jω ) Frecuencia [rad] (d) Y(e jω ) Frecuencia [rad] (e) Figura 4.: Modelo de parábola invertida: (a) representación de un ciclo; (b) tren de dos impulsos; (c) señal formada por dos ciclos; (d) amplitud del espectro; (e) fase del espectro. Para representar la periodicidad de una forma simple se convoluciona (4.57) con la secuencia formada por dos impulsos como sigue: x[n] = δ[n] + δ[n N o ] (4.58) donde N o representa la separación entre ambos impulsos, tal como se muestra en la Figura 4.(b). La transformada de Fourier de (4.58) está dada por: X(e jω ) = 2cos(ωN o /2) e jωn o/2 (4.59)

160 32 CAPÍTULO 4. DECONVOLUCIÓN HOMOMÓRFICA.5 Directo Raíces Amplitud t[segundos] (a) Amplitud t[segundos] (b) Amplitud t[segundos] (c) Figura 4.2: CC del modelo de parábola invertida: (a) MD (trazo continuo) y basado en las raíces (trazo discontinuo); (b) CC de dos ciclos del modelo de parábolas calculado con el MD; (c) utilizando rellenado con ceros. La transformada (4.59) presenta ceros espectrales, o valores para los cuales ésta es igual a cero, y el cálculo de su CC se vuelve problemático debido a que la función logaritmo complejo no está definida en esos valores. Por otra parte, no es posible calcular una expresión matemática cerrada para el CC de la señal de variación suave (4.57). La Figura 4.2(a) muestra el CC ȟ[n] definido principalmente para valores de tiempo negativo, lo que indica el predominio de su componente de fase mínima. El CC ȟ[n] se ha calculado utilizando el MD con desenrrollado de fase. También, al igual que en el experimento anterior, se utiliza el MR para validar este resultado. Es importante hacer notar que el MR solo puede utilizarse en los casos donde las raíces de la señal no están en la circunferencia de radio unidad. De esta manera el método no funciona en la señal (4.58), ya que la componente formada por dos impulsos introduce ceros espectrales. La Figura 4.2(b) muestra el CC para dos ciclos de (4.57) calculado con el MD, utilizando una FFT de 24 puntos y una frecuencia de muestreo de 5 Hz. En principio no es posible, desde un punto de vista teórico, calcular el CC de esta señal, ya que la componente ˇx[n] no está determinada. Sin embargo, en este caso particular, es posible calcular el CC de la convolución entre (4.57) y (4.58), ya que debido a las limitaciones de la precisión numérica los ceros

161 4.5. SIMULACIONES Y EXPERIMENTOS 33 Número de ciclos Método MD Sí Sí No Sí Sí Sí No Sí Sí Sí No Sí Sí DS Sí Sí Sí Sí Sí Sí Sí Sí Sí Sí Sí Sí Sí Tabla 4.3: Determinación de la señal dispersa mediante MD y DS a partir de varios ciclos del modelo de señal de parábola invertida. espectrales desaparecen. Pero también, el valor estimado del CC introduce ruido en todo el cepstrum, distorsionando la parte de señal correspondiente a ȟ[n]. Es posible reducir el nivel de ruido utilizando la técnica de rellenado con cerros. En [23] se sugieren valores de rellenado con ceros iguales a veces la longitud de la señal. La Figura 4.2(c) muestra el CC de dos ciclos de señal, utilizando FFT de longitud 65,536. Puede observarse que el rellenado con ceros mejora las prestaciones del MD reduciendo el nivel de ruido, pero sin llegar a realizar una buena estimación de ȟ[n]. El caso más general de una señal periódica modelada a partir de una secuencia de N impulsos se escribe como sigue: x[n] = N k= δ[n kn o ] (4.6) donde N o representa el periodo. La transformada de Fourier de (4.6) está dada por: X(e jω ) = 2j sen(ωn on/2) sen(ωn o /2) ejωno( N)/2 (4.6) la cual presenta ceros espectrales que dificultan el cálculo directo del CC. La Tabla 4.3 identifica los experimentos para los cuales es posible obtener una estimación de la señal dispersa a partir de una señal formada por varios ciclos de la señal (4.57). De los resultados obtenidos con el MD se observa que, debido a las limitaciones en precisión numérica, no siempre se tienen ceros espectrales en una señal periódica. Sin embargo, hay varios casos donde no se puede determinar el CC de una señal periódica. También, de la Tabla 4.3 se observa que el método de DH basado en la transformación al DS permite estimar la señal dispersa en todos los casos. Las Figuras 4.3(a) y 4.3(b) muestran los resultados obtenidos al aplicar la DH en DS a 2 ciclos de la parábola invertida h[n]. Se observa la existencia de un nivel de ruido presente en la estimación de la señal dispersa. Por otro lado, el algoritmo obtiene buenas prestaciones en la estimación de h[n]. En el siguiente experimento se analiza un tipo de señales periódicas limitadas en banda que incluye reflexiones como parte del ciclo de señal fundamental.

162 34 CAPÍTULO 4. DECONVOLUCIÓN HOMOMÓRFICA Amplitud t[segundos] (a) Amplitud t[segundos] (b) Figura 4.3: Señales originales y resultados de la DH en DS: (a) modelo parabólico de señal (trazo continuo) y su estimación (trazo discontinuo); (b) señal dispersa original (círculos) y su estimación (trazo discontinuo) Señales periódicas limitadas en banda con reflexiones En el presente experimento se utiliza el modelo de señal de parábola invertida (4.57) de la sección anterior. Por otra parte, la señal dispersa está formada por un conjunto de impulsos que representan reflexiones. Las reflexiones pueden generar las dos situaciones siguientes: las reflexiones se dan dentro de la longitud del ciclo formando una señal discreta de fase mínima; las reflexiones se dan dentro de la longitud del ciclo formando una señal discreta de fase no mínima. Se descarta la situación donde las reflexiones se extienden más allá de la longitud del ciclo, debido a que las motivaciones prácticas a estudiar en el capítulo siguiente consideran las reflexiones dentro de la longitud del ciclo básico. Estas reflexiones no modifican la periodicidad de la señal de salida sino que solamente modifican la forma de ésta con respecto a la entrada. Las Figuras 4.4(a) y (4.4)(b) muestran el modelo señal basado en una parábola invertida y una secuencia de reflectores periódicos, respectivamente. En la Figura 4.4(c) se observa cómo el efecto de los reflectores se traduce en una distorsión de la señal con modelo parabólico. También, de la Figura 4.4(e) se observa cómo el ángulo de fase presenta muchas discontinuidades, dificultando el cálculo directo de su CC. En el caso más simple, donde las reflexiones tienen una separación uniforme y un valor de amplitud que decae exponencialmente, se tiene que las prestaciones de los algoritmos de DH, basados en el MD y el DS, son bastante similares a las obtenidas en el experimento que solo considera periodicidad. La Tabla 4.4 identifica los experimentos, variando el número de ciclos, para

163 4.5. SIMULACIONES Y EXPERIMENTOS Amplitud t[segundos] (a) Amplitud t[segundos] (b) Amplitud t[segundos] (c) H(e jω ) Frecuencia [rad] (d) H(e jω ) Frecuencia [rad] (e) Figura 4.4: Modelo de parábola invertida: (a) representación de un ciclo; (b) tren de seis impulsos con reflexiones; (c) señal formada por seis ciclos; (d) amplitud del espectro; (e) fase del espectro. los cuales es posible deconvolucionar la señal dispersa periódica del modelo de señal (4.57). Las Figuras 4.5(a) y 4.5(b) muestran la estimación del ciclo básico de señal y la señal dispersa periódica para 6 ciclos de señal, respectivamente. Ninguno de los métodos, MD o MR, permiten deconvolucionar dichas señales.

164 36 CAPÍTULO 4. DECONVOLUCIÓN HOMOMÓRFICA Número de ciclos Método MD Sí Sí Sí Sí No Sí No Sí Sí Sí No Sí Sí DS Sí Sí Sí Sí Sí Sí No Sí Sí Sí Sí Sí Sí Tabla 4.4: Determinación de la señal dispersa con reflexiones mediante MD y DS a partir de varios ciclos del modelo de señal de parábola invertida. Amplitud t[segundos] (a) Amplitud t[segundos] (b) Figura 4.5: Señales originales y resultados de la DH en DS : (a) modelo parabólico de señal (trazo continuo) y su estimación (trazo discontinuo); (b) señal dispersa original (círculos) y su estimación (trazo discontinuo). En la situación donde las reflexiones presentan una separación no uniforme y el valor de su amplitud no decae exponencialmente, el CC corresponde, en general, a una señal de fase no mínima. Esta situación complica la aplicación de la DH, ya que la componente de fase máxima de la secuencia de impulsos se solapa con el CC de la señal de variación suave. En el Apéndice C se desarrolla la expresión general para el CC de una secuencia de impulsos. De este desarrollo teórico se observa que el problema se resuelve multiplicando la señal por una exponencial decreciente [7, 37]. La multiplicación por una exponencial decreciente se interpreta geométricamente como el acercamiento de las raíces de la señal al origen [87]. La selección adecuada del valor constante de la exponencial decreciente permite tener valores de raíces de la señal con módulos menores que uno, convirtiendo la señal en una señal de fase mínima Conclusiones En este capítulo se han descrito los tres principales problemas que presentan los algoritmos de DH, a saber: el desenrrollado de la fase, el cálculo del logaritmo complejo de valores cercanos o iguales a cero y la estimación cuando la fase es no mínima. También, condicionado por un tipo particular

165 4.6. CONCLUSIONES 37 de aplicación, el análisis de soluciones se ha centrado en el cálculo del logaritmo complejo de señales cuya TZ tiene valores cercanos o iguales a cero. La aplicación en cuestión es similar al problema de identificación de sistemas con señales de entrada y salida periódicas y de banda limitada. Utilizando el algoritmo de la DH basado en la construcción de una señal suplemental fue posible superar el problema de la estimación de la respuesta al impulso de este tipo particular de SLIT. También, se escogió un conjunto de pruebas con una secuencia incremental de dificultad. En el primer experimento fue posible obtener buenas prestaciones a partir de la aplicación directa de la definición del CC. En el segundo experimento, se incrementó la dificultad introduciendo como señal de variación suave una señal de banda limitada. En el tercer experimento se mostraron las buenas prestaciones del algoritmo de la DH para señalas, que además de ser limitadas en banda, tienen un comportamiento periódico. Finalmente, en el cuarto experimento se incrementó la dificultad añadiendo reflexiones a la señal periódica de banda limitada. En todos estos experimentos el algoritmo de DH basado en la transformación de banda completa demostró buenas prestaciones.

166 38 CAPÍTULO 4. DECONVOLUCIÓN HOMOMÓRFICA

167 Capítulo 5 Caracterización de la Impedancia de Entrada Aórtica mediante DH Dos enfoques históricos han sido utilizados para modelar la relación existente entre presión y flujo y las propiedades mecánicas del sistema arterial. Por una parte, el modelo Windkessel describe el sistema arterial como un sistema hidráulico donde la presión y el flujo están relacionados con la compliancia arterial y la resistencia periférica total [95]. Por otra parte, el otro enfoque describe al sistema arterial como un tubo de longitud infinita, donde la impedancia de entrada está asociada a la compliancia arterial local y al área de sección transversal [4]. En la práctica, ambos métodos resultan inadecuados en la caracterización del sistema arterial [97]. Las deficiencias más grades de ambos métodos históricos provienen de la no incorporación de la topología espacial del sistema arterial [95]. Ambos modelos históricos fallan en explicar la diferencia entre las formas de ondas de presión y flujo [95]. Algunas investigaciones teóricas han sugerido que esta diferencia se debe a las reflexiones que las ondas sufren a lo largo de su propagación a través del sistema circulatorio [45, 49, 5]. De hecho, a lo largo de los años, ha habido mucho debate sobre si las ondas originadas en el corazón podrían ser reflejadas y, de serlo, en qué sitios del sistema arterial [45]. La principal prueba presentada sobre la existencia de reflexiones ha sido precisamente la diferencia en las formas de onda entre la presión y el flujo, un sistema arterial libre de reflexiones presentaría formas de onda similares. En la década de 96 se introdujo y extendió el concepto de impedancia de entrada aórtica (IEA), la cual describe la relación en régimen permanente entre la presión en el inicio de la aorta ascendente y el flujo de salida valvular simultáneo. El estudio de la IEA muestra que ninguno de los dos modelos históricos es capaz de describir la impedancia de entrada correctamente, y por tanto modelar el sistema arterial [95, 97]. A pesar de que en [88] se establece que la IEA permite identificar la magnitud de las reflexiones y el sitio donde 39

168 4 CAPÍTULO 5. CARACTERIZACIÓN DE LA IMPEDANCIA DE ENTRADA AÓRTICA MEDIANTE DH éstas se producen, apenas se han realizado trabajos que permitan identificar a partir de la IEA la magnitud y posición de los reflectores. En este capítulo se utilizan técnicas de DH para determinar, en caso de que existan, las reflexiones a través del estudio de la IEA. Se aprovecha la ventaja intrínseca en las técnicas homomórficas que convierten cocientes entre espectros de señales en restas de espectros de señales, siendo esto un factor clave en la estimación de la IEA. El capítulo está organizado como a continuación. La primera sección introduce conceptos básicos sobre la fisiología del sistema cardiovascular. La segunda sección estudia las dos escuelas tradicionales sobre propagación de ondas en el sistema arterial. La tercera y última sección es la aportación de esta Tesis, es decir, se aplican algoritmos de DH en la estimación de la IEA, fundamentalmente se estudian las prestaciones del algoritmo basado en la transformación al DS. También, se incluye la aplicación de los algoritmos de DH en la estimación de la admitancia de entrada aórtica (AEA), que ha demostrado ser un problema mejor condicionado. 5.. El Sistema Cardiovascular El sistema cardiovascular está formado por el corazón, los vasos sanguíneos y la sangre. En términos simples, algunas de sus principales funciones son: () la distribución de oxígeno y nutrientes a todos los tejidos corporales; (2) el transporte del dióxido de carbono y los productos metabólicos de desecho desde los tejidos a los pulmones y a los órganos excretores; (3) la distribución del agua, electrólitos y hormonas por todo el organismo. En el sistema cardiovascular, la sangre es impulsada a lo largo de todo el sistema por el corazón, que funciona como una doble bomba muscular. Cada lado del corazón hace las funciones de una bomba que contiene dos cavidades, una aurícula y un ventrículo. Las aurículas contribuyen al llenado de los ventrículos. Los ventrículos al contraerse expulsan la sangre hacia todo el organismo. Por una parte, la contracción del ventrículo derecho impulsa la sangre a través de la válvula pulmonar hacia la arteria pulmonar, que se ramifica progresivamente formando la circulación pulmonar, llegando a los pulmones. Dentro de los pulmones se da el intercambio de CO 2 por O 2, y la sangre oxigenada vuelve al corazón a través de la aurícula izquierda cerrando el circuito de la circulación pulmonar. Por otro parte, la contracción del ventrículo izquierdo empuja la sangre a través de la válvula aórtica hacia la aorta, primera y mayor arteria de la circulación sistémica. La sangre fluye desde la aorta a las arterias principales, cada una de las cuales irriga un órgano o una región corporal. La sangre vuelve al corazón a través del sistema venoso, cerrando la circulación sistémica.

169 5.. EL SISTEMA CARDIOVASCULAR Ciclo cardíaco El ciclo cardíaco es la secuencia de hechos eléctricos y mecánicos que se producen durante un único latido cardíaco. Se compone de dos fases: la diastólica o de relajación y la sistólica o de contracción. En la Figura 5. se observa la representación gráfica de las curvas de presión y flujo en el ventrículo izquierdo y en la entrada de la aorta. Sobre estas curvas se ha representado la señal del electrocardiograma, que permite identificar las fases que se dan durante cada ciclo. A continuación se explicara el ciclo cardíaco a través de cada una de los eventos que ocurren durante los periodos sistólicos y diastólicos. Sístole auricular. La despolarización auricular producida por la despolarización del nodo sinusal (SA) conduce a la contracción auricular. Dicha contracción completa el llenado ventricular (las aurículas contribuyen en menos del 2 % al volumen ventricular final). El volumen ventricular después del llenado se conoce como volumen del fin de diástole o volumen telediastólico (VTD). Sístole ventricular. La contracción ventricular produce un aumento brusco de la presión ventricular y las válvulas auriculoventriculares (AV) se cierran una vez que dicha presión supera la auricular. Durante la fase inicial de la contracción ventricular, la presión es menor que la de la arteria pulmonar y de la aorta, con lo que las válvulas de salida permanecen cerradas. En este tiempo de generación de presiones no hay entrada ni salida de sangre del ventrículo, y se conoce como contracción isométrica o isovolumétrica, dado que el volumen ventricular no cambia. Eyección. Las válvulas de salida se abren cuando la presión en el ventrículo excede la de la arteria. El flujo dentro de las arterias es inicialmente muy rápido pero a medida que la contracción se reduce, la eyección se reduce. La contracción activa cesa durante la segunda mitad de la eyección, y el músculo se repolariza. La presión ventricular durante la fase de eyección reducida es ligeramente menor que la de la arteria, pero la sangre sigue saliendo a causa del momento. Finalmente el flujo se invierte brevemente, provocando el cierre de la válvula de salida y un pequeño incremento en la presión aórtica, la muesca dicrótica. Diástole: relajación y llenado. En el periodo siguiente al cierre de las válvulas de salida, los ventrículos se relajan rápidamente. Sin embargo, la presión ventricular es aún mayor que la presión auricular, y las válvulas AV permanecen cerradas. Esto se conoce como relajación isovolumétrica. Durante las últimas dos terceras partes de la sístole, la presión auricular se

170 42 CAPÍTULO 5. CARACTERIZACIÓN DE LA IMPEDANCIA DE ENTRADA AÓRTICA MEDIANTE DH R Milivoltios 2 P Q Válvula Semilunar Abierta S T Válvula Semilunar Cerrada Presión (mmhg) Válvula AV Cerrada Válvula AV Abierta 6 Volumén (ml) 2 8 Periodo de contracción isométrica Periodo de eyección Periodo de relajación isométrica Periodo de contracción isométrica Periodo de eyección Periodo de relajación isométrica Sonido (Hz) Primer Sonido Segundo Sonido Tercer Sonido Sístole Diástole Sístole Figura 5.: Curvas de presión y flujo que ilustran la secuencia de hechos mecánicos que se producen durante un único ciclo cardíaco (figura adaptada de []). incrementa como resultado del llenado desde las venas. Cuando la presión ventricular desciende por debajo de la presión auricular, las válvulas AV se abren, y la presión auricular disminuye a medida que los ventrículos se llenan rápidamente. A este hecho contribuye el reflujo elástico de las paredes ventriculares, esencialmente aspirando la sangre. Cuando los ventrículos se relajan

171 5.. EL SISTEMA CARDIOVASCULAR 43 completamente, el llenado se ralentiza. Esto continúa sobre las dos terceras partes de la diástole a causa del flujo venoso. En reposo la diástole tiene una longitud dos veces mayor que la sístole, pero decrece proporcionalmente durante el ejercicio y a medida que la frecuencia cardíaca se incrementa. Pulso arterial. Hay muchas cuestiones que todavía no están claras sobre la génesis del pulso arterial [88]. Como definición general, el pulso arterial es cualquier fluctuación periódica originada por el corazón con una frecuencia igual al ritmo cardíaco. La eyección de la sangre en cada latido cardíaco se convierte en pulsaciones de flujo y presión. Las pulsaciones de flujo son los movimientos longitudinales de sangre a lo largo del sistema arterial. El flujo puede medirse a través de medidores de flujo colocados alrededor de la arteria o insertados a través de cateterismo; también, puede medirse a través de técnicas no invasivas como las técnicas basadas en ultrasonidos. La forma de onda del flujo en la aorta es la de un pulso triangular durante el periodo sitólico que se reduce a un valor constante y casi nulo durante el periodo diastólico, véase Figura 5.2(b). El flujo alcanza valores máximos en la entrada de la aorta (hasta 3 % el valor medio de flujo arterial) y decrece progresivamente hacia la periferia. Esta atenuación se debe a: () las características mecánicas de las arterias; y (2) el incremento en el área total efectiva de la sección transversal de todas las arterias de la periferia con relación al área de la aorta. Las pulsaciones de presión tienen su origen en la sangre que sale del ventrículo izquierdo a través de la aorta, que genera una onda de presión. La fluctuación pulsátil se incrementa desde la aorta hacia la periferia, véase Figura 5.2(b), alcanzando valores de 4 % a 8 % mayores que la presión arterial media. El estudio del pulso arterial tiene su historia ligada al estudio de la medicina y al esfuerzo por entender el funcionamiento del sistema cardiovascular [88]. No fue sino hasta el siglo XVII que se empezó a tener una concepción del funcionamiento del sistema cardiovascular, y hubo que esperar al desarrollo de teorías en dinámica de fluidos y propagación de ondas para tener una mejor concepción del sistema cardiovascular Hemodinámica La hemodinámica es el estudio de las relaciones entre la presión (P), la resistencia (R) y el flujo sanguíneo (Q) en el sistema cardiovascular. Aunque las propiedades de este flujo son enormemente complejas, pueden deducirse en gran manera de una visión más simple de las leyes físicas que gobiernan el flujo de los líquidos a través de tuberías. Así, el flujo sanguíneo a través Aunque el término pulso se refiere a cualquier tipo de pulsación, normalmente los clínicos lo entienden como pulso de presión en una arteria grande y accesible.

172 44 CAPÍTULO 5. CARACTERIZACIÓN DE LA IMPEDANCIA DE ENTRADA AÓRTICA MEDIANTE DH La aorta se distiende durante la sístole Presión (mmhg) 8 6 Aorta ascendente Aorta torácica Aorta abdominal Aorta Femoral Función Windkessel La aorta recobra su forma durante la diástole para propulsar la sangre hacia adelante. 4 Volumen (cm/seg) 6 La onda de presión refleja desde puntos de bifurcación. 2-2 (a) (b) Figura 5.2: Ondas de presión y velocidad: (a) la onda de presión es causa del ensanchamiento de la pared arterial; (b) ondas de presión y velocidad en diferentes puntos de la aorta (Figuras adaptadas de []). de los vasos está determinado por la diferencia de presión en el vaso y por la resistencia de la sangre, matemáticamente se escribe como sigue: Q = P 2 P R = P R (5.) donde P y P 2 corresponden a los valores de presión en los puntos y 2, respectivamente, y P corresponde al diferencial de presión entre ambos puntos. La resistencia al flujo está originada por las fuerzas de fricción dentro del líquido, y depende de la viscosidad de éste y de las dimensiones de la arteria descritas por la ley de Poiseuille: R = 8νL πr 4 (5.2) donde ν es la viscosidad de la sangre, L es la longitud del vaso, r radio del vaso. De (5.2) se observa cómo pequeños cambios en el radio producen grandes variaciones en el valor de flujo. Considerando el sistema cardiovascular en conjunto, los distintos tipos de vasos sanguíneos (arterias, arteriolas, capilares) se hallan dispuestos en serie, con lo cual la resistencia de todo el sistema es igual a la suma de todas las resistencias ejercidas por cada vaso, siendo las arteriolas las que ofrece mayor resistencia [].

Sitemap