Repositorio Institucional de Documentos

Trabajo Fin de M´aster Ingenier´ıa Biom´edica Curso 2010-2011 An´alisis de Estimulaci´on El´ectrica Funcional en Se˜nales de EEG Ana Isabel Cerrada Baena Septiembre de 2011 Director: Javier M´ınguez Zafra Departamento de Inform´atica e Ingenier´ıa de Sistemas Centro Polit´ecnico Superior Universidad de Zaragoza A mis padres, por la dedicaci´on, generosidad y apoyo incondicional que siempre me han prestado, por los consejos pronunciados con sabidur´ıa, por transmitirnos fortaleza y esperanza. Mi camino siempre ir´a junto al vuestro, incluso cuando parezca alejarse, os reencontrar´a. Gracias por haber hecho de mi la mujer que ahora soy. Por supuesto, esto tambi´en es para mis hermanos y mi hermana, que por muchas diferencias que hayamos tenido, han sabido escucharme y apoyarme en los momentos dif´ıciles, as´ı como hacerme sonre´ır y conseguir que los momentos buenos sean los inolvidables. Tambi´en me gustar´ıa dedicar este trabajo al resto de mi familia, en especial, a mis abuelos, por los jazmines en flor y las tizas de colores en los escalones, por la m´usica, las excursiones con caracoles y setas, las historias, los despertares... Por estar siempre con nosotros. Unas palabras muy especiales van dirigidas a t´ı Michael, por tu apoyo siempre incondicional, por escucharme, entenderme y aconsejarme, por ense˜narme que hay otra forma de percibir lo que nos rodea y de afrontar los problemas. Tus ´animos de las ´ultimas semanas y la ilusi´on de empezar nuestra nueva etapa han contribuido mucho en la finalizaci´on de este trabajo. Por ´ultimo, a los que me hab´eis acompa˜nado este a˜no en Zaragoza, haciendo el a˜no un poco... “m´a fasi”. ii Resumen En los sistemas gobernados por interfaces cerebro-ordenador (Brain Computer Interfaces o BCI) se requiere un alto nivel de procesamiento de se˜nal para extraer la informaci´on que resulta de inter´es, por lo que es muy importante el filtrado de artefactos, es decir, de todas aquellas se˜nales que se registran junto con las se˜nales electroencefalogr´aficas (EEG o electroencefalograma), pero que no resultan de inter´es para el investigador. En este trabajo se estudia un artefacto muy concreto, la estimulaci´on el´ectrica funcional (Functional Electrical Stimulation o FES), ya que actualmente la Universidad de Zaragoza coordina un proyecto de investigaci´on CYCIT en el que se intenta controlar una ´ortesis robotizada con una BCI haciendo uso de la electroestimulaci´on. El objetivo de este proyecto es, por tanto, conocer c´omo afecta una se˜nal de FES a las se˜nales cerebrales y conseguir eliminarla de la se˜nal EEG. Este proyecto consta fundamentalmente de cuatro fases adem´as de la parte pr´actica. La primera de ellas consisti´o en la b´usqueda de bibliograf´ıa relacionada con las se˜nales de EEG ([1]) y artefactos ([2]), fundamentos de las BCI y t´ecnicas de separaci´on ciega de fuentes, como son los an´alisis de componentes principales (Principal Components Analysis oPCA) y de componentes independientes (Independent Component Analysis oICA). Posteriormente se procedi´o a un estudio frecuencial de las se˜nales de EEG contaminadas con FES que fueron proporcionadas por la empresa Fatronik. Posteriormente se aplicaron a esas se˜nales las dos t´ecnicas de separaci´on de fuentes estudiadas para tratar de extraer la se˜nal de artefacto FES del EEG. Ni PCA ni la t´ecnica estudiada inicialmente de ICA consiguieron resultados satisfactorios, por lo que se volvi´o a la fase de documentaci´on. En la tercera parte del proyecto se aplicaron nuevos algoritmos de ICA, recurriendo finalmente a la utilizaci´on de tranformadas wavelets y algoritmos de ICA mejorados con esta transformada para tratar de conseguir el objetivo. Todas las t´ecnicas y algoritmos estudiados se han aplicado a tres conjuntos de datos adquiridos bajo condiciones diferentes en distintos individuos, obteniendo la misma falta de resultados positivos en todos ellos. Por ´ultimo, y ante la imposibilidad de encontrar un algoritmo que extrayese el artefacto FES de forma eficiente, no se aplicaron de forma pr´actica los algoritmos utilizados. La parte pr´actica de este proyecto se ha realizado en la spin-off Bit&Brain Technologies, donde se ha contribuido de forma activa en experimentos de BCI. Se ha participado en la preparaci´on del equipo y material necesario para proceder a la adquisici´on de se˜nales de EEG, procediendo a la colocaci´on de los electrodos y visualizaci´on de las se˜nales de EEG en un ordeiv nador personal. As´ı mismo, he participado como sujeto en un experimento de BCI siguiendo un protolo cuya finalidad es extraer caracter´ısticas de la se˜nal EEG que permita el reconocimiento de las emociones de un individuo. v Abstract Electroencephalographic (EEG) recordings usually contain a lot of artifacts besides the signal representing the neural activity of the brain. Especially when using the EEG signal as an input for brain-computer interfaces (BCI) the removal of those artifacts is crucial. This project had a focus on finding methods in order to remove artifacts generated by functional electrical stimulation (FES) from EEG recordings. Information gathered during this project will serve for a research project with the objective of controlling a robotic prothesis by a BCI using FES at the University of Zaragoza. This project basically consisted of four parts. The first step was a literature research regarding the following topics: EEG ([1]) and artifacts ([2]), BCI basics and Blind Source Separation (BSS) such as Principal Component Analysis (PCA) and Independent Component Analysis (ICA).In the second part frequency analysis of several EEG recordings containing FES artifacts were performed. Several BSS methods were applied in order to extract and remove the FES artifacts from the EEG signals. Neither PCA nor ICA achieved satisfying results, so new methods had to be searched. In the third part of the project new ICA algorithms as well as their combination with wavelet transformation ([3]) were investigated. Every method was applied for three different EEG datasets without obtaining any satisfying results. Due to the fact that no satisfying algorithm was found to remove the FES artifact from the EEG recording the practical part of the project consisted of doing experiments with brain-computer interfaces at Bit&Brain Technologies, a spin-off company of the University of Zaragoza. In order to do BCI experiments the electrodes had to be placed on the scalp of the subject as well as the preparation of the EEG acquisition equipment had to be done for monitoring the EEG signal in real time. Further I took part as one of the subjects in an experiment concerning emotion recognition. New ICA algorithms and wavelet transforms, as well as ICA algorithms enhanced with wavelet transforms ([3]) were applied at this part of the project in order to achieve the main goal. Every technique was applied over three different EEG datasets and any successful result was obtained. vi ´ Indice general 1. Introducci´on 1 1.1. Motivaci´on ....................................... 1 1.2. Planteamiento del Problema . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.3. AlcancedelTrabajo .................................. 4 2. La Se˜nal de Estimulaci´on El´ectrica Funcional 5 2.1. RegistrosdeEEG.................................... 5 2.2. An´alisis en el Dominio de la Frecuencia . . . . . . . . . . . . . . . . . . . . . . . 7 2.2.1. M´etodos..................................... 7 2.2.2. Resultados ................................... 8 3. Pre-procesado de la Se˜nal EEG 10 3.1. TransformadaWavelet................................. 10 3.1.1. M´etodos..................................... 10 3.1.2. Resultados ................................... 11 3.2. An´alisis de Componentes Principales . . . . . . . . . . . . . . . . . . . . . . . . . 13 3.2.1. M´etodos..................................... 13 3.2.2. Resultados ................................... 14 4. Procesado de Se˜nales EEG con An´alisis de Componentes Independientes 16 4.1. M´etodos......................................... 16 4.1.1. FastIcaeIcasso................................. 17 4.1.2. ICA para Distribuciones Supergaussianas . . . . . . . . . . . . . . . . . . 18 vii 4.1.3. ICA para Se˜nales con Estructura Temporal . . . . . . . . . . . . . . . . . 18 4.1.4. WaveletICA .................................. 19 4.2. Resultados........................................ 20 4.2.1. Aplicaci´on de Fastica e Icasso . . . . . . . . . . . . . . . . . . . . . . . . . 20 4.2.2. Aplicaci´on de runica .............................. 23 4.2.3. Aplicaci´on de amuse .............................. 25 4.2.4. Aplicaci´on de wavelet-ICA . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 5. Conclusiones y Trabajo Futuro 26 Bibliograf´ıa 27 viii 1. Introducci´on 1.2 Planteamiento del Problema Figura 1.2: A la izquierda se muestra el espectro en frecuencia de la se˜nal x(t). A la derecha, dicho espectro en bajas frecuencias, donde se observan los arm´onicos de la se˜nal La Figura 1.2 muestra el espectro de frecuencias de la se˜nal x(t), con m´as detalle en la gr´afica derecha, donde se observa que existen componentes a baja frecuencia (en 6 Hz) a pesar de que la se˜nal senoidal tiene una frecuencia de 100 Hz y el tren de pulsos de 35 Hz. Esto demuestra la mezcla de frecuencias que se produce en bajas frecuencias. En este proyecto se analiza la se˜nal de FES, se˜nal peri´odica a una frecuencia de 35 Hz y que utiliza un pulso bif´asico de muy corta duraci´on. Esta situaci´on es, por tanto, an´aloga a la situaci´on descrita. En los registros de EEG puede apreciarse que las se˜nales de EEG y FES se encuentran mezcladas, como muestra la Figura 1.3. Se observan los 16 canales de EEG registrados en un trial, donde son f´acilmente reconocibles a simple vista los picos de FES entre t=3yt=6segundos. La escala de amplitud est´a fijada en esta gr´afica a 50µV , y puede comprobarse tambi´en mediante inspecci´on visual que la frecuencia del FES es 35 Hz. Adem´as puede decirse que la transmisi´on del FES resulta no lineal, puesto que las amplitudes de los picos no son constantes. Figura 1.3: A la izquierda se muestran los 16 canales de un trial, donde se aprecia el artefacto de FES. A la derecha se observa un detalle del canal C1, uno de los situados m´as cerca de la corteza motora. 3 1. Introducci´on 1.3 Alcance del Trabajo En el an´alisis frecuencial de la se˜nal EEG se distinguen distintas bandas de frecuencias asociadas a las distintas actividades neuronales, y todas ellas se ven afectadas por la convoluci´on con el espectro del FES. Las principales son las bandas Delta (1-4 Hz), Theta (4-8 Hz), Alpha (8-12 Hz) y Beta (12-30 Hz) ([1]). Los ritmos Mu aparecen en torno a 10 Hz y est´a relacionado con el ´area motora, disminuyendo su potencia cuando existe movimiento en las extemidades. En las bajas frecuencias de la banda Beta, entre (12-15) Hz, tambi´en se registra actividad de la corteza motora cerebral ([5]). Cuando las ´areas motoras cerebrales se activan, la potencia de esta banda disminuye respecto al nivel que tiene cuando no hay movimiento. Para que la informaci´on neurol´ogica de inter´es de la se˜nal EEG que reside en estas bandas de frecuencias no se vea alterada, es necesario eliminar los artefactos de FES. Esta tarea no resulta eficaz empleando ´unicamente t´ecnicas de filtrado convencionales, sino que es necesario aplicar t´ecnicas que nos permitan separar las se˜nales de distinta naturaleza. Es por ello que se ha orientado la realizaci´on del trabajo a la conocimiento y utilizaci´on de t´ecnicas de separaci´on de fuentes como PCA e ICA. 1.3. Alcance del Trabajo El presente Trabajo Fin de M´aster requiere profundos conocimientos acerca de se˜nales de electroencefalograf´ıa, de interfaces cerebro-ordenador, as´ı como de estad´ıstica; por lo que una parte muy importante en el desarollo de este proyecto ha sido la documentaci´on. Los pasos seguidos en la consecuci´on del objetivo de este proyecto han sido los siguientes: B´usqueda y estudio de documentaci´on que me permitiera obtener el conocimiento suficiente sobre se˜nales EEG ([1]), interfaces cerebro-ordenador y t´ecnicas de separaci´on ciega de fuentes. Aplicaci´on de las t´ecnicas estudiadas a diferentes conjuntos de datos de EEG con FES para lograr la extracci´on del artefacto. Estudio y aplicaci´on de t´ecnicas diferentes a las ya utilizadas para tratar de conseguir separar ambas se˜nales. Participaci´on en experimentos de BCI realizados en la empresa Bit&Brain Technologies, con la finalidad de completar los conocimientos te´oricos aprendidos con los conceptos que requiere el uso pr´actico de esta tecnolog´ıa La propuesta del Trabajo Fin de M´aster inclu´ıa un ´ultimo apartado en el que se estudia la viabilidad de la aplicaci´on de este proyecto a nivel comercial. Aunque no se hayan obtenido resultados satisfactorios en la extracci´on de la se˜nal de FES, se han sembrado las bases para la consecuci´on de este objetivo. 4 2. La Se˜nal de Estimulaci´on El´ectrica Funcional En este cap´ıtulo se presentan los diferentes registros de EEG con artefactos FES con los que se cuenta. Se extraer´an las caracter´ısticas de las se˜nales EEG en el dominio temporal, procedi´endose posteriormente al an´alisis frecuencial de las se˜nales. 2.1. Registros de EEG La adquisici´on de se˜nales de EEG se realiza mediante un gorro en el cual se encuentran situados los electrodos. Dependiendo de la actividad cerebral que quiera registrarse, los sensores deber´an estar situados en distintas localizaciones sobre el cuero cabelludo. De las bandas frecuenciales citadas anteriormente, resultar´an de inter´es aquellas que reflejan la actividad de la corteza motora, es decir, los ritmos Mu y las bajas frecuencias de la banda Beta. Las se˜nales de EEG con artefactos FES utilizadas fueron adquiridas por Fatronik, empresa participante en el proyecto CYCIT. Todos los datos proporcionados se han adquirido con un amplificador comercial g.USBamp de g.Tec, al que se hab´ıan conectado 16 electrodos situados sobre el cuero cabelludo en el ´area de la corteza motora, ya que es ah´ı donde deben aparecer los efectos neuronales de la eletroestimulaci´on motora. Los electrodos fueron colocados en: FC3, FC1, FCz, FC2, FC4, C3, C1, Cz, C2, C4, CP3, CP1, CPz, CP2, CP4 and Pz, y su disposici´on se observa en la Figura 2.1. Los electrodos de tierra y referencia estaban en los mastoides izquierdo y derecho. 5 2. La Se˜nal de Estimulaci´on El´ectrica Funcional 2.1 Registros de EEG Figura 2.1: Localizaci´on de los electrodos sobre el cuero cabelludo. En la realizaci´on de este trabajo se ha contado con dos conjuntos de datos que fueron adquiridos en experimentos configurados con distintos par´ametros. En ambos experimentos se situ´o al sujeto frente a una pantalla que le indicaba lo que deb´ıa hacer en cada momento, utilizando protocolos de actuaci´on similares y teniendo ambos la misma estructura de trials que se muestra en la Figura 2.2. En ella se observa que tanto al principio como al final del trial se mostraba la pantalla completamente negra para que el sujeto tratase de no pensar en nada, mientras que entre los segundos 3 y 6 de cada trial se le indic´o al sujeto la tarea a realizar, que tambi´en vari´o dependiendo del experimento. Figura 2.2: Estructura temporal de un trial correspondiente al registro de EEG En el primero de ellos, las tareas a realizar fueron imaginar movimiento en el antebrazo izquierdo, en el derecho o la aplicaci´on de FES mediante unos electrodos situados en la mu˜neca derecha. Estas tres tareas se realizaron el mismo n´umero de veces en una secuencia totalmente aleatorizada hasta registrar 120 trials. La se˜nal fue adquirida con una frecuencia de muestreo de 19.2 KHz, siendo la frecuencia del FES de 25 Hz. En el segundo experimento, se registraron 40 trials donde la ´unica tarea existente entre los segundos 3 y 6 fue la aplicaci´on de FES. En este caso, la se˜nal se muestre´o a una frecuencia de 2.4 KHz y el FES se aplic´o con una frecuencia de 35 Hz. En ninguno de los dos experimentos se filtr´o durante la adquisici´on el ruido de la red el´ectrica, por lo que se ha eliminado mediante software con un filtro notch a 50 Hz. 6 2. La Se˜nal de Estimulaci´on El´ectrica Funcional 2.2 An´alisis en el Dominio de la Frecuencia El primer experimento proporcion´o dos conjuntos de registros de EEG con artefactos, ya que se realiz´o el experimento con amplitudes de FES peque˜nas y posteriormente con amplitudes de artefacto muy grandes. Aunque todos los an´alisis presentes en este proyecto se han realizado con los tres conjuntos de datos de EEG proporcionados, se presentan s´olo las gr´aficas y conclusiones alcanzadas para las se˜nales EEG registradas en el experimento 2, ya que con todos los datos se obtuvieron resultados similares. Normalmente las se˜nales de EEG se muestrean a una frecuencia de 256 ´o 512 Hz, por lo que en esta ocasi´on contamos con una frecuencia de muestreo muy elevada, lo que provoca que en esta se˜nal haya bastante ruido de alta frecuencia. Este ruido resulta ser menor que en el caso de las se˜nales muestreadas a 19.2 KHz. 2.2. An´alisis en el Dominio de la Frecuencia En esta secci´on se presentan las distintas t´ecnicas empleadas para realizar el an´alisis espectral de las se˜nales de EEG. Este an´alisis nos permitir´a alcanzar un conocimiento m´as profundo de la se˜nal de FES, puesto que el comportamiento frecuencial de las se˜nales EEG es ya conocido. Dada la periodicidad de la se˜nal de FES, y seg´un lo visto en el ejemplo de la introducci´on, su espectro en frecuencia no mostrar´a ´unicamente una componente en 35 Hz y el resto de arm´onicos en todas las frecuencias m´ultiplos de la frecuencia del FES. Dado que el FES estimular´a las fibras motoras, y mediante ´esta la corteza motora cerebral, deber´ıa producirse un descenso de la potencia en torno a 10 Hz (ritmo Mu) y en la banda baja de Beta. Resulta entonces de importancia conocer c´omo se distibuye la potencia de la se˜nal EEG con y sin FES en las distintas bandas de frecuencias, para determinar cu´ales son los efectos de este artefacto. 2.2.1. M´etodos Para calcular la densidad espectral de potencia o PSD (Power Spectral Density) de la se˜nal EEG se utilizar´a el M´etodo de Welch, que consiste en una mejora del M´etodo de Bartlett. Ambos est´an basados en la estimaci´on de la PSD de la se˜nal, ˆ Sx(ejω), mediante el periodograma, definido como la transformada de Fourier de la estimaci´on de la funci´on autocorrelaci´on de la se˜nal (Ecuaci´on 2.1). Calculando la PSD con el periodograma se obtiene un estimador sesgado, aunque el sesgo disminuye al incrementar el n´umero de muestras que tenga la se˜nal. Sin embargo, resulta un estimador inconsistente ya que su valor no tiende al valor real conforme aumenta el n´umero de muestras, es decir, su varianza no disminuye cuando lo hace el sesgo [6]. La resoluci´on frecuencial que se obtiene con el periodograma est´a limitada por el n´umero de muestras ˆ Sx(ejω) = N−1 X k=−N+1 ˆrx(k)e−jωk (2.1) El m´etodo de Bartlett consiste en un promediado de periodogramas enventanados, es decir, se divide la se˜nal original de N muestras en L segmentos de M muestras, y se realiza el periodograma 7 2. La Se˜nal de Estimulaci´on El´ectrica Funcional 2.2 An´alisis en el Dominio de la Frecuencia a cada segmento de M muestras para promediarlos todos despu´es. De esta forma se consigue disminuir la varianza a cambio de una disminuci´on de la resoluci´on frecuencial. En el m´etodo de Welch, la se˜nal se divide en segmentos no consecutivos, es decir, existe solapamiento entre las muestras de dos segmentos. El m´etodo de Welch se ha aplicado a diferentes secciones de la se˜nal EEG dependiendo de si ´esta est´a contaminada o no por el artefacto de FES. De esta forma se obtendr´an representaciones frecuenciales en ambas situaciones y podr´an compararse los resultados. Es importante decir que el c´alculo de la PSD se realiza para cada uno de los canales que regitran el EEG, es decir, se calcular´an 16 PDSs. Dado que los segmentos sin FES de cada trial van de t= 0 a t= 3segundos y de t= 6 a t= 10segundos, se han concatenado los segmentos de t= 0 a t= 3segundos de varios trials hasta formar otro de 9 segundos de duraci´on. De igual forma, se han formado segmentos de 9 segundos de EEG con artefacto FES, realizando esta tarea para cada uno de los canales. A cada uno de esos segmentos de 9 segundos nos referiremos a partir de ahora como segmentos artificiales La ventana empleada ha sido una Hanning con una longitud de un segundo, con un 50 % de muestras solapadas en los segmentos. Adem´as hab´ıa que seleccionar los puntos con lo que se realizar´a la FFT, siendo ´este un par´ametro importante para obtener la resoluci´on frecuencial (4f) adecuada. El n´umero de puntos con los que realizar la FFT debe ser al menos el n´umero de puntos de la se˜nal, aunque se ha elegido un mayor n´umero de puntos para sobremuestrear la se˜nal y evitar la p´erdida de informaci´on. La ecuaci´on ( 2.2) indica la resoluci´on frecuencial obtenida. 4f=fs NF F T =2400Hz 214 (2.2) Tras el c´alculo de la densidad espectral de potecia con los segmentos artificiales de cada canal, se realiz´o la misma estimaci´on para otros segmentos artificiales que se calcularon para poder promediar los resultados por canales. 2.2.2. Resultados Tras la estimaci´on de la PSD para cada canal promediando varios segmentos artificiales, solo se muestran aqu´ı los resultados de aquellos sensores colocados en posiciones m´as cercanas a la corteza motora cerebral. Esto es, se muestran los resultados obtenidos para los canales C1, C2, C3 y C4. La Figura 2.3 muestra el resultado de aplicar el m´etodo de Welch a cuatro segmentos artificiales y promediarlos despu´es. A la izquierda se muestra la PSD estmada para el canal C1 hasta 120 Hz de frecuencia, donde la linea roja representa el EEG sin artefactos FES, y la linea azul el EEG con FES. Es en el espetro de esta ´ultima se˜nal donde se observan tres componentes de potencia importante situados en 35, 70 y 105 Hz. La primera componente se debe, como se dijo anteriormente, al artefacto de FES que se aplica al paciente con esa frecuencia. Las componentes presentes en 70 y 105 Hz son los dos primeros arm´onicos del artefacto FES, es decir, el FES no aparece en una ´unica frecuencia como se esperaba al inicio, sino que en todos los m´ultiplos de su frecuencia aparecen los arm´onicos. De esta forma, la potencia de la se˜nal de FES est´a repartida entre todas estas frecuencias. 8 2. La Se˜nal de Estimulaci´on El´ectrica Funcional 2.2 An´alisis en el Dominio de la Frecuencia Figura 2.4: Detalle de la distribuci´on de potencia en el promedio realizado en los canales C2 (izquierda) y C4 (derecha). Figura 2.3: Densidad espectral de potencia de los canales C1 (izquierda) y C2 (derecha) tras promediar cuatro segmentos artificales. La Figura 2.4 muestra la distribuci´on de potencias en el rango de inter´es de las se˜nales de EEG. Ante esta gr´afica cabe destacar dos resultados importantes relacionados con las bandas frecuenciales alpha y beta. Como se dijo anteriormente, los ritmos Mu tienen lugar en la banda alpha (8-12 Hz), y su actividad se ve atenuada en entorno a una frecuencia de 10 Hz cuando existe movimiento muscular ([5]). Dado que se est´a produciendo la contracci´on del m´usculo, este es el efecto que se observa en ambas gr´aficas. Adem´as se observa que hay mayor actividad en la parte inferior de la banda beta, comportamiento que se corresponde con la contracci´on sostenida de los m´usculos. Recordemos que esta distribuci´on de potencia se ha realizado sobre el segmento artificial. En la realidad, este pico podr´ıa ser menos acusado, ya que el FES se aplica s´olo durante 3 segundos, no durante los 9 segundos analizados. Vemos por tanto, que el FES no estimula ´unicamente las unidades motoras sobre las que se est´a aplicando, si no que estos efectos tambi´en se ven reflejados sobre la corteza motora. 9 3. Pre-procesado de la Se˜nal EEG En este cap´ıtulo se expone la teor´ıa en que se fundamentan dos procedimientos aplicados a la se˜nal de EEG con FES. Ambos se utilizan como paso previo en los an´alisis realizados a las se˜nales en el siguiente cap´ıtulo. 3.1. Transformada Wavelet 3.1.1. M´etodos La transformada Wavelet es una potente herramienta matem´atica con la que pueden analizarse se˜nales no estacionarias, es decir, se˜nales cuyas propiedades estad´ısticas var´ıan de forma considerable con el tiempo ([7], [8]). En el ´ambito del procesamiento de se˜nal, la transformada wavelet tiene muchas aplicaciones, como puede ser la eliminaci´on de se˜nales de electromiograf´ıa superficial ([9]), eliminaci´on de artefactos oculares ([10]), preprocesado de se˜nal para el reconocimiento de artefactos ([11]), eliminaci´on de ruido ([12]) o junto con ICA ([3], [12], [13]). La utilizaci´on de la transformada Wavelet viene motivada por la insuficiente informaci´on que proporciona el an´alisis de Fourier en determinadas condiciones. Cuando queremos realizar un an´alisis espectral de se˜nales estacionarias, el an´alisis de Fourier proporciona resultados satisfactorios. Sin embargo, cuando la se˜nal no es estacionaria, su espectro en frecuencia cambia con el tiempo, por lo que el an´alisis de Fourier resulta insuficiente. Dada esta situaci´on se recurre al enventanado de la se˜nal para estudiar el espectro en frecuencia de estas se˜nales. Sin embargo, el an´alisis Wavelet es el que permite realizar mejor esta tarea, ya que utiliza ventanas temporales cuya duraci´on depende de las caracter´ısticas de la se˜nal. Esto quiere decir, cuando se requiere una mayor resoluci´on a bajas frecuencias, la ventana temporal es m´as amplia; mientras que cuando se requiere una mayor resoluci´on de alta frecuencia se toman ventanas m´as estrechas. Por todo ello se conoce este an´alisis como multiresoluci´on [7]. Por lo tanto, puede entenderse este an´alisis como una serie de filtros paso-bajo y paso-alto que separan la se˜nal en diferentes bandas de frecuencias. De esta forma, se separa la parte de se˜nal de m´as baja frecuencia de la de m´as alta frecuencia, repiti´endose este proceso sobre las diferentes partes de la se˜nal. La Figura 3.1 muestra un esquema donde pueden verse los pasos sucesivos en este proceso, en el cual el ancho de banda de la se˜nal resultante se va viendo disminuido a la mitad en cada iteraci´on debido a que la se˜nal se va subsampleando. 10 3. Pre-procesado de la Se˜nal EEG 3.1 Transformada Wavelet Figura 3.1: Esquema de la descomposici´on wavelet El resultado de la descomposici´on wavelet se resume en los dos coeficientes que aparecen en este gr´afico: cAycD. El primero de ellos recibe el nombre de Aproximaci´on y contiene la parte de la se˜nal con frecuencias m´as bajas; mientras que el segundo coeficiente, cD, se denomina Detalle y contiene la informaci´on de m´as alta frecuencia [7]. Cuando se quiere realizar la descomposici´on wavelet de una se˜nal de ancho de banda BHz, si se emplea un ´unico nivel de descomposici´on, ´esta quedar´a dividida en dos se˜nales con ancho de banda de B/2 Hz cada una. En tal caso, cA abarcar´ıa un rango de (0-B/2) Hz; mientras que cDser´ıa del rango (B/2−B) Hz. La gr´afica 3.1 muestra m´as niveles de descomposici´on, donde son siempre los coeficientes cAlos que vuelven a ser divididos en niveles superiores. Esta t´ecnica es utilizada con asiduidad en la eliminaci´on de ruido en se˜nales EEG. Para ello, se descompone la se˜nal con las transformadas wavelets en un n´umero determinado de niveles, seg´un el ancho de banda de se˜nal y cu´anto ruido de alta frecuencia queramos eliminar. En este apartado se procedi´o a la descomposici´on de la se˜nal en cinco niveles para reconstruir s´olo hasta el nivel 3 posteriormente. De esta forma se han descartado los coeficientes cD4ycD5 que contienen las componentes de m´as alta frecuencia de la se˜nal. En el siguiente apartadose muestran los resultados obtenidos. 3.1.2. Resultados Tras la reconstrucci´on de la se˜nal, la Figura 3.2 muestra la comparaci´on de las se˜nales original (EEG con artefacto FES, l´ınea roja) y la se˜nal tras la eliminaci´on del ruido mediante la descomposici´on y reconstrucci´on wavelet (EEG sin ruido de alta frecuencia, l´ınea azul). En la gr´afica de la izquierda se observa perfectamente que el artefacto FES no consigue ser eliminado completamente. A la derecha, al observar en mayor detalle las dos se˜nales, puede apreciarse que adem´as la aproximaci´on de la transformada wavelet no es buena en los instantes donde tienen lugar los artefactos de FES. Sin embargo, este resultado no es el mejor de los obtenidos con esta t´ecnica en este proyecto. Al aplicar la transformada wavelet para la eliminaci´on de ruido, se ha constatado a trav´es de los diferentes conjuntos de datos, que la amplitud del artefacto de FES puede resultar determinante para lograr una mejor eliminaci´on. Aplicando este procedimiento al primer conjunto de datos, 11 3. Pre-procesado de la Se˜nal EEG 3.1 Transformada Wavelet Figura 3.2: Eliminaci´on de ruido en la se˜nal EEG con FES. A la izquierda se observa un trial completo, y a la derecha un detalle de ´este, donde se aprecia mejor la eliminacion de ruido de alta frecuencia. en el cual la amplitud de los artefactos de FES es del mismo orden que la se˜nal de EEG, los resultados obtenidos fueron mejores. El resultado de aplicar el algoritmo a este subconjunto de datos se muestra en la Figura 3.3, donde a la se˜nal de EEG con artefactos de FES de peque˜na amplitud (l´ınea roja) est´a superpuesta la se˜nal filtrada con wavelets. Al calcular el coeficiente de correlaci´on entre ambas, se obtuvo un valor de 0’862, confirmando con este dato el buen resultado que se observa gr´aficamente. Figura 3.3: Eliminaci´on de FES de peque˜na amplitud. 12 4. Procesado de Se˜nales EEG con An´alisis de Componentes Independientes 4.1 M´etodos 4.1.4. Wavelet ICA La transformada wavelet se utiliza en el tratamiento de se˜nales biol´ogicas con distintas finalidades como ya se ha comentado anteriormente. En esta secci´on se utiliza la transformada wavelet para descomponer la se˜nal en diferentes rangos de frecuencia. De esta forma, cuando se conoce la frecuencia a la que afecta el artefacto (en este caso el FES), se puede realizar el an´alisis de ICA s´olo sobre aquellos rangos de frecuencia en los que afectan los artefactos [13]. Una vez estimadas las componentes independientes, se identifican aquellas que representan al artefacto para poder eliminarlo. Al utilizar ICA s´olo en aquellas frecuencias contaminadas por el artefacto, se minimizan las p´erdidas de la se˜nal EEG que reside en la componente independiente del artefacto que debe eliminarse para suprimirlo de la se˜nal. En la parte izquierda de la Figura 4.1 se observa el diagrama del algoritmo propuesto en [13], y que aplicaremos al conjunto de se˜nales de EEG con artefactos de FES. El procedimiento es el siguiente: Realizar la tranformada wavelet con el n´umero de niveles adecuado para separar la se˜nal en los rangos de frecuencia que resultan de nuestro inter´es. Seleccionar, mediante inspecci´on visual, las componentes wavelets (rangos de frecuencia) que est´an afectadas por el artefacto y ralizar ICA con esas componentes. Seleccionar las componentes independientes correspondientes al artefacto y eliminarlas, para deshacer posteriormente ICA y obtener las componentes wavelets modificadas, sin artefacto. Por ´ultimo, reconstruir las componentes wavelets para recuperar las se˜nales de EEG en el dominio temporal, esta vez ya sin FES. Figura 4.1: Algoritmo del an´alisis wavelet-ICA. Dado que el espectro en frecuencia de las se˜nales EEG que nos ocupan est´an en el rango de frecuencias de 0 a 1200 Hz, y el FES tiene una frecuencia de 35 Hz, la descomposici´on wavelet necesaria, si se quisiera aislar el artefacto, es de 5 niveles. Sin embargo, la tarea de aislar el artefacto FES no resulta trivial, ya que como se ha explicado, se encuentra mezclado con la se˜nal de EEG. Aun as´ı, se procede a la descomposici´on para tratar de aplicar el algoritmo descrito en las bandas frecuenciales donde se detecten artefactos. 19 4. Procesado de Se˜nales EEG con An´alisis de Componentes Independientes 4.2 Resultados 4.2. Resultados 4.2.1. Aplicaci´on de Fastica e Icasso En este apartado se muestran los resultados obtenidos al aplicar el software Icasso con determinados par´ametros de entrada. Se han realizado dos bater´ıas de pruebas con Icasso, cada una de ellas con una de las funciones objetivo (derivadas de las funciones de contraste G) definidas en (4.4) y (4.5). Se ha configurado Icasso para ejecutar 20 veces el algoritmo fastica con la intenci´on de estudiar la convergencia del mismo en cada prueba y determinar las componentes independientes. As´ı mismo, se ha configurado el n´umero m´aximo de iteraciones para alcanzar la convergencia en 500 iteraciones. La Figura 4.2 muestra los resultados obtenidos al ejecutar Icasso con la funci´on definida en la Ecuaci´on (4.4). A la izquierda se observa el origen de cada una de las componentes, es decir, en qu´e sensores se han registrado. Esta gr´afica, tambi´en conocida como topoplot, permite hacer una primera aproximaci´on acerca de qu´e componentes son artefactos y qu´e componentes no lo son. Observando ´unicamente estos mapas podr´ıa decirse que las componentes 1 y 5 son las que representan el FES. Esta suposici´on se explica al observarse que la actividad el´ectrica es muy similar en todos los electrodos. Sin embargo, tambi´en podr´ıa decirse lo mismo de las componentes 7 y 11, aunque ´estas se descartan al observar la gr´afica que aparece a la derecha en esa misma figura. Figura 4.2: Componentes Independientes estimadas por Icasso con la funci´on definida en la Ecuaci´on (4.4) Al observar la gr´afica situada a la derecha se observa que las componentes est´an ordenadas seg´un el´ındice Iqo´ındice de similaridad del cluster. Tras obtener las componentes independientes (ICs) en cada ejecuci´on del algoritmo fastica, ´esta son projectadas en un espacio 2D para hacer posible su visualizaci´on. Una vez que se han projectado las 20 soluciones de ICs, se procede a agruparlas en clusters, siendo el ´ındice Iqel que eval´ua c´omo de compacto y aislado es cada cluster. Es decir, es una forma cuantitativa de medir el parecido de las ICs estimadas y agrupadas en el mismo cluster, e interesa que sea lo m´as cercano posible a la unidad. Esto significar´a, que las estimaciones agrupadas en un determinado conjunto representan la misma componente independiente. 20 4. Procesado de Se˜nales EEG con An´alisis de Componentes Independientes 4.2 Resultados Como se observa en la parte izquierda de la Figura 4.3, comienza a haber m´as dispersi´on en el cl´uster a partir de cluster n´umero 6, ya que el ´ındice de calidad comienza a disminuir notablemente a partir de dicho cluster. En la parte derecha se contabiliza el n´umero de veces que se ha estimado cada componente tras la ejecuci´on completa de icasso, habi´endose estimado las componentes pertenecientes a los clusters m´as compactos en cada ejecuci´on del algoritmo fastica. Figura 4.3: ´ Indice de calidad de los clusters (izquierda) y similaridades entre las estimaciones de cada cluster (derecha). A la derecha de esa misma figura se observan las similaridades entre las estimaciones de cada cluster, siendo las componentes 1 a la componente 6 las que mayor ´ındice de similaridad alcanzan. Este concepto se aprecia con mayor facilidad en la Figura 4.4, donde se observa la proyecci´on de las componentes al espacio 2D y los clusters realizados. Figura 4.4: Clusters de las ICs estimadas por icasso. 21 4. Procesado de Se˜nales EEG con An´alisis de Componentes Independientes 4.2 Resultados Analizando los resultados del an´alisis ICA, podemos decir que la componente de FES est´a representada por las ICs 1 y 5, que adem´as siempre son estimadas y con un alto ´ındice Iq, aunque no se ha conseguido separar completamente el FES de la se˜nal de EEG en una ´unica IC. La explicaci´on de este comportamiento es que el algoritmo, cuando converge, no lo hace a un m´ınimo global, sino que converge a m´ınimos locales ([17]). Cuando se tienen funciones no cuadr´aticas, como es este caso, existen muchos m´ınimos y m´aximos, por lo que encontrar el m´ınimo global tiene mayor dificultad. En estas situaciones es importante el punto inicial donde fastica comienza a iterar para lograr encontrar el m´ınimo global. La segunda bater´ıa de pruebas se ha realizado con la funci´on de coste definida en la Ecuaci´on (4.5). En esta ocasi´on, el topoplot mostrado a la zquierda en la Figura 4.5, indica que las ICs que representan el FES son las de las componentes 4 y 5, suposici´on que se confirma al observar la gr´afica derecha de la misma figura. Figura 4.5: Componentes Independientes estimadas por Icasso con la funci´on definida en la Ecuaci´on (4.5) Observando los resultados referentes al ´ındice de calidad en la Figura 4.6, en esta ocasi´on la dispersi´on dentro de cada cluster resulta m´as evidente a partir de la quinta agrupaci´on, es decir, que el valor de Iqdisminuye antes. Se observa tambi´en a la derecha en esta figura los clusters formados, siendo en este caso m´as dispersos, es decir, con menor ´ındice de calidad, que en el caso anterior. 22 4. Procesado de Se˜nales EEG con An´alisis de Componentes Independientes 4.2 Resultados Figura 4.6: ´ Indice de calidad (izquierda) y agrupamiento realizado por Icasso con la funci´on definida en la Ecuaci´on (4.5) Tras la realizaci´on de esta prueba y por lo que se observa en el gr´afico de ICs estimadas, Icasso tampoco ha conseguido en esta ocasi´on recuperar la componente de FES en una se˜nal ´unica, sino que tambi´en la estima en dos componentes. Sin embargo, en esta ocasi´on las ICs estimadas est´an ordenadas de distinta forma, ahora se encuentran en las componentes 4 y 5 an lugar de 1 y 4 como ocurr´ıa en la primera prueba de Icasso. Esto se debe a que las componentes est´an ordenadas seg´un el Iq, distinto a lo que ocurr´ıa en PCA, donde las componentes estaban ordenadas de mayor a menor varianza. 4.2.2. Aplicaci´on de runica Esta secci´on muestra los resultados obtenidos al aplicar el algoritmo runica sobre uno de los trials en su formato original. Se han realizado varias ejecuciones del algoritmo, que se ha realizado con una condici´on de parada de 1 ·10−5en lugar de la condici´on establecida por defecto (1 ·10−6). Este cambio ha venido motivado por la falta de convergencia con la condici´on configurada por defecto. Tras aplicar este algoritmo sobre uno de los trial en su formato original de EEG contaminado con FES, tampoco se logr´o una recuperaci´on de la se˜nal de FES en una ´unica componente, quedando dividida de nuevo en dos componentes. La Figura 4.7 muestra el EEG generado por las componentes 1 y 11. Con este algoritmo consigue eliminarse mejor el FES de una de las componentes, ya que en la gr´afica de la derecha se ve que no afecta a todos los canales. Adem´as la componente 1 tiene una amplitud mucho mayor que la componente 11. 23 4. Procesado de Se˜nales EEG con An´alisis de Componentes Independientes 4.2 Resultados Figura 4.7: EEG generado por las Componentes Independientes 1 y 11 con el algoritmo runica. Realizando una nueva prueba en la que, ahora s´ı, se configura el par´ametro extended para que considere distribuciones supergaussianas de las ICs, se obtienen los resultados se muestran en la Figura 4.8. En ella se observa el EEG generado por dos de las componentes estimadas, que como puede apreciarse, ambas contienen parte del FES a pesar de haber considerado la supergaussianidad. En esta ocasi´on las componentes de FES estimadas resultan ser las componentes 1 y 3. Figura 4.8: EEG generado por las Componentes Independientes 1 y 3 con el algoritmo runica-extended. El motivo por el cual este algoritmo no ha funcionado es porque la se˜nal de FES resulta ser una se˜nal super supergaussiana y con una desviaci´on est´andar muy peque˜na, lo que hace muy dif´ıcil la estimaci´on de la componente. 24 4. Procesado de Se˜nales EEG con An´alisis de Componentes Independientes 4.2 Resultados 4.2.3. Aplicaci´on de amuse La aplicaci´on del algoritmo amuse requiere la elecci´on del tama˜no de la ventana a considerar en las autocovarianzas. En este caso se han estudiado distintos tama˜nos de ventana, habi´endose elegido finalmente una ventana de 87’5 ms, lo suficientemente grande como para considerar 2-3 picos de FES. Los resultados de aplicar este algoritmo se muestran en la Figura 4.9, que refleja las se˜nales de EEG generadas por dos de las componentes estimadas. Existe tambi´en una tercera componente donde quedan residuos de la se˜nal de FES. De nuevo la t´ecnica empleada no separa en una ´unica componente independiente la se˜nal de FES. Figura 4.9: EEG generado por las ICs estimadas por el algoritmo amuse. A la izquierda el producido por la componente 15, de mayor amplitud; y a la derecha el generado por la componente 16. 4.2.4. Aplicaci´on de wavelet-ICA Tras haber aplicado el algoritmo explicado con un nivel de descomposici´on 5, los resultados obtenidos fueron muy similares a los obtenidos con fastica y runica. Con esta t´ecnica, el artefacto tampoco puede recuperarse ya que no se identifican todas las bandas de frecuencias que est´an afectadas por el artefacto. 25 5. Conclusiones y Trabajo Futuro En este proyecto se han empleado diversas t´ecnicas de eliminaci´on de artefactos con las que anteriormente se hab´ıan obtenido buenos resultados en la extracci´on o eliminaci´on de este tipo de se˜nales; aunque en el presente trabajo no se han alcanzado resultados completamente satisfactorios. Sin embargo, las diferencias existentes entre aqu´ellos y los artefactos de FES que han ocupado este estudio son relevantes. Los artefactos causados por la se˜nal de electroestimulaci´on tienen una frecuencia m´as elevada y con cierta periodicidad temporal que imposibilita el filtrado debido a la mezcla frecuencia de se˜nales. Respecto al an´alisis de componentes se ha llegado a la conclusi´on de que PCA es una t´ecnica insuficiente para resolver el problema aqu´ı planteado. As´ı mismo hay que destacar que el estudio de las t´ecnicas aqu´ı empleadas ha permitido conocer la carencia de los algorimos ICA actuales con respecto a la estimaci´on de componentes independientes cuyas distribuci´on sea super supergaussiana. Con respecto a la convergencia de los algoritmos, para minimizar el problema de los m´ınimos locales, se propone iniciar la b´usqueda del m´ınimo global en un punto cercano a ´este, para evitar que el algoritmo se quede en uno de los m´ınimos locales y no se alcance la convergencia en el m´ınimo global. La utilizaci´on de transformadas wavelet ha permitido, mediante la eliminaci´on de ruido de alta frecuencia, atenuar los efectos del FES. Sin embargo, debido al problema de mezclado existente, un filtrado convencional resulta ineficiente en la eliminaci´on de la se˜nal de FES. Se propone continuar la documentaci´on sobre algoritmos ICA que permitan estimar las componentes independientes de distribuciones supergaussianas y super supergaussianas, as´ı como profundizar en la herramienta que supone runica en este sentido. Adem´as se ha comprobado que la actividad neuronal responde ante estos est´ımulos de forma similar a como lo har´ıa cuando la corteza motora se estimula de forma fisiol´ogicamente natural. 26 Bibliograf´ıa [1] Paul L. Nunez and Ramesh Srinivasan. Electric Fields of the Brain. Oxford University Press, 2006. [2] D. Corydon Hammond and Jay Gunkelman. The Art of Artifacting. International Society for Neurofeedback & Research, 2001. [3] Valeri A. Makarov Nazareth P. Castellanos. Recovering eeg brain signals: Artifact suppression with wavelet enhanced independet component analysis. In Journal of Neuroscience Methods, pages 300–312. Elsevier, 2006. [4] SALUD CENETEC. Estimulador nervioso el´ectrico transcut´aneo (tens), Septiembre 2005. [5] A. James Rowan and Eugene Tolunsky. Conceptos b´asicos de EEG. Elsevier, 2004. [6] Luis Vicente Borruel. Estimaci´on espectral en aplicaciones biom´edicas. http://www.masterib.es/, Marzo 2008. [7] Samir Kouro R. and Rodrigo Musalem M. Tutorial introductorio a la teor´ıa de wavelet. T´ecnicas Modernas en Autom´atica, Julio 2002. [8] Introducci´on a la transformada wavelet. descomposici´on de se˜nales. http://www.exa.unicen.edu.ar/escuelapav/cursos/wavelets/apunte.pdf, 2006. [9] E La Foresta B. Azzerboni, M. Carpentien and E C. Morabito. Neural-ica and wavelet transform for artifacts removal in surface emg. IEEE, pages 3223–3228, 2004. [10] K. Sivakumar P. Senthil Kumar, R. Arumuganathan and C. Vimal. Removal of ocular artifacts in the eeg through wavelet transform without using an eog reference channel. Int. J. Open Problems Compt. Math., 1(3):188–200, December 2008. [11] Piotr Durka W.Szelenberger Rafal Ksiezyk, Katarzyna Blinowska and W. Androsiuk. Neural networks with wavelet preprocessing in eeg artifact recognition. [12] Janett Walters-Williams and Yan Li. A new approach to denoising eeg signals - merger of translation invariant wavelet and ica. International Journal of Biometrics and Bioinformatics, 5:130–148, 2011. [13] Nadia Mammone Giuseppina Inuso, Fabio La Foresta and Francesco Carlo Morabito. Wavelet-ica methodology for efficient artifact removal from electroencephalographic recordings. In Proceedings of International Joint Conference on Neural Networks, pages 12–17. IEEE, 2007. 27 [14] Juha Karhunen Aapo Hyv¨arinen and Erkki Oja. Independent Component Analysis. Wiley Interscience, 2001. [15] Luis Vicente Borruel. Separaci´on ciega de fuentes: Pca e ica. http://www.masterib.es/, Marzo 2011. [16] Aapo Hyv¨arinen. Fast and robust fixed-point algorithms for independent component analysis. In IEEE Trans. on Neural Networks, pages 626–634, 1999. [17] Johan Himberg and Aapo Hyv¨arinen. Icasso: Software for investigating the reliability of ica estimates by clustering and visualization. In IEEE Workshop on Neural Networks for Signal Processing, pages 259–268, 2003. [18] Aapo Hyv¨arinen Johan Himberg and Fabrizio Esposito. Validating the independet components of neuroimaging time-series via clustering and visualization. NeuroImage, 2004. [19] Aapo Hyv¨arinen. http://www.cis.hut.fi/projects/ica/fastica/. [20] Anthony J. Bell and Terrence J. Sejnowski. Blind separation and blind deconvolution: An information-theoretic approach. In IEEE International Conference on Acoustics, Speech and Signal Processing, pages 3415–3418, 1995. [21] Mark Girolami Te-Won Lee and Terrence J. Sejnowski. Independent component analysis using an extended infomax algorithm for mixed subgaussian and supergaussian sources. In Neural Computation. Massachusetts Institute of Technology, pages 417–441, 1999. 28