CLASIFICACIÓN DE SEÑALES CEREBRALES RELACIONADAS CON LA IMAGINACIÓN DE MOVIMIENTOS PARA APLICACIONES DE BCI LAURA MARCELA SALAZAR CRUZ


Save this PDF as:
 WORD  PNG  TXT  JPG

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

Download "CLASIFICACIÓN DE SEÑALES CEREBRALES RELACIONADAS CON LA IMAGINACIÓN DE MOVIMIENTOS PARA APLICACIONES DE BCI LAURA MARCELA SALAZAR CRUZ"

Transcripción

1 CLASIFICACIÓN DE SEÑALES CEREBRALES RELACIONADAS CON LA IMAGINACIÓN DE MOVIMIENTOS PARA APLICACIONES DE BCI JORGE ANDRÉS ZACCARO VALVERDE JAVIER HERNÁN CÁRDENAS DÍAZ LAURA MARCELA SALAZAR CRUZ TRABAJO DE GRADO PARA OPTAR POR EL TÍTULO DE INGENIERO ELECTRÓNICO DIRECTORES: ING. JULIÁN ARMANDO QUIROGA SEPÚLVEDA, M.Sc. ING. DIEGO ALEJANDRO PATIÑO GUEVARA, Ph.D. PONTIFICIA UNIVERSIDAD JAVERIANA FACULTAD DE INGENIERÍA DEPARTAMENTO DE ELECTRÓNICA BOGOTÁ DC, 2011

2 PONTIFICIA UNIVERSIDAD JAVERIANA FACULTAD DE INGENIERÍA CARRERA DE INGENIERÍA ELECTRÓNICA RECTOR MAGNÍFICO: PADRE JOAQUIN EMILIO SÁNCHEZ GARCÍA S.J DECANO ACADÉMICO: ING. FRANCISCO JAVIER REBOLLEDO MUÑOZ DECANO DEL MEDIO UNIVERSITARIO: PADRE SERGIO BERNAL RESTREPO S.J DIRECTOR DE CARRERA: ING. JUAN MANUEL CRUZ BOHÓRQUEZ DIRECTORES DE PROYECTO: ING. JULIÁN ARMANDO QUIROGA SEPÚLVEDA. M.Sc., ING. DIEGO ALEJANDRO PATIÑO GUEVARA, Ph.D.

3 NOTA DE ADVERTENCIA Artículo 23 de la Resolución N 13 de Julio de 1946 La Universidad no se hace responsable por los conceptos emitidos por sus alumnos en sus trabajos de tesis. Solo velará por que no se publique nada contrario al dogma y a la moral católica y por que las tesis no contengan ataques personales contra persona alguna, antes bien se vea en ellas el anhelo de buscar la verdad y la justicia. 3

4 Índice NOTA DE ADVERTENCIA 3 1. INTRODUCCIÓN 7 2. MARCO TEÓRICO SEÑALES CEREBRALES Adquisición ANÁLISIS DE SEÑALES Análisis en Frecuencia Análisis no Lineal en Tiempo Métodos de Análisis no Lineal-Exponentes de Lyapunov Análisis Estadístico y Transformaciones Lineales RECONOCIMIENTO ESTADÍSTICO DE PATRONES Espacio de características Métodos paramétricos Métodos no paramétricos Aprendizaje no supervisado Aprendizaje supervisado ESPECIFICACIONES Diagrama de Bloques MEL: Máximo Exponente de Lyapunov DESARROLLOS Descripición de las señales Métodos de extracción de características Discrete Wavelet Transform (DWT), Transformada Discreta Wavelet Continuous Wavelet Transform (CWT), Transformada Continua Wavelet Common Spatial Patterns (CSP), Patrones Espaciales Comunes Máximo exponente de Lyapunov Definición del vector de características PCA e ICA Métodos de Clasificación Linear Discriminant Analysis (LDA), Análisis Lineal Discriminante Support Vector Machines (SVM), Máquinas de Vectores de Soporte Artificial Neural Networks (ANN), Redes Neuronales Artificiales ANÁLISIS DE RESULTADOS Cuadros de resultados Discusión Resultados anteriores CONCLUSIONES 37 4

5 Índice de figuras 1. Esquema descripción del problema. Tomada [37],[9], Matlab Help, [38], [27] Sistema Internacional ubicación de electrodos. Tomado de [37] Algoritmo de Wolf Máquina de Vectores de Soporte Diagrama de bloques del sistema Diagrama de bloques del sistema Configuración en tiempo de cada ensayo. Tomada de [9] Descomposición Wavelet de 3 niveles. Tomada de Matlab Help Transformada Continua Wavelet de un ensayo CWT del canal C3, Ensayo CWT del canal C4, Ensayo Mejores Resultados con PCA Mejores Resultados 4 Subconjuntos

6 Índice de cuadros 1. Bandas generadas por la descomposición utilizada Vector de características Parámetros SVM Parámetros ANN Mejores Resultados Resultados SVM Resultados LDA Resultados ANN Mejores Resultados 4 Subconjuntos

7 1. INTRODUCCIÓN Los avances tecnológicos de las últimas décadas han llevado al ser humano a rodearse de una gran variedad de máquinas y sistemas computacionales, cuyas aplicaciones van desde los medios de transporte y el manejo de información hasta el entretenimiento y las comunicaciones. A medida que diferentes productos han sido refinados en funcionalidad, la interactividad ha pasado a tomar cada vez mayor importancia para el usuario, por lo que han surgido fuertes tendencias enfocadas en que la comunicación entre el hombre y las máquinas sea cada vez más espontánea. Una de ellas es conocida como Brain-Computer Interfaces (BCI), en la que el objetivo es crear sistemas que procesen las señales cerebrales de un individuo y las conviertan en comandos específicos para un sistema de control dado (e.g. computador). Esto es posible gracias a que la actividad cerebral se manifiesta como potenciales eléctricos que pueden ser medidos tanto de forma invasiva (implantes de electrodos) como no invasiva (EEG, MEG, fmri), y a que dicha actividad está estrechamente relacionada con procesos motrices, cognitivos y emocionales del cerebro. Entre las principales aplicaciones de BCI se encuentran la predicción de epilepsia, la rehabilitación de la movilidad, las comunicaciones militares y la interacción con entornos virtuales. Esta última es posible gracias a que la actividad cerebral relacionada con la imaginación de movimientos puede ser registrada, procesada y convertida en comandos de acción coherentes para el contexto, como por ejemplo, en movimientos del puntero de un computador [1], de un personaje en un videojuego o de objetos virtuales en general. Existen varios laboratorios alrededor del mundo que se encuentran desarrollando proyectos relacionados con esta aplicación, entre los cuales se encuentran el Laboratory of Brain-Computer Intefaces en Graz (Austria) [14], el Bernstein Center for Computational Neuroscience en Berlín (Alemania) [15] y el Institute of Neural Engineering de la Universidad Tsinghua en Beijing (China) [16]. Investigadores de estos y otros laboratorios han participado en la Competencia BCI organizada por el grupo Berlin BCI, en la que bases de datos de señales cerebrales son puestas a disposición del público para el entrenamiento y prueba de clasificadores propuestos por los participantes. La realización de este proyecto servirá de preparación para considerar la participación en alguna de las próximas competencias. En cuanto a las proyecciones comerciales, existen compañías como Emotiv Systems y Neurosky que se encuentran desarrollando productos de software y hardware para BCI, y que están interesadas en resolver problemas como el que se tratará en este proyecto. Además, compañías fuertemente comprometidas con la interactividad como Apple y Google, y líderes del mercado de los videojuegos como Microsoft, Nintendo y Sony, podrían estar dispuestas a incorporar aplicaciones de BCI en su portafolio de servicios. En el contexto colombiano, este proyecto marcará un punto de referencia de investigación en esta tendencia, y servirá de invitación para el desarrollo de investigaciones similares, de modo que la Pontificia Universidad Javeriana e investigadores de nuestro país continúen haciendo aportes a esta prometedora aplicación. 7

8 En la Figura 1 se muestra el diagrama de bloques que representa uno de los procedimiento a seguir en problemas de clasificación de señales cerebrales: adquisición, extracción de características, mejoramiento del espacio de características (opcional) y clasificación. Figura 1: Esquema descripción del problema. Tomada [37],[9], Matlab Help, [38], [27]. Dado que la actividad cerebral relacionada con la imaginación de movimientos (MI) se presenta principalmente en las bandas de frecuencia µ (8-12 Hz) y β (18-25 Hz), y que las aplicaciones finales de BCI precisan de la implementación de los clasificadores en tiempo real, se requiere tanto buena resolución en tiempo como en frecuencia. Entre las soluciones más efectivas para este problema se encuentran los métodos basados en wavelets, pues a diferencia de los enfoques por Transformada de Fourier, gozan de la propiedad de multi-resolución. Publicaciones como [9], [10], [11], [12] han comprobado el buen desempeño en la clasificación MI al realizar la extracción de características utilizando la Transformada Wavelet. Por otro lado, también resulta conveniente usar métodos para extraer información temporal de las señales. Uno de los más utilizados para la predicción y detección de eventos es el análisis de caoticidad por Lyapunov, en el que las series de tiempo son representadas como trayectorias en el espacio de fase, y su caoticidad puede ser cuantificada por la distancia entre las mismas. En cuanto al vector de características, como no existen criterios fijos para su construcción, algunos de sus componentes pueden resultar correlacionados. La aplicación de transformaciones como PCA (Principal Component Analysis) o ICA (Independent Component Analysis) podrían llegar a mejorar el espacio de características al reducir su dimensión o la dependencia entre sus componentes. Algunas publicaciones se dedican a comparar el desempeño de ciertos clasificadores para exactamente el mismo problema. Entre las comparaciones que sirvieron de orientación se encuentran: Lineal vs NN-HMM [3], LDA vs SVM [4], [5], [6], LDA-QDA-Fisher y SVM-MLP-kNN. Teniendo en cuenta los resultados presentados en publicaciones como estas y en [11] y [13], se escogieron los clasificadores que presentaron los mejores resultados: Máquinas de Vectores de Soporte (SVM), Análisis Lineal Discriminante (LDA) y Redes Neuronales Artificiales (ANN). El propósito de este trabajo es comparar el desempeño de tres estrategias de clasificación para el reconocimiento de patrones de señales cerebrales asociadas a la imaginación de movimientos de la mano izquierda (clase 1) y la mano derecha (clase 2), y evaluar el desempeño al aplicar PCA o ICA para mejorar el espacio de características o reducir la redundancia de las señales originales. 8

9 2. MARCO TEÓRICO 2.1. SEÑALES CEREBRALES La actividad eléctrica cerebral puede ser registrada e interpretada como una serie de tiempo estocástica no estacionaria, en la que la información asociada a los procesos mentales que la generan (e.g. imaginación de movimientos) no se puede observar de forma directa, sino que debe ser caracterizada por medio de herramientas como la Transformada Wavelet, y clasificada usando métodos como las Máquinas de Vectores de Soporte (SVM) o las Redes Neuronales Artificiales (ANN) Adquisición La adquisición de la actividad cerebral es un procedimiento delicado tanto en el caso invasivo como en el no invasivo. En el primero, se debe ser extremadamente cuidadoso para que los instrumentos utilizados no causen lesiones en el tejido nervioso manipulado. En el segundo, la calidad del procedimiento está determinada por factores como la relación señal a ruido y la sensibilidad de los sensores empleados. La gran mayoría de dispositivos BCI son diseñados para operar de forma no invasiva, por lo que el método de adquisición más utilizado es el Electroencefalograma (EEG). Este consiste en ubicar un conjunto de electrodos (sensores) en posiciones predefinidas de la corteza cerebral para capturar las señales eléctricas generadas sobre regiones determinadas del cerebro. En la Figura 2 se muestra el Sistema Internacional de Ubicación de Electrodos. Figura 2: Sistema Internacional ubicación de electrodos. Tomado de [37]. En aplicaciones BCI, se utilizan términos específicos para hacer referencia a distintos elementos del proceso de adquisición y registro de la actividad cerebral: Definición 2.1. Un canal es el resultado de ubicar un electrodo en la posición indicada por un sistema como el sistema De este modo, al hacer referencia al canal C3 en realidad se está hablando del electrodo ubicado en la posición C3 de un sistema dado. Definición 2.2. Una muestra es un número que representa la actividad eléctrica cerebral en un instante de tiempo. Definición 2.3. Un ensayo es un conjunto de muestras registradas consecutivamente durante un intervalo de tiempo dado. A continuación se presentarán herramientas que se utilizarán a lo largo de todo el documento. 9

10 2.2. ANÁLISIS DE SEÑALES Análisis en Frecuencia La trasformada de Fourier brinda información sobre las componentes de frecuencia de una señal. Sin embargo, esta información no se encuentra relacionada de forma explícita con la información temporal. Por ejemplo, para una señal no estacionaria no se tiene información sobre el tiempo en el cual se presenta o no una determinada componente de frecuencia [17]. Definición 2.4. Una ventana temporal g(t) es una señal de energía finita que se encuentra concentrada en un intervalo acotado de tiempo [17]. Definición 2.5. Un átomo tiempo-frecuencia g τ,ξ (t) se define como g τ,ξ (t) = g τ (t)e jξt, donde g τ (t) es una ventana temporal con centro en τ, y ξ R [17]. Para analizar señales no estacionarias es necesario utilizar átomos tiempo-frecuencia con diferentes soportes en tiempo de acuerdo con la frecuencia que se vaya a analizar. Con este método se puede encontrar la transformada de Fourier de tiempo corto de una señal x(t) L 2. Sin embargo, dado que con esta no es posible tener una buena resolución en tiempo y en frecuencia simultáneamente resulta conveniente utilizar otras técnicas. Definición 2.6. Una wavelet es una función ψ L 2 (R) (energía finita) con media cero que satisface [18]: ψ(t)dt = 0 (1) Es normalizada ψ = 1, y centrada en la vencindad de t = 0. Una familia de wavelets se obtiene escalando ψ por s y trasladándola en u: ψ u,s (t) = 1 ( ) t u ψ (2) s s Estos átomos permanecen normalizados: ψ(u, s) = 1. La principal característica de las wavelets es su localización tiempo-frecuencia, por lo que la mayoría de su energía se encuentra concentrada en un intervalo finito de tiempo. Esta característica hace posible que el análisis wavelet tenga buena resolución en frecuencias bajas (ventanas grandes en tiempo) y en frecuencias altas (ventanas pequeñas en tiempo) simultáneamente [19]. La transformada wavelet descompone señales en wavelets comprimidas y expandidas para analizar sus características tiempo-frecuencia, y goza de la propiedad de multiresolución. Definición 2.7. La transformada wavelet de una señal x(t) L 2 (R) en el tiempo u y escala s está dada por [18]: W x (u, s) = x(t), ψ u,s (t) = 1 ( ) t u x(t)ψ dt (3) s s Para cada par de valores (u, s) R R +, W x (u, s) almacena el resultado de la comparación entre x(t) y ψ u,s (t) [17]. La descomposición de la señal da lugar a un conjunto de coeficientes llamados coeficientes wavelet. De este modo, la señal se puede reconstruir como una combinación lineal de funciones wavelet, teniendo en cuenta que para obtener una reconstrucción exacta de la señal se debe calcular el número adecuado de coeficientes. 10

11 Definición 2.8. Sea f(t) una señal de tiempo continuo que es muestreada uniformemente y f[n] = f(n) la señal discreta de tamaño N. Sea ψ(t) una wavelet cuyo soporte se encuentra en [ K/2, K/2]. Para 2 a j NK 1 una wavelet discreta escalada por a j se define [18]: ψ j [n] = 1 a j ψ ( n a j ) (4) La wavelet discreta tiene Ka j valores diferentes de cero en [ N/2, N/2]. Para evitar problemas de frontera, se consideran f[n] y ψ j [n] como señales periódicas de periodo N. Definición 2.9. La transformada wavelet discreta puede ser escrita como una convolución circular ψ j [n] = ψj [ n] [18]: Análisis no Lineal en Tiempo Serie de Tiempo W f [n, a j ] = N 1 m=0 f[m]ψ j [m n] (5) Una serie de tiempo es una sucesión temporal ordenada de observaciones tomadas a lo largo del tiempo. Estacionariedad Si se tiene una serie de tiempo de valores discretos {x 1, x 2,..., x N }, estacionariedad significa que la función de distribución de probabilidad conjunta f 12 (x 1, x 2 ) depende únicamente de la diferencia temporal t 1 t 2 y no de los valores absolutos t 1 y t 2. De forma más clara, estacionariedad significa que las características de la serie de tiempo, tal como media, varianza o espectro de potencia, no cambian con el tiempo [20], [21] Métodos de Análisis no Lineal-Exponentes de Lyapunov Espacio de fase El espacio de fase es un espacio en el que se puede describir el comportamiento dinámico de un sistema, donde las coordenadas espaciales representan las variables de estado, y la evolución del sistema se puede analizar a partir de trayectorias que dependen de las condiciones iniciales [21]. Uno de los enfoques para construir el espacio de fase es el método de los retardos, el cual consiste en calcular el retardo τ y la dimensión de embebimiento m que permiten construir la matriz [22] Exponentes de Lyapunov X(t) = {x(t), x(t + τ),..., x[t + (m 1)τ]} T (6) Definición Dado un sistema dinámico continuo en un espacio de fase de n dimensiones, se observa la evolución a largo plazo de una n-esfera infinitesimal de condiciones iniciales; la esfera se convertirá en una n-elipsoide debido a la naturaleza de deformación del flujo. El i-ésimo exponente de Lyapunov de una dimensión λ i se define en términos de la longitud del eje principal de la elipsoide p i (t) [23] 1 λ i = lím t t log p i (t) 2 p i (0) (7) 11

12 Los λ i se encuentran ordenados de mayor a menor. Como la orientación de una elipsoide cambia continuamente a medida que esta va evolucionando, las direcciones asociadas con un exponente dado varían alrededor del atractor, por lo que no se puede hablar de una dirección bien definida asociada a un exponente dado [24]. Esto implica que para un sistema dinámico, la sensitividad a las condiciones iniciales está cuantificada por los exponentes de Lyapunov. Por ejemplo, si se consideran dos trayectorias con condiciones iniciales cercanas a un atractor múltiple, cuando el atractor es caótico las trayectorias divergen en promedio a una tasa exponencial caracterizada por el exponente de Lyapunov más grande. Este concepto se puede generalizar para el espectro de los exponentes de Lyapunov, λ i (i = 1,..., n), al considerar una pequeña esfera de condiciones iniciales de n dimensiones, donde n es el número de ecuaciones (ó el número de variables de estado) utilizadas para describir el sistema. Como se mencionó anteriormente, a medida que el sistema evoluciona, la esfera se convierte en un elipsoide cuyos ejes principales se expanden o se contraen en tasas dadas por los exponentes de Lyapunov. La presencia de un exponente positivo es suficiente para diagnosticar caos y representa la inestabilidad local en una dirección particular. Nótese que para la existencia de un atractor, la dinámica global debe ser disipativa, es decir, globalmente estable, y la tasa total de contracción debe sobrepesar la tasa total de expansión. Luego, incluso cuando hay muchos exponentes de Lyapunov positivos, la suma sobre todo el espectro es negativa. Las magnitudes de los exponentes de Lyapunov cuantifican la dinámica de un atractor en términos de teoría de la información. Los exponentes miden la tasa a la cual el sistema crea o destruye información; por lo que los exponentes se expresan en bits de información por segundo, ó bits por órbita para un sistema continuo y bits por iteración para un sistema discreto. La tasa promedio a la cual la información contenida en transientes se pierde puede ser determinada por los exponentes negativos. El decremento asintótico de una perturbación del atractor es gobernada por el exponente menos negativo, razón por la cual es el exponente más fácil de estimar [24]. Métodos para Estimar los Exponentes de Lyapunov Algoritmo de Wolf et al. [24]: Dada una serie de tiempo x(t), se reconstruye el espacio de estados por medio de coordenadas con retardos. Se localiza el vecino más cercano (por medio de la distancia euclideana) del punto X(t) = {x(t), x(t + τ),..., x[t + (m 1)τ]} T (8) y se denota la distancia entre estos dos puntos como L(t 0 ). En un tiempo posterior t 1, la longitud inicial evoluciona a la longitud L (t 1 ) y el elemento de longitud se propaga a través del atractor para un tiempo lo suficientemente corto para que solo una pequeña parte de la estructura del atractor sea examinada. Si la evolución de tiempo es muy grande, se puede observar que L se comprime a medida que las dos trayectorias que lo definen pasan a través de una región de plegamiento del atractor, lo cual llevaría a una mala estimación de λ 1. Ahora se busca por un nuevo punto que satisfaga dos criterios: que su separación L(t 1 ) desde el punto de evolución de referencia sea pequeña, y que la separación angular entre los elementos evolucionados y de reemplazo sea pequeña. (Ver Figura 3). Si no se puede encontrar un punto de reemplazo adecuado, se retienen los puntos que se han utilizado. Este procedimiento se repite hasta que las trayectorias hayan atravesado todo el conjunto de datos, y en este momento se estima 1 λ 1 = t M t 0 M k=1 donde M corresponde al número total de reemplazos [24]. log 2 L (t k ) L ( t k 1 ) (9) 12

13 Figura 3: Representación esquemática de la evolución y procedimiento de reemplazo utilizado para estimar el mayor exponente de Lyapunov a partir de una serie de tiempo. Tomada de [24] Análisis Estadístico y Transformaciones Lineales Common Spatial Patterns (CSP), Patrones Espaciales Comunes Partiendo del hecho que la generación de la actividad cerebral relacionada con el movimiento es contralateral, resulta conveniente utilizar una estrategia que aproveche tanto la información en frecuencia como la localización espacial de las señales EEG. Una de las mayores ventajas del electroencefalograma reside en la información espacial que se agrega a los datos adquiridos con configuraciones multicanal. La existencia de información a priori que relaciona ciertas áreas de la corteza cerebral con acciones específicas como la imaginación de movimientos justifica la utilización de dicha información espacial para el procesamiento de las señales EEG. Esta se puede incluir en formas que van desde criterios de prioridad para la selección de canales hasta representaciones matriciales que consideran todas las series de tiempo como un conjunto. La representación matricial ofrece numerosas ventajas para el procesamiento de señales EEG multicanal, entre las cuales se encuentra la posibilidad de encontrar transformaciones lineales que resalten características de las señales observadas, como la varianza temporal y la independencia estadística. Basándose en los cambios de frecuencia y de energía de un canal EEG durante la imaginación de un movimiento, es posible suponer que dichos cambios a su vez implican cambios en la varianza de la serie de tiempo correspondiente, por lo que dicha medida estadística podría tenerse en cuenta como una característica para clasificación. Una de las técnicas que aprovecha la representación matricial de las señales EEG multicanal es la técnica de Common Spatial Patterns (CSP), pues al aplicar ciertas transformaciones lineales busca generar una matriz cuyas varianzas sean más apropiadas que las de la matriz original para la discriminación de las clases del experimento [25]. 13

14 Para encontrar la transformación lineal deseada, se deben calcular C filtros a partir de las matrices de covarianza R C, donde C es el número de clases y n c es el número de canales. R c = i c R ci n c (10) La matriz de covarianza R ci de cada ensayo X i, se calcula como [26]: R ci = X i X i T traza(x i X i T ), i = 1,..., n c (11) Luego de encontrar la matriz de whitening W a partir de la matriz de covarianza compuesta R, y las descomposiciones propias de las matrices de covarianza de cada clase (valores propios y vectores propios U 0 ), se calculan los filtros espaciales SF c como: SF c = U c T W (12) Donde, W = 1/2 U 0 T, R = C i=1 R c con descomposición propia R = U 0 U 0 T, y U c es la matriz formada por los vectores propios de WR c W T. Principal Component Analysis (PCA), Análisis de Componentes Principales El análisis de componentes es un enfoque no supervisado para encontrar características que sean estadísticamente más apropiadas. El análisis de componentes principales (PCA) es un procedimiento matemático que usa una transformación ortogonal para convertir un conjunto de variables correlacionadas a un conjunto de variables no correlacionadas llamadas componentes principales. Esta transformación de los datos a un nuevo sistema de coordenadas contiene, al proyectar cada uno de los datos a los componentes principales, las varianzas más grandes. El número de los componentes principales es menor o igual al número de variables originales, haciendo posible reducir la dimensión del espacio [27], [28]. En reconocimiento de patrones, los datos se representan con una matriz X en donde cada columna es un ensayo. La tarea consiste en encontrar una matriz ortonormal P tal que Y = PX y la matriz de covarianza de Y sea diagonal. Las filas de P corresponden a los componentes principales, es decir, una base alternativa para los datos y Y son los datos expresados en términos de la base alternativa. Si la matriz de covarianza de Y es diagonal, significa que no existe redundancia entre las diferentes dimensiones de Y. Convenientemente, la cantidad de varianza mostrada por el i-ésimo componente principal es igual a la i-ésima entrada diagonal de la matriz de covarianza de Y. Ordenando los componentes principales por la varianza, se obtiene una forma de clasificarlos con respecto a su importancia [29]. Independent Component Analysis (ICA), Análisis de Componentes Independientes El método de ICA es útil cuando se desea estimar un vector de variables aleatorias latentes S (que no pueden ser observadas directamente) a partir de un vector de variables aleatorias disponibles X, que se pueden expresar por medio de una transformación lineal dada por: X = AS (13) donde A es una matriz de mezcla cuadrada y los componentes de S se asumen no gaussianos. El objetivo es estimar una matriz W = A 1 cuyos componentes sean maximamente independientes, tal que [30]: S = WX (14) 14

15 Cuando se aplica ICA en análisis de señales cerebrales, se asume que los componentes relacionados con la actividad son independientes de los componentes no relacionados y están mezclados linealmente para formar las señales observadas. El número de componentes es menor o igual al número de canales. Si existe alguna correlación entre los datos, debe realizarse preprocesamiento llamado whitening para remover esta correlación (transformando los componentes linealmente para alcanzar varianza unitaria) antes de aplicar ICA [31]. Definición Dos variables aleatorias y 1 y y 2 son independientes si la información en los valores de y 1 no brinda información alguna sobre los valores de y 2 y viceversa. Técnicamente la independencia puede ser definida por las densidades de probabilidad. Sea p(y 1, y 2 ) la función de densidad de probabilidad conjunta de y 1 y y 2, luego la función de probabilidad de densidad marginal de y 1 puede ser considerada como [30]: p 1 (y 1 ) = p(y 1, y 2 )dy 2 (15) Igualmente para y 2. Entonces se define que y 1 y y 2 son independientes si y solo si la función de densidad de probabilidad conjunta es factorizable de la forma: p(y 1, y 2 ) = p(y 1 )p(y 2 ) (16) Esta definición naturalmente se extiende para cualquier número de n variables aleatorias, en cuyo caso la función de densidad de probabilidad conjunta será la multiplicación de n términos. Por otro lado, se puede derivar una propiedad importante de las variables aleatorias independientes. Dadas 2 funciones, h 1 y h 2 siempre se tendrá: E{h 1 (y 1 )h 2 (y 2 )} = E{h 1 (y 1 )}E{h 2 (y 2 )} (17) La restricción de la no gaussianidad en los componentes de la matriz S se debe a que una función de densidad conjunta de vectores con distribuciones gaussianas no ofrece información sobre las direcciones de la matriz de mezcla A RECONOCIMIENTO ESTADÍSTICO DE PATRONES Espacio de características El espacio de características es un espacio en donde se encuentran los vectores que representan a cada uno de los ensayos de un problema dado. La calidad de un espacio de características radica en la separabilidad entre sus clases, y no en su dimensión Métodos paramétricos En los métodos paramétricos se asume que la forma de las funciones de densidad son conocidas [27] Métodos no paramétricos En la mayoría de aplicaciones de reconocimiento de patrones no se puede asumir que se conoce la forma de las funciones de densidad, ya que las formas paramétricas comunes no siempre se acomodan a las densidades encontradas en la práctica. Por lo tanto, resulta necesario utilizar técnicas de extracción de características que permitan describir los datos de forma compacta [27]. 15

16 Aprendizaje no supervisado El aprendizaje no supervisado busca extraer información de muestras no etiquetadas. Si la distribución proviene de una mezcla de densidades de probabilidad descritas por un conjunto de parámetros desconocidos Θ, estos parámetros pueden ser estimados por medio de métodos Bayesianos o de máxima verosimilitud. Un enfoque más general se tiene al definir alguna medida de similitud entre dos clusters, así como un criterio global tal como suma de error cuadrático o traza de una matriz de dispersión [27] Aprendizaje supervisado En el aprendizaje supervisado las muestras de entrenamiento utilizadas para diseñar el clasificador están etiquetadas con respecto a la clase a la cual pertenecen. Linear Discriminant Analysis (LDA), Análisis Lineal Discriminante Una función discriminante que es una combinación lineal de los componentes de x = {x 1, x 2,..., x n } puede ser escrita como [27]: g(x) = w T x + w 0 (18) donde w = {w 1, w 2,..., w n } es el vector de pesos y w 0 es el offset. La función discriminante g(x) da una medida algebraica de la distancia entre x y el hiperplano. Un clasificador lineal de dos categorías implementa la siguiente regla de decisión: Decide w 1 si g(x) > 0 y w 2 si g(x) < 0. Luego, x es asignado a w 1 si el producto interno w T x excede el umbral w 0 y w 2 en caso contrario. Si g(x) = 0, x puede ser asignado a cualquier clase. La ecuación g(x) = 0 define la superficie de decisión que separa los puntos asignados a w 1 de los puntos asignados a w 2. Cuando g(x) es lineal, la superficie de decisión es un hiperplano. En general, el hiperplano H divide el espacio de características en dos semiespacios: la región de decisión R 1 para w 1 y la región R 2 para w 2. Resumiendo, una función lineal discriminante divide el espacio de características por medio de una superficie de decisión que en este caso es un hiperplano. La orientación de la superficie es determinada por el vector normal w y su localización es determinada por el bias w 0. Supóngase que se tiene un conjunto de n muestras y 1, y 2,..., y n, algunas etiquetadas con w 1 y otras con w 2. Se debe determinar los pesos a en una función discriminante lineal g(x) = a T y a partir de estas muestras, es decir, se debe buscar un vector de pesos que clasifique todas las muestras correctamente. Si tal vector de pesos existe, se dice que las muestras son linealmente separables. Una muestra y i es clasificada correctamente si a T y i > 0 y y i está etiquetada con ω 1, o si a T y i < 0 y y i está etiquetada con w 2. Esto sugiere una normalización que simplifica el tratamiento para el caso de dos clases al reemplazar todas la muestras etiquetadas con ω 2 por sus negativos. De esta manera, se pueden olvidar las etiquetas y buscar un vector de pesos a tal que a T y i > 0 para todas las muestras. Tal vector de pesos es llamado vector de separación o vector solución. Support Vector Machines (SVM), Máquinas de Vectores de Soporte Las máquinas de vectores de soporte se encargan de preprocesar los datos para representar los patrones en un espacio de mayor dimensión con respecto al espacio original de caracterísiticas. Con un mapeo no lineal 16

17 apropiado ϕ(.) hacia una dimensión lo suficientemente grande, los datos provenientes de dos clases ω 1 y ω 2 pueden ser separados por un hiperplano. Se asume que cada patrón x k ha sido transformado a y k = ϕ(x k ). Para cada uno de los n patrones, k = 1,..., n se tiene z k = ±1, dependiendo si k ω 1 ó k ω 2. Un discriminante lineal g(y) en el espacio de mayor dimensión y está dado por [27]: g(y) = a T y (19) En donde el vector de pesos y el vector de patrones se encuentran en el espacio y, con a 0 ω 0 y y 0 = 1, respectivamente. Luego, el hiperplano separador asegura: z k g(y k ) 1 k = 1,..., n (20) El objetivo de entrenar una máquina de vectores de soporte es encontrar el hiperplano separador con el margen más grande; esperando que entre más grande el margen, mayor la generalización del clasificador. La distancia del hiperplano a cualquier patrón y es g(y) / a. Asumiendo que existe un margen positivo b, se tiene: z k g(y k ) a b k = 1,..., n (21) Se debe encontrar el vector de pesos a que maximiza b. Como el vector solución puede ser escalado arbitrariamente preservando el hiperplano, se imponen las condiciones que b a = 1 y que a 2 sea máximo para asegurar unicidad. Los vectores de soporte representan los patrones para los cuales z k g(y k ) = 1. Estos corresponden a las muestras de entrenamiento que definen el hiperplano óptimo de separación y son los patrones más difíciles de clasificar. Figura 4: Entrenar una máquina de vectores de soporte consiste en encontrar el hiperplano óptimo, es decir, aquel con la máxima distancia desde los patrones de entrenamiento más cercanos. Los vectores de soporte corresponden a estos patrones más cercanos, a una distancia b del hiperplano. Se muestran tres vectores de soporte en la figura. Tomada de [27]. Artificial Neural Networks (ANN), Redes Neuronales Artificiales Una red neuronal artificial está compuesta por una capa de entrada, capas ocultas y una capa de salida, interconectadas por medio de pesos modificables y compuestas por neuronas. En reconocimiento de patrones, las neuronas de la capa de entrada reciben los componentes de un vector de características (que va a ser utilizado en entrenamiento o en clasificación) y las señales emitidas por las neuronas de la capa de salida 17

18 son funciones discriminantes usadas para clasificación [27]. Cada neurona j de la capa oculta computa su salida en función de una suma ponderada net j de sus entradas, que puede ser expresada como el producto punto de su vector de pesos ω j con la entrada x: net j = d x i ω ji + ω j0 = i=1 d x i ω ji w T j x (22) i=0 donde d es el tamaño de la capa de entrada. Luego, cada salida y j es encontrada al aplicar net j a la función no lineal de activación f(.): y j = f(net j ) (23) De forma similar, cada neurona k de la capa de salida computa su resultado a partir de su vector de pesos ω k y del vector y proveniente de la capa oculta como: net k = n H j=1 y i ω kj + ω k0 = n H j=0 y i ω kj = w k T y (24) donde n H denota el número de neuronas de la capa oculta. Luego, cada salida z k es encontrada al aplicar net k a la función no lineal de activación g(.): z k = g(net k ) (25) Finalmente, la desición de clasificación consiste en asignar el patrón de entrada a la clase k que obtuvo el mayor valor de z k. 18

19 3. ESPECIFICACIONES 3.1. Diagrama de Bloques El diagrama de bloques del sistema de reconocimiento de patrones implementado en este trabajo y su correspondiente explicación se muestra a continuación. Figura 5: Diagrama de bloques del sistema. 19

20 Señales: Corresponde al conjunto de señales EEG que serán pre-procesadas o que pasarán directamente al bloque de selección para el proceso de extracción de características. ICA: Pre-procesa el conjunto de señales con la técnica ICA. Selección 1: Escoge el conjunto de señales que será utilizado (pre-procesado o no) y los métodos de extracción de características. DWT: Método de extracción que calcula la potencia promedio y la energía de las señales de las descomposiciones DWT. CWT: Método de extracción que calcula el promedio de los coeficientes CWT sobre 8 ventanas de tiempo. CSP: Método de extracción que transforma los datos por medio de filtros espaciales. MEL: Método de extracción que calcula el máximo exponente de Lyapunov. Vector de características: Forma el vector de características agrupando las medidas extraídas por cada uno de los métodos. Este puede ser pre-procesado o pasar directamente al bloque de selección para el proceso de clasificación. PCA: Pre-procesa el vector de características con la técnica de PCA. Selección 2: Escoge el vector de características que será utilizado (pre-procesado o no) y los métodos de clasificación. SVM: Método de clasificación que crea, entrena y prueba la SVM a utilizar. ANN: Método de clasificación que crea, entrena y prueba la ANN a utilizar. LDA: Método de clasificación que implementa LDA. 20

21 3.2. MEL: Máximo Exponente de Lyapunov El diagrama de bloques de la implementación propuesta para encontrar el máximo exponente de Lyapunov y su correspondiente explicación se muestra a continuación. Figura 6: Diagrama de bloques del sistema. Señales: Corresponde a las señales a las que se les estimará el máximo exponente de Lyapunov. Retardo τ: Calcula el retardo de cada pareja ensayo-canal a partir de la función de autocorrelación de la señal. Dimensión de embebimiento: Calcula la dimensión de embebimiento utilizando la relación planteada por D. Kugiumtzis [33]. Máximo exponente de Lyapunov: Usando el retardo y la dimensión de embebimiento reconstruye el espacio de fase y estima el máximo exponente de Lyapunov para cada canal. 21

22 4. DESARROLLOS 4.1. Descripición de las señales Este proyecto se realizó con un conjunto de señales EEG descargadas del sitio web del grupo Berlin Brain-Computer Interface en Alemania, destinadas a la comunidad participante en la Competencia BCI 2003 (Data set III) [2]. En este se dispone de 140 ensayos de tres variables (canales C3, CZ, C4) con sus respectivas etiquetas, de los que se utilizaron 112 para entrenamiento y 28 para prueba, en 5 combinaciones de validación cruzada. La actividad cerebral de los voluntarios fue registrada durante 9 segundos a una frecuencia de muestreo de 128 Hz (1152 muestras). La Figura 7 muestra la configuración en tiempo de cada ensayo: dos segundos después de iniciar el registro se presentó un estímulo acústico para indicar el inicio del ensayo, y se mostró una cruz ( + ) durante 1 segundo; finalmente, entre el segundo 3 y el 9 se mostró una flecha hacia la izquierda (clase 1) o hacia la derecha (clase 2), mientras el usuario debía imaginar el movimiento de la mano correspondiente a la dirección indicada. Figura 7: Configuración en tiempo de cada ensayo. Tomada de [9] Métodos de extracción de características Discrete Wavelet Transform (DWT), Transformada Discreta Wavelet La DWT es una herramienta de gran utilidad para la caracterización espectral de una señal, pues al utilizar bancos de filtros y hacer descomposiciones en bandas de frecuencia ortogonales, se puede analizar el contenido de energía de una señal en cada una de ellas por separado. Figura 8: Descomposición Wavelet de 3 niveles. Tomada de Matlab Help La Figura 8 muestra el esquema de descomposición wavelet de 3 niveles utilizado en este trabajo, del cual resultan un vector ca3 con la señal de aproximación, y tres vectores cd1, cd2 y cd3 con las señales de detalle. 22

23 Para implementar este esquema de descomposición se utilizó la función wavedec de Matlab con 3 niveles de profundidad y la wavelet Daubechies10. Las bandas generadas por esta descomposición para la frecuencia de muestreo trabajada (128 Hz) se muestran en el Cuadro 1. Descomposición Rango de Frecuencia (Hz) Nivel D D D A Cuadro 1: Bandas generadas por la descomposición utilizada La función wavedec retorna un vector Coef con las señales descompuestas en cada una de las bandas, y un vector Lint con sus respectivas longitudes. [Coef, Lint] = wavedec(eegdata DW T (:, c), 3, db10 ) (26) Algunas medidas de la señal en cada sub-banda de frecuencia pueden encontrarse a partir de las señales resultantes de la descomposición wavelet [9]. Las medidas seleccionadas para incluir en el vector de características de este trabajo fueron la potencia promedio y la energía. Potencia promedio de las descomposiciones DWT La potencia promedio de cada descomposición en cada canal se calculó como la norma al cuadrado del vector que las contenía, dividida por la longitud de dicho vector. P x = 1 L (x i ) 2 = 1 L L (norm(x))2 (27) i=1 Teniendo en cuenta que las bandas de interés para el problema son µ (8-13 Hz) y β (18-25 Hz), y que estas se encuentran contenidas en las descomposiciones D3 y D2, respectivamente, las cantidades incluidas en el vector de características fueron: C3 µd3: Potencia promedio de la descomposición D3 del canal C3. C3 βd2: Potencia promedio de la descomposición D2 del canal C3. C4 µd3: Potencia promedio de la descomposición D3 del canal C4. C4 βd2: Potencia promedio de la descomposición D2 del canal C4. Energía de las descomposiciones DWT La energía de las descomposiciones de aproximación se encontró utilizando la función wenergy de Matlab. Los parámetros de la función wenergy son el vector de descomposiciones Coef y el vector de longitudes Lint que retorna la función wavedec. [Ea3, Eb3] = wenergy(coef C3, L C3) (28) [Ea4, Eb4] = wenergy(coef C4, L C4) (29) Siguiendo un razonamiento similar al caso de la potencia promedio,las cantidades incluidas en el vector de características fueron: C3 en: Energía de las descomposiciones de aproximacion del canal C3. C4 en: Energía de las descomposiciones de aproximacion del canal C4. 23

24 Continuous Wavelet Transform (CWT), Transformada Continua Wavelet La Transformada Continua Wavelet, al calcular la participación de distintas frecuencias en cada instante de tiempo de una señal dada, permite observar simultáneamente las características tiempo-frecuencia de la misma. La mencionada propiedad de multi-resolución se puede observar al graficar las matrices de coeficientes resultantes (e.g. Figura 9), en las que las escalas bajas (frecuencias altas) tienen mayor resolución temporal y las escalas altas mayor resolución frecuencial. Figura 9: Transformada Continua Wavelet de un ensayo Para implementar la CWT se utilizó la función cwt de Matlab, que genera una matriz con los coeficientes CWT, la cual se puede graficar en tres dimensiones para visualizar una medida proporcional a la energía de una señal en cada pareja tiempo-escala del rango especificado. A diferencia de la función wavedec, la función cwt permite observar simultáneamente la información de tiempo y escala (frecuencia), por lo que es posible determinar visualmente las regiones del plano (u, t) donde se encuentra la mayor cantidad de energía de una señal dada. Las Figuras 10 y 11 muestran la CWT para uno de los ensayos con etiqueta 1 (imaginación mano izquierda). Figura 10: CWT del canal C3, Ensayo 81 24

25 Figura 11: CWT del canal C4, Ensayo 81 Luego de probar esta función con distintos intervalos de escalas s, se determinó que el intervalo [1, 150] era suficiente para analizar el contenido espectral de la señal, y que efectivamente existía información discriminante en la CWT. Suma de los coeficientes CWT sobre regiones muestra/escala Como una medida de la cantidad de energía presente en una región tiempo-escala, se sumaron los coeficientes CWT de cada canal y cada ensayo desde 3 hasta 9 segundos y desde la escala 50 a la 150. Teniendo en cuenta la etiqueta de cada ensayo, se comparó la cantidad de cada canal con el canal contrario, esperando que fuera mayor la del canal contrario al lado de la imaginación del movimiento. Promedio de los coeficientes sobre 8 ventanas de tiempo Publicaciones como [7] y [32] proponen el cálculo de cantidades sobre varias ventanas de tiempo en vez de una sola ventana para todo el ensayo. En particular, [32] propone utilizar el promedio de los coeficientes de la CWT, y sostiene que el número óptimo de ventanas para una SVM es 8. En este caso, las medidas escogidas para incluir en el vector de características fueron los máximos en escalas s sobre las medias del valor absoluto de las muestras en tiempo t j de las ventanas W j para cada canal Ci. C3 W3: Máximo de la ventana 3 del canal C3. C3 W4: Máximo de la ventana 4 del canal C3. C3 W5: Máximo de la ventana 5 del canal C3. C3 W6: Máximo de la ventana 6 del canal C3. C4 W3: Máximo de la ventana 3 del canal C4. C4 W4: Máximo de la ventana 4 del canal C4. C4 W5: Máximo de la ventana 5 del canal C4. C4 W6: Máximo de la ventana 6 del canal C4. donde Ci W j = max s (mean t ( Coefficients CW T trial(s, t j, C i ) )), t j W j (30) 25

26 Common Spatial Patterns (CSP), Patrones Espaciales Comunes Para implementar la técnica de CSP, se creó una función que construye los filtros espaciales a partir de las muestras de entrenamiento EEGdataCSP y las etiquetas EEGlabels: [SpatialF ilter 1, SpatialF ilter 2] = Z F unction F ilters CSP (EEGdata CSP, EEGlabels) (31) Luego se creó una función que aplica las transformaciones a los datos originales X por medio de los filtros espaciales SF c (X c = SF c X) y se construye el vector de características a partir de la varianza de los canales. F eatures CSP trial = Z F unction F eatures CSP (SpatialF ilter 1, SpatialF ilter 2, EEGdata CSP (:, :, t)) (32) donde las características para cada filtro c y canal i del ensayo X son encontradas como: ( ) var(x ci F eaturescsp trial(x ci ) = log ) K k=1 var(x ck ) (33) Finalmente, las medidas incluidas en el vector de características fueron: C3 SF1: Varianza del canal C3 transformado por el filtro espacial 1. C4 SF2: Varianza del canal C4 transformado por el filtro espacial Máximo exponente de Lyapunov Dado que las señales trabajadas provienen de un sistema inherentemente caótico (cerebro), estas deben ser analizadas con herramientas que permitan cuantificar dicho comportamiento. Una de estas es el máximo exponente de Lyapunov. El primer paso para hallar el máximo exponente de Lyapunov es encontrar el retardo y la dimensión de embebimiento para reconstruir el espacio de fase [22]: X(t) = {x(t), x(t + τ),..., x[t + (m 1)τ]} T (34) Para encontrar la matriz de retardos τ, se creó una función que calcula el retardo de cada ensayo-canal a partir de la función de autocorrelación de la señal, asignando a τ el índice del primer cruce por cero. Para calcular la dimensión de embebimiento se utilizó la relación planteada por D. Kugiumtzis [33] τ w = (m 1)τ (35) donde τ w es igual a la media de tiempo entre picos (mtbp) [34], que se encontró como el inverso de la frecuencia fundamental de la serie de tiempo. Finalmente, haciendo uso de la función lyapvec tomada de [21], modificada para suministrar la matriz de retardos y la dimensión de embebimiento, se encontró el máximo exponente de Lyapunov para cada canal, y se incluyeron en el vector de características: C3 MaxExp: Máximo exponente de Lyapunov del canal C3. C4 MaxExp: Máximo exponente de Lyapunov del canal C4. 26

27 4.3. Definición del vector de características Tomando las cantidades presentadas en la sección 4.2, se procedió a construir el vector de características como se muestra en el Cuadro 2. Posición Canal Medida Método 1 C3 W3 2 C3 W4 3 C4 W3 Medias de los coeficientes 4 C4 W4 CWT sobre 8 ventanas 5 C3 W5 de tiempo w i 6 C3 W6 7 C4 W5 8 C4 W6 9 C3 en Energía de los 10 C4 en coeficientes DWT 11 C3 µd3 12 C3 βd2 Potencia promedio de 13 C4 µd3 los coeficientes DWT 14 C4 βd2 15 C3 SF1 Varianzas 16 C4 SF2 por CSP 17 C3 MaxExp Lyapunov 18 C4 MaxExp Cuadro 2: Vector de características 4.4. PCA e ICA Independent Component Analysis (ICA), Análisis de Componentes Independientes: Algoritmo FastICA La regla de aprendizaje en FastICA encuentra una dirección, en este caso un vector unidad w cuya proyección w T x maximiza la no gaussianidad de las variables latentes s i. La no gaussianidad puede ser medida por la aproximación de la negentropía J(w T x) o por derivaciones iterativas de aproximaciones de Newton para una función g dada. La forma básica del algoritmo es: 1. Escoger un vector inicial de pesos w. 2. Hacer w + = E{xg(w T x)} E{g (w T x)}w. 3. Hacer w = w+ w Repetir desde 2 hasta que la diferencia entre los antiguos y nuevos valores en un punto de w converja [30]. Para aplicar este algoritmo se usó la función fastica de Matlab (versión 2.5). Para el preprocesamiento de las señales cerebrales se creó la función [ICA data] =Z F unction P reprocessing ICA(EEGdata ICApre) (36) FastICA recibe la matriz de datos cuyas filas corresponden a los ensayos y retorna una matriz en cuyas filas se encuentran los componentes independientes. Los dos parámetros más importantes son el número máximo de iteraciones y la función de aproximación g (cúbica, tangente hiperbólica, gaussiana y cuadrática). 27

28 Principal Component Analysis (PCA): Análisis de Componentes Principales Para implementar la técnica de preprocesamiento del espacio de características PCA se utilizó la función princomp de Matlab. Princomp recibe la matriz de características con los ensayos en las filas y retorna una matriz con los componentes principales. [Coef f F eatures All] = princomp(eegdata P CApre F eatures) (37) 4.5. Métodos de Clasificación Linear Discriminant Analysis (LDA), Análisis Lineal Discriminante Para implementar el clasificador LDA, se separaron tanto las muestras como las etiquetas de entrenamiento de las de prueba, y se utilizaron las primeras para encontrar los parámetros del clasificador conforme muestran las siguientes ecuaciones: Primero se encontró la media general de los datos a partir de las medias de las clases: µ = N izqµ izq + N der µ der N izq + N der (38) donde µ izq es la media de los vectores del espacio de características pertenecientes a la clase 1, µ der la media de los pertenecientes a la clase 2, y N izq y N der el número de vectores de cada clase. Luego se calculó la matriz de covarianza C de los datos utilizando la función cov de Matlab, y se procedió a calcular el vector de pesos w y el bias w o como: w =C 1 (µ der µ izq ) T (39) w 0 = µw (40) Después, se creó la función discriminante a partir de las cantidades anteriores para generar el hiperplano de separación en el espacio de características: D(x) =w T x + w 0 (41) Finalmente, se utilizó dicha función discriminante con las muestras de prueba para formar el vector de salidas del clasificador LDA, y luego compararlo con las etiquetas de prueba y evaluar el desempeño del mismo Support Vector Machines (SVM), Máquinas de Vectores de Soporte Para crear y entrenar la SVM se utilizó la función svmtrain de Matlab, a la que se le proporcionaron la matriz de características (ensayos por características) y las etiquetas de entrenamiento. SV Mgraz =svmtrain(data T rain SV M, Labels T rain SV M, KernelF unction,...); (42) Para probar el desempeño del clasificador se utilizó la función svmclassify de Matlab, ingresando como parámetros la SVM creada anteriormente y los datos de prueba. SV Moutputlabels = svmclassify(sv Mgraz, DataT est SV M); (43) 28

29 Artificial Neural Networks (ANN), Redes Neuronales Artificiales Para crear la red neuronal feed-forward back-propagation se utilizó la función newff de Matlab, ingresando como parámetros el conjunto completo de características junto a sus etiquetas. ANNfunction = newff(minmax(f eatures All train ), [3 1], { purelin, purelin }); (44) Esta red se entrenó utilizando la función train de Matlab limitando el número de épocas a 100, y se simuló utilizando la función sim de Matlab, ingresando como parámetros la red neuronal creada y los datos de prueba. ANNfunction = train(annfunction, F eatures All train, LabelsT rain ); (45) ANNoutputs = sim(annfunction, F eatures All test ); (46) Debido a que el resultado de una red neuronal puede variar para diferentes ejecuciones de un mismo conjunto de características, se optó por encontrar el resultado final como el promedio de varias ejecuciones. 29

30 5. ANÁLISIS DE RESULTADOS Para analizar la calidad del espacio de características construido y el desempeño de los clasificadores propuestos se abordaron los siguientes enfoques: 1. Combinación de técnicas de extracción: Se probaron todas las combinaciones posibles del vector de características agrupado por técnicas de extracción, para un total de 15 combinaciones. 2. K-folds: Para estudiar la generalización de los clasificadores propuestos, se implementó una validación cruzada de orden 5, en donde todos los subconjuntos tienen el mismo tamaño. Durante el desarrollo del proyecto se probaron distintas combinaciones de los parámetros de los clasificadores con el fin de encontrar los que mejores desempeños presentaran para este problema. Para SVM se probaron el kernel lineal, el polinomial y el RBF; para ANN se probaron arquitecturas de 3, 6 y 12 neuronas en la capa oculta y funciones de activación purelin, logsig y tansig. Los siguientes cuadros muestran los porcentajes de clasificación para el subconjunto de prueba 3 y las características DWT con diferentes combinaciones de parámetros. Kernel Parámetro Valor lineal 92,86 Orden=2 78,57 Polinomial Orden=3 64,29 Orden=4 60,71 σ=0.1 46,43 RBF σ= 1 82,14 σ=10 85,71 Cuadro 3: Parámetros SVM Neurona CO tansig-purelin logsig - purelin purelin - purelin Cuadro 4: Parámetros ANN Basados en los resultados de los cuadros 3 y 4, se escogió el kernel lineal para SVM y una arquitectura de 1 capa oculta con 3 neuronas para ANN, con funciones de activación purelin tanto para la capa oculta como para la de salida. Los mejores resultados obtenidos se resumen en el cuadro 5: 30

31 Grupos de Resultados Clasificación ( %) SVM LDA ANN Promedio Cuadro 5: Mejores Resultados En la siguiente gráfica de barras se muestran los mejores resultados obtenidos para cada clasificador. Figura 12: Mejores Resultados con PCA 5.1. Cuadros de resultados Cada una de los cuadros de resultados se estructura de la siguiente manera: En la columna de características se indican con una x las técnicas de extracción que participan en cada clasificación. En la columna de PCA se indica con una x si se realizó preprocesamiento con esta técnica. En la columna de partición de prueba se muestran los resultados de clasificación para cada subconjunto de prueba. En la última columna se muestra el promedio de los resultados obtenidos para cada combinación de características. A continuación, se muestra un cuadro de resultados para cada uno de los clasificadores propuestos. 31

32 Resultados SVM: Características Subconjuntos de Prueba DWT CWT CSP Lyapunov PCA Promedio x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x Promedio x Cuadro 6: Resultados SVM El máximo porcentaje de clasificación obtenido por SVM fue 92,86 % para el subconjunto de prueba 3, utilizando las características de DWT y la combinación DWT-Lyapunov. El mejor porcentaje promedio de clasificación fue 75,002 %, utilizando las características de DWT y aplicando PCA para su preprocesamiento. 32

33 Resultados LDA: Características Subconjuntos de Prueba DWT CWT CSP Lyapunov PCA Promedio x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x Promedio x Cuadro 7: Resultados LDA El máximo porcentaje de clasificación obtenido por LDA fue 92,86 % para el subconjunto de prueba 3, utilizando las características de DWT y la combinación DWT-Lyapunov. El mejor porcentaje promedio de clasificación fue 72,858 %, utilizando las características de DWT y la combinación DWT-CWT. 33

34 Resultados ANN: Características Subconjuntos de Prueba DWT CWT CSP Lyapunov PCA Promedio x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x Promedio x Cuadro 8: Resultados ANN El máximo porcentaje de clasificación fue obtenido por ANN 92,86 % para el subconjunto de prueba 3, utilizando las características de DWT. El mejor porcentaje promedio de clasificación fue 73,68 %, utilizando las características de DWT y aplicando PCA para su preprocesamiento. 34

35 5.2. Discusión Para estudiar la generalización de los clasificadores propuestos con respecto al número de subconjuntos, se realizó el mismo procedimiento de clasificación con una validación cruzada de orden 4 con subconjuntos de igual tamaño. Los mejores grupos de resultados para este caso se muestran en el Cuadro 9. Grupos de Resultados Clasificación ( %) SVM LDA ANN Promedio Cuadro 9: Mejores Resultados 4 Subconjuntos La siguiente gráfica de barras muestra los resultados anteriores obtenidos para cada clasificador. Figura 13: Mejores Resultados 4 Subconjuntos Generalización Los clasificadores propuestos no presentaron una buena generalización para los distintos subconjuntos de prueba. La única combinación que presentó un buen promedio y poca dispersión fue la de DWT- CWT con SVM y ANN. Esto se debe principalmente a que en la mayoría de los casos el peor porcentaje de clasificación se presenta en el segundo subconjunto de prueba, sesgando el promedio para los tres clasificadores. El mismo fenómeno se pudo evidenciar con el primer subconjunto al hacer validación cruzada de orden 4, lo cual podría significar que las condiciones experimentales de algunos de estos ensayos fueron bastante distintas a las de los demás. Preprocesamiento A pesar de que la diferencia no fue tan significativa al aplicar preprocesamiento con PCA a las características, se pudo observar que en la mayoría de los subconjuntos de prueba se igualó ó se obtuvo un aumento en el porcentaje de clasificación. 35

Sitemap