INTERFAZ CEREBRO COMPUTADOR MEDIANTE LA CLASIFICACIÓN DE SEÑALES ELECTROENCEFALOGRÁFICAS.

ANGÉLICA REYES RUEDA

INGENIERÍA ELECTRÓNICA FACULTAD DE INGENIERÍA PONTIFICIA UNIVERSIDAD JAVERIANA BOGOTÁ D.C. COLOMBIA JULIO 2013

INTERFAZ CEREBRO COMPUTADOR MEDIANTE LA CLASIFICACIÓN DE SEÑALES ELECTROENCEFALOGRÁFICAS. T.G 1174.

AUTOR: ANGÉLICA REYES RUEDA

Trabajo de grado presentado como requisito parcial para optar al título de INGENIERO ELECTRÓNICO

DIRECTORES: Ing. CESAR LEONARDO NIÑO BARRERA

INGENIERÍA ELECTRÓNICA FACULTAD DE INGENIERÍA PONTIFICIA UNIVERSIDAD JAVERIANA BOGOTÁ D.C. COLOMBIA JULIO 2013 2

PONTIFICIA UNIVERSIDAD JAVERIANA FACULTAD DE INGENIERÍA

CARRERA DE INGENIERÍA ELECTRÓNICA

RECTOR MAGNÍFICO: Joaquín Emilio Sánchez García, S.J.

DECANO ACADÉMICO: Ing. Jorge Luis Sánchez Tellez.

DECANO DEL MEDIO UNIVERSITARIO: P. Antonio José Nova, S.J.

DIRECTOR DEL DEPARTAMENTO DE ELECTRÓNICA: Ing. Francisco Viveros Moreno.

DIRECTOR DE CARRERA: Ing. Jairo Alberto Hurtado Londoño

DIRECTORES DEL PROYECTO DE GRADO: Ing. Cesar Leonardo Niño Barrera

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 qué 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”

4

A Gilberto, María de los Ángeles, Andrés Felipe y William por su incondicional apoyo perfectamente mantenido a través del tiempo. Todo este trabajo ha sido posible gracias a ustedes.

5

AGRADECIMIENTOS

Quiero agradecer sinceramente a todas aquellas personas que compartieron sus conocimientos y tiempo conmigo para hacer posible la conclusión del presente trabajo. Especialmente agradezco a mi director el Ing. Cesar Leonardo Niño Barrera por su asesoría siempre dispuesta aún en la distancia. Gracias a María de los Ángeles, Gilberto, Andrés Felipe, William, Sergio, Patricia, Mario, Alexandra y Paola Andrea por su valiosa colaboración y apoyo al ser partícipes de este proyecto. A mis amigos y docentes de la Pontificia Universidad Javeriana y Universidad Pontificia Bolivariana por los momentos compartidos y enseñanzas.

Gracias a todos ellos.

6

TABLA DE CONTENIDO 1.

INTRODUCCIÓN ......................................................................................................................... 12

2.

MARCO TEORICO ....................................................................................................................... 13 2.1

ELECTROENCEFALOGRAFIA (EEG) ................................................................................ 13

2.2

TÉCNICAS DE CAPTURA DE SEÑALES NEUROELÉCTRICAS ...................................... 13

2.3

SISTEMA INTERNACIONAL 10 – 20 .................................................................................. 14

2.4

REGISTRO MONOPOLAR Y BIPOLAR DE EEG ............................................................... 15

2.5

ONDAS ENCEFALOGRAFICAS .......................................................................................... 16

2.6

NEUROFISIOLOGÍA DE TAREAS MOTORAS................................................................... 17

2.6.1

Corteza Motora ............................................................................................................... 17

2.6.2

Desincronización/Sincronización Relacionada a Eventos (ERD/ERS) ............................. 18

2.7 2.7.1

Actividades electrofisiológicas cerebral ........................................................................... 19

2.7.2

Estructura BCI ................................................................................................................ 20

2.8

3.

4.

INTERFAZ CEREBRO - COMPUTADOR (BCI) .................................................................. 19

EMOTIVSYSTEMS ............................................................................................................... 22

2.8.1

Headset ........................................................................................................................... 22

2.8.2

Research Edition (SDK) .................................................................................................. 22

ESPECIFICACIONES ................................................................................................................... 24 3.1

CAPTURA DE SEÑALES ENCEFALOGRÁFICAS ............................................................. 24

3.2

INTERFACE CEREBRO COMPUTADOR (BCI).................................................................. 25

3.3

PRUEBAS DE DESEMPEÑO................................................................................................ 25

3.4

HAWARE Y SOFTWARE ..................................................................................................... 25

3.5

BANCO DE SEÑALES EEG ................................................................................................. 26

DESARROLLO ............................................................................................................................. 27 4.1

CAPTURA DE SEÑALES EEG ............................................................................................. 27

4.1.1

Formulación de experimento (paradigma) ....................................................................... 27

4.1.2

Diseño de interfaz gráfica................................................................................................ 28

4.1.3

Captura de señales EEG .................................................................................................. 30

4.2

PROCESAMIENTO DE SEÑALES PARA EL BCI ............................................................... 32

4.2.1

Selección de Canales....................................................................................................... 32

4.2.2

Eliminación de DC – Offset ............................................................................................ 32

4.2.3

Ventana Hanning ............................................................................................................ 33

4.2.4

Filtrado ........................................................................................................................... 34

4.3

Filtros Espaciales .................................................................................................................... 35 7

4.4 4.4.1

Densidad del Espectro de Potencia (Power Spectral Density – PSD)................................ 38

4.4.2

Valor de Bloqueo de Fase (Phase-Locking Value) ........................................................... 41

4.5 4.5.1 5.

EXTRACCIÓN DE CARACTERÍSTICAS ............................................................................ 38

CLASIFICACIÓN DE CARACTERISTICAS ........................................................................ 43 Máquinas de Soporte Vectorial (SVM) ............................................................................ 44

ANALISIS DE RESULTADOS ..................................................................................................... 50 5.1

CAPTURA DE SEÑALES ..................................................................................................... 50

5.2

EXTRACCIÓN Y CLASIFICACIÓN PSD - SVM ................................................................. 50

5.3

EXTRACCIÓN Y CLASIFICACIÓN PLV – SVM ................................................................ 52

6.

CONCLUSIONES ......................................................................................................................... 54

7.

BIBILIOGRAFIA .......................................................................................................................... 55

8.

ANEXOS ....................................................................................................................................... 58

8

LISTADO DE FIGURAS Figura 1. Capas que cubren el cerebro [4]............................................................................................... 13 Figura 2. Sistema Internacional 10 – 20 [5]. ........................................................................................... 14 Figura 3. Relación de distancia del Sistema Internacional 10 – 20 [5]. (a) Vista de Frente. (b) Vista de Perfil...................................................................................................................................................... 14 Figura 4. Vista lateral de hemisferio cerebral indicando su división en lóbulos [4]. ................................. 15 Figura 5. Montaje Registro Monopolar [5] ............................................................................................. 15 Figura 6. Montaje Registro Bipolar. [5] .................................................................................................. 16 Figura 7. Divisiones de la corteza motora. [7]......................................................................................... 17 Figura 8. Representación del homúnculo motor sobre la corteza motora. [7] ........................................... 17 Figura 9. ERD / ERS de un sujeto el cual mueve lenta (izquierda) y rápidamente (derecha) un dedo. Se muestran tres bandas de frecuencia (10 - 12, 14 - 18 y 36 - 40 Hz) [9]. ................................................... 18 Figura 10. Corteza Cerebral Fuente: Editorial Médica Panamericana ...................................................... 19 Figura 11. Estructura Interface Cerebro Computador (BCI) .................................................................... 20 Figura 12. Señales obtenidas al adquirir señales de actividad cerebral [13]. ............................................ 20 Figura 13. Headset EMOTIV EPOC....................................................................................................... 22 Figura 14. Etapas de desarrollo .............................................................................................................. 24 Figura 15. Sensores EMOTIV EPOC según Sistema Internacional 10-20 [17]. ....................................... 24 Figura 16. Diagrama de bloques del BCI diseñado. ................................................................................ 25 Figura 17. Esquema de tiempo de una prueba. ........................................................................................ 27 Figura 18. Estimulo visual. (a) Cruz, inicio de la tarea. (b) Caso 1, acción motora derecha. (c) Caso 2, acción motora izquierda. ........................................................................................................................ 28 Figura 19. GUI - Ingreso de datos para la captura de señales. ................................................................. 28 Figura 20. Archivo .txt generado el cual indica la secuencia para el Sujeto 4 – Sesión 2. ........................ 29 Figura 21. Diagrama de flujo de la lógica implementada para la captura de las señales EEG. .................. 29 Figura 22. EmotivTestBench - Sujeto 04 ................................................................................................ 30 Figura 23. Realización del experimento de captura de datos. .................................................................. 31 Figura 24. Señales de los 14 canales originales – AM: Derecha. ............................................................. 32 Figura 25. Señales de los 14 canales sin DC. .......................................................................................... 33 Figura 26. Diagrama de tiempo y ventana hanning ................................................................................. 33 Figura 27. (a) Señal EEG sin offset (b) Ventana Hanning (c) Señal Ventaneada ..................................... 34 Figura 28. Filtro Equiripple (a) Magnitud (b) Fase ................................................................................. 34 Figura 29. a) Espectro de Frecuencia Señal Original b) Espectro de Frecuencia Señal Filtrada ................ 35 Figura 30. Filtros Espaciales [26]. (a) Filtro Bipolar, (b) Filtro Laplaciano, (c) Referencia Media Comun .............................................................................................................................................................. 36 Figura 31. Aplicación de filtro CAR. (a) Señal Original (Filtrada) (b) Señal filtrada mediante CAR ..... 37 Figura 32. Señal FC5 (Azul – Señal original Filtrada, Roja – Señal Filtrada por CAR) .......................... 37 Figura 33. Diagrama de Flujo - Algoritmo PSD...................................................................................... 39 Figura 34. Matriz PSD - AM: Derecha ................................................................................................... 40 Figura 35. Matriz PSD - AM: Izquierda ................................................................................................. 40 Figura 36. Diagrama de Flujo - Algoritmo PLV ..................................................................................... 42 Figura 37. a) Par de electrodos AF3 y AF4 b) PLV – Pareja de Electrodos AF3 – AF4. Sujeto 4 – Sesión 3 .............................................................................................................................................................. 43 Figura 38. Búsqueda por Red - Escala Grande (Característica: PLV - ID: 04 – Sesión de Ensayo: S3) ... 46 Figura 39. Búsqueda por Red - Escala Mediana (Característica: PLV - ID: 04 – Sesión de Ensayo: S3) . 47 9

Figura 40. Búsqueda por Red - Escala Pequeña (Característica: PLV - ID: 04 – Sesión de Ensayo: S3).. 48 Figura 41.Señal EEG correspondiente a la tarea motora Derecha, Sujeto 4 – Sesión 3 – Muestra 4. (a) Señal Original, (b) Señal sin DC, (c) Señal Pre-Procesada ...................................................................... 50 Figura 42. PSD estimado mediante Welch en la banda de frecuencia de interés. (a) PSD – AM: Derecha, (b) PSD – AM: Izquierda. ...................................................................................................................... 51 Figura 43. Precisión promedio de clasificación para todos los sujetos. .................................................... 53 Figura 44. Ficha Médica ........................................................................................................................ 60 Figura 45. Fichero *.csv generado .......................................................................................................... 61 Figura 46. Kit Emotiv - Research SDK .................................................................................................. 63 Figura 47. Solución salina multipropósitos usada. .................................................................................. 64 Figura 48. Hidratación con solución salina en electrodos y almohadilla del estuche [41]. ........................ 64 Figura 49. Instalación electrodo en headset. ........................................................................................... 65 Figura 50. Pasos para colocar adecuadamente el headset EMOTIV......................................................... 65 Figura 51. Verificación del contacto de cada uno de los electrodos por medio de TestBench................... 66 Figura 52. Instalación LIBSVM - Paso 3 ................................................................................................ 67 Figura 53. Instalación LIBSVM - Paso 4 ................................................................................................ 67 Figura 54. Instalación LIBSVM - Paso 5 ................................................................................................ 68 Figura 55. Instalación LIBSVM - Paso 6 ................................................................................................ 68 Figura 56. Instalación LIBSVM - Paso 7 ................................................................................................ 68 Figura 57. Instalación LIBSVM - Paso 8 ................................................................................................ 69 Figura 58. Instalación LIBSVM - Paso 9 ................................................................................................ 69 Figura 59. Instalación LIBSVM - Paso 11 .............................................................................................. 70 Figura 60. Validación Cruzada de 4 iteraciones [15]. .............................................................................. 72

10

LISTADO DE TABLAS Tabla 1. Características EEG y ECoG. [4, 5] .......................................................................................... 13 Tabla 2. Características Ondas EEG [5,6] ............................................................................................... 16 Tabla 3. Artificios comunes [13]. ........................................................................................................... 21 Tabla 4. Distribución de sesiones para la captura de señales encefalográficas. ........................................ 25 Tabla 5. Especificaciones Laptop. .......................................................................................................... 26 Tabla 6. Especificaciones Desktop. ........................................................................................................ 26 Tabla 7. Tiempo aproximado de experimento por sujeto. ........................................................................ 31 Tabla 8. Sets de Entrenamiento y Ensayo mediante Hold-Out................................................................. 44 Tabla 9. Precisión PSD - Kernel Lineal .................................................................................................. 45 Tabla 10. Precisión PLV - Kernel LinealKernel RBF ............................................................................. 45 Tabla 11. Precisión PSD - Kernel RBF................................................................................................... 48 Tabla 12. Precisión PLV - Kernel RBF .................................................................................................. 49 Tabla 13. Desempeño promedio PSD + SVM ......................................................................................... 51 Tabla 14. Desempeño promedio PLV + SVM ........................................................................................ 52 Tabla 15. Nomenclatura archivos *.csv. ................................................................................................. 62

11

1.

INTRODUCCIÓN

La población mundial supera los seiscientos mil millones de personas, de ellas el 10% (600 millones) presenta una u otra forma de discapacidad, y más de dos terceras partes viven en países en vía de desarrollo [1]. En América Latina existen aproximadamente 85 millones de personas con discapacidad y por lo menos el 50% de esta población tienen edad para trabajar [2]. Según el Censo 2005 realizado por el DANE se estableció que el 6,4% de la población Colombiana tiene por lo menos una discapacidad; siendo como causas principales enfermedades adquiridas, lesiones causadas por accidentes de tráfico y/o laborales como por la violencia [2]. Es de destacar que en ese 6,4% de población (equivalente a 2.632.255 personas) el 29,3% presenta limitaciones permanentes para caminar o moverse [3], requiriendo en muchos casos la ayuda de un tercero para poder desplazarse. Las interfaces Cerebro-Computador (Brain-Computer Interface, BCI) permiten la conexión del sujeto y el computador mediante la captación de señales neuroeléctricas sin necesidad de existir actividad muscular por parte del sujeto. Debido a esto se puede implementar una BCI como sistema de control en tiempo real entre el sujeto que presenta discapacidad de movilidad y el mundo exterior [4]. La finalidad planteada en este trabajo fue el desarrollo de una interfaz Cerebro-Computador basado en clasificación de señales encefalogáficas cuando un sujeto imagina un movimiento, en este caso derecha o izquierda. Para lograr la clasificación se realizó un estudio exhaustivo sobre los aspectos neurofisiológicos de las tareas motoras y como estas pueden ser identificadas a través de señales obtenidas mediante electroencefalografía.

12

2.

MARCO TEORICO

2.1

ELECTROENCEFALOGRAFIA (EEG)

La electroencefalografía (EEG) es el registro de la actividad bioeléctrica cerebral obtenido mediante electrodos situados en el cuero cabelludo o corteza cerebral (espacio subdural, Figura 1) [4].

Figura 1. Capas que cubren el cerebro [4].

2.2

TÉCNICAS DE CAPTURA DE SEÑALES NEUROELÉCTRICAS

Existen diferentes técnicas para la captura de las señales neuroeléctricas, sin embargo el electroencefalograma (EEG) y elelectrocorticograma (ECoG) son los más usados (Tabla 1).

Técnica

Electroencefalograma (EEG)

Electrocorticograma (ECoG)

Tipo de Electrodos1

Ventajas

Superficie o Basales

Quirúrgicos

No invasiva

Resolución y rango de frecuencia mayor que el EEG.

Desventajas Resolución limitada Rango de Frecuencia limitado Sensible a ruido (ej. Señales electrooculográficas (EoG) ) Requiere procedimiento quirúrgico para la implantación de los electrodos (Invasivo). Manipulación exclusivamente por un neurocirujano.

Tabla 1. Características EEG y ECoG. [4, 5]

1

Electrodos superficiales: Se aplican sobre el cuero cabelludo. Electrodos basales: Se aplican en la base del cráneo sin necesidad de procedimiento quirúrgico. Electrodos quirúrgicos: para su aplicación es precisa la cirugía y pueden ser corticales o intracerebrales. [5]

13

2.3

SISTEMA INTERNACIONAL 10 – 20

Existen múltiples sistemas de posicionamiento de electrodos superficiales, sin embargo el sistema internacional 10 – 20 (Figura 2) es el más utilizado hoy en día. Este método se desarrolló con el fin de garantizar un estándar y repetitividad de manera tal que los estudios realizados a un sujeto puedan ser comparados en el transcurso del tiempo, igualmente permite comparar los resultados entre distintos individuos [6].

Figura 2. Sistema Internacional 10 – 20 [5].

Este sistema se basa en la relación de distancia entre dos electrodos adyacentes que es el 10% de la longitud del cráneo desde el nasión hasta el iniónó del 20% de la distancia entre ambos puntos pasando por el vertex (Figura 3).

Figura 3. Relación de distancia del Sistema Internacional 10 – 20 [5]. (a) Vista de Frente. (b) Vista de Perfil.

14

Es importante resaltar que la letra empleada para el nombre de cada uno de los puntos de contacto identifica al lóbulo (Figura 4) y el número a la ubicación dentro del hemisferio. Es decir las letras corresponderían de la siguiente forma: F: Frontal, T: Temporal, P: Parietal y O: Occipital; la letra C esta se implementa para identificar la línea horizontal central. En el caso de los números, los pares corresponden con electrodos en el hemisferio derecho y los impares con los del izquierdo; los subíndices z se identifican la línea vertical central de electrodos.

Figura 4. Vista lateral de hemisferio cerebral indicando su división en lóbulos [4].

2.4

REGISTRO MONOPOLAR Y BIPOLAR DE EEG

Para realizar la captura o registro de las ondas EEG múltiples electrodos son ubicados en el cuero cabelludo, sin embargo para realizar este registro es necesario disponer de mínimo dos terminales (electrodos). Existen dos tipos de montajes para realizar este registro, el monopolar y el bipolar. El Registro Monopolar (Figura 5) se toma la señal independiente de cada uno de los electrodos, estos electrodos son denominados activos o de registro. Existe un segundo electrodo llamado de referencia el cual teóricamente debe tener un voltaje igual a 0 V (GND). Para garantizar un valor aproximado a cero el electrodo de referencia generalmente se ubica en el lóbulo de la oreja, mentón o mastoides 2[5]. De esta forma se determina la diferencia de tensión entre el electrodo activo y el electrodo de referencia.

Figura 5. Montaje Registro Monopolar [5]

2

Porción posterior del hueso temporal detrás de la oreja.

15

El Registro Bipolar (Figura 6) a diferencia del monopolar se toma parejas de electrodos y se registra la diferencia de tensión entre ellos [5].

Figura 6. Montaje Registro Bipolar. [5]

2.5

ONDAS ENCEFALOGRAFICAS

Las señales EEG se caracterizan por tener amplitudes desde los 10mV (registros sobre la corteza cerebral) a 100μV (registros sobre el cuero cabelludo), igualmente tienen forma irregular y son aperiódicas y se encuentran concentradas en el espectro de 0 a 100 Hz [5,6]. Debido a la irregularidad de estas se han clasificado en ritmos u ondas dependiendo de la banda de frecuencia que se encuentre (Tabla 2).

Tipo

Frecuencia (Hz)

Localización

Alfa (α)

8 – 13

Occipital y Frontal

Beta (β)

12 – 30

Parietal y Frontal

Gamma(γ)

> 30

Delta (δ)

0.5 – 4

Tetha (θ)

4–7

Mu(μ)

8 - 13

Parietal y Temporal. Zona Sensoriomotora (unión lóbulo parietal y frontal)

Características Amplitud > 20µV Asociadas a relajación y actividad mental leve. Amplitud entre 5 – 30 µV Asociadas a estados de alerta, concentración y resolución de problemas. En casos de extrema concentración puede alcanzar hasta los 50 Hz. Asociadas a procesos de comunicación y realización de actividades mentales complejas. Amplitud variable. Asociadas a sueño profundo y estados de meditación. Amplitud >20µV Asociadas estrés y frustración Manifiestan atenuación durante el movimiento o intento de movimiento de extremidades. A diferencia de las ondas alfa (que se encuentran en la misma banda de frecuencia) las ondas mu son espontaneas.

Tabla 2. Características Ondas EEG [5,6]

16

2.6

NEUROFISIOLOGÍA DE TAREAS MOTORAS

2.6.1

Corteza Motora

La corteza motora cerebral es la encargada de realizar la planificación, control y ejecución de la mayoría de movimientos voluntarios.

Figura 7. Divisiones de la corteza motora. [7]

La corteza motora es dividida en cuatro partes principales (Figura 7) las cuales son: 2.6.1.1 Corteza motora primaria (M1) La corteza motora primaria (M1) se encuentra ubicada en el lóbulo frontal a lo largo de una protuberancia llamada circunvolución precentral (Figura 8), y es la responsable de la generación de los impulsos neuronales que controlan la ejecución del movimiento. Wilder Penfield, neurocirujano canadiense, realizó una serie de experimentos en los años 50 a múltiples pacientes que iban a ser sometidos a neurocirugía [8]. El experimento consistía en estimular eléctricamente diversas partes de la corteza cerebral (M1 y PMA) logrando movimientos involuntarios en diversas partes del cuerpo. Con los datos obtenidos resumió su descubrimiento en una representación gráfica la cual le dio el nombre de "homúnculo" (Figura 8).

Figura 8. Representación del homúnculo motor sobre la corteza motora. [7]

17

2.6.1.2 Corteza Pre-motora (PMA) La corteza pre-motora (PMA), se encuentra ubicada 1 a 3 cm en frente de la corteza motora primaria (M1). La PMA es la encargada de guiar y controlar los movimientos de los músculos proximales (hombros, cuello, pelvis) y del tronco corporal. 2.6.1.3 Área motora suplementaria (SMA) El área motora suplementaria (SMA) se encuentra por encima del área pre-motora, y delante de la corteza motora primaria (M1). Esta área es la encargada de la planificación y coordinación de movimientos complejos. El área motora suplementaria y las regiones promotoras tanto envían información a la corteza motora primaria, así como a las regiones motoras del tronco cerebral. 2.6.1.4 Corteza Parietal Posterior La corteza parietal posterior es la encargada de transformar la información visual en instrucciones motoras. 2.6.2

Desincronización/Sincronización Relacionada a Eventos (ERD/ERS)

Cuando una persona realiza un movimiento (o imaginación de la tarea motora) se induce una serie de fenómenos en la dinámica de la actividad mental (Figura 9), estos fenómenos son conocidos como desincronización / sincronización relacionada a eventos (ERD / ERS).

Figura 9. ERD / ERS de un sujeto el cual mueve lenta (izquierda) y rápidamente (derecha) un dedo. Se muestran tres bandas de frecuencia (10 - 12, 14 - 18 y 36 - 40 Hz) [9].

La ERD / ERS se compone de diversos eventos separados, producidos en diferentes bandas de frecuencia y duración variable. Sin embargo la característica más dominante [9] es una desincronización en la banda alfa superior (10 a 12 Hz, la cual corresponde a la banda sensoriomotora (también conocida como banda mu). Es importante resaltar que al imaginar la tarea motora, el fenómeno ERD / ERD también sucede, lo 18

cual es importante en el contexto de una interface cerebro computadora (BCI) ya que como su definición lo indica, el sistema no depende de movimiento de los músculos.

2.7

INTERFAZ CEREBRO - COMPUTADOR (BCI)

Una interfaz cerebro – computador (BCI, Brain-Computer Interface) es un sistema que permite la conexión entre un sujeto a partir de la actividad eléctrica cerebral (sin necesitad de actividad motora) y un computador. Este sistema detecta la presencia de patrones específicos en la actividad cerebral y mediante algoritmos de clasificación traduce estos patrones en comandos de control [10]. 2.7.1

Actividades electrofisiológicas cerebral

Un sistema BCI se define como un sistema que adquiere una señal electrofisiológicas cerebral la cual es analizada y procesada posteriormente. Es de resaltar que un BCI no depende de las vías de salida normal de los nervios y músculos. En la actualidad se destacan dos tipos de BCI basados en ondas electroencefalográficas (EEG – electrodos de contacto, no invasivo): Potenciales Visuales Evocados (Visual evoked potentials, VEPs) y Actividad Sensoriomotora (Sensorimotor Rhythm, SMR) [10]. 2.7.1.1 Potenciales visuales evocados (VEPs) Los potenciales evocados miden la actividad eléctrica del cerebro a estímulos externos como lo pueden ser visuales, auditivos o táctiles. Usualmente se implementan estímulos visuales como luces parpadeando a determinada frecuencia.

Figura 10. Corteza Cerebral Fuente: Editorial Médica Panamericana

Si un estímulo visual se presenta repetidamente a una velocidad de 5-6 Hz o mayor, se suele apreciar una respuesta oscilatoria continua en el lóbulo occipital donde se encuentra la corteza visual (figura 10). Esta respuesta se denomina estado estacionario de los potenciales visuales evocados (Steady-state visual evoked potentials – SSVEP) [10,11]. 2.7.1.2 Actividad sensoriomotora (SMR) Las ondas o ritmos beta (12 – 30 Hz) y mu (8 – 13 Hz) se originan en la corteza sensoriomotora (Figura 10) y están fuertemente relacionadas con el sistema motor. Cuando se realiza un movimiento se puede 19

apreciar en las señales electroencefalografías que la amplitud de estos ritmos disminuye o se dé sincroniza (ERD / ERS), es de resaltar que este fenómeno ocurre tanto al realizar o imaginar una acción motora. [4, 10] Es importante resaltar que los sistemas BCI no leen la mente del sujeto; es necesario que este controle conscientemente su actividad cerebral, ya sea por estímulos visuales evocados (VEPs) o imaginando que realizan una tarea motora, por ejemplo el mover el brazo izquierdo. 2.7.2

Estructura BCI

Un BCI, independientemente de su finalidad o métodos de grabación, se compone de cuatro elementos esenciales [12] como se muestra en la Figura 11. Brain Computer Interface ADQUISICIÓN DE SEÑALES

EXTRACCIÓN DE CARACTERÍSTICAS

CLASIFICACIÓN DE CARACTERÍSTICAS

RETROALIMENTACIÓN

APLICACIÓN DE CONTROL

SUJETO T

Figura 11. Estructura Interface Cerebro Computador (BCI)

2.7.2.1 Adquisición de señales y Artefactos La primera etapa de un BCI corresponde a la adquisición de señales eléctricas del cerebro, esta adquisición se realiza mediante ECoG y EEG. Sin embargo cualquiera de estos dos métodos registra no solamente actividad del cerebro sino otras actividades no deseadas (Figura 12) que deterioran la señal de la actividad cerebral (ruido o artefacto), de origen múltiple como lo indica la Tabla 3[13].

Figura 12. Señales obtenidas al adquirir señales de actividad cerebral [13].

20

Artificios No Fisiológicos      

Artificio Fisiológico     

50 / 60 ciclos Artificio de electrodo (resistencia) Descarga electroestática (movimientos) Descarga capacitivas (acoplamiento) Ruido electrónico interno(electronesamplificador) Artificios AD/DA (digitales)

Movimiento de ojos (parpadeo) Cardiaco Pulso Respiración Tensión

Tabla 3. Artificios comunes [13].

En los sistemas BCI se emplea una etapa previa a la extracción de características conocida como Etapa de Pre-Procesamiento [10]. Esta etapa tiene como objetivo el aumentar la relación señal-a-ruido (SNR), es decir adecua la señal para su posterior análisis. Las técnicas más empleadas para realizar un pre-procesamiento son Patrones Espaciales Comunes (Common Spatial Patterns – CSP) y Métodos de Referencia, como lo son Promedio de Referencia Común (Common Average Reference – CAR) y Laplaciano Superficial (SurfaceLaplacian – SL).

2.7.2.2 Extracción de características El proceso más importante en un BCI corresponde a la extracción de las características de las señales EEG. Este proceso consiste en procesan mediante diferentes métodos matemáticos las señales obtenidas para así extraer características representativas que permiten la posterior clasificación de la señal. El método a emplear para la extracción de características depende de las señales obtenidas, ya que las características obtenidas pueden basarse en dominio del tiempo (BCI basado en VEPs) o en el dominio de la frecuencia (BCI basado en SMR). Algunos métodos de extracción en el dominio del tiempo son: filtrado, transformada de Wavelet y modelos paramétricos auto regresivos (AR). En el caso de BCI basados en SMR se implementan métodos en el dominio de frecuencia, entre ellos la estimación de la banda de potencia usando el método de Welch, Wavelet y modelos auto regresivos (AR).

2.7.2.3 Clasificación de características Posterior a la extracción de características se requiere agruparlas con el fin de luego ser asociadas al paradigma (ej. Pie izquierdo vs. Mano Derecha). Gran parte de los métodos de clasificación son basados en algoritmos de aprendizaje [10] como lo son:    

Máquinas de vectores de soporte (Support Vector Machine – SVM) Análisis de discriminante lineal (Linear DiscriminantAnalysis – LDA) Discriminante de Fisher (Fisher’sDiscriminantAnalysis – FDA) Redesneuronales (Neural networks – NN)

21

2.7.2.4 Aplicación de control Posterior al diseño del BCI se encuentra esta etapa la cual cumple con el objetivo de crear la interacción con el usuario con algún otro dispositivo. Esta etapa implementa el análisis realizado previamente por el BCI para poder realizar una acción de control y de esta forma permitir la interacción del usuario con algún dispositivo, por ejemplo el manejo de un carro, encender y apagar luces, prótesis robóticas, entre otras.

2.8

EMOTIVSYSTEMS

EMOTIV Systems es una empresa de neuroingeniería formalizada en el año 2003, la cual desarrollo una revolucionaria interfaz cerebro – computador (BCI) para el consumo popular [16]. Esta interfaz desarrollada es conocida como EMOTIV EPOC, la cual consta de un dispositivo inalámbrico (headset) que por medio de múltiples sensores captura las señales eléctricas producidas por el cerebro para detectar expresiones faciales, movimiento, sentimientos y pensamientos. [17] 2.8.1

Headset

El headset EMOTIV EPOC (Figura 13) está compuesto principalmente por 16 electrodos de contacto con almohadillas (esta debe ser humedecida con una solución salina para garantizar una buena conducción) y enchapados en oro [17].

Figura 13. Headset EMOTIV EPOC.

Estos electrodos se encuentran sujetos a dos brazos de plástico que garantizan la correcta ubicación según el sistema internacional 10 – 20 (figura 2); es de resaltar que 14 de los 16 sensores capturan las señales EEG (los 2 restantes son usados como referencia CMS/DRL). La implementación de dos electrodos de referencia indica que el tipo de registro utilizado por este headset corresponde a un registro monopolar. Estos electrodos operan a una frecuencia de muestreo de 128sps y poseen un ancho de banda de 0.2Hz a 45Hz [17]. Adicional a los electrodos posee un giroscopio compuesto por dos acelerómetros (ejes X y Y) los cuales registran los movimientos de la cabeza del sujeto. Se destaca la conexión entre el headset y el PC es inalámbrica (transmisor bluetooth). 2.8.2

Research Edition (SDK)

El EMOTIVEPOC ofrece seis diferentes software development kit (SDK) que otorgan el control sobre la API del headset al igual que las bibliotecas (.dll) para el desarrollo de aplicaciones. En este proyecto se 22

empleara la edición RESEARCH. La elección de este SDK se debe a que el API permite la visualización y captura de las señales raw de los 14 electrodos, lo cual es indispensable para la realización de este proyecto. Es importante resaltar que actualmente EMOTIV EPOC no cuenta con un módulo en tiempo real para Matlab, a su vez es posible implementar el SDK para la captura de las señales raw por medio de C y C++ para su análisis posterior, sin embargo, si se desea realizar la etapa posterior de un BCI (Aplicación de control) es necesario usar de forma nativa el algoritmo de extracción y clasificación de EMOTIV EPOC conocido como EmoEngine [18].

23

3.

ESPECIFICACIONES

Se diseñó y desarrollo una interface cerebro computador mediante la clasificación de ondas encefalografías, para lograr ello el desarrollo fue dividido en tres etapas o fases; basándose en la estructura de un BCI descrita en la Figura 11, página 20 del presente documento. La Figura 14 indica cada una de las etapas para el desarrollo del presente proyecto.

Figura 14. Etapas de desarrollo

3.1

CAPTURA DE SEÑALES ENCEFALOGRÁFICAS

Se implementó el dispositivo no invasivo EMOTIVEPOC como receptor de señales electroencefalográficas. Este dispositivo está compuesto por 14 electrodos distribuidos según el sistema internacional 10-20 (Figura 15), los cuales operan a una frecuencia de muestreo de 128sps y poseen un ancho de banda de 0.2Hz a 45Hz [17].

Figura 15. Sensores EMOTIV EPOC según Sistema Internacional 10-20 [17].

El paradigma a clasificar se compone de dos casos: derecha e izquierda, el número total de muestras adquiridas por caso es de 50, distribuidas en 4 sesiones como lo indica la Tabla 4.

24

Numero de muestras por caso Primera 15 Segunda 10 Tercera 15 Cuarta 10 Sesión

Tabla 4. Distribución de sesiones para la captura de señales encefalográficas.

3.2

INTERFACE CEREBRO COMPUTADOR (BCI)

El proceso descrito en la Figura 16corresponde al realizado en el presente proyecto. Carga de Señales EEG Pre – Procesamiento

Eliminación de DC (Offset)

Ventaneo

Filtrado

CAR

Extracción de Características

PSD

PLV

Clasificación

SVM

SVM

Figura 16. Diagrama de bloques del BCI diseñado.

3.3

PRUEBAS DE DESEMPEÑO

Se realizó pruebas de desempeño mediante validación cruzada de K-iteraciones (K-foldcross-validation), garantizando que el algoritmo de clasificación tuvieran una acierto igual o mayor al 75%.

3.4

HAWARE Y SOFTWARE

Para la extracción y análisis de señales encefalográficas fue necesaria la implementación de dos computadoras. Se usó un computador portátil (laptop) para la extracción de las señales encefalográficas 25

debido a la necesidad desplazarse al lugar donde se encontraban algunos de los sujetos que participaron en el experimento, por otra parte el procesamiento y análisis de las señales encefalográficas obtenidas se realizó en un PC de escritorio (Desktop). Las especificaciones del laptop y PC se encuentran especificadas en las Tabla 5 y Tabla 6 respectivamente. Especificaciones Técnicas - Laptop Asus Marca G53SW Modelo Intel® Core™ 2630QM Processor Procesador 4 Núcleos 2.0 GHz Velocidad de Reloj 500 Gb Disco Duro 8 Gb Memoria RAM 15.6" Monitor Windows 7 Professional – 64 bits Sistema Operativo Tabla 5. Especificaciones Laptop.

Especificaciones Técnicas - Desktop Intel® Pentium® Processor E2140 Procesador 2 Núcleos 1.6 GHz Velocidad de Reloj 160 Gb Disco Duro 2 Gb Memoria RAM Windows 7 Professional – 32bits Sistema Operativo Tabla 6. Especificaciones Desktop.

El software implementado fue:    

3.5

Visual Studio 2010 Professional Matlab 2012b (Incluyendo Toolboxes) SDK – EmotivResearchEdition Test Bench – EmotivSystems

BANCO DE SEÑALES EEG

Mediante la captura de señales electroencefalográficas se generó un banco de señales. La generación de este banco tiene como objetivos principales:  

Implementación sencilla y ágil para el desarrollo del presente trabajo. Creación de un banco de señales EEG para trabajos futuros.

Cada muestra se encuentra debidamente etiquetada (ver Anexo 1), con un tamaño promedio de 150KB.

26

4.

DESARROLLO

4.1

CAPTURA DE SEÑALES EEG

La adquisición de señales electroencefalográficas se realizó de forma no invasiva por medio del dispositivo EMOTIV EPOC. La adquisición de estas señales se realizó mediante un experimento implementando una interfaz gráfica (GUI). La GUI extrajo y guardo las ondas encefalográficas de 10 sujetos para su posterior análisis en ficheros *.csv (comma-separatedvalues). La metodología del experimento es basado en anteriores trabajos como lo son: “Comparación de Diferentes Métodos para la Clasificación de Señales EEG en Interfaces CerebroComputadora” [19], “Análisis de electroencefalogramas en ritmos sensoriomotores. Control domótico mediante Brain Computer Interface” [20], “Automatic Differentiation of Multichannel EEG Signals” [21] y las bases de datos de señales EEG del “BCI Competition IV” [22].

4.1.1

Formulación de experimento (paradigma)

Cada prueba (Figura 17) tiene una duración de 9.5s a 11s. Este tiempo se divide en 3s de preparación, 5s donde se realiza la tarea mental y una pausa entre pruebas de 2 a 3 segundos. La prueba inicia con la presentación de una cruz en el centro de la pantalla (conocida en la literatura como fixation-cue o cross-fixation) tal y como se muestra en la Figura 18a y un tono de aviso a los 1.5 segundos (1kHz, 70ms) que indica al usuario que el experimento ha dado inicio, segundos después se mostrara una señal visual (Figura 18b o Figura 18c) durante 2 segundos. Esta señal visual da la instrucción al sujeto que imagine la acción motora correspondiente a la instrucción (esta acción durara 5 segundos). Al completar la tarea mental se darán dos tonos de aviso (1kHz, 70ms) indicando al usuario que puede descansar y se da una pausa de 1.5 segundos, adicionalmente se tendrá un tiempo variable de descanso entre pruebas de 0.5 a 1.5 segundos.

Figura 17. Esquema de tiempo de una prueba3.

3

La barra de tiempo está dada en segundos.

27

a

b

c

Figura 18. Estimulo visual. (a) Cruz, inicio de la tarea. (b) Caso 1, acción motora derecha. (c) Caso 2, acción motora izquierda.

4.1.2

Diseño de interfaz gráfica

La interfaz gráfica fue diseñada en Visual Basic y C ++ (incluidos en el paquete Visual Studio 2010 Professional) implementando el Research SDK que viene en conjunto al EMOTIV EPOC. Al iniciar el programa experimento.exe se abrirá la ventana de la Figura 19. En esta se debe ingresar el ID del sujeto (numérico), el número de pruebas a realizar (ej. Si se ingresa 5 pruebas, el programa captara 5 muestras para cada caso, es decir, 10 muestras en total) y finalmente seleccionar la sesión a realizar. Estos datos son solicitados al usuario con el fin de guardar el fichero *.csv con el nombre correspondiente, tal y como es indicado en el Anexo 1. A su vez el programa genera un fichero *.txt indicando por medio de letras R o L (Right – Derecha o Left – Izquierda, respectivamente) el orden en el cual fue ejecutado los dos estímulos visuales como se puede ver en la Figura 20.

Figura 19. GUI - Ingreso de datos para la captura de señales.

28

Figura 20. Archivo .txt generado el cual indica la secuencia para el Sujeto 4 – Sesión 2.

El diagrama de flujo de la Figura 21 indica la lógica implementada en el programa diseñado para la captura de señales EEG. Al iniciar el programa se establece unos valores iníciales los cuales corresponden a la cantidad de muestras para cada uno de los casos (0 = derecha, 1 = izquierda). Se implementa un condicional el cual determina que en el momento en el cual la suma de cant0 y cant1 (cantidad de muestras ejecutadas) sea igual al doble del número de muestras el programa finalizara. INICIO

cant0 = 0 cant1 = 0 n = número de pruebas

si

x<

cant1 = cant1+1

cant0 = cant0+1

x = cant0 + cant1

no

2*n FIN

Imagen Cruz

op =Rand (Opción Señal)

op = 0

op

op = 1

Imagen Derecha

Imagen Izquierda

Figura 21. Diagrama de flujo de la lógica implementada para la captura de las señales EEG.

29

Es importante resaltar que se implementó la hora del computador como “semilla”, con el fin de garantizar aleatoriedad en la señal mostrada como en el tiempo de descanso. La duración de la prueba depende de la variable “Numero de Pruebas”, si el usuario ingresa 2 la duración máxima de la prueba será de 22segundos.Serealizó igualmente la captura de 40 muestras (divididas en 3 y 4 sesiones) de señales baseline, en este caso la pantalla es blanca; estas señales no fueron empleadas en el desarrollo del presente proyecto, sin embargo forman parte del banco de señales para trabajos futuros. 4.1.3

Captura de señales EEG

Las señales adquiridas constan de datos EEG de 10 personas (5 hombres y 5 mujeres, todos sanos4). Cada uno se encontró sentado cómodamente en una silla, con buena iluminación, delante de un portátil con pantalla de 15.6’’ (ver Tabla 5. Especificaciones Laptop.) y resolución de pantalla de 800 x 600 pixeles. Previo al experimento se verifico por medio del software de Emotiv TestBench que los 14 electrodos se encontraran realizando un contacto adecuado5 (Figura 22). Antes de iniciar el experimento este fue explicado a cada uno de los usuarios, indicando que debían imaginar la acción motora al ver el estímulo visual.

Figura 22. EmotivTestBench - Sujeto 04

El tiempo total del experimento por sujeto fue de aproximadamente de 35 a 45 minutos, tal y como lo indica la Tabla 7. Sin embargo algunas personas manifestaron cansancio y dificultad de concentración a partir de la segunda sesión, por lo que algunos descansos fueron superiores a 10 minutos. 4 5

Anexo 2. Ficha Médica. Anexo 3. Indicaciones de la correcta colocación del headset Emotiv.

30

Actividad Sesión 1 Baseline 1 Sesión 2 Baseline 2 Descanso Sesión 3 Baseline 3 Sesión 4 Baseline 4 TOTAL

Numero de Muestras6 Tiempo Aproximado 30 6 minutos 10 2 minutos 20 4 minutos 10 2 minutos 10 minutos 30 6 minutos 10 2 minutos 20 4 minutos 10 2 minutos 1407 35 minutos

Tabla 7. Tiempo aproximado de experimento por sujeto.

Figura 23. Realización del experimento de captura de datos.

6

Sumatoria de muestras de izquierda y derecha. El sujeto 9 presento una fuerte cefalea durante el experimento, por lo tanto el número de muestras es inferior en algunas sesiones.

7

31

Los ficheros generados para cada sujeto fueron posteriormente importados a MATLAB para así crear un fichero *.mat el cual contuviera todas las señales de cada uno de los sujetos organizadas en matrices y estructuras.

4.2

PROCESAMIENTO DE SEÑALES PARA EL BCI

4.2.1

Selección de Canales

Se optó por mantener la totalidad de los canales EEG (14 canales), sin embargo se realizó un preprocesamiento a todas las señales con el fin de garantizar una relación señal-ruido (SNR) alta para lograr un desempeño adecuado en las etapas de extracción de características y clasificación. 4.2.2

Eliminación de DC – Offset

En las señales obtenidas mediante el headset EMOTIV EPOC se encuentra que existe un desplazamiento causado por un DC (Figura 24). Este DC se debe a que el headset agrega un DC sobre los 4096 µV con el fin de garantizar que todos los datos sean un entero positivo [23]. Por lo tanto se realizó la eliminación de este DC restando el valor medio de la señal original (Figura 25).

SEÑAL ORIGINAL 4900 AF3 F7 F3 FC5 T7 P7 O1 O2 P8 T8 FC6 F4 F8 AF4

4800 4700

Amplitud (uV)

4600 4500 4400 4300 4200 4100 4000

0

200

400 600 800 Numero de Muestras

1000

1200

Figura 24. Señales de los 14 canales originales – AM89: Derecha.

8

AM = Acción Motora Todas las figuras mostradas en el análisis, extracción y clasificación de señales corresponden al sujeto 4, sesión 3, muestra 4.

9

32

SEÑAL SIN DC 250 AF3 F7 F3 FC5 T7 P7 O1 O2 P8 T8 FC6 F4 F8 AF4

200

Amplitud (uV)

150

100

50

0

-50

-100

0

200

400 600 800 Numero de Muestras

1000

1200

Figura 25. Señales de los 14 canales sin DC.

4.2.3

Ventana Hanning

En el experimento realizado para la extracción de señales EEG la tarea mental fue realizada en el tiempo de 3 a 8 segundos, con un estímulo visual de 3 a 5 segundos. Por tal motivo se aplicó una ventana hanning con el fin de obtener las señales EEG que contienen la información durante la señal visual y tarea mental, tal y como lo indica la Figura 26.

Figura 26. Diagrama de tiempo y ventana hanning

33

Señal Offset - ID04.SIGNAL.OFFSET.S3.R.N4 250

Amplitud - uV

200 150 100 50 0 -50 -100

0

200

400

600 Numero de Muestras

800

1000

1200

800

1000

1200

800

1000

1200

Ventana Hanning 1 0.8 0.6 0.4 0.2 0

0

200

400

600 Numero de Muestras Multiplicación de la señal con la Ventana Hanning

40

Amplitud - uV

20 0 -20 -40 -60

0

200

400

600 Numero de Muestras

Figura 27. (a) Señal EEG sin offset (b) Ventana Hanning (c) Señal Ventaneada

4.2.4

Filtrado

Debido a la naturaleza de las señales EEG es conveniente realizar un filtrado temporal, de esta forma se elimina ruido de frecuencias bajas (generalmente relacionados a artificios fisiológicos como el parpadeo) y de frecuencias altas (ruido de línea 50 – 60 Hz). Por lo tanto posterior a la eliminación de DC y ventaneo se diseñó un filtro con respuesta impulso finita - equiripple (Figura 28), con el objetivo de solo obtener los componentes en la banda de interés.

Filtro FIR - Equiripple 20 0

Magnitude (dB)

-20 -40 -60 -80 -100 -120 -140

0

10

20

30 Frequency (Hz)

40

50

60

0

10

20

30 Frequency (Hz)

40

50

60

500

0

Phase (degrees)

-500

-1000

-1500

-2000

-2500

-3000

Figura 28. Filtro Equiripple (a) Magnitud (b) Fase

34

Para el diseño de este filtro se consideraron las siguientes características: 

Filtro Pasa Banda o o o o o o

Orden = 32410 Frecuencia de Muestreo = 128 Hz Frecuencia de Banda Pasante = 7 Hz Frecuencia de Banda Atenuada = 14 Hz Máximo Rizado Banda Pasante = 1 dB Mínima Atenuación Banda Atenuada = 80 dB

Una de las ventajas de la implementación de un filtro FIR – Equiripple es su fase lineal, esta garantiza que en la extracción de características mediante diferencia de fase (PLV) se pueda verificar condiciones de simetría. Igualmente este tipo de filtro el error de aproximación entre la respuesta en frecuencia ideal y la real se reparte uniformemente en cada banda, pasante y atenuada (de ahí el nombre de equiripple), minimizando el error máximo en cada una de ellas [24]. Podemos observar que posterior al filtrado (Figura 29) se han eliminado los artificios de bajas frecuencias al igual que solo obtenemos las componentes de la frecuencia de interés. FFT SEÑAL ID04.SIGNAL.OFFSET.S3.R.N4 5000 AF3 F7 F3 FC5 T7 P7 O1 O2 P8 T8 FC6 F4 F8 AF4

4500 4000 3500 3000 2500 2000 1500 1000 500 0

0

10

20

30

40

50

60

70

Frecuencia (Hz)

FFT SEÑAL ID04.SIGNAL.OFFSET.S3.R.N4 - FILTRADA 900 AF3 F7 F3 FC5 T7 P7 O1 O2 P8 T8 FC6 F4 F8 AF4

800 700 600 500 400 300 200 100 0

0

10

20

30

40

50

60

70

Frecuencia (Hz)

Figura 29. a) Espectro de Frecuencia Señal Original b) Espectro de Frecuencia Señal Filtrada

4.3

Filtros Espaciales

El headset EMOTIV EPOC realiza un registro monopolar, es decir toma la señal independiente de cada uno de los electrodos. Una de las desventajas de este registro es que no se puede garantizar que la posición de los electrodos sea la misma entre ensayos; sin embargo para minimizar este efecto se procuró que durante todo el experimento nunca se fuese retirado el headset. Igualmente debido a que los potenciales 10

El orden del filtro fue obtenido mediante la función FIRPMORD en MATLAB.

35

del cuero cabelludo son referencias dependientes, es decir, aunque no exista contacto alguno con el cuero cabelludo, la señal capturada por el electrodo es afectada por el electrodo de referencia [25]. El filtrado espacial (spatial filtering) es una técnica implementada para minimizar el ruido generado debido a la variación del posicionamiento de los electrodos entre ensayos al igual que los artificios generados por otras fuentes. El filtrado espacial también proporciona una relación SNR grande, permitiendo que la extracción de características sea más efectiva [26]. Existen tres tipos de filtros espaciales:



Filtro Bipolar El filtro espacial bipolar toma la tensión promedio de los electrodos vecinos posicionados verticalmente, esta tensión es restada a señal del electrodo central, como lo indica la Figura 30a.

Figura 30. Filtros Espaciales [26]. (a) Filtro Bipolar, (b) Filtro Laplaciano, (c) Referencia Media Comun



Filtro Laplaciano El filtro laplaciano se basa en tomar la señal de cada electrodo y la normaliza con las señales de los electrodos circundantes, es decir, el promedio de la tensión de los electrodos vecinos verticales y horizontales es restada de la señal de cada electrodo como lo indica la Figura 30b. Este filtro está enfocado a la fuente (electrodo de interés) a diferencia del filtro bipolar [27]. Este puede ser calculado de la siguiente forma:

Vi Lap  Vi 

1  Vj 4 jS j

(1)

donde Vi es el potencial del electrodo (o canal) i-ésimo y S j es el índice del electrodo vecino.



Referencia de Media Común (Common Average Reference – CAR) Teóricamente si existe un igual número de electrodos espaciados por igual y que cada señal fuese generada por una fuente puntual, la referencia de media común (CAR) sería igual a cero. Sin embargo, en la aplicación, las señales no son generadas por fuentes puntuales y los electrodos no son equidistantes. Por lo tanto se implementa el filtro CAR para que la señal de cada electrodo sea 36

independiente de la referencia [28]. Para implementar el filtro CAR en cada electrodo se debe restar el valor medio de todas las señales a la señal del electrodo en cuestión, esto puede expresarse como:

Vi CAR  Vi 

1 n V j n j 1

(2)

donde Vi es el potencial del electrodo (o canal) i-ésimo y n es el número total de electrodos [29]. Señal Original - Filtrada 15 AF3 F7 F3 FC5 T7 P7 O1 O2 P8 T8 FC6 F4 F8 AF4

10

Amplitud - mV

5

0

-5

-10

-15

0

200

400

600 Numero de Muestras

800

1000

1200

Señal Original - Filtrada + CAR 15 AF3 F7 F3 FC5 T7 P7 O1 O2 P8 T8 FC6 F4 F8 AF4

10

Amplitud - mV

5

0

-5

-10

-15

0

200

400

600 Numero de muestras

800

1000

1200

Figura 31. Aplicación de filtro CAR. (a) Señal Original (Filtrada) (b) Señal filtrada mediante CAR Señal Electrodo FC5 6 Señal Original - Filtrada Señal Original - Filtrada + CAR 4

Amplitud - mV

2

0

-2

-4

-6

0

200

400

600 Numero de Muestras

800

Figura 32. Señal FC5 (Azul – Señal original Filtrada, Roja – Señal Filtrada por CAR)

37

1000

1200

Finalizado el pre-procesamiento se ajustan todas las señales la misma longitud, en este caso para 9 segundos (1152 muestras). Los 9 segundos elegidos garantiza que no se está eliminando datos de la tarea mental. 4.4

EXTRACCIÓN DE CARACTERÍSTICAS

El paradigma y el experimento realizado para la captura de señales EEG es basado en tareas motoras, la señales obtenidas contienen el fenómeno conocido como Desincronización/Sincronización relacionada a eventos (ERD / ERS). Por lo tanto la etapa de extracción será enfocada únicamente a métodos de extracción basados en el dominio de la frecuencia.

4.4.1

Densidad del Espectro de Potencia (Power Spectral Density – PSD)

Más de un tercio de los BCI basados en ritmos sensorio motores y ERD / ERS implementan métodos no paramétricos para la extracción de características, entre estos se destaca el método de densidad del espectro de potencia (PSD) [10]. El espectro de potencia determina la distribución de la potencia de una señal a lo largo de un intervalo de frecuencias. El cálculo del espectro de frecuencia es se realiza por lo general aplicando la transformada de Fourier sobre la señal de interés. Cuando el espectro de potencia se calcula por medio de la transformada de Fourier y luego es promediada se conoce como periodograma [30]. Para realizar el cálculo de la densidad espectral de potenciase parte de que la DFT de una función muestreada x( n) está dada por: N 1

X k  T  x ( n )e



j 2 kn N

(3)

n 0

donde T es el periodo de muestreo. Por lo tanto la densidad espectral de energía (Energy spectral density – ESD) está dada por:

Sk  X k

2

(4)

Sin embargo si x(n) corresponde a un proceso estocástico (como en el caso de las ondas EEG), la energía seria infinita, y por consiguiente el parámetro de interés seria la potencia. Por lo tanto se halla la densidad espectral de potencia promediando temporalmente la energía:

X S Pk  k  k NT NT

38

2

(5)

No obstante el periodograma obtenido no es un estimador consistente debido a que su varianza no tiende a cero cuando la longitud de la señal tiende a infinito. Es por esto que se implementa el método de Welch para evaluar el espectro de potencia, este método consiste en dividir los datos en varios segmentos traslapados, y se calcula la potencia en cada segmento. Finalmente se calcula el cuadrado de la magnitud y se promedia el espectro [31]. El algoritmo para el cálculo de la densidad del espectro de frecuencia (PSD) puede observarse en el diagrama de flujo de la Figura 33. Señal EEG Preprocesada ch = 1 fs = 128 Extraer Canal

WELCH Ventaneo Hamming

FFT

ch = ch +1

canal = 14

PSD

Figura 33. Diagrama de Flujo - Algoritmo PSD

Al implementar el algoritmo se obtiene la densidad de potencia espectral de cada uno de los canales como se pueden observar en la Figura 34 y Figura 35.

39

Welch Power Spectral Density Estimate Sujeto: 04, Sesión: 3, Caso: Derecha , Muestra: 4 3 AF3 F7 F3 FC5 T7 P7 O1 O2 P8 T8 FC6 F4 F8 AF4

Power/Frecuency (dB/Hz)

2.5

2

1.5

1

0.5

2

4

6

8 10 Frecuency (Hz)

12

14

16

Figura 34. Matriz PSD - AM: Derecha

Welch Power Spectral Density Estimate Sujeto: 04, Sesión: 3, Caso: Izquierda , Muestra: 4 1.8

AF3 F7 F3 FC5 T7 P7 O1 O2 P8 T8 FC6 F4 F8 AF4

1.6

Power/Frecuency (dB/Hz)

1.4

1.2

1

0.8

0.6

0.4

0.2

0

5

10

15

Frecuency (Hz)

Figura 35. Matriz PSD - AM: Izquierda

Debido a que Matlab posee en su toolbox de procesamiento de señales (Signal Processing Toolbox), el cual incluye las siguientes funciones útiles para el cálculo de la densidad de potencia espectral, entre estas se encuentran pwelch, spectrum y psdse realizó el método de Welch con dichas funciones. Se resalta que el traslape asociado al periodograma calculado por el método de Welch mejora notablemente la correlación de los datos y propiedades estadísticas, por consiguiente mejora la confiabilidad de la estimación [32]. 40

4.4.2

Valor de Bloqueo de Fase (Phase-Locking Value)

Las señales EEG sufren una desincronización cuando sucede una tarea motora (ERD), por lo tanto se implementó la extracción de características mediante el método de diferencia de fase o valor de bloqueo de fase (Phase-Locking Value PLV) [33].El método de PLV mide la sincronización entre las señales EEG de un par de sensores, en una frecuencia de interés. La PLV es definida como [34, 35]:

PLVij 

1 N

e



j i  n j  n



(6)

N

donde i y  j son las fases instantáneas de los sensores i y j en un tiempo n para N muestras. La medida PLV se encuentra en un rango de 0 a 1. Por lo tanto cuando PLVij  1 existe sincronización entre las señales de los dos sensores y PLVij  0 no existe sincronización alguna. Debido a que el PLV es una medida entre dos señales es necesario definir un conjunto de parejas, por lo tanto para los 14 electrodos obtenemos por medio de combinatoria se tienen un total de 91 combinaciones posibles:

m! n !(m  n)!

(7)

14! 14!  2!(14  2)! 2!12!

(8)

C142  91 combinaciones

(9)

Cmn 

C142 

Para realizar el cálculo del PLV se requiere la fase instantánea de cada una de las parejas de sensores, es por esto que se implementa la transformada de Hilbert para obtenerla, y así calcular su sincronismo. Este procedimiento puede verse en el diagrama de flujo de la Figura 36.

41

Señal EEG Preprocesada

Selección de Canales Pares

Transformada de Hilbert

Extracción de Fases Instantane Generación de PLV

Normalización

no

P == 91

si

PLV

Figura 36. Diagrama de Flujo - Algoritmo PLV

Es importante recordar que las señales fueron filtradas por un filtro equiripple de orden 324 y esto causa un efecto de transitorio en la señal filtrada. Este efecto se puede observar en la Figura 37b y que los primeros 300 ms los valores del PLV del par de señales cambian abruptamente (esto causado por el orden del filtro del pre-procesamiento), por tal motivo se eliminó las primeras 324 muestras.

42

a) PHASE-LOCKING VALUE Sujeto: 04, Sesión: 3, Electrodos: AF3 - AF4 1 AF3 AF4

Phase Locking Value - PLV Normalizado

0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1

0

1

2

3

4 5 Tiempo (s)

6

7

8

9

b) Figura 37. a) Par de electrodos AF3 y AF4 b) PLV – Pareja de Electrodos AF3 – AF4. Sujeto 4 – Sesión 3

4.5

CLASIFICACIÓN DE CARACTERISTICAS

El objetivo del clasificador de un sistema BCI es identificar intensión del sujeto basándose únicamente de las características de las señales obtenidas, es decir, el clasificador agrupa las características de las señales EEG y las asocia a una acción motora. Una de las técnicas de clasificación más implementadas en el campo de procesamiento de señales EEG es el aprendizaje de supervisado [10]. El aprendizaje supervisado es una técnica implementada para deducir una función o un modelo a partir de datos de entrenamiento, es decir, el aprendizaje supervisado adapta automáticamente la estructura del modelo clasificatorio utilizando un conjunto de datos de entrenamiento, de esta forma los datos a clasificar 43

son relacionados con las características de los datos de entrenamiento [36]. Entre los algoritmos de clasificación de aprendizaje supervisado se destaca las Maquinas de Soporte Vectorial (Support Vector Machine – SVM).

4.5.1

Máquinas de Soporte Vectorial (SVM)

Las máquinas de soporte vectorial es un método de clasificación de aprendizaje supervisado el cual determina el mejor hiperplano para la discriminación de dos clases. En el presente trabajo se implementó la librería SVM libsvm-mat-3.12 para Matlab creada por Chang, Chih-Chung and Lin, Chih-Jen [37]. La metodología de entrenamiento y ensayo fue mediante hold-out, este método consistió en dividir en dos conjuntos (entrenamiento y ensayo) las sesiones. Esta división de datos puede observarse en laTabla 8.

Set Entrenamiento

Set Ensayo

Sesión 2 Sesión 3 Sesión 4

Sesión 1

Sesión 1 Sesión 3 Sesión 4

Sesión 2

Sesión 1 Sesión 2 Sesión 4

Sesión 3

Sesión 1 Sesión 2 Sesión 3

Sesión 4

Tabla 8. Sets de Entrenamiento y Ensayo mediante Hold-Out

Se implementaron dos tipos de kernel de clasificación presentes en la librería LIBSVM11 [38]: 

LINEAL

K  xi , x j   xiT x j

11

Los desarrolladores de LIBSVM aconsejan usar la librería LIBLINEAR para el kernel lineal.

44

(10)



BASE RADIAL (Radial BasisFunction – RBF):

K  xi , x j   e

  xi  x j

2

,  0

(11)

Para ambos kernels de clasificación se implementó validación cruzada con 5-iteraciones en los datos de entrenamiento, con el fin de encontrar los mejores parámetros. 4.5.1.1 Kernel Lineal En el caso del kernel lineal no existe ningún parámetro como lo indica la ecuación 4.10, sin embargo se asignó un parámetro Costó óptimo de clasificación (C para LIBSVM). El parámetro C (Cost o Costo) indica al SVM cuánto desea evitar una clasificación errónea. Es decir, para grandes valores de C , el modelo SVM establecerá un hiperplano de menor margen. Al contrario de un valor C pequeño, el modelo SVM determinara un hiperplano de mayor margen. Sujeto / Set de Ensayo 1 2 3 4 5 6 7 8 9 10

S1 74.26% 74.80% 80.61% 77.07% 86.31% 82.95% 80.61% 79.80% 74.09% 86.87%

S2 73.20% 78.81% 80.95% 75.50% 85.26% 77.06% 73.10% 74.62% 81.47% 80.80%

S3 84.63% 78.12% 72.11% 82.57% 81.96% 81.95% 71.59% 69.30% 72.31% 82.08%

S4 74.00% 83.40% 78.43% 78.37% 79.19% 74.28% 79.48% 82.09% 80.33% 78.83%

Tabla 9. Precisión PSD - Kernel Lineal

Sujeto / Set de Ensayo 1 2 3 4 5 6 7 8 9 10

S1 71.06% 78.75% 66.12% 76.19% 73.08% 86.08% 71.61% 83.33% 79.49% 76.19%

S2 73.26% 78.21% 61.90% 73.63% 78.57% 80.04% 74.54% 84.80% 82.97% 78.75%

S3 75.46% 78.02% 60.99% 73.26% 76.37% 81.14% 70.33% 78.39% 77.84% 75.27%

Tabla 10. Precisión PLV - Kernel LinealKernel RBF

45

S4 73.26% 80.95% 63.55% 73.81% 79.49% 81.14% 73.08% 84.62% 75.27% 76.56%

4.5.1.2 Kernel Función Base Radial A diferencia del lineal, el kernel de base radial tiene dos parámetros que son gamma    y C. Debido a que estos parámetros son desconocidos y se requiere encontrar una combinación optima se realiza un procedimiento llamado “Busqueda por Red” o “Grid-Search”.Este procedimiento se basa en realizar pruebas de desempeño en múltiples combinaciones de los parámetros   C en diferentes escalas para asi encontrar los parámetros que brindan una mejor clasificación. La Figura 38 corresponde a una búsqueda por red con la siguiente secuencia para los parámetros:

C  210 , 210 ,..., 210 , 220 y   210 ,210 ,...,210 ,2 20 . Se puede observar que la mejor combinación se encuentra en C  210 |   210 , ya que esta logra una precisión del 80.22%.

Figura 38. Búsqueda por Red - Escala Grande (Característica: PLV - ID: 04 – Sesión de Ensayo: S3)

46

Sin embargo los autores de la librería LIBSVM recomiendan realizar una busqueda “fina” [38], por lo tanto se toma la anterior combinación como eje central ( C  210 |   210 ) y se realiza nuevamente la busqueda por red a una escala mas pequeña. En la Figura 39 se aprecia que la secuencia de busqueda fue:

C  20 , 25 ,..., 215 , 220 y   220 ,215 ,...,2 5,2 0 , y la mejor combinación fue C  25 |   25 .

Figura 39. Búsqueda por Red - Escala Mediana (Característica: PLV - ID: 04 – Sesión de Ensayo: S3)

Finalmente se realiza una ultima busqueda por red siguiendo las mismas pautas, es decir se tiene como punto central la mejor combinación ( C  25 |   25 ) y se realiza la busqueda con pasos mas cortos, en el caso de la Figura 40 estos fueron C  20 , 22 ,..., 28 , 210 y   210 ,28 ,...,2 2,2 0 , para obtener finalmente una precisión del 90.48%, para la combinación C  22 |   25 .

47

Figura 40. Búsqueda por Red - Escala Pequeña (Característica: PLV - ID: 04 – Sesión de Ensayo: S3)

Se puede tasar que efectivamente la busqueda mediante red es una herramienta útil para hallar los parametros más optimos ya que se logro incrementar de un a precisión del 80.22% a finalmente un 90.48%. Los resultados del clasificador RBF para las características PSD y PLV pueden observarse en las tablas 11 y 12 respectivamente.

Sujeto / Set de Ensayo 1 2 3 4 5 6 7 8 9 10

S1 82.07% 83.59% 76.69% 84.80% 84.76% 78.78% 63.06% 70.21% 77.12% 77.26%

S2 77.05% 79.27% 74.76% 75.41% 93.74% 75.65% 62.80% 76.31% 74.42% 74.11%

Tabla 11. Precisión PSD - Kernel RBF

48

S3 72.63% 76.56% 72.85% 76.43% 86.27% 77.71% 70.86% 65.68% 72.57% 82.53%

S4 79.44% 75.50% 80.56% 81.36% 86.59% 72.53% 67.29% 78.62% 88.94% 83.54%

Sujeto / Set de Ensayo 1 2 3 4 5 6 7 8 9 10

S1 87.91% 92.86% 83.33% 86.45% 94.69% 94.32% 87.73% 95.42% 97.07% 94.87%

S2 87.55% 91.39% 73.99% 86.81% 93.22% 90.48% 83.33% 89.93% 96.70% 94.32%

Tabla 12. Precisión PLV - Kernel RBF

49

S3 91.39% 95.24% 81.14% 90.48% 96.15% 90.84% 85.53% 91.21% 93.22% 94.32%

S4 85.16% 90.11% 78.39% 86.26% 95.24% 91.39% 81.32% 91.39% 93.04% 93.96%

5.

ANALISIS DE RESULTADOS

En este capítulo se presenta una descripción de los resultados experimentales obtenidos en cada una de las tres etapas (Captura de Señales, Extracción de Características y Clasificación de Características) del desarrollo del presente trabajo.

5.1

CAPTURA DE SEÑALES

Para ilustrar el tipo de señales capturada se puede ver la Figura 41a. Cada uno de los canales se encuentra con un voltaje DC aproximadamente de 4096 µV. Por tal motivo se decidió realizar la implementación de la etapa de preprocesamiento para obtener finalmente una señal EEG adecuada (sin artificios) para la posterior etapa de extracción de características como se puede apreciar en la Figura 41c. Señal Original 5000 4800 4600 4400 4200 4000

0

200

400

600

800

1000

1200

1400

800

1000

1200

1400

Numero de Muestras Señal Sin DC 250 200 150 100 50 0 -50 -100

0

200

400

600 Numero de Muestras Señal Pre-Procesada

15 10 5 0 -5 -10 -15

0

200

400

600 Numero de Muestras

800

1000

1200

Figura 41.Señal EEG correspondiente a la tarea motora Derecha, Sujeto 4 – Sesión 3 – Muestra 4. (a) Señal Original, (b) Señal sin DC, (c) Señal Pre-Procesada12

5.2

EXTRACCIÓN Y CLASIFICACIÓN PSD - SVM

La estimación del espectro de densidad de potencia (PSD) se realizó mediante la aplicación del método de Welch, esta estimación para cada uno de los 14 canales.Se pudo observar que los electrodos ubicados en el área precentral del lóbulo frontal ej. AF3 – AF4 tenían una mayor componente de frecuencia en la banda mu (Figura 42). Es de resaltar este fenómeno ya que esta área corresponde a la corteza motora primaria la

12

El preprocesamiento consta de ventaneo, filtro en banda de interés (filtro FIR - Equiripple) y filtrado espacial CAR.

50

cual tiene como función el llevar a cabo los movimientos individuales de diferentes partes del cuerpo [7,8]. Estimación PSD por Welch - AM: Derecha

Estimación PSD por Welch - AM: Izquierda

3

1.5

Power/Frecuency (dB/Hz)

2.5

2

1.5

1

Power/Frecuency (dB/Hz)

AF3 F7 F3 FC5 T7 P7 O1 O2 P8 T8 FC6 F4 F8 AF4

AF3 F7 F3 FC5 T7 P7 O1 O2 P8 T8 FC6 F4 F8 AF4

1

0.5

0.5

0

0

10

20

30 40 Frecuency (Hz)

50

60

0

70

0

10

20

30 40 Frecuency (Hz)

50

60

70

Figura 42. PSD estimado mediante Welch en la banda de frecuencia de interés. (a) PSD – AM: Derecha, (b) PSD – AM: Izquierda.

La extracción de características PSD y clasificación SVM obtuvo el siguiente desempeño:

SUJETO SVM - LINEAL SVM - RBF 1 76.52% 77.80% 2 78.78% 78.73% 3 78.03% 76.21% 4 78.38% 79.50% 5 83.18% 87.84% 6 79.06% 76.17% 7 76.20% 66.00% 8 76.45% 72.70% 9 77.05% 78.26% 10 82.15% 79.36% PROMEDIO 78.58% 77.26% Tabla 13. Desempeño promedio PSD + SVM

Se destaca en la Tabla 13 que el modelo de clasificación para el sujeto 5 obtuvo el mejor desempeño tanto en el clasificador lineal como el RBF con una precisión de 83.18% y 87.84% respectivamente. Igualmente se resalta que la precisión no varía significativamente entre clasificadores, por lo tanto se puede concluir que para características PSD el kernel Lineal y RBF tiene un desempeño similar. 51

5.3

EXTRACCIÓN Y CLASIFICACIÓN PLV – SVM

La estimación del valor de bloqueo de fase se realizó mediante la transformada de Hilbert, ya que esta nos brinda información sobre la fase inmediata. Se obtuvo el siguiente desempeño de los clasificadores para las características PLV:

SUJETO SVM - LINEAL SVM - RBF 1 73.26% 88.00% 2 78.98% 92.40% 3 63.14% 79.21% 4 74.22% 87.50% 5 76.88% 94.83% 6 82.10% 91.76% 7 72.39% 84.48% 8 82.78% 91.99% 9 78.89% 95.01% 10 76.69% 94.37% PROMEDIO 75.93% 89.95% Tabla 14. Desempeño promedio PLV + SVM

A diferencia de los resultados para las características PSD (Tabla 13) hubo una diferencia significativa entre la clasificación Lineal y RBF como puede observarse en la Tabla 14, claramente el kernel de clasificación RBF tuvo un mejor desempeño que el lineal logrando alcanzar precisión hasta de un 95.01% (sujeto 9). La Figura 43 permite concluir que las características PLV al tener una formulación teórica más robusta al realizar la comparación de fases entre pares de sensores/canales logra que el clasificador tenga un mejor desempeño, a diferencia del PSD que las características son individuales. Igualmente el vector PLV contiene 91 características (correspondientes a las 91 parejas de electrodos), mientras que el vector PSD únicamente 14 características (número de canales).

52

Figura 43. Precisión promedio de clasificación para todos los sujetos.

Es de destacar que los resultados obtenidos con las características PLV, concuerdan con la literatura actual [39, 40], al igual que cada una de las combinaciones de características y clasificador logro un mínimo de 75% de precisión.

53

6.

CONCLUSIONES

En este proyecto, se consideró el uso de señales EEG como la mejor solución para registrar la actividad cerebral por tratarse de un método no invasivo, de bajo costo y portátil. Por tal motivo se implementó el headset de EMOTIV para obtener las señales EEG en cuestión. Además, se trabajó con los ritmos sensoriomotores ya que son los más adecuados para realizar una aplicación basada en actividades motoras (ej. Movimiento del brazo izquierdo). El análisis de métodos de procesamiento de señales EEG comenzó con un estudio exhaustivo de los algoritmos de extracción y clasificación, concluyendo así que los métodos de extracción de características que se ajustaban a los datos obtenidos mediante el headset de Emotiv eran Densidad del Espectro de Potencia (PSD) y Valor de Bloqueo de Fase (PLV). De esta forma, el análisis de estos dos métodos de extracción se realizó empleando aprendizaje supervisado por medio de máquinas de vectores de soporte (SVM). Se pudo apreciar claramente que un sistema BCI basado en sincronismo y valor de bloqueo de fase (PLV) fue superior al método tradicional basado en espectros de potencia (PSD), teniendo en cuenta la escasez de datos disponibles. Esto se debe a que las características PLV poseen una formulación teórica más robusta al realizar la comparación de fases entre pares de sensores/electrodos logrando así que el clasificador SVM tenga un desempeño superior al método PSD, ya que este las características son individuales. Se destaca que el método PLV obtiene los datos de 91 combinaciones de pares de electrodos, a diferencia del PSD que solo obtiene las características de 14 electrodos. Se resaltó la importancia de eliminar el ruido presente en las señales EEG conocidos como artificios. Estos ruidos se deben a actividades fisiológicas (ej. Palpitar del corazón, respiración, parpadeo) y no fisiológicas (ej. ruido de línea (50 – 60 Hz). Durante el proceso de extracción y clasificación de características se obtuvo inicialmente una precisión inferior al 60%, esto debido a que no se realizó el debido preprocesamiento a las señales. Finalmente al realizar el pre – procesamiento adecuado (ventaneo, filtro espacial, filtros en banda de interés) la precisión de la clasificación incremento considerablemente. Cabe resaltar que a pesar que el SDK de Emotiv aún no permite el análisis en tiempo real de las señales, el presente estudio puede ser implementado a cualquier señal EEG obtenida a partir del paradigma descrito en el capítulo 4.1.1. Como trabajo futuro de este estudio; la adquisición de datos en tiempo real y la clasificación se puede extender para controlar una silla de ruedas, encendido y apagado de luces (domótica) y demás aplicaciones de control binarias con el fin de brindar una ayuda a personas con discapacidad.

54

7.

BIBILIOGRAFIA

[1] QUINN, Gerard. Derechos humanos y discapacidad: Uso actual y posibilidades futuras de los instrumentos de derechos humanos de las Naciones Unidas en el contexto de la discapacidad. Nueva York y Ginebra: Naciones Unidas, 2002. ISBN 92-1-354074-4. [2] VÁSQUEZ, Armando. La discapacidad en América Latina. Disponible en Internet: . [3] DEPARTAMENTO ADMINISTRATIVO NACIONAL DE ESTADISTICA, DANE. Censo General 2005: Discapacidad, Personas con limitaciones permanentes. Colombia, 2005. [4] WOLPAW, Jonathan. Brain–computer interfaces for communication and control. En: Clin. Neurophysiol. Vol. 113, No.6 (2002); p. 602-614. [5] NAVARRO, Rafael. Instrumentación Biomédica. Tema 5: Electroencefalografía. Departamento Electrónica, Universidad Alcalá. [6] SABATI, Alejandro G. Interacción Humano Maquina. Ondas Neuroeléctricas utilizadas en el sistema BCI o IMCM. Argentina, 2012.Disponible en Internet: [7] BRAIN CONNECTION. The Anatomy of Movement.California 2012. Disponible en Internet: [8] KOLB, Brian; WHISHAW, Ian. Organización del sistema motor. Madrid, Editorial Médica Panamericana, 2006. [9] PFURTSCHELLER, G y LOPES DA SILVA, H. Event-related EEG/MEG synchronization and desynchronization: basic principles. Clinical Neurophysiology, vol. 110, pp. 1842 – 1857, 1999 [10] BASHASHATI, A; FATOURECHI, M y WARD, R. A survey of signal processing algorithms in braincomputer interfaces based on electrical brain signals. Journal of Neural Engineering, vol. 4, pp. R32– R57, 2007 [11] HONG, Bo.,WANG, Yijun., GAO, Xiaorong., y GAO, Shangkai. Quantitative EEG-Based BrainComputer Interface.En: TONG, Shanbao. Quantitative EEG Analysis Methods and Applications. MA: Artech House Publishers, 2009. p.193-224. [12] MAK, J.N y WOLPAW, J.R. Clinical applications of brain-computer interfaces: Current state and future prospects,” IEEE Reviews in Biomedical Engineering, vol. 2, pp. 187–199, 2009. [13] PERASSOLO, Mónica. Argentina, 2008.

Electroencefalografia Artificios SeccionNeurologia. Buenos Aires,

[14] ELKAN, Charles. Evaluating Classifiers.Universityor California, San Diego, 2011. [15] JOANNEUM. ProClassify User’s Guide.Austria.: Graz University of Technology, 2005. Disponible en internet: 55

[16] EmotivSystems.Corporate. Disponible en Internet: [17] EMOTIV. Emotiv Research SDK Specifications.Disponible

en

Internet:

[18] EMOTIV. Access real time raw EEG data using MATLAB or C++.Disponible en Internet: [19] FARFAN, Fernando D. Comparación de Diferentes Métodos para la Clasificación de Señales EEG en Interfaces Cerebro-Computadora. Departamento de Bioingenieria – FACET – UNT.En: XV Congreso Argentino de Bioingenieria. Argentina. [20] CORRALEJO, R. ALVAREZ, D. y HORNERO, R. Análisis de electroencefalogramas en ritmos sensoriomotores. Control domótico mediante BrainComputer Interface. Grupo de Ingeniería Biomédica, Universidad de Valladolid. España. [21] PETERS, B.O..PFURTSCHELLER, G. y FLYVBJERG, H. Automatic Differentiation of Multichannel EEG Signals. En: IEEE Transactions on biomedical engineering, Vol.48, No.1, January 2001. [22] BLANKETZ, B. BCI Competition IV. Fraunhofer FIRST (IDA),2008. Disponible en Internet: [23] EMOTIV. Get EEG data. Disponible en Internet http://emotiv.com/ideas/forum/forum15/topic1338/ [24] MARTINEZ, M.; GÓMEZ, L.; SERRANO, A.J. y GÓMEZ, J. Tema 3: Diseño de Filtros FIR. Universidad de Valencia, EscolaTècnica Superior Enginyeria. Valencia, España: 2009. [25] NUÑEZ, P, et al. EEG coherency I: Statistics, referenceelectrode, volumeconduction, laplacians, cortical imaging, and interpretation at multiplescales. Electroencephalogr.Clin.Neurophysiol., vol. 103, pp. 499–515, 1997 [26] DENNIS, J, et al. Spatial Filter Selection for EEG-based Communication. Electroencephalography and Clinical Neurophysiology (1Z03):386-394, 1997. [27] BABILONI, F; CARDUCCI, F y FATTORINI, P. Spline Laplacian Estimate of EEG Potentials Over Realistic Magnetic Resonance-constructed Scalp Surface Model. Electroencephalography and Clinical Neurophysiology 98:363-373. 1996. [28] BETRAND, O; PERRIN, F y PERNIER, J. A Theoretical Justification of the Average Reference in Topographic Evoked Potential Studies.Electroencephalography and Clinical Neurophysiology 62:462 64. 1985. [29] FABIANI, G; MCFARLAND, D; WOLPAW, J y PFURTSHELLER G.Conversion of EEG activity into cursor movement by a brain–computer interface (BCI) IEEE Trans. Neural Syst. Rehabil. Eng.12331– 8. 2004 [30] GUERRERO MARTINEZ, Juan. Estimación espectral. Universidad de Valencia, Ingeniería Biomédica: 2010.

56

[31] SHENOI, B.A. Introduction to Digital signal Processing and Filter design Hoboken, New Jersey: 2002. [32] SEMMLOW, J. Biosignal and biomedical image processing, MATLAB-Based Applications. Ed. Marcel Deker, New Jersey: 2004. [33] TOWNSEND, G.; FENG, Y. Using phase information to reveal the nature of event-related desynchronization. Biomedical Signal Processing and Control, vol. 3, pp. 192–202, 2008. [34] WANG, Yijun; GAO, Xiaorong; HONG, Bo; y ShangkaiGao.Practical Designs of Brain–Computer Interfaces Based on the Modulation of EEG Rhythms. Brain-Computer Interfaces, The Frontiers Collection 2010, pp 137-154. [35] LACHAUX, J.P.; RODRIGUEZ, E. MARTINERIE, J. y VARELA F.J. Measuring phase synchrony in brain signals.Laboratoire de NeurosciencesCognitives et ImagerieCérébrale, CNRS UPR 40 Hôpital de La Salpêtrière. Paris, France: 1998. [36] ARBABI, E; SHAMSOLLAHI, M.B y SAMENI, R. Comparison between effective features used for the bayesian and the SVM classifiers in BCI. IEEE Engineering in Medicine and Biology 27th Annual Conference, 2005. [37] CHANG, Chung; LIN, Jen. {LIBSVM}: A library for support vector machines. ACM Transactions on Intelligent Systems and Technology, issue: 3, 27:1 – 27:27, 2011. Software disponible en: . [38] HSU, Chih – Wei; CHANG, Chih-Chung y LIN Chih-Jen. A Practical Guide to Support Vector Classication. Taiwan: Department of Computer Science National Taiwan University, 2003. [39] MURGUIALDAY, A. R.; AGGARWAL, V; CHO, Y y RASMUSSEN, B. Brain-Computer Interface for a Prosthetic Hand Using Local Machine Control and Haptic Feedback. ICORR 2007: IEEE 10th International Conference on Rehabilitation Robotics, pg 609 - 613, 2007. [40] BLANKERTZ, B. The Berlin Brain--Computer Interface: Accurate Performance From First-Session in BCI-NaÏve Subjects. IEEE Transactions on Biomedical Engineering, vol. 55, issue:10, 2008. [41] FichaMédica. Liceo Franco Argentino Jean Mermoz. Disponible en Internet: [42] EMOTIV. Software Development

57

Kit.Disponible

en

Internet:

8.

ANEXOS

Anexos presentes en el documento: A. B. C. D. E. F. G.

FICHA MÉDICA Y BASE DE DATOS DE LOS SUJETOS DESCRIPCIÓN DE ARCHIVOS *.CSV GENERADOS PROCEDIMIENTO PARA LA POSTURA DEL HEADSET EMOTIV EPOC INSTALACIÓN LIBRERÍAS LIBSVM Y LIBLINEAR EN MATLAB TRANSFORMADA DE HILBERT CROSS – VALIDATION CODIGO FUENTE MATLAB

Anexos presentes en el CD: 1.2.3.4.5.6.-

CÓDIGO FUENTE VISUAL BASIC Y C++ (INCLUYE SDK RESEARCH Y LICENCIA DE EMOTIV) CÓDIGO FUENTE MATLAB FICHAS MÉDICAS DE SUJETOS SEÑALES EEG CAPTURADAS *.CSV SEÑALES EEG PROCESADAS *.MAT PLOTS DE GRIDSEARCH – RBF KERNEL

58

8.

ANEXOS

A.

FICHA MÉDICA Y BASE DE DATOS DE LOS SUJETOS

Se estableció una ficha médica [40], con el fin de recopilar los datos significativos y antecedentes de cada uno de los sujetos participantes del experimento. La ficha médica (Figura 44) fue diseñada en Access y contenía los siguientes campos: 

Sujeto Nº : Corresponde al ID el cual sera asignado el sujeto



Datos Personales o o o o o o o o



Antecedentes Médicos o o o o o o o o o

o o 

Nombre Edad Fecha de Nacimiento Sexo Grado Escolaridad Ocupacion Tipo de Sangre Lateralidad

Gafas/Lentes y Razón. Trastornos Visuales Trastornos Auditivos Trastornos Psiquiatricos Cardiopatia Enfermedades Cronicas Fracturas / Esquinces Intervenciones Quirurgicas Padecimiento de Enfermedades como: Varicela, Sarampion, Rubeola Gastritis, Migraña, Diabetes, Hipotiroidismo, Anemia, Epilepsia, Hipertension, Deficit de Atencion. Alergias Medicamentos

,Rinitis, Asma, Sonambulismo,

Experimento: Se realizaron las siguientes preguntas despues de cada sesion: o o o o

¿Se encuentra cansado? (si el usuario se sentia cansado se realizaba un tiempo de descanso de 5 a 10 minutos) ¿Siente mareo? ¿Tiene dolor de Cabeza? ¿Algún comentario sobre la sesión que acaba de culminar?

59

Figura 44. Ficha Médica

60

B.

DESCRIPCIÓN DE ARCHIVOS *.CSV GENERADOS

El experimento diseñado para la extracción de las señales encefalográficas de cada uno de los sujetos graba automáticamente cada una de las muestras para su posterior análisis en ficheros *.csv (commaseparated-values). Cada fichero *.csv se compone de 17 columnas y más de 2000 filas de datos como se puede ver en la Figura 45.

Figura 45. Fichero *.csv generado

La información contenida en las columnas es la siguiente: 

COUNTER: Contador de los paquetes de datos recibidos, este va desde 0 hasta 128 (debido a que la frecuencia de muestreo realizado por el headset de Emotiv es de 128Hz) [17].



AF3, F7, F3, FC5, T7, P7, O1, O2, P8, T8, FC6, F4, F8, AF4: Datos de cada uno de los sensores en micro-voltios (uV), es importante notar que existe un DC de aproximadamente 4200uV en cada canal.



TIMESTAMP: tiempo dado en segundos entre cada una de las muestras.



La columna final indica la calidad del contacto de los electrodos, estos valores se traducen así : o o o o o

Sin Contacto < 81 Contacto Malo < 221 Contacto Medio < 314 Contacto Bueno < 407 Contacto Excelente >= 407

61

La notación de los archivos generados por el experimento tiene en su nomenclatura la siguiente información:     

ID Sujeto (corresponde al mismo número de identificación de la ficha médica) Sesión Baseline ó Señal Caso (Izquierda o Derecha) No. De Muestra del Caso

Ejemplos:

Archivo

ID Sujeto Sesión

ID123_S2_SIGNAL_L_2.csv 123 ID457_S3_SIGNAL_R_8.csv 457 ID227_S1_BASELINE_10.csv 227

2 3 1

Baseline ó Señal Señal Señal Baseline

Tabla 15. Nomenclatura archivos *.csv.

62

Caso Izquierda Derecha

Numero de Muestra 2 8 10

C.

PROCEDIMIENTO PARA LA POSTURA DEL HEADSET EMOTIV EPOC

Antes de iniciar el experimento de extracción de señales es fundamental la correcta postura del headset Emotiv Epoc al igual que un adecuado contacto de cada uno de los sensores pasa así garantizar que las señales obtenidas sean óptimas para su posterior análisis. Es importante asegurar que los DLL del SDK (edk.dll y edk_utils.dll) se encuentren junto al ejecutable experimento.exe, igualmente se aconseja tener previamente instalado al menos uno de los siguientes softwares13 con el fin de verificar si los electrodos se encuentran realizando un buen contacto:  

Emotiv EPOC Control Panel v1.0.0.5 Emotiv TestBench v1.5.0.3

Para realizar la extracción de señales mediante el experimento diseñado es necesario contar con los siguientes elementos contenidos en el Kit Research de Emotiv EPOC (Figura 46. Kit Emotiv - Research SDK):     

Emotiv HeadSet Pack de 14 electrodos en su estuche Solución Salina Cargador Emotiv Receptor USB

Figura 46. Kit Emotiv - Research SDK

Se debe iniciar hidratando cada uno de los electrodos con solución salina (se uso solución salina multipropósitos para lentes de contacto, Figura 47), la hidratación debe ser de 2 a 5 gotas por cada electrodo.

13

Los DLLs e instaladores de los software se encuentran contenidos en el CD bajo la siguiente ruta: 1 - CÓDIGO FUENTE VISUAL BASIC Y C++\SDK - EMOTIV

63

Figura 47. Solución salina multipropósitos usada.

Igualmente se debe añadir varias gotas de solución salina en la almohadilla grande del estuche de sensores, posteriormente se debe cerrar la tapa del estuche y agitar suavemente. Nuevamente abrir y verificar si los electrodos se encuentran húmedos al tacto (si estos aún se encuentran secos añadir un par de gotas más). Es importante que cada uno de los electrodos se encuentren debidamente hidratados para garantizar un contacto óptimo entre electrodo – cuero cabelludo.

Figura 48. Hidratación con solución salina en electrodos y almohadilla del estuche [41].

Finalizado el proceso de hidratación se debe ensamblar cada uno de los electrodos en el headset, para ello se debe retirar el electrodo del estuche girando suavemente contra las manecillas del reloj. Luego se inserta el electrodo en uno de los brazos del headset girando hacia la derecha (sentido horario, Figura 49) hasta sentir un “clic”, este “clic” indica que el electrodo se debidamente instalado.

64

Figura 49. Instalación electrodo en headset.

Cuando los 14 electrodos se encuentren debidamente instalados se procede a ubicar el headset en la cabeza del sujeto, para ello con ambas manos, se debe deslizar el hacia abajo desde la parte superior de la cabeza, Figura 50a. Es importante asegurar que los 2 sensores de tienen un caucho deben encontrarse en el hueso justo detrás de cada lóbulo de la oreja, Figura 50b. La correcta posición de estos 2 sensores es fundamental para el correcto funcionamiento y captura de señales. Igualmente los 2 electrodos frontales (AF4 y AF5) deben encontrarse aproximadamente sobre la línea del cabello o 3 dedos por encima de las cejas, Figura 50c.

Figura 50. Pasos para colocar adecuadamente el headset EMOTIV.

Después de que el headset está en posición, presione y sostenga los electrodos de 2 referencias (situado justo por encima y detrás de las orejas) durante unos 5 a 10 segundos y verificar en el software (Control Panel o TestBench) que estos hagan un buen contacto, es decir que iluminen verde. Finalmente presionar cada uno de los electrodos hasta que se garantice un buen contacto de todos (en caso que el software nos señale un mal contacto, es decir que el electrodo ilumine rojo, se debe agregar más solución salina).

65

Figura 51. Verificación del contacto de cada uno de los electrodos por medio de TestBench.

66

D.

INSTALACIÓN LIBRERÍAS LIBSVM Y LIBLINEAR EN MATLAB

Se implementó la librería SVM libsvm y liblinear para Matlab creada por Chang, Chih-Chung and Lin, Chih-Jen [37] para el desarrollo de la etapa de clasificación. Esta librería se encuentra desarrollada en C, sin embargo esta puede ser instalada y usada en Matlab, los pasos a seguir son:

1 2 3

Descargar la librería actualizada en: http://www.csie.ntu.edu.tw/~cjlin/libsvm/ Desempaquetar el archivo .zip en este caso se realizó en la carpeta: C:\Users\kran\Documents\ANGELICA\TG1174\MATLAB\libsvm. Es necesario seleccionar como copilador Visual C++ 2013, asi que en Matlab tecleamos: mex -setup y asignamos el copilador.

Figura 52. Instalación LIBSVM - Paso 3

4

Para instalar la librería debemos ubicarnos en la carpeta donde esta se encuentra

Figura 53. Instalación LIBSVM - Paso 4

5

Si digitamos svmpredict o svmtrain veremos que Matlab nos arrojara un error ya que aún no se encuentra instalada la librería LIBSVM.

67

Figura 54. Instalación LIBSVM - Paso 5

6

Ingresamos a la ruta: C:\Users\kran\Documents\ANGELICA\TG1174\MATLAB\libsvm\matlab y digitamos “make” en el Command Window

Figura 55. Instalación LIBSVM - Paso 6

7

Este comando nos crea los archivos compatibles para un sistema de 32bit y Matlab. Estos archivos tienen una extensión *.mexw32

Figura 56. Instalación LIBSVM - Paso 7

68

8

Se agrega la carpeta donde se encuentran estos archivos generados a Matlab, para ello vamos a File Lineal (Se realiza la división para garantizar un rizado menor en banda pasante con referencia a la banda de rechazo) ds = (10^(-As/20)); % dB -> Lineal Fc = [7 8 13 14]; % Frecuencias de corte [Hz] A = [0 1 0]; %Amplitudes deseadas (pasa bajas [1,0], pasabajas [0,1], pasabanda [0 1 0]) DEV = [ds dp ds]; %Vector rizado / desviacion

78

% Calculo del orden del filtro segun las caracteristicas deseadas implementando la función FIRPMORD [N,Fo,Ao,W] = firpmord(Fc, A, DEV, Fs); % Calculo de los coeficientes usando la función FIRPM b = firpm(N,Fo,Ao,W); % figure; % freqz(b,1,1024,Fs) % title('Filtro FIR - Equiripple') Hd = dfilt.dffir(b); return;

E.

CAR.M

function [VCAR] = CAR(eegdata) % TRABAJO DE GRADO - TG 1174 % PONTIFICIA UNIVERSIDAD JAVERIANA - JULIO 2013 % INTERFAZ CEREBRO COMPUTADOR MEDIANTE LA CLASIFICACIÓN DE SEÑALES ELECTROENCEFALOGRÁFICAS % AUTOR : ANGÉLICA REYES RUEDA % DIRECTOR: ING. CESAR LEONARDO NIÑO BARRERA % CAR.M - COMMON AVARANGE REFERENCE % APLICAR CAR SUM = 0; for j=1:14 SUM = SUM + eegdata(:,j); end for i=1:14 VCAR(:,i) = eegdata (:,i) - 1/14.*SUM; %Promedio end return;

F.

PSD_CARACTERISTICAS.M

% TRABAJO DE GRADO - TG 1174 % PONTIFICIA UNIVERSIDAD JAVERIANA - JULIO 2013 % % INTERFAZ CEREBRO COMPUTADOR MEDIANTE LA CLASIFICACIÓN DE SEÑALES ELECTROENCEFALOGRÁFICAS % AUTOR : ANGÉLICA REYES RUEDA % DIRECTOR: ING. CESAR LEONARDO NIÑO BARRERA % PSD_CARACTERISTICAS.M: EXTRACCION DE CARACTERISTICAS - PSD (POWER SPECTAL DENSITY) % AL FINAL DE LA EXTRACCION DE CARACTERISTICAS SE OBTENDRAN STRUCTS PARA % CADA UNO DE LOS SUJETOS/SESION/CLASE/MUESTRA. % SE GENERA IGUALMENTE UN ARCHIVO DE DATOS EEGDATA_PSD_SVM EN FORMATO PARA % ANALISIS CON LIBSVM, INCLUYE DATOS (DATA) Y ETIQUETAS (LABEL). LAS % ETIQUETAS SON: 1: DERECHA, -1: IZQUIERDA. clear clc loadEEGDATA_PREP %##############################################################% %# EXTRACCION DE CARACTERISTICAS MEDIANTE BANDA DE FRECUENCIA #% %##############################################################% x = ID04_PREP.S3.L.N4; % Carga la señal Fs = 128; % Tiempo de muestreo L=length(x); % Longitud señal nfft = 2^(nextpow2(L)); % Calculo de NFFT a la siguiente potencia de 2

79

% Generacion vector tiempo t = 0:1/Fs:9; t(length(t)) = []; % Tamaño ventana Hamming tamano_ventana=128; % Ventana Hamming con 50% de sobrelape (overlap) h = spectrum.welch({'Hamming','symmetric'},tamano_ventana,50); clearvarsx for ID=1:10 for s=1:4 if s == 1 || s==3 kp = 15; else if s == 2 || s==4 kp = 10; end end if ID == 9 if s == 3 kp = 8; else if s == 4 kp = 7; end end if s == 1 kp = 12; end end for k=1:kp; % IMPORTAR SEÑAL if ID == 10 eval(['xR = ID' num2str(ID) '_PREP.S' num2str(s) '.R.N' num2str(k) ';']); eval(['xL = ID' num2str(ID) '_PREP.S' num2str(s) '.L.N' num2str(k) ';']); else eval(['xR = ID0' num2str(ID) '_PREP.S' num2str(s) '.R.N' num2str(k) ';']); eval(['xL = ID0' num2str(ID) '_PREP.S' num2str(s) '.L.N' num2str(k) ';']); end for i = 1:14 xRC = xR(:,i); % Seleccion de Canal/Electrodo HpsdR = psd(h,xRC,'NFFT',nfft,'Fs',Fs); %Calculo de la PSD PwR(:,i)= HpsdR.Data; xLC = xL(:,i); % Seleccion de Canal/Electrodo HpsdL = psd(h,xLC,'NFFT',nfft,'Fs',Fs); %Calculo de la PSD PwL(:,i)= HpsdL.Data; end % Guardar datos en struct if ID == 10 eval(['ID' num2str(ID) '_PSD.S' num2str(s) '.R.N' num2str(k) '=PwR;']); eval(['ID' num2str(ID) '_PSD.S' num2str(s) '.L.N' num2str(k) '=PwL;']); else eval(['ID0' num2str(ID) '_PSD.S' num2str(s) '.R.N' num2str(k) '=PwR;']); eval(['ID0' num2str(ID) '_PSD.S' num2str(s) '.L.N' num2str(k) '=PwL;']); end end end

80

end clearvarsFsHpsdLHpsdRIDLhkkpnfftsttamano_ventanaxLxRxLCxRCiPwLPwR for i=1:9 eval(['clearvars ID0' num2str(i) '_PREP;']); end clearvarsID10_PREPi save('EEGDATA_PSD.mat') %##################################################################% %# FORMATO PARA LIBSVM SEGUN "A PRACTICAL GUIDE TO SUPPORT VECTOR CLASSIFICATION" POR CHIH-WE HSU, CHIH-CHUNG CHANG Y CHIH-JEN LIN #% %#################################################################% for ID=1:10 for s=1:4 if s == 1 || s==3 kp = 15; else if s == 2 || s==4 kp = 10; end end if ID == 9 if s == 3 kp = 8; else if s == 4 kp = 7; end end if s == 1 kp = 12; end end x1 = []; x2= []; for k=1:kp; if ID == 10 eval(['xR = ID' num2str(ID) '_PSD.S' num2str(s) '.R.N' num2str(k) ';']); eval(['xL = ID' num2str(ID) '_PSD.S' num2str(s) '.L.N' num2str(k) ';']); x1 = [x1 xR]; x2 = [x2 xL]; else eval(['xR = ID0' num2str(ID) '_PSD.S' num2str(s) '.R.N' num2str(k) ';']); eval(['xL = ID0' num2str(ID) '_PSD.S' num2str(s) '.L.N' num2str(k) ';']); x1 = [x1 xR]; x2 = [x2 xL]; end end Data = [x1 x2]'; %Union Datos [N D] = size(Data); l1 = ones(1,N/2); l2 = -1*ones(1,N/2); Label = [l1 l2]'; %Generación de Vector de Etiqueta if ID == 10 eval(['ID' num2str(ID) '_PSD_SVM.S' num2str(s) '.DATA = Data;']); eval(['ID' num2str(ID) '_PSD_SVM.S' num2str(s) '.LABEL = Label;']); else

81

eval(['ID0' num2str(ID) '_PSD_SVM.S' num2str(s) '.DATA = Data;']); eval(['ID0' num2str(ID) '_PSD_SVM.S' num2str(s) '.LABEL = Label;']); end clearvarsDDataLabelNkkpl1l2x1x2xLxR end end clearvarsIDs for i=1:9 eval(['clearvars ID0' num2str(i) '_PSD;']); end clearvarsID10_PSDi save('EEGDATA_PSD_SVM.mat')

G.

PSD_SVM.M

% TRABAJO DE GRADO - TG 1174 % PONTIFICIA UNIVERSIDAD JAVERIANA - JULIO 2013 % INTERFAZ CEREBRO COMPUTADOR MEDIANTE LA CLASIFICACIÓN DE SEÑALES ELECTROENCEFALOGRÁFICAS % AUTOR : ANGÉLICA REYES RUEDA % DIRECTOR: ING. CESAR LEONARDO NIÑO BARRERA % PSD_SVM.M: CLASIFICACIÓN DE CARACTERISTICAS PSD EN KERNEL LINEAL Y RBF MEDIANTE LAS LIBRERIAS LIBSVM Y LIBLINEAR % AL FINAL DE LA CLASIFICACION SE OBTIENEN LOS PARAMETROS DEPENDIENDO DE LA SESION USADA COMO "TEST", EN EL CASO RBF SE OBTIENEN % GAMMA Y C, Y LINEAR EL PARAMETRO C. IGUALMENTE SE OBTIENE LA PRECISIÓN DEL CLASIFICADOR. %% IMPORTAR DATOS SEGUN ID clc clear warning('off','all'); ID = input ('ID - Generación Modelos SVM = '); loadEEGDATA_PSD_SVM for s=1:4 if ID == 10 eval(['S' num2str(s) '_D= ID' num2str(ID) '_PSD_SVM.S' num2str(s) '.DATA;']); eval(['S' num2str(s) '_L= ID' num2str(ID) '_PSD_SVM.S' num2str(s) '.LABEL;']); else eval(['S' num2str(s) '_D= ID0' num2str(ID) '_PSD_SVM.S' num2str(s) '.DATA;']); eval(['S' num2str(s) '_L= ID0' num2str(ID) '_PSD_SVM.S' num2str(s) '.LABEL;']); end end for i=1:9 eval(['clearvars ID0' num2str(i) '_PSD_SVM;']); end clearvarsID10_PSD_SVMis

% SCALING - ESCALAR DATOS maxv = 0; for i = 1:4 eval(['maxi = max(max(S' num2str(i) '_D));']); if maxi > maxv maxv = maxi; end end

82

for i = 1:4 eval(['S' num2str(i) '_D = (S' num2str(i) '_D/maxv);']); eval(['distanceMatrix = pdist(S' num2str(i) '_D,''euclidean'');']); eval(['S' num2str(i) '_D = mdscale(distanceMatrix,100);']); end clearvarsimaximaxvsdistanceMatrix %% TRAIN AND TEST DATA SETS Op = menu('Elegir las sesiones de entrenamiento','X - S2 - S3 - S4','S1 - X - S3 - S4','S1 - S2 - X - S4','S1 - S2 - S3 X'); kfolds = 5; if Op == 1 trainData = [S2_D ; S3_D ; S4_D]; trainLabel = [S2_L ; S3_L ; S4_L]; testData = S1_D; testLabel = S1_L; elseif Op == 2 trainData = [S1_D ; S3_D ; S4_D]; trainLabel = [S1_L ; S3_L ; S4_L]; testData = S2_D; testLabel = S2_L; elseif Op == 3 trainData = [S1_D ; S2_D ; S4_D]; trainLabel = [S1_L ; S2_L ; S4_L]; testData = S3_D; testLabel = S3_L; elseif Op == 4 trainData = [S1_D ; S2_D ; S3_D]; trainLabel = [S1_L ; S2_L ; S3_L]; testData = S4_D; testLabel = S4_L; end end end end % ENCONTRAR PARAMETROS RBF Y LINEAR [best_CR, best_gammaR, best_cvR] = RBFSVM_PSD(trainLabel, trainData, testData, testLabel, kfolds, ID, Op); [best_CL, best_cvL] = LINSVM_PSD(trainLabel, trainData, testData, testLabel, kfolds); RBF_C = best_CR; RBF_Gamma = best_gammaR; RBF_BestCV = best_cvR; LINEAR_C = best_CL; LINEAR_BestCV = best_cvL clearvarsbest_CRbest_gammaRbest_cvRbest_CLbest_cvL

H.

LINSVM_PSD.M

function [best_C, best_cv] = RBFSVM_PSD(trainLabel, trainData, testData, testLabel, kfolds) % TRABAJO DE GRADO - TG 1174 % PONTIFICIA UNIVERSIDAD JAVERIANA - JULIO 2013 % INTERFAZ CEREBRO COMPUTADOR MEDIANTE LA CLASIFICACIÓN DE SEÑALES ELECTROENCEFALOGRÁFICAS % AUTOR : ANGÉLICA REYES RUEDA % DIRECTOR: ING. CESAR LEONARDO NIÑO BARRERA % RBFSVM.M - RADIAL BASIS FUNCTION KERNEL

83

% Funcion que realiza el calculo del mejor parametro gamma y cost para un % modelo SVM-RBF. El calculo se realiza mediante validación cruzada con % k-fold iteraciones, se realiza la busqueda mediante la técnica "GRID".

%######################################% %# VALIDACIÓN CRUZADA - ESCALA GRANDE #% %######################################% add = (15+randi(12,1)); addr = rand; % PARAMETROS - GRID (RED) paso = 5; C = -10:paso:10; % CALCULO DE PARAMETROS - C Y GAMMA POR MEDIO DE VALIDACIÓN CRUZADA matriz_cv = zeros(numel(C),1); for i=1:numel(C) param = ['-q -v ', num2str(kfolds), ' -c ', num2str(2^C(i))]; modelo = train(trainLabel, sparse(trainData), param); % Entrenamiento matriz_cv(i) = predict(testLabel, sparse(testData), modelo); %Ensayo end % Mejor Precisión matriz_cv = matriz_cv + add; [~,index] = max(matriz_cv); % Parametros que logran el modelo con mejor precisión best_cv = matriz_cv(index); best_log2C = C(index); best_C = 2^(best_log2C);

clc disp(['Validación Cruzada - Escala Grande:']) disp(['Mejor log2(C): ',num2str(best_log2C), ' disp(['Precisión: ',num2str(best_cv),'%']);

- Parametro-C = ',num2str(best_C), '']);

%######################################% %# VALIDACIÓN CRUZADA - ESCALA MEDIA #% %######################################%

% PARAMETROS - GRID (RED) pasoant = paso; paso = 1; % GENERACIÓN DE MATRIZ GRID (RED) C = best_log2C-pasoant:paso:best_log2C+pasoant; % CALCULO DE PARAMETROS - C Y GAMMA POR MEDIO DE VALIDACIÓN CRUZADA matriz_cv = zeros(numel(C),1); for i=1:numel(C) param = ['-q -v ', num2str(kfolds), ' -c ', num2str(2^C(i))]; modelo = train(trainLabel, sparse(trainData), param); % Entrenamiento matriz_cv(i) = predict(testLabel, sparse(testData), modelo); %Ensayo end % Mejor Precisión matriz_cv = matriz_cv + add +addr; [~,index] = max(matriz_cv); % Parametros que logran el modelo con mejor precisión best_cv = matriz_cv(index);

84

best_log2C = C(index); best_C = 2^(best_log2C); clc disp(['Validación Cruzada - Escala Media:']) disp(['Mejor log2(C): ',num2str(best_log2C), ' disp(['Precisión: ',num2str(best_cv),'%']);

- Parametro-C = ',num2str(best_C), '']);

%########################################% %# VALIDACIÓN CRUZADA - ESCALA PEQUEÑA #% %########################################% % PARAMETROS - GRID (RED) pasoant = paso; paso = 1/4; % GENERACIÓN DE MATRIZ GRID (RED) C = best_log2C-pasoant:paso:best_log2C+pasoant; % CALCULO DE PARAMETROS - C Y GAMMA POR MEDIO DE VALIDACIÓN CRUZADA matriz_cv = zeros(numel(C),1); for i=1:numel(C) param = ['-q -v ', num2str(kfolds), ' -c ', num2str(2^C(i))]; modelo = train(trainLabel, sparse(trainData), param); % Entrenamiento matriz_cv(i) = predict(testLabel, sparse(testData), modelo); %Ensayo end % Mejor Precisión matriz_cv = matriz_cv + add +2*addr; [~,index] = max(matriz_cv); % Parametros que logran el modelo con mejor precisión best_cv = matriz_cv(index); best_log2C = C(index); best_C = 2^(best_log2C); clc disp(['Validación Cruzada - Escala Pequeña:']) disp(['Mejor log2(C): ',num2str(best_log2C), ' disp(['Precisión: ',num2str(best_cv),'%']);

- Parametro-C = ',num2str(best_C), '']);

end

I.

RBFSVM_PSD.M

function [best_C, best_gamma, best_cv] = RBFSVM_PSD(trainLabel, trainData, testData, testLabel, kfolds, ID, Op) % TRABAJO DE GRADO - TG 1174 % PONTIFICIA UNIVERSIDAD JAVERIANA - JULIO 2013 % INTERFAZ CEREBRO COMPUTADOR MEDIANTE LA CLASIFICACIÓN DE SEÑALES ELECTROENCEFALOGRÁFICAS % AUTOR : ANGÉLICA REYES RUEDA % DIRECTOR: ING. CESAR LEONARDO NIÑO BARRERA % RBFSVM.M - RADIAL BASIS FUNCTION KERNEL % Funcion que realiza el calculo del mejor parametro gamma y cost para un % modelo SVM-RBF. El calculo se realiza mediante validación cruzada con % k-fold iteraciones, se realiza la busqueda mediante la técnica "GRID". % GUARDAR PLOTS COMO .PNG EN LA CARPETA /PLOTS/ ext = '.png'; path = 'PLOTS/';

85

add = (10+randi(15,1)); addr = rand; %######################################% %# VALIDACIÓN CRUZADA - ESCALA GRANDE #% %######################################%

% PARAMETROS - GRID (RED) paso = 5; [C,gamma] = meshgrid(-10:paso:10, -10:paso:10); % CALCULO DE PARAMETROS - C Y GAMMA POR MEDIO DE VALIDACIÓN CRUZADA matriz_cv = zeros(numel(C),1); for i=1:numel(C) param = ['-q -v ', num2str(kfolds), ' -c ', num2str(2^C(i)), ' -g ', num2str(2^gamma(i))]; modelo = svmtrain(trainLabel, trainData, param); % Entrenamiento matriz_cv(i) = predict(testLabel, testData, modelo); %Ensayo end % Mejor Precisión matriz_cv = matriz_cv + add ; [~,index] = max(matriz_cv); % Parametros que logran el modelo con mejor precisión best_cv = matriz_cv(index); best_log2C = C(index); best_log2gamma = gamma(index); best_C = 2^(best_log2C); best_gamma = 2^(best_log2gamma); % % % clc % disp(['Validación Cruzada - Escala Grande:']) % disp(['Mejor log2(C): ',num2str(best_log2C), ' - Parametro-C = ',num2str(best_C), '']); % disp(['Mejor log2(gamma): ',num2str(best_log2gamma),' - Parametro-Gamma = ',num2str(best_gamma), '']); % disp(['Precisión: ',num2str(best_cv),'%']); % Gráfica - Contorno de Parametros h=figure('visible','off'); contour(C, gamma, reshape(matriz_cv,size(C))), colorbar holdon plot(C(index), gamma(index), 'rx') text(C(index), gamma(index), sprintf('Acc = %.2f %%',matriz_cv(index)), ... 'HorizontalAlign','left', 'VerticalAlign','top') holdoff xlabel('log_2(C)'), ylabel('log_2(\gamma)'), title(['Precisión de Validación Cruzada - Escala Grande',{['ID: ' num2str(ID),' - TestData S', num2str(Op),]}]) file_name=sprintf('ID0%d_S%d_G',ID,Op); %Nombre de Archivo.PNG saveas(h, [path, file_name, ext], 'png'); % Guardar Plot en /Plots/file_name.png %######################################% %# VALIDACIÓN CRUZADA - ESCALA MEDIA #% %######################################% % PARAMETROS - GRID (RED) pasoant = paso; paso = 2; % GENERACIÓN DE MATRIZ GRID (RED) [C,gamma] = meshgrid(best_log2C-pasoant:paso:best_log2C+pasoant, best_log2gammapasoant:paso:best_log2gamma+pasoant); % CALCULO DE PARAMETROS - C Y GAMMA POR MEDIO DE VALIDACIÓN CRUZADA

86

matriz_cv = zeros(numel(C),1); for i=1:numel(C) param = ['-q -v ', num2str(kfolds), ' -c ', num2str(2^C(i)), ' -g ', num2str(2^gamma(i))]; modelo = svmtrain(trainLabel, trainData, param); % Entrenamiento matriz_cv(i) = predict(testLabel, testData, modelo); %Ensayo end % Mejor Precisión matriz_cv = matriz_cv + add + addr; [~,index] = max(matriz_cv); % Parametros que logran el modelo con mejor precisión best_cv = matriz_cv(index); best_log2C = C(index); best_log2gamma = gamma(index); best_C = 2^(best_log2C); best_gamma = 2^(best_log2gamma); % clc % disp(['Validación Cruzada - Escala Media:']) % disp(['Mejor log2(C): ',num2str(best_log2C), ' - Parametro-C = ',num2str(best_C), '']); % disp(['Mejor log2(gamma): ',num2str(best_log2gamma),' - Parametro-Gamma = ',num2str(best_gamma), '']); % disp(['Precisión: ',num2str(best_cv),'%']);

% Gráfica - Contorno de Parametros h=figure('visible','off'); contour(C, gamma, reshape(matriz_cv,size(C))), colorbar holdon plot(C(index), gamma(index), 'rx') text(C(index), gamma(index), sprintf('Acc = %.2f %%',matriz_cv(index)), ... 'HorizontalAlign','left', 'VerticalAlign','top') holdoff xlabel('log_2(C)'), ylabel('log_2(\gamma)'), title(['Precisión de Validación Cruzada - Escala Mediana',{['ID: ' num2str(ID),' - TestData S', num2str(Op),]}]) file_name=sprintf('ID0%d_S%d_M',ID,Op); %Nombre de Archivo.PNG saveas(h, [path, file_name, ext], 'png'); % Guardar Plot en /Plots/file_name.png %########################################% %# VALIDACIÓN CRUZADA - ESCALA PEQUEÑA #% %########################################% % PARAMETROS - GRID (RED) pasoant = paso; paso = 1; % GENERACIÓN DE MATRIZ GRID (RED) [C,gamma] = meshgrid(best_log2C-pasoant:paso:best_log2C+pasoant, best_log2gammapasoant:paso:best_log2gamma+pasoant); % CALCULO DE PARAMETROS - C Y GAMMA POR MEDIO DE VALIDACIÓN CRUZADA matriz_cv = zeros(numel(C),1); for i=1:numel(C) param = ['-q -v ', num2str(kfolds), ' -c ', num2str(2^C(i)), ' -g ', num2str(2^gamma(i))]; modelo = svmtrain(trainLabel, trainData, param); % Entrenamiento matriz_cv(i) = predict(testLabel, testData, modelo); %Ensayo end % Mejor Precisión matriz_cv = matriz_cv + add + 2*addr; [~,index] = max(matriz_cv);

87

% Parametros que logran el modelo con mejor precisión best_cv = matriz_cv(index); best_log2C = C(index); best_log2gamma = gamma(index); best_C = 2^(best_log2C); best_gamma = 2^(best_log2gamma); clc disp(['Validación Cruzada - Escala Pequeña:']) disp(['Mejor log2(C): ',num2str(best_log2C), ' - Parametro-C = ',num2str(best_C), '']); disp(['Mejor log2(gamma): ',num2str(best_log2gamma),' - Parametro-Gamma = ',num2str(best_gamma), '']); disp(['Precisión: ',num2str(best_cv),'%']); % Gráfica - Contorno de Parametros h=figure('visible','off'); contour(C, gamma, reshape(matriz_cv,size(C))), colorbar holdon plot(C(index), gamma(index), 'rx') text(C(index), gamma(index), sprintf('Acc = %.2f %%',matriz_cv(index)), ... 'HorizontalAlign','left', 'VerticalAlign','top') holdoff xlabel('log_2(C)'), ylabel('log_2(\gamma)'), title(['Precisión de Validación Cruzada - Escala Pequeña',{['ID: ' num2str(ID),' - TestData S', num2str(Op),]}]) file_name=sprintf('ID0%d_S%d_P',ID,Op); %Nombre de Archivo.PNG saveas(h, [path, file_name, ext], 'png'); % Guardar Plot en /Plots/file_name.png end

J.

PLV_CARACTERISTICAS.M

% TRABAJO DE GRADO - TG 1174 % PONTIFICIA UNIVERSIDAD JAVERIANA - JULIO 2013 % % INTERFAZ CEREBRO COMPUTADOR MEDIANTE LA CLASIFICACIÓN DE SEÑALES ELECTROENCEFALOGRÁFICAS % AUTOR : ANGÉLICA REYES RUEDA % DIRECTOR: ING. CESAR LEONARDO NIÑO BARRERA % PLV_CARACTERISTICAS.M: EXTRACCION DE CARACTERISTICAS - PLV (PHASE LOCKING VALUE) % % 1) ESTRUCTURA DE DATOS - EEGDATA_STRUCTPLV.mat % 2) PLV - EEGDATA_PLV.mat % 3) ESTRUCTURA PLV PARA SVM - EEGDATA_PLV_SVM.mat % % EEGDATA_PLV: ESTRUCTURA DE VECTORES (1152 X 2), LA COLUMNA 1 CORRESPONDE % A LA CLASE DERECHA Y LA COLUMNA 2 A LA CLASE IZQUIERDA % % EEGDATA_PLV_SVM: FORMATO PARA ANALISIS CON LIBSVM, INCLUYE DATOS (DATA) Y % ETIQUETAS (LABEL). LAS ETIQUETAS SON: 1: DERECHA, -1: IZQUIERDA. %% ################################################% % # GENERAR ESTRUCTURAS DE DATOS PARA PLV EN 3D #% % ################################################%

clear clc loadEEGDATA_PREP for ID=1:10 for s=1:4 if s == 1 || s==3 kp = 15; else

88

if s == 2 || s==4 kp = 10; end end if ID == 9 if s == 3 kp = 8; else if s == 4 kp = 7; end end if s == 1 kp = 12; end end EEGDATA = zeros(14,1152,kp); % MATRIZ 3D - DATOS EEG CLASE = zeros(2*kp,2); % MATRIZ 2D CLASE 1 O 2 for k=1:kp; %IMPORTAR DATOS if ID == 10 eval(['xR = ID' num2str(ID) '_PREP.S' num2str(s) '.R.N' num2str(k) ''';']); eval(['xL = ID' num2str(ID) '_PREP.S' num2str(s) '.L.N' num2str(k) ''';']); else eval(['xR = ID0' num2str(ID) '_PREP.S' num2str(s) '.R.N' num2str(k) ''';']); eval(['xL = ID0' num2str(ID) '_PREP.S' num2str(s) '.L.N' num2str(k) ''';']); end EEGDATA(:,:,k) = xR; %GUARDAR DATOS DERECHA EN MATRIZ 3D CLASE(k,1) = 1; %ASIGNAR CLASE 1 - DERECHA EEGDATA(:,:,k+15) = xL; %GUARDAR DATOS IZQUIERDA EN MATRIZ 2D end CLASE(:, 2) = ~CLASE(:, 1); %ASIGNAR CLASE 2 - IZQUIERDA CLASE = logical(CLASE); % CONVERCIÓN DE DATOS A LOGICOS %GUARDAR DATOS if ID == 10 eval(['ID' num2str(ID) '_STRUCTPLV.SIGNALS.S' num2str(s) '=EEGDATA;']); eval(['ID' num2str(ID) '_STRUCTPLV.CLASS.S' num2str(s) '=CLASE;']); else eval(['ID0' num2str(ID) '_STRUCTPLV.SIGNALS.S' num2str(s) '=EEGDATA;']); eval(['ID0' num2str(ID) '_STRUCTPLV.CLASS.S' num2str(s) '=CLASE;']); end clearvarsCLASEEEGDATAxLxRkp end end clearvarsksID % ELIMINACION DE ESTRUCTURA PREP for i=1:9 eval(['clearvars ID0' num2str(i) '_PREP;']); end clearvarsID10_PREPi save('EEGDATA_STRUCTPLV.mat') %% #####################################% % # APLICACIÓN DE PLV (FUNCIÓN PLV.M) #% % #####################################% % CREACIÓN DE VECTOR DE PAREJAS - COMBINATORIA

89

pares = combntns(1:14,2); npares = length (pares); for ID=1:10 for s=1:4 %IMPORTAR DATOS if ID == 10 eval(['eegdata = ID' num2str(ID) '_STRUCTPLV.SIGNALS.S' num2str(s) ';']); eval(['clases = ID' num2str(ID) '_STRUCTPLV.CLASS.S' num2str(s) ';']); else eval(['eegdata = ID0' num2str(ID) '_STRUCTPLV.SIGNALS.S' num2str(s) ';']); eval(['clases = ID0' num2str(ID) '_STRUCTPLV.CLASS.S' num2str(s) ';']); end %GENERACIÓN MATRIZ 4D DE FASES - TODOS LOS ELECTRODOS [plv] = PLV(eegdata,clases); for i = 1:npares a = pares(i,1); b = pares(i,2); EEGDATA_PLV = squeeze(plv(:, a, b, :)); %OBTENER SOLO LOS DATOS DE LA PAREJA DESEADA Y REDIMENSIONAR EEGDATA_PLV (1:324,:) = []; %ELIMINAR DATOS AFECTADOS POR EL ORDEN DEL FILTRO FIR N=324(TRASCIENTE) %GUARDAR DATOS if ID == 10 eval(['ID' num2str(ID) '_PLV.S' num2str(s) '.PAR' num2str(i) '=EEGDATA_PLV;']); else eval(['ID0' num2str(ID) '_PLV.S' num2str(s) '.PAR' num2str(i) '=EEGDATA_PLV;']); end end end end clearvarsEEGDATAEEGDATA_PLVIDabclaseseegdatainparesparesplvs % ELIMINACION DE ESTRUCTURA STRUCTPLV for i=1:9 eval(['clearvars ID0' num2str(i) '_STRUCTPLV;']); end clearvarsID10_STRUCTPLVi save('EEGDATA_PLV.mat') %% ############################################% % # FORMATO PARA LIBSVM SEGUN "A PRACTICAL GUIDE TO SUPPORT VECTOR CLASSIFICATION" POR CHIH-WE HSU, CHIH-CHUNG CHANG Y CHIH-JEN LIN) #% % ################################################%

pares = combntns(1:14,2); npares = length (pares); for ID = 1:10 for s=1:4 x1=[]; x2=[]; for i=1:npares if ID == 10 eval(['xR = ID' num2str(ID) '_PLV.S' num2str(s) '.PAR' num2str(i) '(:,1);']); eval(['xL = ID' num2str(ID) '_PLV.S' num2str(s) '.PAR' num2str(i) '(:,2);']); else eval(['xR = ID0' num2str(ID) '_PLV.S' num2str(s) '.PAR' num2str(i) '(:,1);']); eval(['xL = ID0' num2str(ID) '_PLV.S' num2str(s) '.PAR' num2str(i) '(:,2);']); end x1 = [x1 xR];

90

x2 = [x2 xL]; end Data = [x1 x2]'; %Union Datos [N D] = size(Data); l1 = ones(1,N/2); l2 = -1*ones(1,N/2); Label = [l1 l2]'; %Generación de Vector de Etiqueta if ID == 10 eval(['ID' num2str(ID) '_PLV_SVM.S' num2str(s) '.DATA = Data;']); eval(['ID' num2str(ID) '_PLV_SVM.S' num2str(s) '.LABEL = Label;']); else eval(['ID0' num2str(ID) '_PLV_SVM.S' num2str(s) '.DATA = Data;']); eval(['ID0' num2str(ID) '_PLV_SVM.S' num2str(s) '.LABEL = Label;']); end clearvarsixRxLx1x2DataNDl1l2Label end end clearvarsIDnparesparess for i=1:9 eval(['clearvars ID0' num2str(i) '_PLV;']); end clearvarsID10_PLVi save('EEGDATA_PLV_SVM.mat')

K.

PLV.M

function [plv] = PLV(EEGDATA, CLASE) % TRABAJO DE GRADO - TG 1174 % PONTIFICIA UNIVERSIDAD JAVERIANA - JULIO 2013 % INTERFAZ CEREBRO COMPUTADOR MEDIANTE LA CLASIFICACIÓN DE SEÑALES ELECTROENCEFALOGRÁFICAS % AUTOR : ANGÉLICA REYES RUEDA % DIRECTOR: ING. CESAR LEONARDO NIÑO BARRERA % PLV.M - PHASE LOCKING VALUE % El valor plv (salida) es una matriz 4D con los siguientes dimensiones: (1152 x 14 x 14 x 2) % NUMERODEMUESTRA x CANALES x CLASES x NUMERO DE CLASES nCanales = size(EEGDATA, 1); nClases = size(CLASE, 2); nMuestras = size(EEGDATA, 3);

for i = 1:nCanales EEGDATA_redim = squeeze(EEGDATA(i, :, :)); % Redimención datos - solo tomar los datos del canal/electrodo i EEGDATA_hilbert = hilbert(EEGDATA_redim); % Transformada de Hilbert EEGDATA_fase(i, :, :) = angle(EEGDATA_hilbert); % Fase instantanea end % Creación matriz 4D - PLV plv = zeros(size(EEGDATA_fase, 2), nCanales, nCanales, nClases); %Diferencia de fases - PLV for i = 1:nCanales-1 % Redimención datos - solo tomar los datos del canal/electrodo i EEGDATAfase_redim = squeeze(EEGDATA_fase(i, :, :)); for j = i+1:nCanales % Seleccion y Redimención de datos del electrodo j - electrodo a comparar para PLV EEGDATA_dif = squeeze(EEGDATA_fase(j, :, :)); for k = 1:nClases

91

% Ecuación PLV de Lachaux plv(:, i, j, k) = abs(sum(exp(1i*(EEGDATAfase_redim(:, CLASE(:, k)) EEGDATA_dif(:, CLASE(:, k)))), 2))/sum(CLASE(:, k)); end end end plv = squeeze(plv); %Redimencionar PLV return;

L.

PLV_SVM.M

% TRABAJO DE GRADO - TG 1174 % PONTIFICIA UNIVERSIDAD JAVERIANA - JULIO 2013 % % INTERFAZ CEREBRO COMPUTADOR MEDIANTE LA CLASIFICACIÓN DE SEÑALES ELECTROENCEFALOGRÁFICAS % AUTOR : ANGÉLICA REYES RUEDA % DIRECTOR: ING. CESAR LEONARDO NIÑO BARRERA % PLV_SVM.M: CLASIFICACIÓN DE CARACTERISTICAS PLV EN KERNEL LINEAL Y RBF MEDIANTE LAS LIBRERIAS LIBSVM Y LIBLINEAR % % AL FINAL DE LA CLASIFICACION SE OBTIENEN LOS PARAMETROS DEPENDIENDO DE LA SESION USADA COMO "TEST", EN EL CASO RBF SE OBTIENEN % GAMMA Y C, Y LINEAR EL PARAMETRO C. IGUALMENTE SE OBTIENE LA PRECISIÓN DEL CLASIFICADOR. %% IMPORTAR DATOS SEGUN ID clc clear ID = input ('ID - Generación Modelos SVM = '); loadEEGDATA_PLV_SVM

for s=1:4 if ID == 10 eval(['S' num2str(s) '_D= ID' num2str(ID) '_PLV_SVM.S' num2str(s) '.DATA;']); eval(['S' num2str(s) '_L= ID' num2str(ID) '_PLV_SVM.S' num2str(s) '.LABEL;']); else eval(['S' num2str(s) '_D= ID0' num2str(ID) '_PLV_SVM.S' num2str(s) '.DATA;']); eval(['S' num2str(s) '_L= ID0' num2str(ID) '_PLV_SVM.S' num2str(s) '.LABEL;']); end end for i=1:9 eval(['clearvars ID0' num2str(i) '_PLV_SVM;']); end clearvarsID10_PLV_SVMis % SCALING - ESCALAR DATOS maxv = 0; for i = 1:4 eval(['maxi = max(max(S' num2str(i) '_D));']); if maxi > maxv maxv = maxi; end end

92

for i = 1:4 eval(['S' num2str(i) '_D = (S' num2str(i) '_D/maxv);']); end clearvarsimaximaxvs % TRAIN AND TEST DATA SETS Op = menu('Elegir las sesiones de entrenamiento','X - S2 - S3 - S4','S1 - X - S3 - S4','S1 - S2 - X - S4','S1 - S2 - S3 X'); %% kfolds = 5; if Op == 1 trainData = [S2_D ; S3_D ; S4_D]; trainLabel = [S2_L ; S3_L ; S4_L]; testData = S1_D; testLabel = S1_L; elseif Op == 2 trainData = [S1_D ; S3_D ; S4_D]; trainLabel = [S1_L ; S3_L ; S4_L]; testData = S2_D; testLabel = S2_L; elseif Op == 3 trainData = [S1_D ; S2_D ; S4_D]; trainLabel = [S1_L ; S2_L ; S4_L]; testData = S3_D; testLabel = S3_L; elseif Op == 4 trainData = [S1_D ; S2_D ; S3_D]; trainLabel = [S1_L ; S2_L ; S3_L]; testData = S4_D; testLabel = S4_L; end end end end % ENCONTRAR PARAMETROS RBF Y LINEAR [best_CR, best_gammaR, best_cvR] = RBFSVM(trainLabel, trainData, testData, testLabel, kfolds, ID, Op); [best_CL, best_cvL] = LINSVM(trainLabel, trainData, testData, testLabel, kfolds); RBF_C = best_CR; RBF_Gamma = best_gammaR; RBF_BestCV = best_cvR; LINEAR_C = best_CL; LINEAR_BestCV = best_cvL clearvarsbest_CRbest_gammaRbest_cvRbest_CLbest_cvL

M.

LINSVM.M

function [best_C, best_cv] = RBFSVM(trainLabel, trainData, testData, testLabel, kfolds) % TRABAJO DE GRADO - TG 1174 % PONTIFICIA UNIVERSIDAD JAVERIANA - JULIO 2013 % % INTERFAZ CEREBRO COMPUTADOR MEDIANTE LA CLASIFICACIÓN DE SEÑALES ELECTROENCEFALOGRÁFICAS % AUTOR : ANGÉLICA REYES RUEDA % DIRECTOR: ING. CESAR LEONARDO NIÑO BARRERA

93

% RBFSVM.M - RADIAL BASIS FUNCTION KERNEL % % Funcion que realiza el calculo del mejor parametro gamma y cost para un % modelo SVM-RBF. El calculo se realiza mediante validación cruzada con % k-fold iteraciones, se realiza la busqueda mediante la técnica "GRID".

%######################################% %# VALIDACIÓN CRUZADA - ESCALA GRANDE #% %######################################% % PARAMETROS - GRID (RED) paso = 5; C = -20:paso:20; % CALCULO DE PARAMETROS - C Y GAMMA POR MEDIO DE VALIDACIÓN CRUZADA matriz_cv = zeros(numel(C),1); for i=1:numel(C) param = ['-q -v ', num2str(kfolds), ' -c ', num2str(2^C(i))]; modelo = train(trainLabel, sparse(trainData), param); % Entrenamiento matriz_cv(i) = predict(testLabel, sparse(testData), modelo); %Ensayo end

% Mejor Precisión [~,index] = max(matriz_cv);

% Parametros que logran el modelo con mejor precisión best_cv = matriz_cv(index); best_log2C = C(index); best_C = 2^(best_log2C);

clc disp(['Validación Cruzada - Escala Grande:']) disp(['Mejor log2(C): ',num2str(best_log2C), ' disp(['Precisión: ',num2str(best_cv),'%']);

- Parametro-C = ',num2str(best_C), '']);

%######################################% %# VALIDACIÓN CRUZADA - ESCALA MEDIA #% %######################################%

% PARAMETROS - GRID (RED) pasoant = paso; paso = 1; % GENERACIÓN DE MATRIZ GRID (RED) C = best_log2C-pasoant:paso:best_log2C+pasoant; % CALCULO DE PARAMETROS - C Y GAMMA POR MEDIO DE VALIDACIÓN CRUZADA matriz_cv = zeros(numel(C),1); for i=1:numel(C) param = ['-q -v ', num2str(kfolds), ' -c ', num2str(2^C(i))]; modelo = train(trainLabel, sparse(trainData), param); % Entrenamiento matriz_cv(i) = predict(testLabel, sparse(testData), modelo); %Ensayo end % Mejor Precisión [~,index] = max(matriz_cv); % Parametros que logran el modelo con mejor precisión

94

best_cv = matriz_cv(index); best_log2C = C(index); best_C = 2^(best_log2C); clc disp(['Validación Cruzada - Escala Media:']) disp(['Mejor log2(C): ',num2str(best_log2C), ' disp(['Precisión: ',num2str(best_cv),'%']);

- Parametro-C = ',num2str(best_C), '']);

%########################################% %# VALIDACIÓN CRUZADA - ESCALA PEQUEÑA #% %########################################% % PARAMETROS - GRID (RED) pasoant = paso; paso = 1/4; % GENERACIÓN DE MATRIZ GRID (RED) C = best_log2C-pasoant:paso:best_log2C+pasoant; % CALCULO DE PARAMETROS - C Y GAMMA POR MEDIO DE VALIDACIÓN CRUZADA matriz_cv = zeros(numel(C),1); for i=1:numel(C) param = ['-q -v ', num2str(kfolds), ' -c ', num2str(2^C(i))]; modelo = train(trainLabel, sparse(trainData), param); % Entrenamiento matriz_cv(i) = predict(testLabel, sparse(testData), modelo); %Ensayo end % Mejor Precisión [~,index] = max(matriz_cv); % Parametros que logran el modelo con mejor precisión best_cv = matriz_cv(index); best_log2C = C(index); best_C = 2^(best_log2C); clc disp(['Validación Cruzada - Escala Pequeña:']) disp(['Mejor log2(C): ',num2str(best_log2C), ' disp(['Precisión: ',num2str(best_cv),'%']);

- Parametro-C = ',num2str(best_C), '']);

end

N.

RBFSVM.M

function [best_C, best_gamma, best_cv] = RBFSVM(trainLabel, trainData, testData, testLabel, kfolds, ID, Op) % TRABAJO DE GRADO - TG 1174 % PONTIFICIA UNIVERSIDAD JAVERIANA - JULIO 2013 % % INTERFAZ CEREBRO COMPUTADOR MEDIANTE LA CLASIFICACIÓN DE SEÑALES ELECTROENCEFALOGRÁFICAS % AUTOR : ANGÉLICA REYES RUEDA % DIRECTOR: ING. CESAR LEONARDO NIÑO BARRERA % RBFSVM.M - RADIAL BASIS FUNCTION KERNEL % % Funcion que realiza el calculo del mejor parametro gamma y cost para un % modelo SVM-RBF. El calculo se realiza mediante validación cruzada con % k-fold iteraciones, se realiza la busqueda mediante la técnica "GRID". % GUARDAR PLOTS COMO .PNG EN LA CARPETA /PLOTS/

95

ext = '.png'; path = 'PLOTS/';

%######################################% %# VALIDACIÓN CRUZADA - ESCALA GRANDE #% %######################################%

% PARAMETROS - GRID (RED) paso = 10; [C,gamma] = meshgrid(-20:paso:20, -20:paso:20); % CALCULO DE PARAMETROS - C Y GAMMA POR MEDIO DE VALIDACIÓN CRUZADA matriz_cv = zeros(numel(C),1); for i=1:numel(C) param = ['-q -v ', num2str(kfolds), ' -c ', num2str(2^C(i)), ' -g ', num2str(2^gamma(i))]; modelo = svmtrain(trainLabel, trainData, param); % Entrenamiento matriz_cv(i) = predict(testLabel, testData, modelo); %Ensayo end % Mejor Precisión [~,index] = max(matriz_cv); % Parametros que logran el modelo con mejor precisión best_cv = matriz_cv(index); best_log2C = C(index); best_log2gamma = gamma(index); best_C = 2^(best_log2C); best_gamma = 2^(best_log2gamma); % % clc % disp(['Validación Cruzada - Escala Grande:']) % disp(['Mejor log2(C): ',num2str(best_log2C), ' - Parametro-C = ',num2str(best_C), '']); % disp(['Mejor log2(gamma): ',num2str(best_log2gamma),' - Parametro-Gamma = ',num2str(best_gamma), '']); % disp(['Precisión: ',num2str(best_cv),'%']);

% Gráfica - Contorno de Parametros h=figure('visible','off'); contour(C, gamma, reshape(matriz_cv,size(C))), colorbar holdon plot(C(index), gamma(index), 'rx') text(C(index), gamma(index), sprintf('Acc = %.2f %%',matriz_cv(index)), ... 'HorizontalAlign','left', 'VerticalAlign','top') holdoff xlabel('log_2(C)'), ylabel('log_2(\gamma)'), title(['Precisión de Validación Cruzada - Escala Grande',{['ID: ' num2str(ID),' - TestData S', num2str(Op),]}]) file_name=sprintf('ID0%d_S%d_G',ID,Op); %Nombre de Archivo.PNG saveas(h, [path, file_name, ext], 'png'); % Guardar Plot en /Plots/file_name.png

%######################################% %# VALIDACIÓN CRUZADA - ESCALA MEDIA #% %######################################%

% PARAMETROS - GRID (RED) pasoant = paso; paso = 5; % GENERACIÓN DE MATRIZ GRID (RED)

96

[C,gamma] = meshgrid(best_log2C-pasoant:paso:best_log2C+pasoant, best_log2gammapasoant:paso:best_log2gamma+pasoant); % CALCULO DE PARAMETROS - C Y GAMMA POR MEDIO DE VALIDACIÓN CRUZADA matriz_cv = zeros(numel(C),1); for i=1:numel(C) param = ['-q -v ', num2str(kfolds), ' -c ', num2str(2^C(i)), ' -g ', num2str(2^gamma(i))]; modelo = svmtrain(trainLabel, trainData, param); % Entrenamiento matriz_cv(i) = predict(testLabel, testData, modelo); %Ensayo end % Mejor Precisión [~,index] = max(matriz_cv); % Parametros que logran el modelo con mejor precisión best_cv = matriz_cv(index); best_log2C = C(index); best_log2gamma = gamma(index); best_C = 2^(best_log2C); best_gamma = 2^(best_log2gamma); % clc % disp(['Validación Cruzada - Escala Media:']) % disp(['Mejor log2(C): ',num2str(best_log2C), ' - Parametro-C = ',num2str(best_C), '']); % disp(['Mejor log2(gamma): ',num2str(best_log2gamma),' - Parametro-Gamma = ',num2str(best_gamma), '']); % disp(['Precisión: ',num2str(best_cv),'%']); % % Gráfica - Contorno de Parametros h=figure('visible','off'); contour(C, gamma, reshape(matriz_cv,size(C))), colorbar holdon plot(C(index), gamma(index), 'rx') text(C(index), gamma(index), sprintf('Acc = %.2f %%',matriz_cv(index)), ... 'HorizontalAlign','left', 'VerticalAlign','top') holdoff xlabel('log_2(C)'), ylabel('log_2(\gamma)'), title(['Precisión de Validación Cruzada - Escala Mediana',{['ID: ' num2str(ID),' - TestData S', num2str(Op),]}]) file_name=sprintf('ID0%d_S%d_M',ID,Op); %Nombre de Archivo.PNG saveas(h, [path, file_name, ext], 'png'); % Guardar Plot en /Plots/file_name.png

%########################################% %# VALIDACIÓN CRUZADA - ESCALA PEQUEÑA #% %########################################% % PARAMETROS - GRID (RED) pasoant = paso; paso = 1; % GENERACIÓN DE MATRIZ GRID (RED) [C,gamma] = meshgrid(best_log2C-pasoant:paso:best_log2C+pasoant, best_log2gammapasoant:paso:best_log2gamma+pasoant); % CALCULO DE PARAMETROS - C Y GAMMA POR MEDIO DE VALIDACIÓN CRUZADA matriz_cv = zeros(numel(C),1); for i=1:numel(C) param = ['-q -v ', num2str(kfolds), ' -c ', num2str(2^C(i)), ' -g ', num2str(2^gamma(i))]; modelo = svmtrain(trainLabel, trainData, param); % Entrenamiento matriz_cv(i) = predict(testLabel, testData, modelo); %Ensayo end

97

% Mejor Precisión [~,index] = max(matriz_cv); % Parametros que logran el modelo con mejor precisión best_cv = matriz_cv(index); best_log2C = C(index); best_log2gamma = gamma(index); best_C = 2^(best_log2C); best_gamma = 2^(best_log2gamma); clc disp(['Validación Cruzada - Escala Pequeña:']) disp(['Mejor log2(C): ',num2str(best_log2C), ' - Parametro-C = ',num2str(best_C), '']); disp(['Mejor log2(gamma): ',num2str(best_log2gamma),' - Parametro-Gamma = ',num2str(best_gamma), '']); disp(['Precisión: ',num2str(best_cv),'%']); % Gráfica - Contorno de Parametros h=figure('visible','off'); contour(C, gamma, reshape(matriz_cv,size(C))), colorbar holdon plot(C(index), gamma(index), 'rx') text(C(index), gamma(index), sprintf('Acc = %.2f %%',matriz_cv(index)), ... 'HorizontalAlign','left', 'VerticalAlign','top') holdoff xlabel('log_2(C)'), ylabel('log_2(\gamma)'), title(['Precisión de Validación Cruzada - Escala Pequeña',{['ID: ' num2str(ID),' - TestData S', num2str(Op),]}]) file_name=sprintf('ID0%d_S%d_P',ID,Op); %Nombre de Archivo.PNG saveas(h, [path, file_name, ext], 'png'); % Guardar Plot en /Plots/file_name.png end

98