Axel Hutt

Keywords: Brain-computer interfaces, EEG, motor imagery, robotics control, multilabel approach, Common Spatial Pattern, Riemannian geometry i . . . . . . . . . . . . . . . . . . . . . .

Brain-Computer Interfaces (BCIs) replace the natural nervous system outputs by artificial ones that do not require the use of peripheral nerves, allowing people with severe motor impairments to interact, only by using their brain activity, with different types of applications, such as spellers, neuroprostheses, wheelchairs, or among others robotics devices. A very popular technique to record signals for BCI implementation purposes consists of electroencephalography (EEG), since in contrast with other alternatives, it is noninvasive and inexpensive. In addition, due to the potentiality of Motor Imagery (MI, i.e., brain oscillations that are generated when subjects imagine themselves performing a movement without actually accomplishing it) to generate suitable patterns for scheming self-paced paradigms, such combination has become a common solution for BCI neuroprostheses design. However, even though important progress has been made in the last years, full 3D control is an unaccomplished objective.

In order to explore new solutions for overcoming the existing limitations, we present a multiclass approach that considers the detection of combined motor imageries, (i.e., two or more body parts used at the same time). The proposed paradigm includes the use of the left hand, right hand, and both feet together, from which eight commands are provided to direct a robotic arm comprising fourteen different movements that afford a full 3D control. To this end, an innovative switching-mode scheme that allows managing different actions by using the same command was designed and implemented on the OpenViBE platform. Furthermore, for feature extraction a novel signal processing scheme has been developed based on the specific location of the activity sources that are related to the considered body parts. This insight allows grouping together within a single class those conditions for which the same limb is engaged, in a manner that the original multiclass task is transformed into an equivalent problem involving a series of binary classification models. Such approach allows using the Common Spatial Pattern (CSP) algorithm; which has been shown to be powerful at discriminating sensorimotor rhythms, but has the drawback of being suitable only to differentiate between two classes. Based on this perspective we also have contributed with a new strategy that combines together the CSP algorithm and Riemannian geometry. In which the CSP projected trials are mapped into the Riemannian manifold, from where more discriminative features can be obtained as the distances separating the input data from the considered class means. These strategies were applied on three new classification approaches that have been compared to classical multiclass methods by using the EEG signals from a group of naive healthy subjects, showing that the proposed alternatives not only outperform the existing schema, but also reduce the complexity of the classification task.

Résumé

Les interfaces cerveau-ordinateur (ou BCI en anglais pour Brain-Computer Interfaces) mettent en place depuis le système nerveux central un circuit artificiel secondaire qui remplace l'utilisation des nerfs périphériques, permettant entre autres à des personnes ayant une déficience motrice grave d'interagir, uniquement à l'aide de leur activité cérébrale, avec différents types d'applications, tels qu'un système d'écriture, une neuro-prothèse, un fauteuil roulant motorisé ou un bras robotique. Une technique répandue au sein des BCI pour enregistrer l'activité cérébrale est l'électroencéphalographie (EEG), étant donné que contrairement à d'autres techniques d'imagerie, elle est non invasif et peu coûteuse. En outre, l'imagination motrice (MI), c'est-à-dire les oscillations des neurones du cortex moteur générées lorsque les sujets imaginent effectuer un mouvement sans réellement l'accomplir, est appropriée car détectable dans l'EEG et liée à l'activité motrice pour concevoir des interfaces comme des neuro-prothèses non assujetties à des stimuli. Cependant, même si des progrès importants ont été réalisés au cours des dernières années, un contrôle 3D complet reste un objectif à atteindre. Afin d'explorer de nouvelles solutions pour surmonter les limitations existantes, nous présentons une approche multiclasses qui considère la détection des imaginations motrices combinées. Le paradigme proposé comprend l'utilisation de la main gauche, de la main droite, et des deux pieds ensemble. Ainsi, par combinaison, huit commandes peuvent être fournies pour diriger un bras robotisé comprenant quatorze mouvements différents qui offrent un contrôle 3D complet. À cette fin, un système de commutation entre trois modes (déplacement du bras, du poignet ou des doigts) a été conçu et permet de gérer les différentes actions en utilisant une même commande. Ce système a été mis en oeuvre sur la plate-forme OpenViBE. En outre, pour l'extraction de caractéristiques une nouvelle approche de traitement d'information fournie par les capteurs a été développée sur la base de l'emplacement spécifique des sources d'activité liées aux parties du corps considérées. Cette approche permet de regrouper au sein d'une seule classe les différentes actions pour lesquelles le même membre est engagé, d'une manière que la tâche multiclasses originale se transforme en un problème équivalent impliquant une série de modèles de classification binaires. Cette approche permet d'utiliser l'algorithme de Common Spatial pattern (CSP) dont la capacité à discriminer des rythmes sensorimoteurs a été largement montrée mais qui présente l'inconvénient d'être applicable uniquement pour différencier deux classes. Nous avons donc également contribué à une nouvelle stratégie qui combine un ensemble de CSP et la géométrie riemannienne. Ainsi des caractéristiques plus discriminantes peuvent être obtenues comme les distances séparant les données des centres des classes considérées. Ces stratégies ont été appliquées sur trois nouvelles approches de classification qui ont été comparées à des méthodes de discrimination multiclasses classiques en utilisant les signaux EEG d'un groupe de sujets sains naïfs, montrant ainsi que les alternatives proposées permettent non seulement d'améliorer l'existant, mais aussi de réduire la complexité de la classification.

Mots-clés: Interfaces cerveau-ordinateur, EEG, imagination motrice, contrôle robotique, classification multilabels, Common spatial pattern, géométrie Riemanienne

General Introduction

Brain-Computer Interfaces (BCI) systems replace the natural nervous system outputs by artificial ones that do not demand the use of peripheral nerves, allowing people with severe motor impairments to interact, only by using their brain activity, with different types of applications, such as spellers, neuroprostheses, wheelchairs, or among others robotics devices. [START_REF] Wolpaw | BCI: Principles and Practice[END_REF], [START_REF] Gomez-Gil | Brain computer interfaces -a review[END_REF]. The use of BCIs is not restricted and they have also raised the interest for developing rehabilitation schemes aimed at stroke patients [START_REF] Pfurtscheller | Event-related synchronization (ERS): an electrophysiological correlate of cortical areas at rest[END_REF], [START_REF] Fok | An EEG-based brain computer interface for rehabilitation and restoration of hand control following stroke using ipsilateral cortical physiology[END_REF], as well as interfaces for healthy people (e.g., in gaming applications).

There are several brain signals that can be used for BCI control. One of particular interest, since it has the benefit of favoring self-paced paradigms (i.e., an asynchronous operation mode), is known as Motor Imagery (MI); which can provide an intuitive and natural control to enable people with restricted mobility an ample interaction with their environment. Motor imagery involves the mental practice of a motor task without exhibiting any apparent movement [START_REF] Pfurtscheller | Motor imagery and direct brain-computer communication[END_REF]. The imagination of a movement includes similar brain regions and functions as those engaged in programming and executing a real one. The main difference is that in the former case, movements are blocked at a cortico-spinal level [START_REF] Neuper | Evidence for distinct beta resonance frequencies in human EEG related to specific sensorimotor cortical areas[END_REF].

The EEG activity elicited during motor imageries consists of changes of sensorimotor rhythms that exhibit an oscillatory behavior that is of particular interest within two specific frequency ranges: the alpha (8-13 Hz) and beta [START_REF] Barachant | Multi-class brain computer interface classification by riemannian geometry[END_REF][START_REF] Dornhege | Increase information transfer rates in BCI by CSP extension to multi-class[END_REF][START_REF] Dornhege | Boosting bit rates in noninvasive eeg single-trial classifications by feature combination and multiclass paradigms[END_REF][START_REF] Weibo | Eeg feature comparison and classification of simple and compound limb motor imagery[END_REF][START_REF] Royer | EEG) control of a virtual helicopter in 3-dimensional space using intelligent control strategies[END_REF][START_REF] Zhou | A novel technique for phase synchrony measurement from the complex motor imaginary potential of combined body and limb action[END_REF][START_REF] Scherer | Can synergistic finger movements be represented as the superposition of individual finger movements?[END_REF][START_REF] Blankertz | Optimizing spatial filters for robust EEG single-trial analysis [revealing tricks of the trade[END_REF][START_REF] Lindig-León | Comparison of sensorimotor rhythms in EEG signals during simple and combined motor imageries over the contra and ipsilateral hemispheres[END_REF][START_REF] Lindig-León | A multilabel classification method for detection of combined motor imageries[END_REF][START_REF] Mason | Temporal control paradigms for direct brain interfaces -rethinking the definition of asynchronous and synchronous[END_REF][START_REF] Mason | A comprehensive survey of brain interface technology designs[END_REF][START_REF] Mak | Clinical applications of brain-computer interfaces: current state and future prospects[END_REF]) bands [START_REF] Pfurtscheller | Beta rebound after different types of motor imagery in man[END_REF]- [START_REF] Babiloni | Human cortical electroencephalography (EEG) rhythms during the observation of simple aimless movements: a high-resolution EEG study[END_REF], since they have been shown to be modulated during and following the preparation and execution of imagined movements [START_REF] Avanzini | The dynamics of sensorimotor cortical oscillations during the observation of hand movements: an EEG study[END_REF]. This modulation takes two different forms: an oscillatory power reduction during the motor task known as Event-Related Desynchronization (ERD); and an oscillatory power increase known as Event-Related Synchronization (ERS) that appears just after the end of the motor task over a few hundred milliseconds [START_REF] Pfurtscheller | Event-related EEG/MEG synchronization and desynchronization: basic principles[END_REF]. During the preparation of a motor task certain cortical and subcortical regions are activated; contrarily, in the course of motor imageries most of the activity is found within the primary motor cortex over the corresponding contralateral hemisphere [START_REF] Pfurtscheller | On the existence of different types of central beta rhythms below 30 hz[END_REF].

The topographic representation and band power change during motor imageries are known to be located in specific regions on the primary motor cortex (broadly, left hand activities are observed over electrode C4, right hand activities on C3, feet activities on C z , etc.) [START_REF] Pfurtscheller | Event-related EEG/MEG synchronization and desynchronization: basic principles[END_REF], [START_REF] Barachant | Multi-class brain computer interface classification by riemannian geometry[END_REF].

The number of commands afforded by motor imagery-based BCIs depends on the amount of brain states that the system is able to detect which, in turn, is limited to the number of body parts that users are able of producing suitable patterns with. Moreover, the difficulty of such a task increases as the total of available commands grows, which is worthwhile only if a xxi GENERAL INTRODUCTION suitable amount of effectively separable human brain states can be assigned [START_REF] Dornhege | Increase information transfer rates in BCI by CSP extension to multi-class[END_REF]. There are two approaches for solving multiclass problems: On the one hand some techniques directly forecast the predicted label by using a straightaway method, such as distance-based algorithms; and on the other hand some alternatives apply multiclass extensions to decompose the problem into an ensemble of binary tasks that can be accomplished in different ways, for instance, by using the classical one-versus-one or one-versus-all approaches [START_REF] Dornhege | Boosting bit rates in noninvasive eeg single-trial classifications by feature combination and multiclass paradigms[END_REF].

An alternative multiclass scheme that allows affording users with a plenty repertory of different commands is known as multilabel approach, which addresses the detection of combined motor imageries (i.e., two or more body parts used at the same time) [START_REF] Weibo | Eeg feature comparison and classification of simple and compound limb motor imagery[END_REF], [START_REF] Royer | EEG) control of a virtual helicopter in 3-dimensional space using intelligent control strategies[END_REF]. In comparison with the standard scheme, this approach has the advantage of considerably increasing the number of different brain states while using the same number of body parts (in order of 2 P compared to P , where P is the number of body parts).

In contrast to the standard multiclass paradigms, upon which most of the research has been devoted, less work has been reported about brain oscillatory patterns induced by combined motor imageries. The empirical mode decomposition technique in combination with phase synchronization analysis, has been proposed and applied to evaluate the human brain synchrony of the supplementary motor and the primary motor areas during complex motor imagination of combined body and limb action, from where the classification performance suggested the feasibility for BCI applications in motor nerve reconstruction of lower limbs [START_REF] Zhou | A novel technique for phase synchrony measurement from the complex motor imaginary potential of combined body and limb action[END_REF]. In addition, aiming at overcoming the limited amount of output commands resulting from the scarce number of classes, and for the purpose of achieving continuous 3D control, both hands motor imagery was adopted to complement the set of instructions in a simple limb motor imagery based-BCI [START_REF] Royer | EEG) control of a virtual helicopter in 3-dimensional space using intelligent control strategies[END_REF]. More recently, the characteristics of EEG patterns during simple and combined motor imageries were investigated to confirm the separability of seven different classes for motor imagery based-BCI systems implementation [START_REF] Weibo | Eeg feature comparison and classification of simple and compound limb motor imagery[END_REF]. Another study reported that the movement of several fingers at the same time is not the simple superposition of individual finger movements at the ECoG scale, but rather the activation of additional neural networks. Which implies that for BCI implementation purposes, the generation of control signals encoding complex movements would likely use signal components beyond individual members superposition [START_REF] Scherer | Can synergistic finger movements be represented as the superposition of individual finger movements?[END_REF]. However, to our best knowledge the latter conclusion remains to be confirmed at the EEG scale; which has become the main interest of the present work.

The importance of this question lies on the particular interest of optimizing the Common Spatial Pattern (CSP) algorithm for its application on multilabel classification tasks. This method is a very convenient and widely used algorithm to extract motor imagery features from EEG signals. It generates a series of spatial filters that are applied to decompose multi-dimensional data into a set of uncorrelated components [START_REF] Blankertz | Optimizing spatial filters for robust EEG single-trial analysis [revealing tricks of the trade[END_REF], which aim at extracting elements that simultaneously maximize the variance of one class, while minimizing the variance of the other one. Because of the manner in which the CSP algorithm achieves such a contrast, this approach has the drawback of being suitable only to discriminate between two classes and, in order to solve problems involving xxii many conditions, multiclass extensions are required [START_REF] Gomez-Gil | Brain computer interfaces -a review[END_REF], [START_REF] Dornhege | Boosting bit rates in noninvasive eeg single-trial classifications by feature combination and multiclass paradigms[END_REF]. Classical extensions such as the one-versus-one and one-versus-all approaches present several disadvantages. Perhaps the most critical one is the extensive amount of models that is required as the number of classes increases, which negatively impacts on the computation time and on the reliability of results. In this sense, by confirming that the motor imagery of a combination of two or more body parts involved at the same time corresponds, in fact, to the simple superposition of the individual members at the EEG scale, it is possible to design more convenient approaches for adapting not only the CSP method but any other, since, under this scheme, it becomes possible to use an ensemble of binary classifiers to determine at each activity source whether or not its associated body part is engaged in the task, and thus to infer the final decision from the combination of all subproblems [START_REF] Lindig-León | Comparison of sensorimotor rhythms in EEG signals during simple and combined motor imageries over the contra and ipsilateral hemispheres[END_REF], [START_REF] Lindig-León | A multilabel classification method for detection of combined motor imageries[END_REF].

In the present study we provide some evidence suggesting that this hypothesis is plausible, upon which we have developed three new feature extraction methods for discriminating combined motor imageries, namely, One-step Multilabel (OsM) approach, Hierarchical Multilabel (HM) approach, and One-step Hierarchical Multilabel (OsHM) approach, which to our best knowledge have not been considered before. The motivation was to provide a full 3-D control for a commercial robotic arm that affords users with fourteen predefined motions. The proposed paradigm involves the use of three body parts, namely, right hand, left hand and feet; which in combination with the rest condition, constitutes an 8-class problem. For the implementation, a switching-mode control that allows managing different actions by using the same command was proposed. Six healthy subjects participated in one recording session comprising three runs of ten trials per class, which totals 30 trials per class. Results were slightly superior as those obtained by the classical one-versus-one and one-versus-all approaches, which is of particular interest if we consider that, in order to solve the aforementioned 8-class problem both methods demand, respectively, 28 and 8 feature extraction models; which is by far more than the 3, 7 and 6 models that are required by the OsM, HM and OsHM approaches.

The remainder of this document is organized as follows: The first chapter is dedicated to provide a general basis of the principles and practices of Brain-Computer Interface (BCI) systems. These aspects will be further deepened along the second chapter in the context of motor imagery.

Along these sections we will present the 4-class database, which comprises a series of EEG recordings from five healthy subjects who participated in a motor imagery experiment involving the use of the right hand, left hand and both hands together; which, in addition with the rest condition, leads to a paradigm of four different classes. The 4-class database has been created to carry out a preliminary study to assess the similarities of the activity that is found over the electrodes associated with the use of the left and right hands during simple and combined motor imageries. In this way, by proposing a small number of classes, it is possible to register a considerable amount of trials without causing users fatigue. Upon these data we will present a series of statistical results based on the analysis of the power spectrum observed at each electrode during the four different motor tasks, which supports the hypothesis that EEG data recorded during a combined motor imagery can be characterized as the superposition of the activity that xxiii GENERAL INTRODUCTION is independently generated by each one of the sources associated to the body parts engaged in the task. Subsequently, along the third chapter, we will describe the characteristics of the robotic arm that has motivated the practical solutions developed throughout this project. In the same line, we will describe the operation scheme that has been designed; which consists of a smart implementation based on a switching-mode control, that allows managing different actions by using the same command. As mentioned, this device affords users with fourteen predefined motions, which demands far more commands than those obtained with the 4-class database.

With this in mind, we have generated the 8-class database, which comprises EEG signals from 6 healthy subjects inducing the eight mental states that it is possible to generate by the motor tasks involving the use of the right hand, left hand, both feet together, and rest. Similarly as with the 4-class database, we will present a qualitative analysis of the ERD patterns that are found in the EEG signals during the eight different tasks, together with a statistical comparison of the power spectrum at each electrode among each pair of conditions. These results create the basis upon which we will justify the grouping of classes that has been suggested for the feature extraction methods contributed by the present work. In the fourth chapter the concept of feature extraction is introduced and elaborated in the context of EEG signals during motor imagery. Throughout these pages, the CSP algorithm will be elucidated together with two of its variants, namely, Analytic Common Spatial patterns (ACSP) and Common Spatial Pattern by Joint Approximate Diagonalization (CSP by JAD). Additionally, we will present a strategy called Filter Bank Common Spatial Pattern (FBCSP) that optimizes the algorithm performance by decomposing the EEG signals into a series of frequency bands, which favors the selection of a discriminative subject-specific frequency range. Subsequently, we will present the notion of classification, for which two commonly used algorithms will be described, namely, Linear Discriminant Analysis (LDA) and Support Vector Machines (SVM). Moreover, in order to explore different solutions, the present work considers as well the use of Riemannian geometry; with which classification can be directly achieved within the manifold of symmetric and positive definite matrices that is given by the covariance matrices of the EEG signals. Afterwards, in the last section of this chapter, we will introduce a new classification method for motor imagery-based BCI systems that makes use of the combination of the CSP algorithm and Riemannian geometry, which represents one of the contributions of the present work. In the fifth chapter the concept of multiclass approach will be introduced by describing three classical methods, namely, one-versus-one, one-versus-all, and hierarchical algorithms. At this point we will describe the three new multilabel methods that we have implemented, namely, One-step Multilabel (OsM) approach, Hierarchical Multilabel (HM) approach, and One-step Hierarchical Multilabel (OsHM) approach. The basis of each one of these new methods lies in the results that were presented along chapters 2 and 3, upon which we have proposed a convenient manner to carry out data partition, so that the entire problem is reduced into a series of binary classification tasks. Finally, along the sixth chapter we will present the results, in terms of accuracy, that were obtained by each one of the described methods over both databases. The corresponding discussion is presented in chapter Brain-computer interface (BCI) systems consist of a combination of modules that allow recording the user brain activity in order to analyze it and find specific patterns indicating that the person is executing a particular mental task, which in turn is associated to a specific command for controlling the function of certain devices, such as computers, robotics, wheelchairs, prostheses, orthoses, etc. In this way, BCIs replace the natural nervous system pathways by artificial ones, which allow people with a severe motor impairment to interact with their environment by using only their brain activity. Most of BCI applications aim at improving the quality of life of their users, who have lost certain function(s) due to a degenerative disease or an accident, in which case the system will attempt to restore the missing ability. There are other functions that BCIs undertake according to the user needs and system scopes, which will be presented in the following chapter.

The initial part of this document is dedicated to introduce the notion of BCI systems. To this end, their general architecture will be presented along the first section, where the function of each of their components will be described. Additionally, the concepts underlying the BCI terminology, as well as the different types of brain activity that are suitable for their design, will be detailed.

Finally, we will present the basic applications that BCIs are aimed for based on the nature of the brain activity that is used to design different paradigms. In this sense, we will emphasize on Motor Imagery (MI) since, as it will be further exposed, its nature facilitates robotics control.

BCI definition

Brain-Computer Interfaces consist of a direct communication pathway between the brain and an external device. BCIs are often directed at researching, mapping, assisting, augmenting, or repairing human cognitive or sensory-motor functions. A recent definition involving all usages of BCI systems is stated as follows (see Figure 1.1):

A BCI is a system that measures central nervous system activity and converts it into an artificial output that replaces, restores, enhances, supplements, or improves natural outputs and thereby changes the ongoing interactions between the central nervous system and its external or internal environment [START_REF] Wolpaw | BCI: Principles and Practice[END_REF].

BCI architecture

BCI systems can be presented in different variations, however, from a general point of view there are three main components that every BCI includes in its architecture, namely, a signal acquisition module, a signal processing stage, and the corresponding output, which can presented in many different forms according to the application of the system [START_REF] Wolpaw | BCI: Principles and Practice[END_REF], [START_REF] Mason | Temporal control paradigms for direct brain interfaces -rethinking the definition of asynchronous and synchronous[END_REF]- [START_REF] Mak | Clinical applications of brain-computer interfaces: current state and future prospects[END_REF]. In addition, BCIs might include a feedback module to provide users with some specific information, such as brain signals and/or classifiers outputs. Figure 1.1 shows the general framework of a BCI where it is possible to observe each one of the aforementioned modules.

Signal acquisition

As it will further seen along this chapter, brain signals can be recorded at different scales depending on whether the system is invasive (i.e., direct contact with the cortex) or not (superficial contact onto the skin). In both cases electrodes are located over the regions of interest and the acquired signals are amplified and digitized to be used by the BCI system. In this stage some preprocessing steps can be considered, such as the use of a notch filter to suppress the line noise, and/or the use of band-pass filters to extract the frequencies of interest.

Signal processing

Signal processing is the analysis of the brain signals that is carried out in order to extract and interpret useful information related to the user intention. Basically, there are two main aspects of signal processing in the context of BCI, namely, feature extraction and feature translation.

The first concept consists in determining the characteristics that confer distinctiveness to the brain activity that is elicited during the corresponding task, so that it can be recognized from other states. On the other hand feature translation consists in the classification of the extracted features in conjunction with the assignment of the corresponding action.

Feature extraction

The goal of feature extraction is to select relevant information from the recorded signals in order to fit a classification algorithm (classifier) that is able to categorize such elements according to their position within the classification space, so that the system can proceed to translate the corresponding command.

In most of the cases this process reduces the dimensionality of data by excluding redundant information, and it enhances the attributes of the signal by exploiting the aspects that are related to the physiological phenomena. The derived values are projected into the feature space, where elements from different conditions are expected to be located as separated as possible, so that distinct groups can be easily identified.

Feature classification

The aim of classification is to assign the feature vectors to one of K discrete classes, where K corresponds with the number of different brain states that the system is able to detect and to translate into commands for controlling the interface.

The classification algorithm can be implemented for distinguishing between two classes (binary problems), or several classes (multiclass problems). In this sense, given that the difficulty of the classification task increases with the number of classes, there is a trade-off between providing a natural and intuitive control by affording users with many commands, or ensuring an effective performance with a limited action.

System output

As mentioned, the number of considered mental states (classes) are associated to control actions (commands). These actions can be used to control an unlimited number of different types of applications, the most common ones are based on the use of devices such as computers for communication or control, electric wheelchairs, and robotic limbs to be used in rehabilitation or for accomplishing motion actions. Some applications include movement restoration, in which BCI systems bypass a spinal cord injury and activate muscles with electrical stimulation.

Feedback

In many cases, BCI systems can provide useful information that may be used to improve their control. For instance, o detect and inform users about execution errors. This feature is not always included in BCI applications, even though it can be specially useful during on-line sessions, since errors may be critical and it becomes important to inform users about their progress. Recently, BCI researchers attempt to incorporate feedback as a tool for training and developing users skills [START_REF] Lotte | A:8 teaching the bci skill: feedback and human training approaches[END_REF].

BCI taxonomy

Aside from the general architecture shared by BCI systems, there are other aspects from which specific features can be distinguished. Such as the nature and manner in which the targeted brain activity is induced, the operating protocols and the type of decoding which, among other characteristics, provide a considerable number of different alternatives for BCI applications that can be pursued according to the users' needs. In the following sections we provide a description of some important features for BCI design.

Dependent and independent BCIs

The terms dependent BCI and independent BCI were recently introduced [START_REF] Wolpaw | Memory in neuroscience: rhetoric versus reality[END_REF], according to the basic definition, both use brain signals, but the differ in their dependence on natural outputs induced by the central nervous system. In this sense, a dependent BCI uses brain signals that depend on muscle activity. For instance, the P300 speller demands user to direct their gaze towards the screen position where the symbol they want to select is located. Even though dependent BCIs do not provide the brain with a new output that is independent of natural outputs, they may still be useful [START_REF] Wolpaw | BCI: Principles and Practice[END_REF].

On the contrary, an independent BCI does not depend on natural outputs at all, since muscle activity is not fundamental for inducing the brain signals that are used by the BCI. For instance, in motor imagery-based BCIs, users can modulate their brain signals by imagining a motor task in order to control the system output. For people with severe motor impairments independent BCIs are likely to be more useful. However, it is important to acknowledge that nowadays most of the BCIs are neither completely dependent nor completely independent. Just as it is necessary to control gaze direction to use the P300 speller, most of motor imagery-based BCIs require the user having sufficient visual function (and thus gaze control) to watch the results of the BCIs output commands (e.g., cursor movements).

Exogenous and endogenous BCIs

Brain activity for BCI control can be elicited spontaneously as the result of an internal process, or can be triggered by an external stimulus. The former case receives the name of endogenous brain activity, which can be induced by some mental tasks such as the imagination of body movements.

The second case is known as exogenous activity and can be evoked by either visual, auditory or somatosensory stimuli intended to induce patterns that are easily detectable.

Endogenous BCIs are likely to provide more freedom. For example, an endogenous BCI can allow users to move a cursor towards any direction in a 2-dimensional (or even 3-dimensional) space, whereas an exogenous BCI can allow subjects to make a choice between different symbols located on the same space. However, an important advantage of exogenous BCIs is that they do not require users to be trained, and thus they can be effectively used since the very first time, contrarily to endogenous BCIs, which sometimes require long periods of training, without implying that they would work accurately afterwards.

Passive and active BCIs

Another recent term, passive BCI, is applied to BCI applications that use brain signals associated with aspects of the user's current state, such as level of attention [START_REF] Zander | Towards passive brain-computer interfaces: applying braincomputer interface technology to human-machine systems in general[END_REF]. For instance, a BCI might detect EEG features preceding lapses in attention and produce an output (e.g., a sound) that alerts the user and restores attention. In this sense, the term passive is intended to distinguish between this kind of BCI applications from those that provide communication and control (i.e., active BCIs). However,passive and active are subjective terms that lack clear neuroscientific definitions.

Moreover, continued use of a passive BCI might induce adaptations of the central nervous system that improve its performance, so that the term passive becomes no longer applicable. Thus, it seems preferable to categorize BCI applications simply as shown in Figure 1.1, in which case passive BCIs will generally fit into the enhance or supplement category.

Hybrid BCIs

The concept of hybrid BCI can be applied in two different ways [START_REF] Graiman | Brain-computer interfaces[END_REF]. It may point to a BCI that uses two different kinds of of brain signals (e.g., visual evoked potentials and sensorimotor rhythms) to generates its outputs or, alternatively, it may describe a system that combines a BCI output with a natural muscle-based output. In the latter usage, the BCI output supplements a natural output (e.g., as illustrated in Figure 1.1).

BCI operating protocols

As with any other task, a BCI user may by sometimes actively using the BCI (e.g., using it to communicate with someone), or may be not using it at other times (e.g., is doing something else, sleeping, daydreaming, etc.) [START_REF] Mason | Temporal control paradigms for direct brain interfaces -rethinking the definition of asynchronous and synchronous[END_REF], [START_REF] Mason | A comprehensive survey of brain interface technology designs[END_REF]:

-The state of active BCI use is referred to as intentional control.

-The state of other activities is referred to as no intentional control, or more simply, nocontrol.

Ideally, BCI use should be as natural and convenient as a person's muscle-based actions, so that the BCI is able to operate in a self-paced mode, in which the system is readily available and users can achieve intentional control when desired. To afford such user-paced operation, a BCI must be able to reliably distinguish between the user's no-control state and the user's intentional control state. This is a difficult requirement because a wide range of brain activity signals might correspond to no-control, and the system must be able to recognize all of them as indicating such state. If the BCI is not able of accurately distinguishing between no-control and intentional control, brain activity during no control will sometimes be translated into unintended actions. For many BCI applications, even a small number of such false-positive errors could make BCI use very frustrating and possibly impractical.

Currently, BCI operating protocols address the need of avoiding false-positive errors in one of two ways: Trough synchronous protocols, which limit BCI actions to times when the user can safely be assumed to be in the intentional-control state; and, alternatively, through self-paced (or asynchronous) protocols, which are able to distinguish satisfactorily between the no-control state and the intentional-control state.

Synchronous protocols

As mentioned, synchronous protocols limit the times during which the BCI will actually translate the user's brain activity into an action [START_REF] Mason | Temporal control paradigms for direct brain interfaces -rethinking the definition of asynchronous and synchronous[END_REF], [START_REF] Müller-Putz | Brain-computer interfaces for control of neuroprostheses: from synchronous to asynchronous mode of operation[END_REF]. A common approach for providing synchronous protocols consists in giving the user a specific sensory stimulus and to analyze the brain activity that is induced during a fixed period immediately afterward. In this way, the system converts brain signals into outputs only during that poststimulus period, so that it limits and specifies the times when the BCI will convert brain activity into an action. Alternatively, the protocol might provide cues to inform users of an impending period of expected intentional control.

Self-paced protocols

A self-paced protocol (also referred to as an asynchronous protocol) is able to distinguish between no-control and intentional-control states. This mode of operation allows users to control the BCI whenever they desire. Self-paced protocols can afford more natural and dynamic interactions between the user and the BCI, but they are not yet suitable for general out-of-the-lab use.

Self-paced BCIs can rely on motor imagery as well. For instance, if the system detects the EEG activity associated with one of several specific mental tasks (such as imagining a left-hand movement), then that activity is translated into a particular command.

Continuous and discrete decoding

In general, BCI systems have been used to provide one of two types of decoding: continuous or discrete. Continuous decoding is best described as real-time control of an ongoing movement, such as of a cursor or a robotic limb. Contrarily, discrete decoding consists of classifying specific events, as when a letter or icon is selected from among a set of possible choices.

Brain Signals for BCI

Most of the activity for BCIs implementation purposes consists of brain's electrical signals, which can be recorded at different scale levels, from the centimeter scale of electroencephalography (EEG) through the millimeter scale of electrocorticography (ECoG) to the tens-of-microns scale of neuronal action potentials (spikes). Which recording resolution is more advantageous for which purposes remains as a research topic, since even though the activity recorded from basic units suggests that spikes could provide several degrees of freedom and might thus be the best signals for BCIs use, the evidence that adaptive methods can elicit multiple brain states also from EEG signals, points out that the difference between BCI performance provided by these different scales may not be as nearly important as the difference in their topographical resolutions.

Recently, nonelectrical metabolic signals, such as those generated by functional near-infrared spectroscopy (fNIRS) and functional magnetic resonance imaging (fMRI), have drawn the attention of BCIs developers due to their potentiality of exposing suitable brain activity patterns. On the one hand, fNIRS has the advantages of being noninvasive and inexpensive, whereas on the other hand, despite being much more expensive, fMRI is not invasive either and offers a much higher resolution.

Besides of the suitability of brain signals for BCI purposes, a very important aspect consists of the access to the corresponding sources, which may demand an invasive procedure that constrains the feasibility of certain techniques. In the following sections, some examples of invasive and non-invasive recording methods are provided together with their advantages and disadvantages.

Invasive recording techniques

Invasive recordings are often used in humans only for evaluation purposes in patients suffering from epilepsy, Parkinson's disease or others neurological conditions. In these cases electrodes are implanted temporarily (5-14 days [START_REF] Lal | Methods towards invasive human brain computer interfaces[END_REF]) to evaluate the possibility of undergoing surgery for the extreme cases where patients are resistant to pharmacological treatments [START_REF] Nair | Chronic subdural electrodes in the management of epilepsy[END_REF]]. The main goal of this technique is to identify the brain tissue where seizures are originated. To this end, electrodes are directly implanted on or in the cortex through a surgical procedure, which provides a high temporal and spatial resolution. In general, BCI research by using invasive recordings methods are based on animals, since humans are mostly available only as patients undergoing a surgical procedure who accept to collaborate in the research topic and, in many cases, the implanted electrodes are not located over regions that are of interest for BCI purposes.

Invasive techniques can be subdivided into two major kinds of recordings regarding the type of electrodes that are implanted, namely, electrocorticogram and intracortical recordings.

Electrocorticography

Electrocorticography (ECoG) recordings, also called intracranial EEG (iEEG), use an array of electrodes (on a subdural grid) that is placed in direct contact with the surface of the cortex

BRAIN SIGNALS FOR BCI

without penetrating the brain. ECoG electrodes are usually arranged in a matrix whose size depends on the region that needs to be covered. This technique offers a high spatial resolution and a high signal-to-noise ratio due to the proximity with which neuronal activity is registered.

However, it also presents important inconveniences, such as those related to the required surgery and the following recovery, which also difficult the possibility of a subsequent rearrangement of the electrodes position.

Intracortical Recordings

Microelectrodes can be placed close (i.e., extracellular recordings) or within the neurons (i.e., intracellular recordings). Usually BCI studies based on both types of techniques are restricted to animals. Signals recorded at this scale are action potentials, also referred to as local field potentials in the case of intracellular recordings, or to as spikes for extracellular recordings.

-Neuronal Action Potentials (Spikes): Consist of short electrical pulses emitted by single neurons during information transmission down to the axon. Such recordings can be achieved by placing electrodes close to the cell or to the axon. When only the activity of one neuron is observed, registers are referred to as single-unit recordings. In contrast, when electrodes are in contact with more than one neuron, recordings are referred to as multi-unit. A "spike train" or "temporal coding" corresponds to a series of spikes generated by a neuron either spontaneously, or in response to an external stimulus.

-Local Field Potentials (LFP): Reflect the activity that is generated by a population of neurons, which can be measured through extracellular electrodes. Neurons in the neighborhood generate overlapping spikes, which causes a large amplitude of LFPs, evidencing the synchronous activity between neurons. This is useful to understand their behavior to a specific stimulus.

Non-invasive recording techniques

Non-invasive techniques are distinguished for recording brain activity in a superficial manner (i.e., without penetrating the skin), so that the process is painless and safe for users. Compared to invasive recording methods, non-invasive alternatives offer several advantages; since they do no represent any risk for their users, such as the hemorrhages, infections, or/and bioincompatibility that may be caused by implanted electrodes, without forgetting the hospitalization, recovery, and the trauma induced by the required surgery. Because of this reason, non-invasive techniques are by far more popular for BCI purposes. Some of these methods include functional Magnetic

Resonance Imaging (fMRI), functional Near Infrared (fNIR), magnetoencephalography (MEG), and electroencephalography (EEG).

Functional Magnetic Resonance Imaging (fMRI)

Blood flow and oxygenation changes are closely related to neural activity [START_REF] Roy | On the regulation of the blood-supply of the brain[END_REF]. Based on these principles, fMRI is intended to measure changes in the blood flow for detecting active brain regions.

Basically, when neurons go active, blood flow carries glucose and oxygen. Oxygen is transported by hemoglobin in the form of oxyhemoglobin, which presents magnetic properties. These cells quickly degrade the local oxygen and increase the level of deoxyhemoglobin (hemoglobin without oxygen).

The advantage of this technique is its high spatial resolution (1 mm 2 ), though it presents a low temporal resolution (5-8 seconds), given that inflow of blood is not instantaneous. In addition, this technique is cumbersome and expensive.

Functional Near Infrared

fNIR measures the absorption of near infrared light (NIR) in the brain for monitoring changes in tissue oxygenation. It offers a suitable spatial resolution (1 cm 2 ), but a relatively low temporal resolution, given that determination of changes in the oxygenation of the tissues can take from 5 up to 8 seconds. fNIR is a portable system, since sensors are placed in the prefrontal cortex by using a sort of "headband", which provides users with comfort. It consists of an emerging technology that is similar to fMRI, but it overcomes many of its drawbacks by being affordable, portable and wearable. However, despite these advantages, fMRI has better spatial resolution, and it offers the possibility of obtaining an image of the entire brain [START_REF] Bunce | Functional nearinfrared spectroscopy[END_REF].

Magnetoencephalography (MEG)

It is based on measurements of the magnetic fields generated by the electrical activity produced in the brain [START_REF] Cohen | Magnetoencephalography: detection of the brain's electrical activity with a superconducting magnetometer[END_REF]. This technique favors a high spatial (millimeters) and temporal (milliseconds) resolution, in addition to a relatively high signal-to-noise ratio (SNR). However, the equipment is expensive, cumbersome and imposes special requirements, such as a magnetically shielded room to reduce artifacts generated by electrical currents or ferromagnetic objects. Despite its drawbacks MEG can be a useful tool for brain-computer interfaces.

Electroencephalography (EEG)

Given that this technique is safe, easy to use, affordable, and offers a high temporal resolution (in the order of milliseconds), it has become the most common method for recording brain signals in humans. Roughly speaking, EEG records the electrical activity generated within different regions of the brain. EEG recordings use a set of electrodes that are positioned on specific locations over the scalp. The majority of electrodes require a conductive gel for improving signals transmission, though recently some varieties offer the alternative of "dry electrodes", which facilitates their use and result more comfortable for users. EEG electrodes can be either passive or active. The former case refers to sensors that are connected to an amplifier by using a cable, whereas the latter type include the amplifier in their taxonomy, which prevents from amplifying the additional noise that is present in the environment and/or caused by movements in the cable. Additionally, some brands offer wireless EEG systems that allow complete mobility to their users. Even though the spatial resolution of this kind of systems is not as convenient (at best 1-2 cm), and considering that they present a low signal-to-noise ratio, this solution still offers several advantages for a daily-based use. Along this chapter we will provide a framework of the principles and use of electroencephalographic (EEG) recordings in the context of motor imagery-based BCIs. The first section aims at describing the international 10-20 system of electrode placement since, during the effectuated experiments, sensors location was based on such standard. Subsequently, the different frequency bands that are of interest for the analysis of EEG signals will be introduced, together with the regions over the motor cortex that are activated during a motor task. In this regard, the physiological aspects that take place within the aforementioned areas, such as Event-Related Desynchronization (ERD)

and Event-Related Synchronization (ERS), as well as their lateralization, will be elucidated.

Furthermore, we will develop the concept of oscillatory power by providing its definition and a computation example, which will be useful to visualize the ERD/ERS patterns that are expected during a motor imagery task.

As a consequence of this review, we will point to the fact that the amount of provided commands to control a motor imagery-based BCI, depends on the different brain states that the system is able to detect; which, in turn, is limited to the number of body parts that users are able of producing suitable patterns with. In this sense, we will introduce the concept of multilabel approach; which dramatically increases the quantity of afforded commands by making use of the combination of two or more body parts at the same time, and thus favoring an intuitive and natural control of the system. At this point we will introduce the main hypothesis of the present study; which states that EEG data recorded during a combined motor imagery can be characterized as the superposition of the activity that is independently generated by the sources associated to the body parts engaged in the task. In order to verify this premise, we will introduce a database, to which we will refer as "4-class database", comprising the EEG recordings of five healthy subjects performing three different motor imageries, namely, right hand, left hand and both hands; which, together with the rest condition, sets up four different brain states. Upon these data we will present a series of statistical results based on the analysis of the power spectrum observed at each electrode during the four different motor tasks that supports the aforementioned hypothesis, as well as a qualitative study of the ERD patterns that are generated under the considered conditions.

From the confirmation of this hypothesis, we will develop in the following chapters a new strategy for multilabel classification; which consists in one of the contributions of the present work.

EEG signals

Most of BCI applications use electric or magnetic signals produced by the brain activity in order to transfer the user's intention to the system, as in this way it is possible to obtain recordings noninvasively from sensors on or above the scalp. The two main noninvasive extracranial techniques used for BCI purposes are electroencephalography (EEG) and magnetoencephalography (MEG), from which EEG is the most widely used considering that it is inexpensive, convenient, and adaptable to diverse environments; thereby, has been developed commercially into portable and even wireless models for a broad diversity of purposes.

The underlying phenomena that is observed through EEG signals is caused by the activity of neurons, which are inherent polarized by membrane transport proteins that pump ions out and into the cells. The exchange of ions occurs constantly between neurons and the extracellular milieu, either to maintain the resting potential or to propagate action potentials. Ions with the same charge repel to each other, and when many of them are pushed out at the same time from different neurons, they can push their neighbors, who push their neighbors, and so on, which generates an electric wave. When such a wave reaches the electrodes on the scalp, its ions can push or pull electrons on the metal surface of the sensors. Since metal easily conducts the exchange of electrons, the difference of voltage between any two electrodes can be measured over time, which can be registered as EEG recordings. The electric potential produced by an individual neuron is far too small to be registered; therefore EEG activity always reflects the synchronous activity of thousands or millions of a network of neurons having similar spatial orientation. When neurons are not similarly oriented, their ions do not contribute in creating a measurable wave.

Most of the EEG activity is believed to be produced by pyramidal neurons within the cortex;

since they are well aligned and fire together. However, voltage magnitude falls as the distance between the source and sensors increases due to the volume conduction phenomena, which points out to the loss and spread of currents across the different physiological tissues during their transit to reach the electrodes, so that activity from deep sources is more difficult to detect that the currents generated near the skull. EEG activity shows oscillations at a variety of frequencies, which present characteristic ranges and spatial distributions that are associated with different brain states.

EEG brain rhythms

EEG patterns are typically characterized in terms of rhythmic activity and transients. Rhythmic activity consists in the synchronized oscillation of a population of neurons that fire together at certain frequency within a particular region, which allows making the link between certain physiological processes and their associated brain function, opening the possibility of using such modulations as interpretable commands to control BCI applications. Most of the EEG activity under standard clinical protocols is found within the range from 1 to 30 Hz, components that are observed below or above this range are likely to be artifacts.

In the context of EEG, rhythmic activity is separated into frequency bands wherein the observed oscillations have demonstrated to be associated with certain biological significance. The range and nomenclature of these bandwidths is stated as follows (see Figure 2.1):

-Delta rhythm (< 4 Hz): It is a normal high amplitude wave detectable in the EEG recordings of healthy subjects. In adults it appears within frontal areas whereas, in the case of children, it is present over posterior regions. Delta waves are commonly found during deep sleep, dreamless sleep, trance, deep hypnosis and, rarely, during wakefulness.

In this sense their presence has been reported along some continuous-attention tasks [START_REF] Kisley | Gamma and beta neural activity evoked during a sensory gating paradigm: effects of auditory, somatosensory and cross-modal stimulation[END_REF]. This rhythm is affected by drugs, alcohol, anxiety, depression and other factors that directly or indirectly disturb sleep. Pathological signs might be related with subcortical, midline or diffuse lesions, as well as with metabolic encephalopathy hydrocephalus.

-Theta rhythm (4 -7.5 Hz): It is found within the parietal-temporal area during drowsiness, REM (i.e., Rapid Eye Movement) sleep, and in the transition from sleep to wakefulness, with higher amplitude in young children. It is also elicited by emotional stress, alcohol ingestion and fatigue [START_REF] Matsuoka | Theta rhythms: state of consciousness[END_REF]. Furthermore, it has been associated with inhibition of elicited responses, since it increases in situations when a person is actively trying to repress certain reactions [START_REF] Kisley | Gamma and beta neural activity evoked during a sensory gating paradigm: effects of auditory, somatosensory and cross-modal stimulation[END_REF]. In pathological cases it is associated with focal subcortical lesions, metabolic encephalopathy, deep midline disorders, and some instances of hydrocephalus. The theta rhythm can be characterized into two types: the hippocampal theta rhythm, observed in the hippocampus of mammals (e.g., cats), and the cortical theta rhythm usually recorded on the human scalp.

-Alpha rhythm (8 -13 Hz): First defined by [START_REF] Chatrian | A glossary of terms most commonly used by clinical electroencephalographers[END_REF] as a rhythm found over posterior regions, generally with higher voltage over the occipital areas, during wakefulness. Power within this band is increased with closed eyes in a state of physical relaxation, over the occipital, parietal and posterior temporal regions. Alpha waves are also associated with inhibition control. Given the link between alpha waves and relaxed mental states, their power increase has become a desirable outcome for some types of biofeedback training. EEG can be used to provide the subject with feedback when alpha waves increase, enabling some individuals to consciously increase this activity.

-Beta rhythm (14 -30 Hz): Beta waves are detectable in almost every healthy adult over the parietal and frontal lobes. Their presence at multiple and varying frequencies with a low amplitude is often associated with active concentration [START_REF] Baumeister | Influence of phosphatidylserine on cognitive performance and cortical activity after induced stress[END_REF]. Over the motor cortex, beta activity is suppressed prior to and during muscle contractions of isotonic movements and it is frequently mixed with a mu rhythm (see below) [START_REF] Baker | Oscillatory interactions between sensorimotor cortex and the periphery[END_REF]. Beta activity is increased when movements have to be resisted or voluntarily suppressed [START_REF] Lalo | Phasic increases in cortical beta activity are associated with alterations in sensory processing in the human[END_REF]. In some cases beta activity is also considered within separately ranges as lower beta band (14-20 Hz), and upper beta band (15 -30 Hz).

-Gamma rhythm (> 30 Hz): This rhythm is typically associated with perception and consciousness, thought it has not fully demonstrated and, according to some scientists, their presence may be caused by electromyographic artifacts [START_REF] Whitham | Scalp electrical recording during paralysis: quantitative evidence that EEG frequencies above 20 hz are contaminated by EMG[END_REF].

-Mu rhythm (8 -12 Hz): These waves are present within the same range as the alpha activity but they differ in their topographical location and they are triggered by endogenous responses. Mu activity decreases during movements or preparation of movements and increases with relaxation. These phenomena are known, respectively, as Event-Related Desynchronization (ERD) [43] and Event-Related Synchronization (ERS) [START_REF] Pfurtscheller | Event-related synchronization (ERS): an electrophysiological correlate of cortical areas at rest[END_REF].

International 10-20 system

The way how sensors are positioned along the scalp for EEG acquisition is collectively referred as electrode placement [START_REF] Wolpaw | BCI: Principles and Practice[END_REF]. A standard technique to cover the different regions of interest over the head is known as International 10-20 system [START_REF] Klem | The ten-twenty electrode system of the international federation[END_REF], which owes its name to the proportional distance separating the electrodes between each other. To this end, nasion (e.g., top of the nose)

and inion (e.g., most prominent point of the occipital bone) are taken as references, from which five aligned electrodes are placed at 10%, 20%, 20%, 20% and 10% of the distance separating both points. By using the same proportion, 16 more electrodes are placed to cover the area inside the diameter defined by the aforementioned references. This configuration leads to a total of 21 electrodes, which are represented in black in figure 2.2. The standard nomenclature that is used to identify electrodes by the major skull bones is indicated on each one of the sensors [START_REF]American electroencephalographic society. guideline thirteen: guidelines for standard electrode position nomenclature[END_REF]. This nomenclature uses a letter for referring to the principal lobe where electrodes are located, and identifies all sensors either on the left or right sides by using, respectively, odd and even numbers, while electrodes along the midline are indicated by "z" ("z"=zero).

In addition, there are two other electrode placements corresponding to systematic extensions of the standard clinical 10-20 standard, namely, 10-10 and 10-5 placement systems [START_REF] Oostenveld | The 5% electrode system for high-resolution EEG and ERP measurements[END_REF]. The total of sensors comprised by the 10-10 system includes the 21 electrodes of the 10-20 placement together with 53 additional electrodes, which are represented in gray in figure 2.2. Possible intermediate electrodes comprised in the 10-5 system are indicated by dots or white circles. The 68 white positions (named by extending the standard nomenclature), together with the 74 electrodes included in the 10-10 system, result in a subset of 142 electrodes that affords a more exhaustive and uniform coverage of the head [START_REF] Klem | The ten-twenty electrode system of the international federation[END_REF]. All these techniques are widely, although not universally, used. In practice, electrodes montage varies across research groups, specially when a large amount of channels (> 64) is required, for which other placement systems have been developed in order to obtain a more regularly spaced sampling of scalp potentials [START_REF] Tucker | Spatial sample of head electrical fields: the geodesic network net[END_REF].

Primary motor cortex

The cerebral cortex is an accessible area to electrode probes as well as to scalp recordings where the functionality of executive motor and communication behaviors take place, thus it represents a region of particular interest in BCI research [START_REF] Miller | Brain-computer interfaces, principles and practice[END_REF].

The characteristics of the cortical surface provide convenient landmarks for the identification of particular regions within the brain, which are often used to guide implementation of intracortical electrodes. Back in 1870, Eduard Hitzig and Gustav Fritsch observed the movement of a dog's limb while the surface of its brain was being electrical stimulated on the opposite side [START_REF] Fritsch | Über der elektrische erregbarkeit des grosshirns[END_REF]. This observation was crucial in many respects. On the one hand, it demonstrated that, as muscles, the brain is electrically excitable. And, on the other hand, by making the link between limb movement and the excitation of a particular region, the important interrogation of whether different brain areas, more precisely the cerebral cortex, were engaged to different functions was confirmed. The relation between particular regions and the control of particular body parts was subsequently coined as somatotopic organization. Based on their experiments, Hitzig and Fritsch ultimately described the area located in the frontal lobe of the cerebral cortex that we now know as the primary motor cortex. Nowadays, in order to identify eloquent cortex (i.e., areas where damage will result in paralysis, loss of sensation, or linguistic impairment), neurosurgeons consistently map the brains of awake human patients undergoing surgical procedures for treatment of severe epilepsy or tumor resection by using electrical stimulation. These techniques were first applied by Wilder Penfield, whose work led to the now well-known map of the motor homunculus [START_REF] Penfield | Somatic motor and sensory representation in the cerebra cortex of man as studied by electrical stimulation[END_REF], [START_REF] Penfield | Some mechanisms of consciousness discovered during electrical stimulation of the brain[END_REF]. This map represents the areas of the motor cortex associated with different motor functions. It is a distorted image of the body in which parts requiring more precisely graded control (e.g., the hand, and more specifically the thumb) have disproportionally large representations. additional ones represented by dots and white circles, which gives a total of 142 electrodes that provide a more exhaustive and uniform coverage of the head. The nomenclature for the standard 10-20 and 10-10 and 10-5 extensions is indicated. Note that electrodes on the right side have even numbers and electrodes on the left side have odd numbers, while electrodes along the midline are indicated by "z" ("z"=zero) [START_REF] Oostenveld | The 5% electrode system for high-resolution EEG and ERP measurements[END_REF].

lateral. Although neighboring body parts are typically represented by the homunculus within neighboring areas, they are spatially distorted given that some of them require a more complex control, which demand a bigger surface. For instance, control of the several muscles involved in facial or hand movement is much more complex than it is for the much larger biceps muscle that flexes the elbow. Consequently, a greater amount of cortical area is engaged in the control of the face or hand than to the upper arm. These general principles of somatotopic organization can be observed in sensory as well as motor areas of the cortex.

Fig. 2.3:

The motor homunculus by Wilder Penfield represents the somatotopic organization that maps the areas within the motor cortex associated with different motor functions. Electrodes C3, C z , and C4 have been included to show the location of the different body parts with respect to the international 10-20 system.

Sensorimotor rhythms

Sensorimotor rhythms (SMRs) are oscillations in the electric or magnetic fields registered over sensorimotor cortices (i.e., posterior frontal and anterior parietal areas) [START_REF] Pfurtscheller | Niedermeyer's electroencephalography[END_REF], which are mostly

SENSORIMOTOR RHYTHMS

found between three main frequency bands: mu (8 -12 Hz), beta (14 -30 Hz) and gamma (> 30 Hz). In the case of EEG, recordings are mostly limited to mu, beta and lower-frequency gamma activity, whereas ECoG and MEG signals can be modulated within higher frequencies [START_REF] Wolpaw | BCI: Principles and Practice[END_REF].

Event-related desynchronization

The transitory amplitude decreasing of oscillations in the alpha and beta bands taking place in localized regions and in direct relation to an event is known as event-related desynchronization (ERD) [START_REF] Pfurtscheller | Event-related cortical desynchronization detected by power measurements of scalp EEG[END_REF]. This phenomenon can be observed not only during, but also before, sensory and cognitive processing and motor behavior of the awake brain [START_REF] Van Winsum | The functional significance of event-related desynchronization of alpha rhythms in attentional and activating tasks[END_REF]- [START_REF] Mohl | The role of the right parietal region in a movement time estimation task[END_REF]. Therefore, ERD is not only an electrophysiological correlate of cortical activation related to stimulus processing or motor output, but is also characteristic for cortical areas or neural structures preparing to process sensory information or ready and prepared to execute a motor command [START_REF] Pfurtscheller | Motor imagery and direct brain-computer communication[END_REF].

Event-related synchronization

The opposite phenomenon is known as event-related synchronization (ERS), which consists of a localized increase of power in the alpha and beta bands generated by bursting activity [START_REF] Pfurtscheller | Synergetics of rhythms in biological systems[END_REF]. Both ERS and ERD have a specific topographical distribution depending upon the state of the brain, and they can be either externally (e.g., by stimuli) or internally (e.g., by voluntary behaviour) triggered; which in terms of the BCI terminology correspond, respectively, with synchronous and asynchronous events.

Both cortical and subcortical structures take part as extensive networks responsible for the SMRs observed within the cortex. The use of depth electrodes has revealed mu and beta activity along the thalamus, subthalamic nucleus, and penduculopontine area [START_REF] Androulidakis | Oscillatory activity in the penduculopontine area of patients with parkinson's disease[END_REF]- [START_REF] Williams | Dopamine-dependent changes in the functional conectivity between basal ganglia and cerebral cortex in humans[END_REF]. The relationships among these areas are complex. For instance, during movement, mu ERD occurs within the motor cortex and mu ERS occurs in the subthalamic nucleus, whereas beta ERD occurs uniformly in motor cortex, thalamus, and the subthalamic nucleus [START_REF] Klostermann | Task-related differential dynamics of EEG alphaand beta-band synchronization in cortico-basal motor structures[END_REF].

Lateralization

As mentioned, ERD corresponds with an activity pattern of sensorimotor rhythms occurring during real or imagined movements [START_REF] Pfurtscheller | Motor imagery and direct brain-computer communication[END_REF]. ERD% lateralization, i.e. the difference between the relative power changes of the contra-and ipsilateral hemispheres, constitutes the basis of the majority of binary MI-based BCIs [START_REF] Pfurtscheller | Motor imagery and direct brain-computer communication[END_REF], [START_REF] Zich | Lateralization patterns of covert but not overt movements change with age: an EEG neurofeedback study[END_REF]- [START_REF] Zich | Real-time EEG feedback during simultaneous EEG-fMRI identifies the cortical signature of motor imagery[END_REF].

Previous studies also reported a significant ERD over the contralateral side and a significant ERS over the ipsilateral side of the brain for planned and terminated movements, respectively. For example, alpha and beta ERD begin contralaterally about 1.5 -2 s before self-paced movement, becoming bilateral during movement execution [START_REF] Pfurtscheller | Event-related EEG/MEG synchronization and desynchronization: basic principles[END_REF], [START_REF] Leocani | Event-related coherence and event-related desynchronization/synchronization in the 10 hz and 20 hz EEG during self-paced movements[END_REF]- [START_REF] Serrien | Coordination constraints during bimanual versus unimanual performance conditions[END_REF]. The significant ipsilateral ERS was found to be characteristic of the beta band for both after physical and imagined termination of motor tasks, which has been referred as post-movement rebound [START_REF] Pfurtscheller | Event-related EEG/MEG synchronization and desynchronization: basic principles[END_REF], [START_REF] Pfurtscheller | On the existence of different types of central beta rhythms below 30 hz[END_REF], [START_REF] Stancàk | Desynchronization of cortical rhythms following cutaneous stimulation: effects of stimulus repetition and intensity, and of the size of corpus callosum[END_REF]- [START_REF] Neuper | ERD/ERS patterns reflecting sensorimotor activation and deactivation[END_REF].

Spatial mapping of ERD/ERS

EEG multichannel recordings are frequently recorded against a common reference electrode, thus data are reference-dependent. In order to convert reference-dependent raw data into reference-free data, different techniques are available, such as Common Average Reference (CAR), Hjorth

Laplacian Reference (HR), and Local Average Reference (LR) [START_REF] Pfurtscheller | Event-related synchronization (ERS): an electrophysiological correlate of cortical areas at rest[END_REF], [START_REF] Da | Electroencephalography: basic principles, clinical applications and related fields[END_REF], [START_REF] Pfurtscheller | Source localization using event-related desynchronization (ERD) within the alpha band[END_REF]. Reference-free topographic maps typically provide better defined patterns than referential representations and are especially recommended for sensorimotor rhythms analysis [START_REF] Pfurtscheller | Event-related EEG/MEG synchronization and desynchronization: basic principles[END_REF]. For spatial mapping of ERD/ERS it is possible to use a spherical head model together with an interpolation method (e.g., cubic interpolation).

Event-related potentials

Different events can induce time-locked changes in the activity of neural populations, which are typically known as Event-Related Potentials (ERPs), from which the most remarkable correspond to sensory stimuli responses [START_REF] Pfurtscheller | Event-related EEG/MEG synchronization and desynchronization: basic principles[END_REF]. Averaging techniques are generally applied in order to enhance the morphology of ERPs. The underlying premise is that the patterns of interest have a relatively constant time-delay to the stimulus, while the ongoing activity behaves as additive noise, thus the averaging process improves the signal-to-noise ratio. However, there are some omissions that do not hold to this assumption; on the one hand, it is known that Evoked Potentials (EPs) can 3). Spatial mapping was carried out by using a cubic interpolation method [START_REF] Höhne | BBCI toolbox-offline analysis[END_REF]. Note that each pair of maps was independently normalized. In the same way, it is possible to observe the lateralization of the ERS patterns associated to the motor tasks (around C3, see fig. 2.3). Spatial mapping was carried out by using a cubic interpolation method [START_REF] Höhne | BBCI toolbox-offline analysis[END_REF]. Note that each pair of maps was independently normalized.

be considered to result from a reorganization of the phases of the ongoing EEG signals [START_REF] Sayers | The mechanism of auditory evoked eeg responses[END_REF].

Furthermore, it was also demonstrated that visual stimuli can reduce the amplitude of the ongoing EEG activity [START_REF] Vijn | Visual stimulation reduced eeg activity in man[END_REF], showing that the model assuming that an ERP can be represented by a signal added to uncorrelated noise does not hold in general. And on the other hand, considering that the changes evoked by certain events that block or desynchronize the ongoing alpha activity are timelocked to the event but not phase-locked, they cannot be extracted by a simple linear method, such as averaging, but may be detected by frequency analysis. In summary, event-related phenomena cause frequency specific changes of the EEG activity which, in general terms, consist either of decreases or of increases of power in particular frequency bands. This may be considered to be caused, respectively, by a decrease (ERD) or an increase (ERS) in synchrony of the underlying neural populations [START_REF] Pfurtscheller | Event-related synchronization (ERS): an electrophysiological correlate of cortical areas at rest[END_REF], [START_REF] Pfurtscheller | Event-related cortical desynchronization detected by power measurements of scalp EEG[END_REF], [START_REF] Pfurtscheller | Graphical display and statistical evaluation of event-related desynchronization (ERD)[END_REF].

Contrary to the classical ERPs, which can be characterized as a series of transient postsynaptic responses of main pyramidal neurons responding to a particular stimulus, ERD/ERS phenomena can be considered as the result of changes in one or more variables that control oscillations in neural networks. The properties of EEG oscillations are determined by [START_REF] Pfurtscheller | Event-related EEG/MEG synchronization and desynchronization: basic principles[END_REF], [START_REF] Da Silva | Neural mechanisms underlying brain waves: from neural membranes to networks[END_REF],

[84]:

1. the intrinsic membrane properties of the neurons and the dynamics of synaptic processes;

2. the strength and extent of the interconnections between the network elements, most often formed by feedback loops. Different kinds of feedback loops can be distinguished: involving thalamo-cortical or cortico-cortical either at short or at long distances;

3. the modulating influences from general or local neuro-transmitter systems.

Time course of ERD/ERS

Typically, ERD/ERS measurements embraced in particular frequency bands are represented relative to the power of the same EEG channels recorded during a baseline period a few seconds before the event occurs, as a percentage magnitude. Event related changes in ongoing EEG need time to thrive and to recover, especially when slow alpha band rhythms are involved, therefore the interval between two consecutive events must last at least some seconds. For voluntary limb movement studies, it is recommended to use an inter-interval of at least 10 s and, in the case of a foreperiod reaction time task, this interval must be even longer.

The classical method to compute the time course of ERD/ERS includes the following steps [START_REF] Pfurtscheller | Event-related EEG/MEG synchronization and desynchronization: basic principles[END_REF]:

1. bandpass filtering of all event-related trials within the frequency band of interest;

2. squaring of the amplitude samples to obtain power samples;

3. averaging of power samples across all trials;

4. averaging over time samples to smooth data and reduce variability.

To obtain ERD/ERS values in terms of percentage, the power within the frequency band of interest following the event is given by x, whereas that of the preceding baseline or reference period is given by BL. ERD/ERS values are computed as the percentage of power decrease or increase, respectively, according to the following expression:

ERD/ERS% = x 2 -BL 2 BL 2 × 100, (2.1)
where x 2 is the mean of the squared EEG signal, BL 2 is the mean of a baseline sliding window taken at the beginning of the corresponding trial, and ERD/ERS% is the percentage of the oscillatory power estimated at each step of the sliding window. Oscillations greater than the baseline during the MI task produce positive percentage values that are assumed as a modulation of ERS, whereas the negative values generated by smaller oscillations are considered as a modulation of ERD. observe that this modulation (highlighted in blue) takes place approximately 0.5 s before the motor imagery starts, and that it lasts for some seconds following its termination, after which it begins to recover slowly to the reference line and going beyond to positive values (highlighted in red), which correspond to the ERS modulation characteristic of the post-movement rebound.

Combined movements

In contrast to simple motor tasks, combined imagined movements consider the use of two or more body parts at the same time. In this way it is possible to dramatically increase the number of commands afforded by a BCI system. However, as presented along the introduction of this document, most of the research has been centred on the EEG patterns induced by simple motor imageries, and less work has been reported regarding the study of compound tasks. In this sense, to our best knowledge, the question whether the EEG activity induced by combined imagined movements can be characterized as the superposition of the activity that is independently generated by the sources associated to the body parts engaged in the task. To first approach this question, we have a created a 4-class database comprising the rest condition and motor imageries of left hand, right hand and both hands together. With which a statistical analysis was performed in order to determine whether or not the activity in the associated areas (around electrodes C3 and C4, see fig. 2.3) behaves the same in the case of the combined task with respect to each one of the simple motor imageries.

Paradigm and time scheme

The EEG cap was fitted with 26 electrodes, namely, Fp1; Fp z ; Fp2; F z ; FC5; FC3; FC1; FC z ; FC2; FC4; FC6; C5; C3; C1; C z ; C2; C4; C6; CP5; CP3;CP1; CP z ; CP2; CP4; CP6, and P z , re-referenced with respect to the common average reference across all channels and located over the international 10-10 system positions to cover the primary sensorimotor cortex (see fig. ). The rest condition was also considered and it was indicated by the absence of arrows. After second 6, the task cue disappeared and the crosses were remaining for the next 6 seconds indicating the pause period before the next trial started.

ERD/ERS% analysis

Each trial was considered 2 seconds before the task cue appeared on the screen and 6 seconds after it disappeared, so that each segment corresponds to a 14s-length signal over 26 electrodes. All segments were filtered using a 5th-order Butterworth filter within the mu and beta bands, in conformance to the nature of SMRs described in section 2.5 in page 24.

The time course of ERD/ERS% of each segment was computed by squaring its samples and smoothing the resulting signal through a 2-second-sliding window with a 125 ms shifting step according to equation 2.1. Figure 2.9 presents a topographic representation of the grand average across subjects of ERD/ERS% values, which were computed in terms of equation 2.1 over all electrodes within the mu range during the 6-second period of the four motor tasks. Over the contralateral hemisphere it can be observed that, during simple motor imageries, ERD% values are similar to those generated in the corresponding side throughout the combined task, which considering the entire set of electrodes, present the minimum values around C3 and C4 during the course of the three motor tasks. In the case of left hand motor imageries, the minimum ERD% values are induced around electrode C4, contrastingly to the activity observed during right hand motor imageries, which presents the lowest ERD around electrode C3, whereas for the combined task, the minimum power is found simultaneously over electrodes C3 and C4 (see fig. 2.3). On the other hand, over the ipsilateral hemisphere, simple motor imageries present similar ERS% values to the observed ones during the rest condition. Furthermore, in the ipsilateral hemisphere, simple motor imageries produce similar ERS% values to the ones estimated on both sides for the rest condition, which is represented by the green trace.

The comparison between the rest condition together with the three motor imageries once the motor tasks were finished is shown in figure 2.11. Which shows, during the 6 seconds following the motor tasks, the corresponding spatial maps computed over all electrodes of the grand average across subjects of the ERD/ERS% values estimated in terms of equation 2.1 within the beta range. It can be observed that the activity is distributed in a similar way as during the motor tasks (i.e., around electrodes C3 and C4, see fig. 2.9), but with much higher values; which indicates that the corresponding regions present an ERS modulation characteristic of the termination of motor tasks (i.e., post-movement rebound). In the case of left hand motor imageries, the maximum ERS% values are induced around electrode C4, contrarily to the activity observed after right hand motor imageries, which presents the highest ERS around electrode C3, whereas for the combined task, the maximum power is found simultaneously over electrodes C3 and C4. On the other hand, over the ipsilateral hemisphere, simple motor imageries present similar ERS% values to the observed ones during the rest condition.

The grand average across subjects of the time course of ERD/ERS within the beta band for all motor tasks over electrodes C3 and C4 is shown, respectively, in figure 2.12-A and figure 2.12-B.

Similar to figure 2.10, though with smaller magnitude, there are significant ERD% values over the contralateral side during the course of simple motor imageries, which present a comparable morphology to the one induced by the combined task on both sides. However, it worth noting that, once the motor tasks have finished, traces go back to positive values much faster than in the mu band, phenomena that has been already reported by [START_REF] Avanzini | The dynamics of sensorimotor cortical oscillations during the observation of hand movements: an EEG study[END_REF]. Furthermore, the ERS% values observed after the simple motor tasks termination, which correspond to the so-called post-movement rebound, present the same lateralization that the ERD% values observed during the course of motor imageries, so that in the case of left hand motor imageries, the highest ERS% values take place over electrode C4 (right hemisphere), whereas for the right hand motor imageries, the highest ERS% values are induced in electrode C3 (left hemisphere). In the case of the combined task, the post-movement rebound is appreciable on both sides although, in comparison with the one generated by simple motor imageries, it presents significantly smaller values. On the other hand, over the ipsilateral hemisphere after simple motor imageries, there is no post-movement rebound and the time course of the corresponding traces looks similar to the one generated by the rest condition.

Statistical analysis

From figure 2.10 and figure 2.12 it can be observed that the lowest ERD% values are reached within the interval [1-2 s]. For this reason, in order to better evaluate the differences between conditions that take place during the course of the considered motor tasks, all segments were analyzed only during such a period. Furthermore, since the ERD% values characteristic of a motor task execution are well observed over both mu and beta bands, all trials were band-pass filtered in the range from 8 to 25 Hz (e.g., mu-band + beta-band) by using a 5th order Butterworth filter, within which the power spectrum was computed across all channels. Subsequently, the aforementioned range was scanned by taking subintervals of 2 Hz each to compute, in terms of R 2 , the correlation between the power spectrum magnitude over all pairs of conditions at each electrode [START_REF] Fok | An EEG-based brain computer interface for rehabilitation and restoration of hand control following stroke using ipsilateral cortical physiology[END_REF].

High R 2 values indicate that there is a significant difference between the considered motor Fig. 2.12: ERD/ERS% time series within the beta range for all motor tasks over electrodes C3 and C4. ERD/ERS values were estimated according to equation 2.1 over the grand average across subjects (n=5). The oscillatory power on electrodes C3 and C4 is presented, respectively, in (A) and (B).

tasks at the given frequency and over the corresponding positions. On the contrary, small R 2 values show that the conditions under analysis do not present a significant difference, and therefore are similar between each other [START_REF] Lindig-León | Comparison of sensorimotor rhythms in EEG signals during simple and combined motor imageries over the contra and ipsilateral hemispheres[END_REF]. The situation is similar as the one described for the left hand motor task but with the highest R 2 values found in the other way around, since the ERD/ERS% values occur in the same manner but with respect to the right hemisphere (electrode C3). Figure 2.13-E shows the comparison between left hand and right hand motor imageries. In this case, the highest R 2 values are also at 12 Hz and simultaneously distributed around electrodes C3 and C4. This is because the brain activity within both hemispheres presents an opposite behavior during the course of these two types of motor imageries and, therefore, there is always a significant difference between each other.

In the same way, the comparison between both hands motor imageries and the rest condition presented in figure 2.13-F, shows high R 2 values in both hemispheres; since the both hands motor task generates ERD% values in both sides, whereas the rest condition produces ERS%. From figure 2.13 it can be observed that, for all comparisons, the most significant R 2 values occur around 12 Hz. Therefore, in order to visualize the correlation between electrodes in a more intuitive way, R 2 values within this frequency were used to compute the topographic maps that are shown in figure 2.14, where the comparisons for each pair of conditions are presented.

Along this section, it has been reported that during the course of the considered motor imageries, brain activity has important components around 12 Hz, and that it is distributed over the primary sensorimotor cortex according to the body part that is engaged at each task, especially around electrodes C3 and C4 (see fig. 2.3). With this in mind, it is possible to estimate the level of significance between the considered conditions. To this end, the power spectrum of the one-second segment (i.e., from 1 to 2 s) of each filtered trial within the frequency range from 8 to 25 Hz over electrodes C3 and C4 was computed for each subject. Subsequently, all values at 12 Hz belonging to the same condition at each electrode were gathered together and, an unpaired t-test was applied to determine with a 5% level of significance, whether or not two different conditions present comparable activity over the left (electrode C3) and right (electrode C4) hemispheres. subjects, except for subject 2, which exhibits no significance difference in most of the cases, that simple motor imageries in their contralateral hemisphere are comparable with the combined task in the corresponding side; and that simple motor imageries in their ipsilateral hemisphere are comparable to the rest condition in the corresponding side. with what has been already observed along section 1.6.2; where we have showed, in a qualitative way, the topographic similarities of the combined motor imageries involving both hands with each type of simple hand motor imageries at their contralateral side. In this way, the statistical analysis can be used to establish that EEG data recorded during a combined motor imagery can be characterized as the superposition of the activity that is independently generated by the sources associated to the body parts engaged in the task, which confirms the hypothesis stated at the beginning of this chapter, and upon which we will develop the multilabel classification models contributed by the present work.

3

Robotic arm control CP6, and P z , re-referenced with respect to the common average reference across all channels and located over the extended international 10-20-system positions to cover the primary sensorimotor cortex (see Figure 2.7).

All experiments were carried out with the informed consent agreement of each participant and following the statements of the WMA declaration of Helsinki on ethical principles for medical research involving human subjects.

Paradigm and time scheme

Subjects were seated in a comfortable chair with the arms at their sides in front of a screen

showing the task cue to be performed, which consisted of one of the eight mental states that it is possible to generate with all combinations including the use of the right hand, left hand, both feet together, and rest (i.e., right hand; left hand; feet; both hands together; left hand in combination with feet; right hand in combination with feet; both hands together with feet, and rest). The whole session consisted of 3 runs, containing each one 10 trials per task, which totals

Oscillatory power analysis

Similarly to the analysis presented along section 2.6.2 in page 27, trials were considered 2 seconds before the task cue appeared on the screen and 6 seconds after it disappeared, so that each segment corresponds to a 14s-length signal segment over 26 electrodes. All trials were filtered using a 5th-order Butterworth filter within the mu frequency band, in conformance to the nature of SMRs described in the same section.

The time course of ERD/ERS% of each segment was computed by squaring its samples and smoothing the resulting signal through a 2-second-sliding window with a 125 ms shifting step If we take a closer look over the time course of ERD/ERS% over electrode C3 (i.e, the activity source associated with the use of the right hand) presented in Figure 3. and both hands). However, in comparison to the contrast observed over electrode C3 between the tasks involving the use of the right hand and those that do not involve it, differences are not as significant. In fact, we can observe that the use of both left and right hands induce some desynchronization over the central region, possible because the associated sources are quite close to electrode C z , which might cause it to incorporate some of this activity.

Finally, the same comparison is illustrated in Figure 3.5 for the time course of ERD/ERS% generated over electrode C4 (i.e, the activity source linked to the use of the left hand). In this case, similarly to the analysis over electrode C3, it is possible to observe a significant contrast between the low ERD/ERS% values obtained from the activity induced by the motor imageries involving the use of the left hand (i.e., left hand, left hand in combination with feet, both hands, and both hands in combination with feet) from the high ERD/ERS% values associated to the tasks that do not involve the use of the left hand (i.e., rest, feet, right hand, and right hand in combination with feet).

Statistical analysis

From the previous analysis and in correspondence to what we have reported in section 2.6.2 for the case of the 4-class database, it can be observed that the oscillatory power around the main activity sources associated to the body parts included in this paradigm presents similar behaviour when the corresponding limbs are engaged in a task, independently from the fact of whether the motor imagery is simple or combined. In order to formalize this analysis in a quantitative manner, we have performed the same procedure described in section 2.6.3 for evaluating the correlation between electrodes for all comparisons among each pair of motor imageries. To this end, all trials were band-pass filtered by using a 5th order Butterworth filter within the frequency range from 8 to 25 Hz (e.g., mu-band + beta-band, since ERD% values appear over both ranges). Subsequently, all segments were selected within the time interval from [2-5 s], given that for all subjects the lowest ERD% values were found over this period, within which the power spectrum was computed across all channels. Finally, the mu-beta range was scanned by taking subintervals of 2 Hz each to compute, in terms of R 2 , the correlation between the power spectrum magnitude over all pairs of conditions at each electrode [START_REF] Fok | An EEG-based brain computer interface for rehabilitation and restoration of hand control following stroke using ipsilateral cortical physiology[END_REF]. Note that for the analysis over the trials contained in the 8-class database we have selected a wider time window than the one that was used for the 4-class database. This is due to the morphology that was observed for this set of subjects, which in general presented longer desynchronization periods.

Similarly to what we have found in the case of the 4-class database, for all subjects the largest R 2 values were located around 12 Hz. Figure 3.6 shows the topographic maps generated for one of the subjects throughout all comparisons between the power spectrum magnitudes at this frequency for each pair of conditions.

High R 2 values indicate that there is a significant difference between the considered motor tasks over the corresponding positions. On the contrary, small R 2 values show that the conditions under analysis do not present a significant difference, and therefore are similar between each other [START_REF] Lindig-León | Comparison of sensorimotor rhythms in EEG signals during simple and combined motor imageries over the contra and ipsilateral hemispheres[END_REF]. It can be observed that the highest R 2 values occur when the motor imageries differ on the limbs that are engaged during their execution. For instance, if we analyze the first row of Figure 3.6, which corresponds to the comparison between the power spectrum calculated during rest and each one of the values generated by the other motor tasks, we notice that the most significant differences occur just around the activity sources associated to the limbs that are engaged during each one of the motor imageries. This occurs because during rest none of the sources induces ERD modulation, whereas for the other tasks this activity appears over the associated sources, which causes such differentiation. The same behaviour follows for all pair of comparisons, which Note that the motor imageries involving the use of the feet (i.e., feet, left hand in combination with feet, right hand in combination with feet, and both hands in combination with feet) present lower ERD/ERS% values than those that do no it include it (i.e., rest, left hand, right hand, and both hands). However, this contrast is not as significant as the one observed on electrode C3, possibly due to the fact that the regions associated to both the left and right hands are quite close to electrode C z , which might cause it to incorporate some of this activity.

coincides with what it was observed over the ERD/ERS% time series from the previous section.

Additionally to the correlation analysis between conditions, we present a series of box plots for comparing the power spectrum magnitudes at 12 Hz between all motor tasks over the regions associated to the main sources. Again, it is possible to observe that motor imageries involving the use of the limb associated to each electrode present similar low power values at these positions, whether those that are not engaged generate higher values. This is particularly visible from figure 3.7, which presents the corresponding analysis over electrode C3 (i.e., the region associated to the right hand). Note that the distribution of the mean values confirm that there is no significant difference between motor imageries including the use of the right hand over this region, independently from the fact of whether the motor imagery is simple or combined. In contrast there is a significant difference with all the other tasks that do not include the use of the right hand which, in turn, present a uniform distribution between them. The same analysis is presented for electrode C z (i.e., the region associated to the feet) and electrode C4 (i.e., the region associated to the left hand) in figures 3.8 and 3.9 respectively. The distribution of data shows that there is no significant difference between motor imageries including the use of the right hand over this region, independently from the fact of whether the motor imagery is simple or combined. In contrast there is a significant difference with all the other tasks that do not include the use of the right hand which, in turn, present a uniform distribution between them. The distribution of data shows that there is no significant difference between motor imageries including the use of the left hand over this region, independently from the fact of whether the motor imagery is simple or combined. In contrast there is a significant difference with all the other tasks that do not include the use of the left hand which, in turn, present a uniform distribution between them.

4

Signal Processing As previously mentioned, a BCI system aims at recognizing determined patterns in the user's brain activity that are related to a mental task that the person performs in order to transmit the associated command. To this end two central stages are required, namely, feature extraction and feature classification. The first concept consists in determining the characteristics that confer distinctiveness to the brain activity that is elicited during the corresponding task, so that it can be recognized from other states. In most of the cases this process reduces the dimensionality of data by excluding redundant information, and it enhances the attributes of the signal by exploiting the aspects that are related to the physiological phenomena. The derived values are projected into the feature space, where elements from different conditions are expected to be located as separated as possible, so that distinct groups can be easily identified. Once features have been extracted from the original input data, the next step consists in creating a classification algorithm (classifier) that is able to categorize such elements according to their position within the classification space, so that the system can proceed to translate the corresponding command.

In which favors the selection of a discriminative subject-specific frequency range. In this way, the dimensionality of the feature space is drastically increased and some redundant information is generated. At this point we will introduce the concept of f eature selection, which allows choosing the most descriptive elements, so that computational costs are optimized and models are prevented of considering information that might interfere with discriminative patterns. In this sense we will present two mutual-information-based methods for feature selection, namely the Mutual Information-based Best Individual Feature (MIBIF) approach, and the Information Theoretic Feature Extraction (ITFE) algorithm. The first one is applied together with the FBCSP method, which represents a powerful solution that has drawn the attention of the community after winning the BCI competition IV [START_REF] Tangermann | Review of the bci competition iv[END_REF]. The second one is presented in combination with the CSP by JAD algorithm as an non-heuristic strategy to select the filters that result from the joint diagonalizaton of multiple covariance matrices.

Since the principle of the CSP algorithm is formulated upon the information contained within the covariance matrices of the EEG signals, we will first provide an insight into their definition and interpretation. Subsequently, along the same path, we will formulate the eigenvalue problem that solve the simultaneous diagonalization of the covariance matrices related to two conditions, which represents the foundation of the CSP algorithm and provides the contrast that is useful for the further classification.

The second section of this chapter is dedicated to present the notion of classification, for which two commonly used algorithms will be described, namely, Linear Discriminant Analysis (LDA) and Support Vector Machines (SVM). These classification algorithms are based on a discriminant function, which implies that they formulate a relation that maps each feature vector directly onto a class label. There are other classification methods that do not use discriminant functions to categorize input elements, such as neural networks, decision trees and, among others, distance-based algorithms. The latter approach has drawn special attention due to the possibility of applying it in the context of Riemannian geometry, which has also demonstrated to be potent at differentiating sensorimotor rhythms. In this regard and, in order to explore other options, the present work considers as well this solution and provides a brief outline as an alternative to the CSP algorithm; which endows a classification that can be directly achieved within the manifold of symmetric and positive definite matrices given by the covariance matrices of the EEG signals that are considered as descriptors, from where it is possible to extract the spatial information without using spatial filtering, and thus simplifying the classification process.

Finally, the last section of the chapter introduces a new classification method for motor imagerybased BCI systems that makes use of the combination of the CSP algorithm and Riemannian geometry, which represents one of the contributions of the present work. The motivation of this approach lies in the fact that the covariance matrices of the projected signals preserve the differences attained through the spatial filtering, making it advantageous to apply the CSP algorithm before mapping them into the Riemannian space, so that more discriminative features are generated.

Feature extraction

Covariance matrix

Let x ∈ R n be a random vector with multiple dimensions containing a scalar random variable within each element, then its covariance matrix is defined by the matrix whose element in the (i, j) position corresponds to the covariance between the i th and j th elements of x. In an intuitive manner, the covariance matrix generalizes the concept of covariance to multiple dimensions. For instance, if the elements of the column vector

x =       x 1 x 2 . . . x n      
are random variables, each with finite variance, then the covariance matrix Σ is the matrix whose (i,j) entry is given by

Σ i,j = cov(x i , x j ) = E[(x i -µ i )(x j -µ j ) ]
where

µ i = E[x i ]
is the expected value of the i th element of vector x. In other words,

Σ =       E[(x 1 -µ 1 )(x 1 -µ 1 )] E[(x 1 -µ 1 )(x 2 -µ 2 )] . . . E[(x 1 -µ 1 )(x n -µ n )] E[(x 2 -µ 2 )(x 1 -µ 1 )] E[(x 2 -µ 2 )(x 2 -µ 2 )] . . . E[(x 2 -µ 2 )(x n -µ n )] . . . . . . . . . . . . E[(x n -µ n )(x 1 -µ 1 )] E[(x n -µ n )(x 2 -µ 2 )] . . . E[(x n -µ n )(x n -µ n )]      
Since the covariance of the i th element with itself is its variance, each element on the (i,i) position within the covariance matrix diagonal corresponds to the variance of the i th element.

Moreover, because cov(x i , x j ) = cov(x j , x i ) every covariance matrix is symmetric. In addition, every covariance matrix is positive semi-definite, which implies that their eigenvalues are nonnegative [START_REF] Taboga | [END_REF].

The distribution of data is reflected in their covariance matrix (see Figure 4.1). The diagonal spread is related to the covariance (off-diagonal of the matrix), whereas the dispersion in the horizontal and vertical axis is associated to the variance (diagonal of the matrix).

Normal distribution and eigenvalue decomposition

The normal distribution, also referred as Gaussian, is a widely used model for the distribution of continuous variables, which appears in several different contexts and can be motivated from many perspectives [START_REF] Bishop | Pattern Recognition and Machine Learning[END_REF]. For a single variable x, the normal distribution can be written in the following terms

N (x|µ, σ 2 ) = 1 (2πσ 2 ) 1/2 exp - 1 (2σ 2 ) (x -µ) 2 , (4.1)
where µ and σ 2 correspond, respectively, with the mean and variance of x. In the case of a D-dimensional vector x, the multivariate normal distribution takes the form

N (x|µ µ µ, Σ Σ Σ) = 1 (2π) D/2 1 |Σ Σ Σ| 1/2 exp - 1 2 (x -µ µ µ) Σ Σ Σ -1 (x -µ µ µ) , (4.2)
where µ µ µ corresponds with the D-dimensional mean vector of x and Σ Σ Σ with its D × D covariance matrix, whose determinant is represented by |Σ Σ Σ|.

Now, the functional dependence of the normal distribution on x, takes the quadratic form of the Mahalanobis distance

∆ 2 = (x -µ µ µ) Σ Σ Σ -1 (x -µ µ µ), (4.3)
which is found within the exponent of equation 4.2. The Mahalanobis distance is a measure of the distance between a point x i and the distribution N , and it represents a multidimensional generalization of the idea of measuring how many standard deviations away x i is from µ µ µ [START_REF] Mahalanobis | On the generalised distance in statistics[END_REF].

If x i is at µ µ µ, then the distance is zero, and it grows as x i moves away. Along each principal component axis, it measures the number of standard deviations from x i to µ µ µ. If each of these axes is rescaled to have unit variance, (i.e., when Σ Σ Σ is the identity matrix), then the Mahalanobis distance reduces itself to the Euclidean distance in the transformed space.

We can consider Σ Σ Σ as a symmetric matrix without loss of generality, since any antisymmetric component would disappear from the exponent. Now consider the eigenvalue decomposition for the covariance matrix

Σ Σ Σu i = λu i , (4.4)
where i = 1, ..., D. Since Σ Σ Σ is a real and symmetric matrix, its eigenvalues will be real and, its eigenvectors, can be chosen to form an orthonormal set, so that

u i u j = I ij , (4.5)
where I ij is the i, j element of the identity matrix satisfying

I ij =    1, if i = j 0, otherwise. (4.6)
The covariance matrix Σ Σ Σ can be represented as an expansion in terms of its eigenvectors, so that

Σ Σ Σ = D i=1 λ i u i u i , (4.7)
then as well, the inverse covariance matrix Σ Σ Σ -1 , can be rewritten in the form

Σ Σ Σ -1 = D i=1 1 λ i u i u i . ( 4
∆ 2 = D i=1 y i 2 λ i , (4.9)
where y i has been defined as

y i = u i (x -µ µ µ) (4.10)
The variables {y_i} can be interpreted as a new coordinate system defined by the orthonormal vectors u i that are shifted and rotated with respect to the original x_i coordinates. By forming the vector y = (y i , ..., y D ) we obtain

y = U(x -µ µ µ) (4.11)
where U is a matrix whose rows are given by u i . From equation 4.5 it follows that U is an orthogonal matrix (i.e., it satisfies UU = I, and also U U = I), where I is the identity matrix.

The quadratic form, and hence the normal density, will be constant on surfaces for which equation 4.10 is constant. If all of the eigenvalues λ i are positive, then these surfaces represent ellipsoids, with their centres at µ µ µ and their axes oriented along u i , and with scaling factors in the directions of the axes given by λ 1/2 , as illustrated in figure 4.2.

For the normal distribution to be well defined, it is necessary for all of the eigenvalues λ i of the covariance matrix to be strictly positive, otherwise the distribution cannot be properly normalized. A matrix whose eigenvalues are strictly positive is said to be positive definite.

Common Spatial Patterns (CSP)

A very convenient and widely used algorithm to extract motor imagery features from EEG signals is known as Common Spatial Pattern (CSP); which generates a series of spatial filters that are applied to decompose multi-dimensional data into a set of uncorrelated components [START_REF] Blankertz | Optimizing spatial filters for robust EEG single-trial analysis [revealing tricks of the trade[END_REF]. These filters aim at extracting elements that simultaneously maximize the variance of one class, while minimizing the variance of the other one.

The first step of the CSP method involves the computation of the normalized covariance

matrix Σ k Σ k = E k E k trace(E k E k ) , (4.12)
where

E k ∈ R C×T , k ∈ {1,
Σ 1 W = λΣ 2 W (4.13)
with λ denoting the diagonal matrix containing the eigenvalues of Σ 1 . The spatial filtered signal Y k can be obtained from the EEG trials E k as

Y k = W E k . (4.14)
Each column vector of W represents a spatial filter, and there are as many as the number of channels in the EEG signal. Note that the CSP filters are paired and that not all of them are relevant for the subsequent classification; therefore, after being sorted according to the values of λ (either in an ascending or descending order instinctively), only the first m and the last m columns of W are selected (In the present work, for all methods and subjects m = 2 pairs of filters were considered [START_REF] Ang | Filter bank common spatial pattern algorithm on BCI competition IV datasets 2a and 2b[END_REF]).

The selected elements generate the spatial filters coefficient matrix W, from which the m pairs of CSP features of the i th trial for the band-pass filtered EEG measurements can be computed as In fact, as described at the beginning of this section, the CSP algorithm considers the whole set of electrodes that are placed during the recording session, which normally includes those covering the sensorimotor cortex. The number of resulting filters equals the number of EEG channels and, in most of the cases, at least two pairs (i.e., four filters) are selected. With this scenario the performance of the method is even better than the one observed in the previous example. Figure 4.5 shows the projections of two consecutive single trials for right hand and left hand motor imageries over the first two and the last two CSP filters (from up to down) computed with equation 4.14. In order to obtain the spatial filters coefficient matrix, all trials were considered over the whole set including the 26 electrodes used in the setup (see figure 2.7 in section 2.6.1), and by using a 4-second window starting 0.5 s after the task cue. All trials were previously and independently filtered between the frequency range from 8 to 30 Hz by using a 5th order Butterworth filter. We can observe that the signals variance is significantly increased during the task period for one type of motor imagery while it is decreased for the other one when using one pair of CSP filters and, the same behavior, but in the other way around, when using the other pair.

v i = log diag( W E i E i W) trace[ W E i E i W] (4.
Because of the manner in which the CSP algorithm achieves such a contrast by the simultaneous increasing/decreasing of the signals variance, this approach is only suitable to discriminate between two classes and, in order to solve problems involving many classes, multiclass extensions are required [START_REF] Gomez-Gil | Brain computer interfaces -a review[END_REF], [START_REF] Dornhege | Boosting bit rates in noninvasive eeg single-trial classifications by feature combination and multiclass paradigms[END_REF]. Along this chapter we will introduce a multiclass approach based on the joint approximate diagonalization (JAD) of the mean covariance matrices of multiple classes involved

Analytical Common Spatial Patterns (ACSP)

In contrast to the standard CSP algorithm, ACSP considers data in its complex form, which can provide more information that when considering only the real-value counterpart [START_REF] Falzon | The analytic common spatial patterns method for EEG-based BCI data[END_REF]. Without this, both methods aim at discriminating between two classes of data by finding a set of spatial filters that simultaneously maximize the variance of one class, while minimizing the variance of the other one. The ACSP method starts by converting the band-passed EEG signals s(t) into their complex-valued counterparts s(t), which can be achieved through the computation of the Hilbert transform in the following terms

s(t) = 1 π p.v. s(τ ) t -τ dτ, (4.16)
where p.v. represents the Cauchy principal value. Then, the analytic representation for a real-valued signal s(t) is given by

z(t) = s(t) + j s(t) (4.17)
The second step in the ACSP method involves the computation of the normalized complexvalued covariance matrix [START_REF] Nicolas-Alonso | Analytic common spatial pattern and adaptive classification for multiclass motor imagery-based BCI[END_REF]

Σ k = z k z * k trace(z k z * k ) (4.18)
where z k ∈ C n×T , k ∈ {1,2} denotes the analytic representation of a single EEG trial belonging to class k over N number of channels and T number of samples, and with * representing the hermitian transposition. Then, the means of the covariance matrices Σ 1 and Σ 2 are computed by averaging for each class the covariance matrices of the successive training trials. So that, similarly to the standard CSP approach, the complex-valued spatial filters W can be computed by solving the following generalized eigenvalue decomposition problem [START_REF] Blankertz | Optimizing spatial filters for robust EEG single-trial analysis [revealing tricks of the trade[END_REF] Σ

1 W = λΣ 2 W (4.19)
with λ denoting the diagonal matrix containing the eigenvalues of Σ 1 . The spatial filtered signal Y k can be obtained from the analytic representation of the EEG trials Z k as

Y k = WZ k . (4.20)
Finally, in order to generate the ACSP features equation 4.15 can be applied as discussed in the previous section. Figure 4.7 shows a toy example where the performance of the ACSP method can be appreciated in comparison to the standard version. Finally, in the last column, the obtained features after applying equation 4.15 over the projected trials by using the ACSP algorithm (G) and by using the CSP method (H) are presented. Note that separation in the case of the ACSP approach is enhanced in comparison to the standard algorithm.

Feature selection

During feature extraction many elements may be generated from several EEG channels, different time segments, multiple frequency bands, and even from many classes (e.g., under a one-versusone approach or under a one-versus-all approach, see sections 5.1.1 and 5.1.2), before being concatenated into a single feature vector. From this situation arises one of the major difficulties in BCI design, which consists in choosing relevant elements from the broad number of possible features [START_REF] Gomez-Gil | Brain computer interfaces -a review[END_REF]. High dimensional feature spaces are not desirable due to the curse of dimensionality, which implies that the required amount of training data to obtain reliable results increases exponentially with the dimensionality of the feature space [START_REF] Jain | Statistical pattern recognition: a review[END_REF]. Considering the long periods of time that the training process demands, in combination with the fatigue that it causes to the users, the availability of training data is usually small in BCI research, therefore, such requirements are difficult to meet and it becomes crucial to discard part of the data, specially because in many cases, some may not be even discriminative.

The process of choosing only the most discriminative elements is known as feature selection, which can be accomplished by examining all the possible subsets of features. However, the number of possibilities grows exponentially and, even for a moderate number of features, an exhaustive search may result impractical. Some more efficient algorithms can be applied with the aim of minimizing the number of features while maximizing the classification performance. A widely used approach consists in estimating the mutual inf ormation between features and the existing classes which, as explained in the next section, provides a measure of how descriptive are the involved terms to facilitate the classification decision.

Mutual information

The mutual information of two random variables is a measure of the mutual dependence between them. More precisely, it estimates the amount of information that can be obtained from one random variable through the other one. In this sense, mutual information is closely related to the concept of entropy, which defines the amount of information contained in a discrete random variable x ∈ R n in the following form

H(x) = - n i=1 p(x i )log b P (x i ), (4.21)
where p(x i ) represents the probability density function of x and b represents the base of the logarithm. A common choice is b = 2, for which equation 4.21 is referred as Shannon entropy, and whose units are expressed in bits. Considering a joint distribution p(x, y), if a value of x is already known, then the additional information needed to specify the corresponding value of y ∈ R m is given by

H(y|x) = - n i=1 p(x i ) m j=1 p(y j |x i )log b p(y j |x i ), (4.22)
which is called the conditional entropy of y given x. Now, mutual information is a measure of the inherent dependence expressed in the joint distribution of x and y under the assumption of independence, and it is defined as

I(x, y) = H(x) -H(y|x). (4.23)
Thus mutual information can be understood as the reduction in the uncertainty about x after being told the value of y (or vice versa). From a Bayesian point of view, we can consider p(x) as the prior distribution for x and p(x|y) as the posterior distribution after we have observed new data y. The mutual information therefore represents the reduction in uncertainty about x as a consequence of the new observation y [START_REF] Bishop | Pattern Recognition and Machine Learning[END_REF].

In the context of feature selection, a measure of mutual information reveals how much information the presence/absence of a term contributes to making the correct classification decision.

Filter Bank Common Spatial Pattern (FBCSP)

The effective use of the CSP algorithm discussed in section 4.1.3 requires the assignment of certain parameters, namely, the frequency range for the band-pass filtering of the EEG signals, the time interval relative to the task cue for extracting the trials, and the number of CSP filters pairs to be used [START_REF] Blankertz | Optimizing spatial filters for robust EEG single-trial analysis [revealing tricks of the trade[END_REF]. These parameters are generally fixed by a standard systematization of the method, so that the frequency range involves the alpha and beta bands (e.g., within the range from 8 to 30

Hz), the time segment is taken after 0.5 or 1 seconds after the stimuli, and 2 or 3 pairs of filters are considered. However, the performance of the CSP algorithm can be potentially enhanced if subject-specific parameters are taken into account [START_REF] Blankertz | The non-invasive berlin brain-computer interface: fast acquisition of effective performance in untrained subjects[END_REF]. In this sense, the Filter Bank Common Spatial Pattern (FBCSP) algorithm enhances the performance of the original CSP method by automatically selecting the most discriminative subject-specific frequency range for band-pass filtering of EEG measurements [START_REF] Ang | Filter bank common spatial pattern (FBCSP) in brain-computer interface[END_REF].

The FBCSP includes four progressive stages of signal processing and machine learning that are applied on the EEG data:

v b,i = log diag( W b E b,i E b,i W b ) trace[ W b E b,i E b,i W b ]
where v b,i ∈ R 2m ; W b represents the first m and the last m columns of W b ; diag(•) returns the diagonal elements of the square matrix; and trace[•] returns the sum of the diagonal elements in the square matrix. (Note that m = 2). The resulting feature vector for the ith trial is then arranged as follows

v i = [v 1,i , v 2,i , ..., v 9,i , ],
where v i ∈ R 1×(9 * 2) , i = 1, 2, ..., n; with n denoting the total number of EEG trials.

Feature selection to extract the most discriminative CSP features from

V = [v 1 , v 2 , ..., v n ].
To this end several algorithms can be used, however, in the present work we have decided to apply the Mutual Information-based Best Individual Feature (MIBIF) [START_REF] Ang | Filter bank common spatial pattern algorithm on BCI competition IV datasets 2a and 2b[END_REF], [START_REF] Ang | Filter bank common spatial pattern (FBCSP) in brain-computer interface[END_REF], which is described in the following section.

4. Classification to translate the generated and selected features. For this task there are also many possible choices, in the present work we will explore different solutions by using linear discriminant analysis and support vectors machines. Both techniques are discussed later along this chapter.

Mutual Information-based Best Individual Feature (MIBIF) algorithm

In brief, this algorithm consists in computing the mutual information of the considered features and sorting the results in a descending order, so that only k features are selected. For a 2-class problem, the MIBIF starts by computing the mutual information in the following way

H(f j ; w) = H(w) -H(w|f j ), (4.24)
with w = 1, 2 ∈ C representing the class labels, and the terms H(w) and H(w|f j ) consisting respectively of the entropy and conditional entropy, both defined as

H(w) = - 2 w=1 p(w)log 2 p(w), (4.25)
H(w|f j ) = - where the feature value of the ith trial from f j is represented by f j,i . Now, the probability p(w|f j,i ) can be computed by using the Bayes theorem as

p(w|f j,i ) = p(f j,i |w)P (w) p(f j,i ) , (4.28)
where p(w|f j,i ) is the conditional probability of class w given f j,i ; p(f j,i |w) is the conditional probability of f j,i given class w; P (w) is the prior probability of class w; and p(f j,i ) is given by p(f j,i ) = The conditional probability p(f j,i |w) can be estimated by using a Parzen window [START_REF] Parzen | On estimation of a probability density function and mode[END_REF] in the following terms

p(f j,i |w) = 1 n w k∈Iw φ(f j,i -f j,k , h), (4.30)
where n w is the number of data samples belonging to class w; I w is the set of indexes of the training data trials belonging to class w; f j,k is the feature value of the kth trial from f j and φ is a smoothing kernel function with a smoothing parameter h given by

φ(y, h) = 1 √ 2π exp - y 2 2h 2 , (4.31)
where

h opt = 4 3n 1 5 σ. (4.32)
By substituting equations 4.30 and 4.31 into 4.28 we can proceed to compute the mutual information associated to each feature, from which the k elements presenting the highest values will be selected. For the implementation of the FBCSP algorithm we have decided to use k = 4.

Note that, since the CSP features are paired, this number fluctuates between 4 components (if the selected elements come from the same 2 pairs) to 8 (if each one of the elements belongs to a different pair).

Common Spatial Pattern by Joint Approximate Diagonalization (CSP by JAD)

Classical CSP extensions for solving multiclass problems are carried out by performing a series of binary CSP (e.g., by computing CSP under a one-vs-one approach or under a one-vs all approach).

Another possible technique consists in computing the joint approximate diagonalization (JAD) over the mean covariance matrices of the involved classes, which is the same principle of the eigenvalue decomposition performed by the standard CSP method for two classes [START_REF] Dornhege | Boosting bit rates in noninvasive eeg single-trial classifications by feature combination and multiclass paradigms[END_REF].

Given EEG data recorded over N number of channels and T number of time samples from K different classes, the strategy of the CSP by JAD method consists in finding a transformation W ∈ R N ×N that diagonalizes the covariance matrices R x |c i , so that

WR x|c i R x|c i W = D c i , i = 1, ..., K, (4.33)
with D c i ∈ R N ×N diagonal matrices, and with W representing the spatial filters coefficient matrix. An effective approach to generate the transformation matrix W consists of applying the Fast Frobenius Diagonalization (FFDIAG) [START_REF] Ziehe | A fast algorithm for joint diagonalization with non-orthogonal transformations and its application to blind source separation[END_REF]. Now, each column vector of W represents a spatial filter, however, not all of them are relevant for the subsequent classification, so that only the most discriminative columns should be selected.

To this end, a preprocessing scheme based on information theoretic feature extraction (ITFE) can be implemented, from which only the L elements generating the highest mutual information values are considered as the optimal filters [START_REF] Grosse-Wentrup | Multi-class common spatial patterns and information theoretic feature extraction[END_REF]. The general framework of ITFE and its application to solve the aforementioned problem is presented in the following section.

Information theoretic feature extraction (ITFE)

In the context of non-invasive BCIs, the goal is to find a dimension-reduced representation of EEG data that maximizes the mutual information of class labels and extracted EEG components.

We thus have X = R T ×N , and choose X = R T for the dimension-reduced feature space. In this sense, EEG components are sequentially extracted. Now, let x ∈ X be a random variable (e.g., the recorded EEG data), from which c ∈ C (e.g., the class label to be predicted). Then, the general idea of ITFE, is to find a transformation f *: x → X that maps the original feature space X into a discrete set X while preserving information on the class labels c in x Now, note that the mutual information of c and x = w x can be written as (see section 4.3.1) we assume σ 2 x = 1, which does not imply loss of generality given that w can always be scaled to meet this assumption, we obtain an estimate of the mutual information of c and x in the following terms [START_REF] Grosse-Wentrup | Multi-class common spatial patterns and information theoretic feature extraction[END_REF] I(c, w j x) ≈ -

f * = argmax f ∈F {I(c, f (x))} (4.
I(c, w x) = H(w x) -H((w x)|c) = H(x) -H(x|c) = H(x) - K i=1 P (c i )H(x|c).
K i=1 P (c i )log w j R x|c i w j - 3 16
K i=1 P (c i )((w j R x|c i w j ) 2 -1) 2 (4.37)
The complete procedure of the multiclass ITFE is summarized in algorithm 1. Note that the output of the method corresponds with the L columns of w generating the highest mutual information. In [START_REF] Grosse-Wentrup | Multi-class common spatial patterns and information theoretic feature extraction[END_REF] this number was arbitrarily set to L = 8, from which reliable results where obtained. This parameter can also be chosen by using automatic methods [START_REF] Yang | Automatic selection of the number of spatial filters for motor-imagery BCI[END_REF]. However, this aspect is excluded from the scope of the present work. Thus, in order to evaluate the performance of the CSP by JAD algorithm in comparison to the efficiency attained by the other considered methods, we will use the recommended value of L = 8. Once the L filters have been selected, feature extraction can be accomplished in terms of equation 4.15, so that discriminative elements reflecting the contrast between the variance of the projected EEG trials are generated. (see section 2.6.1). EEG signals from the entire set of channels were band-pass filtered within the frequency range from 8 to 30 Hz, and features were extracted by using equation 4.15 over the selected w j filters. For visualization purposes, we show only the EEG features calculated from the 2 filters generating the highest mutual information, so that the resulting vectors are defined within a 2-dimensional space. Note that elements belonging to different conditions are well separated, however, by considering a higher number of filters their separation can be improved.

Algorithm 1 Multiclass Information Theoretic Feature Extraction (ITFE)

Input: Covariance matrices R x |c i , i = 1, ..., K. 1) Compute joint approximate diagonalization, so that wR x|c i R x|c i w = D c i , i = 1, ..., K (e.g., with the FFDiag-algorithm [START_REF] Ziehe | A fast algorithm for joint diagonalization with non-orthogonal transformations and its application to blind source separation[END_REF]). 2) For each column w j , j = 1, ..., N of w scale w j , so that w j R x w j = 1, and estimate mutual information according to

I(c, w j x) ≈ - K i=1 P (c i )log w j R x|c i w j - 3 16
K i=1 P (c i )((w j R x|c i w j ) 2 -1) 2
3) Select the L columns of w with highest mutual information. Output: Spatial filter coefficient matrix w ∈ R N ×L . Fig. 4.8: Classification space generated by the CSP by JAD algorithm for a 3-class problem including left hand, right hand, and both hands motor imageries. EEG signals from 26 channels were band-pass filtered within the frequency range from 8 to 30 Hz, and a 4-second window starting 0.5 s after the task cue initiation was considered for the feature extraction procedure. Note that only 2 filters were selected in order to generate a 2-dimensional representation and, even though classes are already well-separated, performance of the CSP by JAD algorithm is expected to improve with a higher number of selected filters.

Classification

The aim of classification is to assign a given input vector x to one of K discrete classes C k , where k = 1,..., K. In an ideal case, classes are considered to be disjoint, so that each input is assigned to one and only one class. The input space, or classification space, is then divided into different regions (e.g, each region associated to each one of the K classes), where boundaries are called decision boundaries or decision surf aces [START_REF] Bishop | Pattern Recognition and Machine Learning[END_REF].

Discriminant functions

In the context of classification, a discriminant function defines a decision boundary that separates the input space into K different regions associated to K classes denoted C k . The evaluation of an input vector x reflects its position with respect to the decision boundary, so that it can be associated to a specific region, and hence to the corresponding class. In this section we address attention to linear discriminants, namely those for which such a function defines an hyperplane in the following terms y(x) = w x + w 0 .

(4.38)

The evaluation of y(x) generates a real number, however, the goal is to predict discrete class labels. To this end, we can consider the sign generated after its evaluation, which determines the position of the input vector with respect to the decision boundary, given that any point generating a negative value will be found below the hyperplane, whereas a positive one will be situated above of it. With this in mind, we can consider a generalization of this model in which the linear function is transformed by using f (•) so that

y(x) = f (w x + w 0 ). (4.39)
In this context f (•) is known as the activation f unction, and can be defined as

f (y(x)) =    C 1 , if y(x) < 0 C 2 , if y(x) > 0 , (4.40)
where the two regions represented by the decision surface have been associated with classes C 1 and C 2 , so that any point located below the decision boundary belongs to C 1 and, any point located above it, belongs to C 2 . Note that the corresponding decision boundary is defined by the relation y(x) = 0, the assignment of any vector situated on the decision surface to a given class depends on the context of the classification problem.

Even if the function f (•) is nonlinear, the decision surfaces are linear functions of x, therefore, the class of models described by equation 4.40 are called generalized linear models [START_REF] Mccullagh | Generalized linear models[END_REF]. When classes can be completely separated by linear decision boundaries, the corresponding data set is said to be linearly separable.

To facilitate the discussion, we start with the case of two classes and then we consider the extension to K>2.

Two classes

As mentioned, the simplest model of a linear discriminant function is given by a linear function of the input vector so that

y(x) = w x + w 0 ,
where the parameters w and w 0 represent, respectively, the weight vector and bias of the resulting hyperplane. If the evaluation of y(x) < 0, the input vector x is assigned to class C 1 , otherwise, to class C 2 . The corresponding decision boundary is therefore defined by the relation y(x) = 0, which defines a (D -1)-dimensional hyperplane within the D-dimensional input space. Let us consider two points, x 1 and x 2 , both of them lying on the decision surface. Given that y(x 1 ) = y(x 2 ) = 0, it follows that w (x 1 -x 2 ) = 0 and hence the vector w is orthogonal to every vector located on the decision surface, therefore, w determines the orientation of the hyperplane. In the same way, if x is a point lying within the decision surface, then y(x) = 0, so that the normal distance from the origin to the decision surface corresponds to

w x w = - w 0 w . (4.41)
Based on this, we can see that w 0 defines the location of the decision boundary. These properties are represented in figure 4.9 for a 2-dimensional case.

Moreover, the evaluation of y(x) generates a signed measure of the perpendicular distance r of the point x from the decision boundary. To observe this, consider an arbitrary point x and let

x ⊥ be its orthogonal projection onto the decision surface, so that

x = x ⊥ + r w w . ( 4

Multiple classes

Let us consider the extension of linear discriminants for cases involving more than two classes (K > 2). Intuitively, we might think about building a K-class discriminant by combining a number of two-class discriminant functions. For instance, consider the use of K(K -1)/2 binary discriminant functions, each one trained in order to separate every pair of classes. This approach is known as a one-versus-one classifier (see section 5.1.1), and consists in assigning the class label C k to each input vector according to a majority vote amongst the discriminant functions.

However, this approach leads to ambiguous regions, as exhibited in figure 4.10-A, where the class assignment is compromised.

We can overcome these difficulties by considering a single K-class discriminant including K linear functions of the form

y k (x) = w k x + w k0 , (4.44)
so that a point x is assigned to class C k if y k (x) > y j (x), ∀j = k, or to class C j otherwise. The decision boundary between class C k and class C j is therefore given by y k (x) = y j (x), which can be defined as

(w k -w j ) x + (w k0 -w j0 ) = 0. (4.45)
This model corresponds to a (D -1)-dimensional hyperplane whose decision regions are always singly connected and convex (see fig. 4.10-B). Equation 4.45 has the same form as the one discussed in section 4.6.1 for finding the decision boundary in a two-class case, and hence analogous geometrical properties apply. In fact, this formalism is also pertinent for solving a binary problem, from which the decision boundary can be defined based on two discriminant functions y 1 (x) and y 2 (x). However, the formulation based on a single discriminant function y(x) suggests a simpler approach.

In the next two sections we will explore two approaches to learning the parameters of linear discriminant functions, namely, Linear Discriminant Analysis (LDA) and Support Vector Machines (SVMs).

Linear Discriminant Analysis (LDA)

The linear discriminant analysis for classification can be motivated by using a Bayesian approach, from which the probability of an input vector x belonging to class C k can be estimated in the following terms

p(C k |x) = p(x|C k )p(C k ) p(x) . (4.46)
Now, if we want to determine to which class x belongs to, it is possible to apply the Bayes decision rule, which implies that Decide

C 1 if p(C 1 |x) > p(C 2 |x); otherwise decide C 2 .
The generalization of such decision rule to multiple classes is straightforward and it can be rewritten as (-

g(x) = max k=1,...,K (p(x|C k )p(C k )) = max k=1,...,K (log(p(x|C k )p(C k ))).
1 2 (x -µ k ) Σ -1 (x -µ k ) + log(p(C k ) + δ)). (4.48)
The first term in equation 4.48 is the square form of the Mahalanobis distance (see section 4.1.2), which correspond with the separation from the input data to the corresponding class mean µ k , p(C k ) is the prior probability of class k, and δ is a constant involving all terms that do not depend on k, and hence can be removed, which leads to

g(x) = max k=1,...,K (x Σ -1 µ k - 1 2 µ k Σ -1 + log(p(C k ))). (4.49)
Now, we can define the linear discriminant function as

y k (x) = x Σ -1 µ k - 1 2 µ k Σ -1 + log(p(C k )), (4.50)
then,

g(x) = max k=1,...,K (y k (x)). (4.51)
For a two class problem, the decision boundary between k = 1 and k = 2 corresponds with

y 1 (x) = y 2 (x), (4.52)
or equivalently the following holds

x Σ -1 (µ 1 -µ 2 ) - 1 2 (µ 1 + µ 2 ) Σ -1 (µ 1 -µ 2 ) + log p(C k=1 ) p(C k=2 ) = 0. ( 4
   w = Σ -1 (µ 1 -µ 2 ) w 0 = 1 2 (µ 1 + µ 2 ) Σ -1 (µ 1 -µ 2 ) + log p(C k=1 ) p(C k=2 ) (4.54)
The parameters Σ, µ and p(k) are obtained from the training data set. If we consider a specific class covariance Σ k for each class, then the classifier is no longer linear and is referred to as quadratic discriminant analysis (QDA).

Support Vector Machines (SVM)

A support vector machine is a supervised classification algorithm based on the statistical learning theory developed by Vladimir Vapnik back in the 1970s [START_REF] Vapnik | Estimation of Dependences Based on Empirical Data: Springer Series in Statistics[END_REF]. Their functionality is based on a binary task, for which they aim at finding the hyperplane that maximizes the separation between the closest data belonging to different conditions (i.e., support vectors). This separation is called the margin, and it can be enhanced by mapping the corresponding features into a new space of higher dimension. To this end, they apply a kernel function that transforms the input space into a new one, where discrimination of classes from different conditions is reinforced.

To start with, let us consider the function defining an hyperplane given by equation 4.38

y(x) = w x + w 0 ,
Now, the goal of SVMs is to find the hyperplane that maximizes the margin between the selected support vectors. Such a margin is set so that the hyperplane lies equidistant from the support vectors belonging to different classes with a geometric distance of 1, as it is shown in figure 4.11. Then, the geometric distance for any input vector will satisfy

d(x, w, b) = |w x + b| w ≥ 1 w , (4.55)
so that |w x + b|= 1 for the support vectors and |w x + b|> 1 for the other training points (non-support vectors). This implies that the margin of the hyperplane ρ can be defined as twice the distance separating the support vectors from the hyperplane, so that

ρ = 2 w . (4.56)
In this way, the classification algorithm is able to find the optimal solution. Although, there is an infinite number of hyperplanes that solve a linearly separable problem, there is only one that maximizes the margin ρ in such a way that the parallel hyperplanes to the decision boundary that are drawn on the support vectors satisfy w x + b=1 and w x + b=-1 (see fig. to maximize the margin distance ρ, which can be achieved by minimizing w , or in order to facilitate the task without loss of generality, by minimizing 1 2 w 2 in the following terms

J(w) = min 1 2 w 2 , s.t., |w x + b|≥ 1, ∀i ∈ X , (4.57)
with X representing the training data set indexes. In most of the cases, problems are not completely linearly separable, so that it is not possible to meet this criterion. In order to accept a more flexible solution, the constraint in equation 4.57 is replaced by

w x + b ≥ 1 -ξ i , ∀i ∈ X , (4.58)
where the term ξ i represents the so-called slack variables, which allow support vectors to be in the margin if the evaluation of equation 4.58 complies with 0 ≤ ξ i ≤ 1. To avoid having a large number of misclassified data, the penalty term P is added to J(w), so that

J(w) = min 1 2 w 2 + P i∈X ξ i , s.t., y i (w x i + b) ≥ 1 -ξ i , ξ i ≥ 0, ∀i ∈ X (4.59)
The compromise between maximizing the margin and minimizing the amount of slack is

given by the constant P , also referred to as box constraint. There are two different approaches for solving this problem, on the one hand, J(w) can be directly addressed upon the primal optimization problem [START_REF] Chapelle | Training a support vector machine in the primal[END_REF]. And, on the other hand, it is possible to use Lagrange multipliers for solving the dual optimization problem, which actually is the most common technique [START_REF] Sch | Learning with Kernels : Support Vector Machines, Regularization, Optimization, and Beyond[END_REF]. By using the second alternative, the optimization problem is expressed in terms of α i in the following terms

D(α) = max i∈X α i - 1 2 i∈X j∈X y i y j α i α j x i x j , s.t., i∈X α i y i = 0 0 ≥ α i ≤ P, ∀i ∈ X . (4.60)
Once α has been found, the optimal parameters are given by

y(x) = i∈X α i y i x i + b. (4.61)
Only the support vectors x i , which define the limits of the margin and the separating hyperplane, meet the condition α i > 0.

As mentioned, when the elements from different classes are not linearly separable, it is possible to extend the SVM algorithm to a non-linear approach by implicitly mapping data into a higher-dimensional feature space, where classes are linearly separable. This solution is achieved by applying the so-called kernel trick [START_REF] Aizerman | Theoretical foundations of the potential function method in pattern recognition learning[END_REF], which uses a kernel function K = ϕ(x i )ϕ(x j ) to define the new feature space, with ϕ(•) representing the transformation over the input data. When such a kernel function is applied, equation 4.60 takes the following form

D(α) = max i∈X α i - 1 2 i∈X j∈X y i y j α i α j K(x i x j ) , s.t., i∈X α i y i = 0 0 ≤ α i ≤ P, ∀i ∈ X . (4.62)
The most common kernel functions are linear, polynomial and radial basis, which are defined as

-Linear: K(x i , x j ) = x i x j -Polynomial: K(x i , x j ) = (γx i x j + r) d , γ > 0 -Radial basis function: K(x i , x j ) = e (-γ x i -x j 2 ) , γ > 0
The original SVM model is not directly applicable to classification problems involving more than two classes. However, there are several approaches for its multiclass extension, such as the one-versus-one and one-versus-all techniques which will be both explained in chapter 5.

Distance-based classification

In order to explore different solutions to the CSP algorithm, we will present a small outlook to an alternative approach that is known as distance-based classification. The idea of these approaches is to assign the class label of an unknown element according to the shortest distance that separates it from the class measures defined during the training phase. This approach has drawn special attention due to the possibility of applying it in the context of Riemannian geometry, which will be elucidated in the following section.

Riemannian geometry

Square real symmetric positive-definite (SPD) matrices have the property of being always diagonalizable with strictly real positive eigenvalues, whose space represent a differentiable Riemannian manifold M. Using this manifold offers several advantages that can be easily applied in the context of signal processing [START_REF] Barachant | Multi-class brain computer interface classification by riemannian geometry[END_REF], [START_REF] Barachant | Covariance toolbox[END_REF]. Note that the covariance matrices are SPD, so that they belong to a manifold representing a Riemannian space from where their classification can be directly achieved.

Let denote by S(n) = {S ∈ M (n), S = S} the space of all n × n symmetric matrices in the space of square real matrices M (n); and P (n) = {P ∈ S(n), u Pu > 0, ∀ u ∈ R n } the set of all n × n SPD matrices. The derivatives at matrix P on the manifold lie in a vector space T p , which is the tangent space at that point (Figure 4.12). Now, let Γ Γ Γ(t): → P (N ) be any differentiable path from Γ(0) = P 1 to Γ(0) = P 2 . The length of Γ Γ Γ(t) is given by

L(Γ Γ Γ(t)) = 1 0 Γ Γ Γ Γ Γ Γ(t) , dt. (4.63)
The minimum length curve connecting two points in the manifold is called the geodesic, and it gives the Riemannian distance between the two points δ R (P 1 , P 2 ) = log(P -1

1 P 2 ) F = n i=1 log 2 λ i 1/2 , (4.64)
where λ i , i=1,...,n are the real eigenvalues of P -1 1 P 2 . For each point P ∈ P (n), the set of tangent vectors at P defines a tangent space. Each tangent vector S i can be seen as the derivative at t = 0 of the geodesic Γ Γ Γ i (t) between P and the exponential mapping P i = Exp P (S i ), given by Exp p (S i ) = P i = P

1 2 exp P -1 2 S i P -1 2 P 1 2 . (4.65)
The inverse mapping is given by the logarithmic mapping defined as Log p (P i ) = S i = P Fig. 4.12: Tangent space at point P, S i a tangent vector at P and Γ Γ Γ(t), the geodesic between P and P i defined in the Riemannian manifold.

By using the Riemannian geodesic distance, the Riemannian mean of I ≥ 1 SPD matrices is given by

G(P 1 , ..., P I ) = argmin P∈P (n) I i=1 δ 2 R (P, P i ) (4.67)
This mean is also referred to as the geometric mean. For a manifold of non-positive sectional curvature like P (n), such local minimum exists and is unique [START_REF] Karcher | Riemannian center of mass and mollifier smoothing[END_REF]. However, there is no closedform expression to compute the mean and optimization algorithms must be applied. An efficient iterative algorithm to compute the Riemannian mean of I SPD matrices is given in [START_REF] Fletcher | Principal geodesic analysis on symmetric spaces: statistics of diffusion tensors[END_REF]. the right hand. It can be observed that separation of classes is well achieved, though there are several outliers that difficult the task. In particular, it is possible to observe for both classes a cluster of data that is shifted with respect to the rest. These elements correspond to first of the four runs that were conducted during the recording session. Presumably, there were changes in the electrodes impedance as the session advanced, and/or in the subject mental state.

CSP and Riemannian geometry

As mentioned, the MDRM algorithm uses the band-pass filtered EEG trials covariance matrices as descriptors, which are subsequently considered within the Riemannian manifold, and from which it is possible to extract the spatial information without using spatial filtering. Nevertheless, the spatial filters generated over two classes with the CSP algorithm extract information that is relative to both conditions, since they attempt to enhance uncorrelated information. Since the covariance matrices of the CSP projected trials preserve the differences attained through the spatial filtering, we have considered that it is advantageous to apply the CSP algorithm before mapping them into the Riemannian space in order to generate more discriminative features. In other words, the spatial information that is obtained by only considering the band-pass filtered EEG covariance matrices within the Riemannian space, as it is done by the MDRM method, offers independent information of each class. Whereas the CSP projected covariance matrices provide uncorrelated information from both classes.

5

Multiclass and multilabel approaches Along this chapter we will introduce the grouping methods for classification that have been developed in order to solve the multilabel classification problem that arises from the consideration of the combined motor imageries included in the 8-class database, in fact, two of these methods are also pertinent for the motor tasks comprised in the 4-class database. To this end, before going further, it is important to remember that, on the one hand, it is necessary to extract the features presenting the most discriminative information for contrasting two or more conditions, so that redundancies are eliminated and specific patterns enhanced. And on the other hand, once such features have been extracted, the predicted labels can be inferred by applying a classification algorithm that separates the features according to their position within the classification space (see sections 4.1 and 4.6 for further details).

In the first section we will introduce the concept of multiclass approach in the context of classification by describing three classical methods, namely, one-versus-one, one-versus-all, and hierarchical algorithms. These methods will be subsequently used as benchmarks to compare the performance of the algorithms contributed by the present work. Such contributions will be presented in the second section as multilabel approaches for sensorimotor rhythms discrimination.

The basis of each one of these new methods lies in the results that were presented along chapters 2 and 3, where we have suggested that EEG data recorded during a combined motor imagery can be modelled as the superposition of the activity that is independently generated by the sources associated to the body parts engaged in the task. Upon this generalization we have proposed a convenient manner to carry out data partition, so that the entire problem is reduced into a series of binary classification tasks. This strategy has been applied in three new different methods that vary in the manner of how the features extracted from each binary problem are classified, namely, One-step Multilabel (OsM) approach, Hierarchical Multilabel (HM) approach, and One-step Hiearchical Multilabel (OsHM) approach.

An important aspect of this way of grouping classes is that it optimizes the use of the CSP method described in section 4.1.3, which as mentioned consists of a convenient approach to differentiate sensorimotor rhythms that has the drawback of being suitable for separating only two classes. In this regard, in contrast to the one-versus-one and one-versus-all approaches, the proposed methods generate features throughout the same training process for all classes, so that the distances separating elements from the decision boundaries within the classification space present an appropriate scale that favors their subsequent classification.

Multiclass approaches

As mentioned in section 4.6, the problem of classifying an element in order to assign the corresponding class label consists in determining on which side of the decision surface it is located.

For 2-class problems this is a relative simple task that can be solved by assessing the sign returned by the evaluation of the function that defines the decision boundary over the input element.

Nevertheless, when the classification task involves more than 2 classes, this problem results more complicated. In fact, in this case it is necessary to generate multiple decision surfaces for separating all the considered classes, and to assign the corresponding labels according to the relative position of the input vectors to all boundaries within the classification space.

There are many possible ways for separating classes in order to set down the decision boundaries.

Each approach varies in the number of functions that, according to the amount of classes, are required to divide the classification space into the associated regions. In this regard, the best approach is the one that satisfactorily adapts itself to the distribution of data while demanding a reasonable number of computations. In the following sections we will describe three of the most common methods, namely, one-versus-one, one-versus-all, and hierarchical approaches.

MULTICLASS APPROACHES

One-versus-one approach

As the name suggests, the one-versus-one approach consists in generating a model to independently separate each one of the groups. To this end, it is necessary to train K(K -1)/2 different models, where K represents the number of classes, on all possible pairs of groups and then, in the case of classification, to assign the labels of the testing set according to the class having the highest number of "votes" [START_REF] Bishop | Pattern Recognition and Machine Learning[END_REF]. This technique has the drawback of leading to ambiguities in the resulting classification. However, this issue does not represent its biggest inconvenience; since it can be tackled by generating convex regions, as described in section 4.6.1.1. In fact, the main problem is that, for a large K, such solution requires a considerable amount of time for setting all the required models. Similarly, for the evaluation of testing data, significantly computation is required.

One-versus-all approach

Another common approach consists in constructing K separate models, where K represents the number of classes, in which the k th model is generated by grouping within one group all elements belonging to class C k and, in the second class, all elements from the remaining K -1 classes. This is known as the one-versus-all approach. However, using the decisions of the individual classifiers can generate inconsistent results; in which an input is assigned to multiple classes simultaneously. This problem is sometimes addressed by choosing the label associated to the largest distance separating the input data from the decision boundary among all models [START_REF] Vapnik | Statistical learning theory[END_REF]. Unfortunately, this heuristic approach presents the inconvenience that the different models were trained on different tasks, and there is no guarantee that the distances generated by all models will have appropriate scales. Another problem that has the one-versus-all approach is that the training

Hierarchical approach

Hierarchical methods allow solving multiclass problems by using a tree of binary classifiers, whose root discriminates between two groups; each one containing a half of the classes. Each succeeding node includes again only one half of the classes from the selected group, and the process is recursively repeated until each node contains a single class, from which the final decision can be inferred [START_REF] Lindig-León | Hierarchical classification using riemannian geometry for motor imagery based bci systems[END_REF].

The hierarchical approach can be a very convenient method, since even when several models are required to be generated during the training stage, only a few of them are used during the validation stage. However, there are two important drawbacks to be considered; first, if there is a single erroneous prediction throughout the classification chain, the final result will be wrong. And second, given that from the second stage all models are selected according to the previous decision, all the subsequent stages have to be applied consecutively, which can lead to long computations.

Multilabel approaches

Conventional single-label classification involves learning from a set of examples that are associated with a single label l from a set of disjoint labels L, |L|> 1. As mentioned, if |L|= 2, the learning task is called binary classification, while if |L|> 2, then it is called multiclass classification. In multilabel problems the examples are associated with a set of labels Y ⊆ L [START_REF] Tsoumakas | Multilabel classification: an overview[END_REF]. In the past, multilabel classification was mainly motivated by medical diagnosis tasks. For example, it can be the case that a patient suffers from both hypertension and atrial fibrillation at the same time.

Nowadays multilabel classification methods are increasingly required by modern applications, such as protein function classification, and among others, music categorization [START_REF] Tsoumakas | Mulan: a java library for multi-label learning[END_REF]. Essentially, there are two ways for addressing multilabel classification; problem transf ormation methods, which as its name indicates consists in transforming the multilabel task either into one or more single-label classifications problems, and algorithm adaptation methods; which extend specific learning algorithms in order to handle multilabel data directly.

In the case of motor imagery-based BCI systems, the multilabel approach arises when combined movements are included in the paradigm; which as mentioned, allows to drastically increase the number of afforded commands. Now, if we consider the assumption that was justified by the analysis presented in sections 2.6.3 and 3.1.3 for the 4-class and 8-class databases respectively, it turns out that there is a particular way for grouping classes that allows to define a problem transformation method based on the manner of how the activity is distributed within the different sources. In simple terms, this activity appears as an ERD modulation over the regions whose associated limb is engaged in the motor imagery. In the case of combined imaginations, such modulation can be characterized as the superposition of the activity that is independently generated by each source during simple tasks [START_REF] Lindig-León | Comparison of sensorimotor rhythms in EEG signals during simple and combined motor imageries over the contra and ipsilateral hemispheres[END_REF], [START_REF] Lindig-León | A multilabel classification method for detection of combined motor imageries[END_REF]. With this in mind we have proposed to group EEG data separately for each activity source; so that, all motor imageries involving the use of the associated body part, are gathered together into one class (hereinafter to be referred to as C ERD ), and the remaining conditions that do not include it in another class (hereinafter to be referred to as C IDLE ). In this way, it is possible to reduce the entire problem into a series of binary tasks; which simplifies the complexity of the classification and boost feature extraction by allowing the use of the CSP algorithm, which is only suitable to discriminate between two classes. By using this type of partition we have implemented three new different methods that vary in the manner of how the features extracted from each binary problem are classified, namely, One-step Multilabel (OsM) approach, Hierarchical Multilabel (HM) approach, and One-step Hiearchical Multilabel (OsHM) approach.

Before introducing the aforementioned problems we present, in order to simplify their description, a convention for labeling each one of the different tasks according to the presence/absence of motor imagery-related activity at each one of the main sources. If the sources are active, a value of "1" is used to represent them; otherwise, the value corresponds to "0". The number assigned to each source is used to generate a P -bit word, with P representing the number of body parts considered by the paradigm (2 for the 4-class database, and 3 for the 8-class one), where the less significant bit corresponds to the area around C4 (right side) and the most significant bit to the area around C3 (left side). Hereafter, in order to describe the proposed methods, we will refer to each class by considering the decimal equivalence of the aforementioned convention (see figure 5.4 for a complete association between motor tasks and labels). Fig. 5.4: Labeling convention. The proposed convention assigns a value of "1" to represent the sources that are active in terms of motor imagery-related activity; otherwise, the value corresponds to "0". The number assigned independently to each source is used to generate a 3-bit word, whose decimal equivalence will be used when referring to each one of the classes considered in the paradigm. Note that the same convention is applied to the 4-class database, for which the sources associated to the left hand and right hand are used to generate a 2-bit word, where the less significant bit is associated to the use of the left hand, and the most significant bit to the use of the right hand.

One-step Multilabel (OsM) approach

The One-step Multilabel (OsM) method considers, in the case of the 4-class problem formulated by the 4-class database, two feature extraction modules to independently determine whether the sources linked to the left hand and right hand are engaged during a motor task (see Figure 5.5). Similarly, in the case of the 8-class problem given by the 8-class database, besides these two models there is a third one to analyse the activity generated by the source associated to the feet (see Figure 5.6). Each module is defined as a binary problem where one class corresponds to the group of motor tasks showing ERD modulation at the corresponding source (i.e., motor imageries involving the use of the limb associated to that source), and the other class to the group of motor tasks that do not generate such activity and present similar patterns to the idle state (i.e., motor imageries that do not involve the use of the limb associated to that source).

The whole procedure involves two feature extraction models in the case of the 4-class problem (compared to the six and four models that are required to solve the same problem by using the one-versus-one and one-versus-all approaches respectively), and three feature extraction models for the 8-class one (compared to the 28 and 8 models that are required to solve the same problem by using the one-versus-one and one-versus-all approaches respectively). In both cases, outputs are concatenated to form a single vector, so that only one classification model is required.

When applied in combination with the CSP method, and by selecting 2 pairs of filters (4 in total), for each one of the two binary problems devoted to solve the 4-class task, the OsM approach generates a total of 8 features per trial. Whereas, in the case of the 8-class problem, the resulting vectors comprise 12 features per trial.

Hierarchical Multilabel (HM) approach

The second method consists of a hierarchical decision process that comprises three feature extraction models distributed in two consecutive stages for solving the 4-class problem given by the 4-class database (see Figure 5.7), and seven feature extraction models placed within three consecutive stages for solving the 8-class problem formulated by the 8-class database (see Figure 5.8). Such stages constitute a tree of binary classifiers whose root discriminates between two groups containing each one a half of the classes. Each succeeding node includes only the half of the classes that was selected in the previous stage, and the process is recursively repeated until the last node contains only two classes, from which the final decision can be inferred. Data partition is carried out in a similar manner as the OsM method; so that the presence/absence of motor imagery-related activity associated to one of the hands is analyzed during the first stage in order to determine whether or not it is engaged in the motor task. To this end, two groups are conformed by gathering together the motor tasks involving the selected hand in the ERD-class, and by putting together the motor tasks that do no involve it in the IDLE-class. Note that the hand selection to start the classification process is arbitrary, and that it is possible to start either with the left or the right hand. Once the decision is inferred, half of the classes are discarded.

After this decision there are only two classes left in the case of the 4-class problem, whereas for the 8-class task there are four classes remaining (i.e., in both cases either the group of motor tasks involving the selected hand or the group that does not include it). Based on this criterion, two feature extraction and classification models are generated in the second stage; from which, according to the previous prediction, only one of them is applied during the validation phase. In this part of the process the aim is to identify whether or not the source linked to the hand that has not been analyzed presents ERD activity, given that the hand selected after the first stage has been already recognized as active or inactive in terms of motor imagery-related activity. In the case of the 4-class problem the second stage corresponds with the last one, given that each one of the two models corresponds to a binary problem including only two classes, from where the final decision can be inferred. One of these models correspond to the one formulated by the engagement of the hand analyzed in the previous stage, in which case the goal is to determine whether or not the other hand is active, so that the predicted label would result in the use of both hands if the second hand is also found to be engaged, or only the use of the hand that was already detected in the previous step if the second hand was found to be inactive. In the same way, the second model is generated under the assumption that the hand analyzed during the first stage was labelled as inactive, so that the goal is to determine whether the other hand is also inactive, so that the prediction would correspond to the rest state, or whether the second hand is active, in which case the predicted label would be assigned according to the manner in which data partition was proposed in the formulation of the problem. In the case of the 8-class problem, the second stage corresponds to an intermediate process that aims at determining the state of the second hand before to start analyzing the activity generated by the feet. Thus, if the inference of the first step has considered the ERD-class, the corresponding model generates a new ERD-class by gathering together the two motor tasks that involve the use of both hands; and in a new IDLE-class the two motor tasks that involve the use of the hand selected in the previous model but that do not involve the use of the hand under the current analysis. Contrarily, if the inference of the previous step has considered the IDLE-class, the corresponding model creates a new ERD-class by putting together the two motor tasks that involve the use of the hand under analysis but that do not include the use of the hand selected by the previous step; and a new IDLE-class by gathering together the motor tasks that do not involve the use of any hand. After the second stage there are only 2 classes remaining, from which the final decision can be inferred in the third and last step. To this end, there are four different instances that are determined by the combination of decisions made by the preceding stages, each one of them defining a binary problem that considers the activity generated by the source linked to the feet. Thus, in the case that any hand was found to generate ERD activity, the classification task consist of distinguishing between rest and feet. Contrarily, if only the left hand was found to be involved in the motor task, the classification aims at discriminating between left hand and left hand in combination with feet. In the same way, if only the right hand has generated ERD activity, the task consists of classifying right hand versus right hand in combination with feet and, finally, if both hands were found to be involved during the motor task, the goal is to distinguish between both hands and both hands in combination with feet.

One-step Hierarchichal Multilabel (OsHM) approach

The third method is a combination of the first two approaches; it implements the same binary problems related to the sources associated to the hands, as it was done with the OsM method, and incorporates the four instances associated with the source linked to the feet that were formulated by the HM algorithm. The whole procedure involves 6 feature extraction models (compared to the 28 and 8 models that are required to solve the same problem by using the one-versus-one and one-versus-all approaches respectively), whose outputs are concatenated to form a single vector, so that only one classification model is required. In this way, by considering all the possibilities involving the use of the feet, more discriminative features are generated in comparison with those obtained by the OsM approach. And, given that all features are considered together, it is not necessary to apply different classification stages as it was done with the HM approach, which speeds up the process and helps to compensate outliers.

When applied in combination with the CSP method, and by selecting 2 pairs of filters (4 in total) for each one of the 6 binary problems, the OsHM approach generates a total of 24 features per trial. Figure 5.9 shows the architecture of the OsHM approach for solving the 8-class problem formulated by the 8-class database. Note that in this case there is no adaptation for the 4-class database since, given that it only includes the use of the hands, such approach is equivalent to the OsM method. mented; one model for the first stage, which is trained to discriminate between the two motor tasks involving the use of the selected hand (denoted by C ERD ) from the two motor tasks that do not involve it (denoted by C IDLE ), and two binary models in the second stage which are trained to discriminate between the motor tasks involving the hand that has not been analyzed (denoted by C ERD ) from those that do not involve it (denoted by C IDLE ), given that the hand selected for the first stage has been already recognized as active or inactive in terms of motor imagery-related activity. During the evaluation phase only two models are needed, which are selected according to the prediction of the previous stages. Note that this scheme illustrates the case where data partition is carried out by analysing first the engagement of the left hand. one model for the first stage, which is trained to discriminate between the four motor tasks involving the use of the selected hand (denoted by C ERD ) from the four motor tasks that do not involve it (denoted by C IDLE ), two binary models in the second stage which are trained to discriminate between the motor tasks involving the hand that has not been analyzed (C ERD ) from those that do not involve it (C IDLE ), given that the hand selected for the first stage has been already recognized as active or inactive in terms of motor imagery-related activity, and finally, four binary models in the third stage to determine whether or not the feet are involved in the motor task, given that the hands have been already identified as active or inactive. During the evaluation phase only three models are needed, which are selected according to the prediction of the previous stages. Note that this scheme illustrates the case where data partition is carried out by analysing first the engagement of the left hand. the recording sessions. Furthermore, we also enumerate the characteristics of the classification functions, as well as the complexity parameters that were estimated for the optimization of the different kernel functions considered for the SVM algorithm applied on the 4-class database, namely, linear, polynomial, and radial basis function. Subsequently, we describe the cross-validation procedures that were implemented on both databases for assessing all methods performance as well as for estimating the aforementioned complexity SVM parameters.

After indicating all experimental specifications we present, in terms of the mean accuracy among subjects, the offline results that were generated from applying each one of the feature extraction methods described along chapters 4 and 5, the detailed results obtained for each subjects, are presented in appendix C.

Experimental parameters

In order to carry out an EEG motor imagery-based session, it is necessary to set a series of parameters that depend on the characteristics of the experiment. For instance, it is necessary to choose the amount and position of the recording electrodes; which as discussed in section 2.2, should be located in such a way that the primary motor cortex is covered, so that the elicited activity generated by the considered activity sources is well conducted. Moreover, if the paradigm is asynchronously implemented, it is required to establish the time window in relation with the task cue for classification, which should include the period where the signal patterns associated to the motor task are prevalent (i.e., not before it starts reflecting the associated activity, nor after it begins to stabilize). Additionally, the frequency range for the band-pass filtering must be selected in a way that all the components containing the patterns of interest are considered, whereas those related to the background activity are suppressed.

The tuning of these parameters represents a research concern, since they have direct influence on the system performance and, given that there is a huge inter-subject variability with respect to the characteristics of the brain signals [START_REF] Blankertz | The non-invasive berlin brain-computer interface: fast acquisition of effective performance in untrained subjects[END_REF], many efforts have been made in order to develop automated algorithms for their optimization [START_REF] Ang | Filter bank common spatial pattern algorithm on BCI competition IV datasets 2a and 2b[END_REF], [START_REF] Suk | Subject and class specific frequency bands selection for multiclass motor imagery classification[END_REF]. However, these aspects are out of the scope of the present work and we have decided to use the same parameters along all subjects and for both databases, which have been selected according to the observations described in chapters 2 and 3. All the selected values for the mentioned parameters that were used during the experimental setup are listed below -Recording electrodes -Fp1; Fp z ; Fp2; F z ; FC5; FC3; FC1; FC z ; FC2; FC4; FC6; C5; C3; C1; C z ; C2; C4; C6; CP5; CP3;CP1; CP z ; CP2; CP4; CP6, and P z , re-referenced with respect to the common average reference across all channels and located over the international 10-10 system positions to cover the primary sensorimotor cortex (see fig. 2.7) -Filtering

EXPERIMENTAL PARAMETERS

-Each trial was independently filtered to better emulate an online setup by using a 5th order Butterworth filter within the frequency range from 8 to 30 Hz.

-Time window size -A 4-second window starting 0.5 s after the task cue was considered for classification.

Classification algorithms

As mentioned in chapter 4, besides feature extraction, it is necessary to implement a classification algorithm that according to the position of the feature vectors situated within the classification space assigns their corresponding class label, so that input data are translated into the commands that allow controlling the interface. On this regard, we have decided to evaluate different classification algorithms to assess the performance of the contributed approaches presented in chapter 5 in comparison with classical methods. To this end, for the 4-class database,we present all results obtained by Linear Discriminant Analysis and Support Vectors Machine classification, both of which were described in sections 4.6.1.3 and 4.6.1.4, respectively. As mentioned, in the case of the 8-class database we restrict the evaluation only to LDA classification, since the SVM algorithm generates poorer performance.

For Linear Discriminant Analysis classification we consider two cases; the premise of equivalent covariance matrices of the involved classes (LDA), and the assumption of their specific covariance matrices (QDA). All feature extraction methods (excepting those using a Riemannian distance-based scheme) were evaluated by using the Linear Discriminant Analysis classification under a one-versus-one and one-versus-all approach (see sections 5.1.1 and 5.1.2, respectively).

On the other hand, for the Support Vectors Machine classification applied on the 4-class database, we consider three different kernel functions, namely, linear, polynomial, and radial basis function. Each one of these models requires a set of complexity parameters. For instance, it is necessary to estimate the box constraint, which controls the trade-off between under-and overfitting. For the linear kernel this is the only parameter that we have considered for the model optimization, whereas for the polynomial and radial basis functions, in addition to the box constraint, we address, respectively, the optimization of the polynomial order and the gamma parameter (i.e., the inverse of the radius of influence of samples selected by the model as support vectors). All parameters have been estimated by using the cross-validation procedure that is described in the following section.

For the 4-class database, all feature extraction methods (excepting the Riemannian distancebased methods, and those using the ACSP algorithm, since it generates complex features requiring a special strategy for applying the SVM approach, which is out of the scope of the present work), were evaluated in combination with the SVM algorithm under a one-versus-one and one-versus-all approach by using the three aforementioned kernel functions.

Cross-validation

Due to the overfitting problem, performance obtained on the training set is not a good indicator of predictive performance on unseen data. When the amount of data is sufficiently large, then one approach simply consists in using some of the available data to train a range of models, or a given model with a range of values for its complexity parameters, and then to compare them on independent data, sometimes called validation set, and select the one generating the best results [START_REF] Bishop | Pattern Recognition and Machine Learning[END_REF]. In order to generate robust models, it is important to use as much of the available data as possible for training. However, in many applications the available amount of data is limited, and hence, if the validation set is small, it will produce a noisy estimate of predictive performance.

One possible solution for this inconvenience is to use cross-validation; which consists of using a proportion (s -1)/s of the available data for training while making use of all of the remaining data to asses performance. When data is particularly limited, it may be appropriate to consider the case S = N , where N is the total number of data points, which gives the leave-one-out technique.

An important drawback of cross-validation is that the number of training runs that must be performed is increased by a factor of S, and this can be problematic for models in which the training is itself computationally expensive. A further problem is that there might exist multiple complexity parameters for a single model. In these cases it becomes necessary to apply a nested cross-validation, which consists of carrying out another cross-validation process over the (s -1)/s training data to evaluate the impact of different values for the involved parameters. Exploring combinations of settings for such parameters could, in the worst case, require a number of training runs that is exponential in the number of parameters.

Cross-validation on the 4-class database

As mentioned in section 2.6.1, the 4-class database contains 4 runs of 20 trials per class for each subject, which totals 80 trials per class. Given that we have observed some intra-subject variability over the distribution of data belonging to different runs (see figure 4.13 as an example), which might be presumably due to mental state changes and/or fluctuations in the electrodes impedance throughout the recording session, we have decided to implement a 4-fold cross-validation using the trials of an entire run (20 trials per class) as the testing set, and the trials of the 3 remaining runs (60 trials per class) as the training set. Such procedure is repeated so that all runs are used to train and test the corresponding models. In this way, we can better emulate the variability that it is expected to be found across different sessions, and to generate results whose performance is independent from a (un)fortunate (un)favorable data partition. Figure 6.1 shows the 4-fold cross validation process that has been implemented to evaluate the classification performance obtained by all the discussed feature extraction algorithms applied in combination with the LDA and QDA methods, as well as for the distance-based approaches. Note that none of these classifiers requires any complexity parameter, and therefore the cross-validation process pursues only to evaluate their performance without needing a nested procedure. Once all possible choices have been implemented, the performance scores for the 4 folds are averaged to obtain the performance estimate for each model. In contrast to the linear discriminant analysis algorithm, classification based on SVM requires different complexity parameters that depend on the selected kernel function. In order to determine the optimal values, one possible solution involves performing a nested cross-validation; which consists in partitioning again the training set generated at each iteration to perform a new cross-validation process for assessing the model performance over a range of values. Figure 6.2 shows the nested procedure that has been implemented on the 4-class database to evaluate the performance of the discussed methods applied in combination with SVM classification. As in the previous case, we have decided to use a 4-fold cross-validation considering, at each iteration, all trials belonging to one run as the testing set, and the 3 remaining ones as the training set.

In addition, towards determining the optimal complexity parameters, another cross-validation procedure is performed over the training set generated at each iteration, for which a range of different values is evaluated over the classification model in order to select the optimal ones. As discussed in section 4.6.1.4, regardless the selected kernel function, the SVM algorithm requires the estimation of the box constraint, which controls the trade-off between under-and overfitting. For the linear kernel this is the only parameter that we have considered for the model optimization, whereas for the polynomial and radial basis functions, in addition to the box constraint, we address, respectively, the optimization of the polynomial order and the gamma parameter (i.e., the inverse of the radius of influence of samples selected by the model as support vectors). For both the box constraint and the gamma parameter we have evaluated a series of values covering the range from 10 -2 to 10 2 [103], [START_REF] Cristianini | An Introduction to Support Vector Machines and other kernel-base learning methods[END_REF], whereas for the polynomial order we have restricted the analysis, in order to avoid overfitting, within a range from 2 to 5. As illustrated, the nested cross-validation involves also 4 iterations, for each one of them, the different values within the proposed range are evaluated. When more than one parameter need to be estimated, then the corresponding range is evaluated for each combination of values. Finally, the optimal parameters are selected from the highest performance score to implement, by using the original training set, the optimal classification model which, in turn, is evaluated on the original testing set.

Cross-validation on the 8-class database

In contrast to the 4-class database, the 8-class one contains rather less trials per class due to the increase in time that it takes to accomplish such a high number of different motor tasks; which causes users fatigue and their loss of concentration. As mentioned in section 3.1, the 8-class database contains 3 runs of 10 trials per class, which sums 30 trials per class for a total of 240 trials recorded in the whole session for each subject. Because of this limited amount of data, we have decided to apply a different cross-validation approach that the one described in the previous section, which consists of a 30-fold procedure that considers for each iteration one trial per class (8 trials in total) as the testing set, and the remaining 29 trials per class (232 trials in total) as the training set. This approach allows using more trials for training than the ones that would be available by following the same scheme proposed for the 4-class database, where the trials belonging to one run were considered as the testing set, and the remaining ones as the training set.

In fact, we could have implemented the leave-one-out technique discussed at the beginning of this section to generate more robust models, however, since it drastically increases the computational expenses, we have decided to use the mentioned procedure illustrated in figure 6.3. As it can be observed, throughout each fold the testing set was constituted by one trial consecutively selected from each class (8 trials in total), and the remaining 29 trials per class (232 trials in total) as the training set. After all folds, all possible choices are implemented and the performance scores for the 30 folds are averaged.

Results: summary

Finally, in the following section, we summarize all the reported results on tables 6.1 and 6.2 for the 4-class and the 8-class databases respectively, where performance is shown in terms of the grand average accuracy across subjects. Note that in cases where some of the classification functions were not applied, either because the features were complex or because the feature method corresponded to a distance-based algorithm, this has been indicated as "NA" (non applicable). In the same way, we have used the terms "OVO" and "OVA" for referring to the one-versus-one and one-versus-all approaches respectively. by the OsM, HM, and OsHM approaches in comparison to the classical one-versus-one and one-versus-all algorithms for each subject. In all cases features were extracted by using the CSP method in combination with LDA classification. In addition, Figure 6.5 shows the mean accuracy with the corresponding standard error of the mean (SEM) among all subjects for each one of the methods. And, in order to observe the performance of each method over the best subjects, Figure 6.6 present the mean accuracy together with the SEM among subject, 3, 10 and 11. Note that the For the purpose of a fair comparison between the three methods, we will use the aforementioned data that were taken from the training set generated during the same partition of the 30-fold cross validation for subject 3. The first CSP patterns computed over these data for each one of the 7 binary problems defined by the MH approach are shown, for the case in which the process starts by analysing the activity associated with the use of the left hand, in Figure 12 and, for the case in which the process starts by analysing the activity linked to the use of the right hand, are presented in Figure 13. As observed in the previous method, there are strong weights distributed around the associated activity sources over the motor cortex confirming that this partition is also plausible. Figure 14 shows, step by step, the first 2 features generated for each class across the 3 stages formulated by the HM approach. It can be observed that for all problems conditions are well separated and, contrarily to the first method, classification problems involving the feet show better discrimination;

which is due to the fact that, for these cases, classes were conformed by using only 2 conditions, and hence more specific patterns are considered. Even though features are well separated after applying the HM approach, classification results do not outperform the OsM algorithm. This is possible related to the fact that the HM approach involves three classification stages and, if there is a misclassified element at any point, the predicted label is misclassified. In order to overcome this problem we have developed the OsHM approach; which, as mentioned, consists of a combination of the OsM and HM methods that allows generating discriminative features while favoring a straightforward classification. In fact, the OsHM approach yielded a better performance than the HM method for almost all subjects, however, it does not outperform the OsM either. In this sense, even when the classes formulated by the four instances involving the use/not-use of the feet according to the last stage of the HM approach are well separated (see Figure 15), data whose classes were not considered during the models training cause some ambiguities. For instance, in the case of a left hand motor imagery when applying the model both hands-versus-both-hands & feet, decision is not easy since the left hand is engaged in both tasks. These ambiguities occur within each one of the four models, which presumably hinders an effective classification.

The distribution of the first 2 CSP features generated by the models defined by OsHM approach for evaluating the activity associated to the use of the feet is shown in Figure 15. Note that only the features associated to the feet are presented, since the ones for the left and right hands correspond with the same ones that were already mentioned for the OsM method (see Figure 11). In the same way, the corresponding CSP patterns are, in the case of the hands, those that were presented for the OsM method (see Figure 10-A-C) and, in the case of the feet, the four instances that were described for the HM approach (see Figure 13-D-E-F-G or Figure 12-D-E-F-G, instinctively).

Additionally, in order to compare the specificity of the one-versus-one and the OsM methods in combination with the ACSP algorithm and LDA classification (which had led to the best performance), a graphical representation of the confusion matrices generated by both approaches is presented for all subjects (see Figures 1617181920212223242526). Note that the surface of each element represents the sum of predictions generating the same label for the trials of each class. The distribution of ambiguities is similar along both methods, however, it is possible to appreciate the slightly superiority of the OsM approach for almost all subjects. Also note that labels are assigned according to the nomenclature discussed in section 5.2 (see Figure 5.4 in page 85). The main hypothesis upon which the proposed grouping of classes has been established states that, EEG data recorded during a combined motor imagery, can be characterized as the superposition of the activity that is independently generated by the sources associated to the body parts engaged in the task. In this way it is possible to generate a series of binary classification problems in which the activity induced by the main sources is considered separately, so that one class is conformed by all mental states inducing ERD activity within the corresponding source (i.e., all motor imageries involving the use of the associated limb), and the other class by the mental states that do not induce ERD activity over the associated source (i.e., all motor imageries that do not involve the use of the associated limb). This hypothesis has been motivated by the distribution of the relative oscillatory power that is observed over the primary motor cortex during the different motor imageries; which, as can be seen from figures 2.9 in page 29 and 5.4 in page 85, appears in the contralateral hemisphere of the body parts that are involved during the different motor tasks. Moreover, we have presented a series of statical analysis to confirm the plausibility of such assumption (see figure 2.15 in page 35).

By using this scheme, we have applied the CSP algorithm to solve the binary problems associated to each activity source. In this way, it is possible to directly apply the CSP algorithm without having to resort to any multiclass extension, such as the one-versus-one or one-versus-all approaches, and hence to optimize its performance. The first method that we have presented, namely, One-step Multilabel (OsM) approach, directly applies this solution, so that P binary problems, where P corresponds to the number of body parts considered by the paradigm, are formulated (2 binary problems for the 4-class database and 3 binary problems for the 8-class database). Subsequently, the extracted features are concatenated and used to train/validate a classification model.

The two other contributed algorithms, namely, Hierarchical Multilabel (HM) and One-step Hierarchical Multilabel (OsHM) methods, consist of two variants of the OsM approach. These two methods were motivated by the poor discrimination that was observed over the features generated by the binary problem associated to the central activity source (i.e., related to the feet) after applying the first approach on the 8-class database. As mentioned, this deficiency to well separate the corresponding elements is presumably caused by the fact that the central activity source is located too close from the sources associated with the hands, which influences its behavior. By considering this aspect, the two remaining methods have been conceived in such a way that the main decisions are based on the activity elicited by the sources related to the hands and, on a complementary level, by the activity that is linked to the feet. Thus, the HM method consists of a hierarchical approach for which the first stage of the classification involves a binary problem comprising in one class all motor imageries involving the use of one hand and, in the second class, the motor imageries that do not include it. Once the prediction is made, half of the classes are discarded and the process continues in the second stage where, according to the group that was previously chosen, one of two models is selected. Each one of these models consist again in a binary problem using the same type of partition: from the remaining classes and given that the hand that was already analyzed was found to be active or inactive in terms of motor imagery-related activity, one class is conformed by the motor tasks involving the use of the second hand and, the other class, by the motor tasks that do not include it. Finally, once the prediction for the second stage is made, there are only 2 classes left; which present the same activity associated with the left and right hands but differ in the activity that is linked to the feet. In this way, by considering only 2 classes, the task associated with the feet is simplified and addressed by a third and last stage. In this regard, it must be noticed that the HM method stops in the second stage when it is applied on the 4-class database since the feet were not included in the paradigm. It should be also noticed that the HM method can be applied by considering in the first stage either the left hand or the right hand. In this sense, there was a slightly improvement when starting with the right side, though it is not really significant. Even when separation was better in the case of the 8-class database for the binary problems involving the feet, results remain below that the ones achieved by the OsM and OsHM methods. This might be presumably due to the fact that, if there is a misclassified element at any stage, the predicted label will be wrong, whereas for the other two methods, classification is carried out by a single model over the concatenation of features.

Finally, the OsHM algorithm results from the combination of the two previous methods; which consists in taking the 2 binary problems that were introduced for each hand by the OsM method, and incorporates the 4 models that were formulated by the HM approach for the analysis of the activity generated by the feet. All features are computed independently and simultaneously, and the resulting elements are concatenated in order to train/validate a classification model. The

OsHM approach was applied only on the 8-class database since in the case of the 4-class one it is equivalent to the OsM method.

The best grand average performance across subjects for the 8-class database was obtained by the OsM method in combination with the ACSP algorithm, which is followed be the OsHM in combination with the CSP algorithm. As mentioned in the discussion, even though the OsHM approach generates more discriminative features than the OsM method by formulating the four instances involving the use/not-use of the feet in relation to the engagement of the hands, some ambiguities are caused by data whose classes were not considered during the models training, which hinders an effective classification. However, both approaches are slightly better than those generated by the classical methods, which is of particular interest if we consider that, in order to solve the same problem, the one-versus-one and one-versus-all approaches require respectively, 28 and 8 feature extraction models; which is by far more than the 3, and 6 that are needed, respectively, by the OsM and OsHM approaches. In the case of the 4-class database the best results were obtained by the HM approach in combination with FBCSP algorithm and linear discriminant analysis for classification. In this sense, since the problem involves only two classification stages, the risk of failure is not as high as in the case of the 8-class database.

Results and the corresponding discussion were presented, respectively, in chapters 6 and 7; where, aiming at comparing the performance yielded by each one of the proposed feature extraction methods, we have also included the accuracy generated by the classical one-versus-one and one-versus-all approaches, as well as the one obtained by the CSP by Joint Diagonalization (CSP by JAD) algorithm. In order to explore different solutions for all cases involving the use of the CSP algorithm, we have presented the results generated by two of its variants, namely, ACSP and FBCSP. Moreover, we have also reported the performance of the Minimum Distance to Riemannian Mean (MDRM) algorithm, which has been drawing the attention of the BCI community due to its capability for discriminating sensorimotor rhythms. In this regard, we have contributed with a new approach that uses the CSP method in order to increase the contrast between two conditions before mapping their covariance matrices into the Riemannian space, which even when it does not outperform the CSP approach, yields better results than the MDRM applied independently, showing that it is advantageous to use the CSP algorithm as a preprocessing method.

Interestingly, for classification results were always better by using linear discriminant analysis than support vector machines; presumably because the latter approach is more likely to get overtrained.

To conclude, we have proposed three new classification methods that are capable of distinguishing an unusually high number of motor imageries to enable the full 3D control of a robotic arm by using EEG. Under the proposed control scheme, the system manages to map 14 different movements by using only three body parts. Results show a varying accuracy depending on the subject, but it must be considered that the amount of data was limited and that the grand average performance of two of the proposed methods exceeded the classical one-versus-one and one-versus-all approaches.

Contributions

The main contribution of the present work is the evidence that EEG data recorded during a combined motor imagery can be modelled as the superposition of the activity that is independently generated by the sources associated to the body parts engaged in the task. This assumption can be exploited in different ways, as we have done by proposing an innovative way of grouping classes to address a multilabel paradigm. From this scheme, we have presented three new methods that not only optimize the feature extraction method by allowing to directly use the CSP algorithm on different binary problems, but they also reduce the complexity of the classification task, since contrarily to the classical feature extraction approaches based on the one-versus-one and one-versus-all algorithms, for which different features are required to train each one of the involved classifiers, the proposed methods use the same over all models, which has a positive impact on 8.2. LIMITATIONS the computation time.

Besides the three feature extraction methods that have been proposed, we have also explored and contributed with different solutions to the standard CSP algorithm. As in the case of the new approach presented in section 4.6.2.2, which uses the spatial filters generated by the CSP algorithm as a preprocessing step before mapping the covariance matrices of the CSP projected trials into the Riemannian manifold, from where the class labels of unseen elements can be assigned by using the Minimum Distance to Riemannian Mean (MDRM) algorithm. Even though the performance of such approach was not better than the one achieved by the standard CSP algorithm, it outperforms in all cases the one that is obtained by directly applying the MDRM method, showing that it is advantageous to merge both approaches. The justification of this behaviour is presumably given by the fact that, the spatial filters generated with the CSP algorithm extract information that enhances the difference between two conditions; which is preserved by the covariance matrices of the projected trials. Since such covariance matrices serve as descriptors, from where it is possible to extract the spatial information, such enhancement is reflected in the method performance.

Furthermore, we have also implemented some variants of well-known methods that, to our best knowledge, have not been reported before; such as the Filter Bank Analytical Common Spatial Pattern (FBACSP) algorithm, and a hierarchical approach based on the Minimum Distance to Riemannian Mean (MDRM) method [START_REF] Lindig-León | Hierarchical classification using riemannian geometry for motor imagery based bci systems[END_REF].

Finally, another important aspect to consider from the implementation design, is the switchingmode control, which allows employing the same commands to carry out different actions. Under this scheme it is possible to manage an even higher number of commands, from which it is possible to afford users which a full 3D control.

Limitations

The principal limitation of the proposed methods lies in their major attribute; which consists of the way how classes are grouped to generate the binary problems associated to each activity source. While this may optimize performance and computation time for users whose EEG patterns present a strong lateralization, it may not be effective for those subjects whose EEG activity lacks this differentiation. However, most of the feature extraction methods for sensorimotor rhythms discrimination are based on such a contrast, so this issue does not really represent an inherent drawback of the proposed approach, but a reality that has to be faced by BCI researchers.

Another limitation related to the operation scheme is the switching-mode function; which has the benefit of drastically increasing the number of available commands, but presents the drawback of slowing down the control of the device. In this regard, it is reasonable to consider imposing a particular assignment of classes that associates the most discriminative ones to the most frequent commands. In this way, misclassification rate of common functions, such as the switching component, would be reduced and thus efficiency enhanced.

123

In general, the proposed implementations share the same complications and problems as any motor imagery based EEG BCI system. For instance, one of the biggest challenges is to understand and solve the problem of BCI-Illiteracy, which implies that a substantial portion of users (estimated within 15% and 30%) do not manage to achieve control on this kind of systems [START_REF] Vidaurre | Towards a cure for bci illiteracy[END_REF]. This condition appears in different levels and, for some users, it is possible to generate a

Future work

The main interest is to implement the proposed solution to enable an online control, which will provide users with a self-paced operation mode affording them with the autonomy to carry out actions when desired (i.e., an asynchronous protocol). To this end there are two aspects to be considered: the online detection of the different commands, and the duration of each motion. For the first consideration, either the OsM approach or the OsHM method can be applied by using a sliding window, and thus to detect at each step the motor task that the user is performing. For the second aspect, we propose to implement the post-movement rebound detection as we have presented in [START_REF] Lindig-León | Alpha rebound improves on-line detection of the end of motor imageries[END_REF]. In this way, users will have more control to direct the device according to their needs; contrary to a predefined motion duration, where the users goal could lie in-between the range of movement. Moreover, regardless the application that is attempted to be controlled by using the proposed implementation, there are some considerations that can help to facilitate the classification task. For instance, in the case of the robotic arm, as it was discussed in section D.1, there is a limited range of movements once the device has reached an extreme position, so that even if the system detects a command demanding a motion that is beyond the possibilities of movement, the corresponding action will not be executed. In this sense, an effective approach will be to reduce the classification task in such a way that only the mental states related to a suitable 124 8.3. FUTURE WORK command are considered, so that less classes are evaluated and thus the classification task is facilitated, which would result in an improvement of the system performance. Another aspect that may improve the system performance consists in using users reaction upon a predicted action as feedback to continuously calibrate the system. It is well-known that a subject's recognition of having committed a response error is accompanied by specific variations that can easily be observed in averaged event-related potentials (ERP) [START_REF] Blankertz | Single trial detection of eeg error potentials: a tool for increasing BCI transmission rates[END_REF], which can be used in many ways either for avoiding false positive actions, or for generating adaptive classification models. As part of the feature work, it remains to explore this solution in order to provide an adjustable system.

As mentioned, a very important aspect that has to be considered in order to develop effective applications consists in conceiving adequate training techniques. While important progress has been made on this field [START_REF] Jeunet | Why standard brain-computer interface (BCI) training protocols should be changed: an experimental study[END_REF], the use of combined movements demands an even higher competence than the required over simple motor imageries; which has to be address if satisfactory performance is to be achieved. In this sense, it is interesting to note that the proposed methods favor a sensorimotor rhythms discrimination that is based on their lateralization, which represents one of the main physiological characteristics of such mental task. In this regard, such scheme for feature extraction might be more suitable than the one-versus-one approach for providing feedback during users' training, given than the resulting improvement would be based on the achievement of this specific differentiation, which is directly related to the metal process.

Finally, as mentioned in section 6.1, there are some subject-specific experimental parameters that can be tuned to increase the system performance, namely, the frequency for band-pass filtering of the EEG signals, the time window of the EEG trials taken relative to the stimuli, and the subset of CSP filters to be used [START_REF] Blankertz | The non-invasive berlin brain-computer interface: fast acquisition of effective performance in untrained subjects[END_REF]. In this sense, it remains to integrate optimization algorithms to find the best possible adjustment of such parameters working in combination with the proposed methods, and thus to boost the general performance of the system. In the same way, regularization of the classification algorithms would be also pertinent.

A Résumé Étendu

Les interfaces cerveau-ordinateur (BCI) substituent aux sorties naturelles du système nerveux des sorties artificielles qui n'exigent pas l'utilisation des nerfs périphériques. Ceci permet à des patients atteints de troubles moteurs sévères d'interagir, en utilisant seulement leur activité cérébrale, avec différents types d'applications: correcteur orthographique, neuroprothèse, fauteuil roulant, ou dispositifs robotiques [START_REF] Wolpaw | BCI: Principles and Practice[END_REF], [START_REF] Gomez-Gil | Brain computer interfaces -a review[END_REF]. L'utilisation des interfaces cerveau ordinateur n'est pas restreinte, et elles ont également suscité de l'intérêt pour le développement de programmes de réadaptation destinés à des patients atteints d'AVC [START_REF] Pfurtscheller | Event-related synchronization (ERS): an electrophysiological correlate of cortical areas at rest[END_REF], [START_REF] Fok | An EEG-based brain computer interface for rehabilitation and restoration of hand control following stroke using ipsilateral cortical physiology[END_REF], ainsi que des interfaces pour les personnes en bonne santé (par exemple, dans les applications de jeu).

Plusieurs types de signaux cérébraux peuvent être utilisés pour le contrôle BCI. L'un d'entre eux, particulièrement intéressant puisqu'il a l'avantage de favoriser les paradigmes auto-stimulés (c'est-à-dire un mode de fonctionnement asynchrone), est connu sous le nom d'imagerie motrice (IM); ceci peut apporter un contrôle intuitif et naturel pour permettre à des personnes à mobilité réduite une interaction riche avec leur environnement. L'imagerie motrice implique la pratique mentale d'une tâche motrice sans effectuer de mouvement apparent [START_REF] Pfurtscheller | Motor imagery and direct brain-computer communication[END_REF]. L'imagination d'un mouvement implique les mêmes régions du cerveau, et des fonctions similaires à celles engagées dans la programmation et l'exécution d'un mouvement réel. La principale différence est que dans le premier cas, les mouvements sont bloqués au niveau cortico-spinal [START_REF] Neuper | Evidence for distinct beta resonance frequencies in human EEG related to specific sensorimotor cortical areas[END_REF].

L'activité EEG obtenue lors de l'imagerie motrice consiste en des changements de rythmes sensori-moteurs qui présentent un comportement oscillatoire particulièrement intéressant dans deux gammes de fréquences spécifiques: les bandes alpha (8-13 Hz) et beta [START_REF] Barachant | Multi-class brain computer interface classification by riemannian geometry[END_REF][START_REF] Dornhege | Increase information transfer rates in BCI by CSP extension to multi-class[END_REF][START_REF] Dornhege | Boosting bit rates in noninvasive eeg single-trial classifications by feature combination and multiclass paradigms[END_REF][START_REF] Weibo | Eeg feature comparison and classification of simple and compound limb motor imagery[END_REF][START_REF] Royer | EEG) control of a virtual helicopter in 3-dimensional space using intelligent control strategies[END_REF][START_REF] Zhou | A novel technique for phase synchrony measurement from the complex motor imaginary potential of combined body and limb action[END_REF][START_REF] Scherer | Can synergistic finger movements be represented as the superposition of individual finger movements?[END_REF][START_REF] Blankertz | Optimizing spatial filters for robust EEG single-trial analysis [revealing tricks of the trade[END_REF][START_REF] Lindig-León | Comparison of sensorimotor rhythms in EEG signals during simple and combined motor imageries over the contra and ipsilateral hemispheres[END_REF][START_REF] Lindig-León | A multilabel classification method for detection of combined motor imageries[END_REF][START_REF] Mason | Temporal control paradigms for direct brain interfaces -rethinking the definition of asynchronous and synchronous[END_REF][START_REF] Mason | A comprehensive survey of brain interface technology designs[END_REF][START_REF] Mak | Clinical applications of brain-computer interfaces: current state and future prospects[END_REF] [7]- [START_REF] Babiloni | Human cortical electroencephalography (EEG) rhythms during the observation of simple aimless movements: a high-resolution EEG study[END_REF], puisqu'il a été montré qu'ils sont modulés pendant et après la préparation et l'exécution des motrice d'un seul membre [START_REF] Royer | EEG) control of a virtual helicopter in 3-dimensional space using intelligent control strategies[END_REF]. Plus récemment, on a étudié les caractéristiques des modèles EEG au cours d'imagerie motrice simple et combinée, afin de confirmer la séparation de sept classes différentes pour la mise en oeuvre des systèmes BCI à base d'imagerie motrice [START_REF] Weibo | Eeg feature comparison and classification of simple and compound limb motor imagery[END_REF]. Une autre étude a rapporté que le mouvement de plusieurs doigts en même temps n'est pas la simple superposition de mouvements de doigt individuels à l'échelle ECoG, mais plutôt l'activation de réseaux de neurones supplémentaires. Ceci implique que pour mettre en oeuvre une interface BCI, la génération de signaux de contrôle codant des mouvements complexes devrait probablement utiliser des composantes de signal d'ordre plus élevé que la superposition des signaux individuels des membres [START_REF] Scherer | Can synergistic finger movements be represented as the superposition of individual finger movements?[END_REF]. Cependant, à notre connaissance, cette dernière conclusion reste à confirmer à l'échelle EEG; ceci est la question principale de notre travail.

B Time series of ERD/ERS% patterns

Along this section a series of complementary images are presented in order to show the ERD/ERS% patterns that are obtained independently for each subject. Images were obtained by applying the same method described in sections 2.5 and 3.1.2.

For each subject such procedure was applied over the band-passed filter trials over the alpha [START_REF] Pfurtscheller | Beta rebound after different types of motor imagery in man[END_REF][START_REF] Toro | Event-related desynchronization and movement-related cortical potentials on the ECoG and EEG[END_REF][START_REF] Babiloni | Human cortical electroencephalography (EEG) rhythms during the observation of simple aimless movements: a high-resolution EEG study[END_REF][START_REF] Avanzini | The dynamics of sensorimotor cortical oscillations during the observation of hand movements: an EEG study[END_REF][START_REF] Pfurtscheller | Event-related EEG/MEG synchronization and desynchronization: basic principles[END_REF][START_REF] Pfurtscheller | On the existence of different types of central beta rhythms below 30 hz[END_REF][START_REF] Barachant | Multi-class brain computer interface classification by riemannian geometry[END_REF][START_REF] Dornhege | Increase information transfer rates in BCI by CSP extension to multi-class[END_REF] and alpha+beta bands, within which time series were computed along the time period from [0.5-4.

C Offline Results

Along this section the detailed description of the offline results presented in chapter 6 will be presented, which consist of the accuracy yielded by each subject for each one of the described methods. The first case corresponds with the results generated by the standard CSP algorithm described in section 4.1.3, together with its analytical version; the ACSP approach, which has been discussed in section 4.2. Subsequently, we report the performance obtained by the FBCSP described in section 4.4, and by its analytical version; the Filter Bank Analytical Common Spatial Pattern (FBACSP) algorithm, which corresponds with one of the contributions of the present work. Afterwards, we present another variant of the CSP method from which feature extraction is directly achieved by considering all classes at once, without having to resort to any complementary technique, namely, the CSP by JAD approach (see section 4.5). Since the aforementioned method consists in recursively approximating a transformation matrix for the joint diagonalization of real covariance matrices, its analytical version is not considered within the scope of the present work.

As discussed in section 4.6.2.1, the MDRM consists of another approach that allows to directly extract features without using any auxiliary technique. This method has been applied over both databases as an alternative technique for sensorimotor rhythms discrimination. Moreover, we have also proposed a convenient variant of the MDRM method, where the CSP algorithm is applied before computing the covariance matrices that define the Riemannian manifold, which are used as the EEG descriptors, so that more discriminative features are generated (see section The idea is to generate more discriminative covariance matrices, which are used as descriptors within the Riemannian manifold, by computing them over the CSP projected trials instead of over the band-pass filtered segments, as it is done by the standard MDRM algorithm. Other than that, the method follows the same principle. For generating the CSP filters we have applied both the one-versus-one and the one-versus-all grouping techniques. C.9.6 Results: HM method by using the CSP + MDRM approach

We have also applied the HM approach in combination with the CSP algorithm preceding the MDRM method, so that the EEG covariance matrices become more discriminative between conditions after the spatial filtering (see section 4.6.2.2).

Tables C. [START_REF] Pfurtscheller | Patterns of cortical activation during planning of voluntary movement[END_REF] and C.56 show, respectively, the performance over the 4-class database when the analysis starts on the activity generated by the left hand, and when it starts on the right hand.

Similarly We have also applied the OsHM approach in combination with the CSP algorithm preceding the MDRM algorithm as discussed in section 4.6.2.2. Which, in short, consists in the same procedure described in the preceding section with the only difference that the covariance matrices, which are used as the EEG descriptors, are computed over the trials projection obtained after applying the CSP algorithm, instead of over the band-pass filtered signal segments (as it is done by the standard MDRM). Note that the CSP filters are computed exactly as it was described in section C.10.1. Table C.64 presents the correspondig results. In this sense, the proposed paradigm adds a greater difficulty in comparison to the standard protocols, since the movement imagination of more than one body part at the same time requires a higher degree of coordination. Furthermore, training sessions can be long and tedious for users, specially if the task is not well accomplished, which might cause frustration and make the training process even more difficult. It must be mentioned that this aspect represents one of the main areas of BCI research, which was beyond of the scope of the present project and, therefore, it remains to be addressed as part of the related feature work. One of the interest of using motor imagery for controlling robotic devices is the benefit of providing users with a self-paced paradigm, which makes sense only if the implementation is built upon an online setup (i.e., commands are translated in real time).

In order to implement an online setup there are several aspects to be addressed; perhaps the most important one consists in achieving an accurate classification that delivers results as quickly as possible, so that users are provided with an almost real-time feedback. In this sense, the present work has focused in the implementation of feature extraction methods that favor classification performance from a multilabel basis, so that more discriminative features are generated and computational costs are optimized. There are other aspects to be considered for building online applications, especially if the goal is to use them on a daily basis in out-of-the-lab environments.

These aspects are to be addressed as part of the feature work of this project. In the two following sections we mention some characteristics that can favor the implementation of an online setup in combination with the proposed feature extraction methods.

D.2 3D control of a robotic arm, implementation example

The following session presents the implementation of the One-step Multilabel method in combination with the CSP algorithm and LDA classification for controlling the aforementioned robotic

3. 4

 4 ERD/ERS% time series within the mu range for all motor tasks included in the 8-class database over electrode C z . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5 ERD/ERS% time series within the mu range for all motor tasks included in the 8-class database over electrode C4 . . . . . . . . . . . . . . . . . . . . . . . . . . 3.6 Topographic representations of the correlation between the conditions included in the 8-class database in terms of R 2 . . . . . . . . . . . . . . . . . . . . . . . . . . 3.7 Box plots comparing the power spectrum magnitude at 12 Hz among all the considered conditions over electrode C3 . . . . . . . . . . . . . . . . . . . . . . . 3.8 Box plots comparing the power spectrum magnitude at 12 Hz among all the considered conditions over electrode C z . . . . . . . . . . . . . . . . . . . . . . . . 3.9 Box plots comparing the power spectrum magnitude at 12 Hz among all the considered conditions over electrode C4 . . . . . . . . . . . . . . . . . . . . . . . 4.1 Distribution of data and their covariance matrices. . . . . . . . . . . . . . . . . . 4.2 Normal distribution. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3 EEG trials during left hand and right hand motor imageries in electrodes C3 and C4. 4.4 EEG trials during left hand and right hand motor imageries in electrodes C3 and C4 after having applied the CSP algorithm. . . . . . . . . . . . . . . . . . . . . . 4.5 CSP projected trials. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.6 CSP features extracted from the changes in variance of the projected trials. . . . 4.7 Analytic Common Spatial Pattern algorithm. . . . . . . . . . . . . . . . . . . . . 4.8 Classification space generated by the CSP by JAD algorithm for a 3-class problem. 4.9 Geometry of a linear discriminant function in two dimensions generated by LDA. 4.10 Discriminant functions for multiple classes. . . . . . . . . . . . . . . . . . . . . . 4.11 Support Vectors Machines classification. . . . . . . . . . . . . . . . . . . . . . . . 4.12 Riemannian manifold . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.13 MDRM classification. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.14 MDRM classification in combination with the CSP algorithm. . . . . . . . . . . . 5.1 One-versus-one classification for a 4-class problem. . . . . . . . . . . . . . . . . . 5.2 One-versus-all classification for a 4-class problem. . . . . . . . . . . . . . . . . . . 5.3 Hierarchical approach for a 4-class problem. . . . . . . . . . . . . . . . . . . . . . 5.4 Labeling convention. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.5 Architecture of the One-step Multilabel Approach (OsM) for solving the 4-class problem. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.6 Architecture of the One-step Multilabel Approach (OsM) for solving the 8-class problem. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.7 Architecture of the Multilabel Hierarchical approach (MH) for solving the 4-class problem. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xii 5.8 Architecture of the Multilabel Hierarchical approach (MH) for solving the 8-class problem. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.9 Architecture of the One-step Multilabel Approach (OsM) for formulated for the 8-class database. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.1 4-fold cross-validation process for assessing performance over the 4-class database by using the discussed feature extraction methods in combination with the LDA and QDA algorithms, as well as for the distance-based approaches. . . . . . . . . 6.2 Nested cross-validation process for assessing performance over the 4-class database by using the discussed feature extraction methods in combination with the SVM algorithm. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.3 30-fold cross-validation process for assessing performance over the 8-class database for all discussed methods in combination with the LDA and QDA algorithms, as well as for the distance-based approaches. . . . . . . . . . . . . . . . . . . . . . . 6.4 Comparison of the accuracy results obtained by using the proposed methods in combination with the CSP algorithm for all subjects of the 8-class database. . . 6.5 Comparison of the mean accuracy results obtained by using the proposed methods in combination with the CSP algorithm for all subjects of the 8-class database. . 6.6 Comparison of the mean accuracy results obtained by using the proposed methods in combination with the ACSP algorithm for all subjects of the 8-class database. 6.7 Comparison of the accuracy results obtained by using the proposed methods in combination with the ACSP algorithm for all subjects of the 8-class database. . . 6.8 Comparison of the mean accuracy results obtained by using the proposed methods in combination with the ACSP algorithm for all subjects of the 8-class database. 6.9 Comparison of the mean accuracy results obtained by using the proposed methods in combination with the ACSP algorithm for all subjects of the 8-class database. 10 Most discriminative CSP patterns obtained for each one of the binary problems formulated by the OsM approach. . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 First two features obtained with the CSP method for each one of the binary problems formulated by the OsM method. . . . . . . . . . . . . . . . . . . . . . . 12 Most discriminative CSP patterns obtained for each one of the binary problems formulated by the HM method starting with the left hand . . . . . . . . . . . . . 13 Most discriminative CSP patterns generated for each one of the binary problems formulated by the MH method starting with the right hand . . . . . . . . . . . . 14 Distribution of the 2 first CSP features generated by the MH approach. . . . . . 15 Distribution of the two first CSP features generated by the OsHM approach for the four instances involving the feet motor task. . . . . . . . . . . . . . . . . . . . 16 Confusion matrices obtained by the one-versus-one and OsM approaches for subject 1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiii 2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 Confusion matrices obtained by the one-versus-one and OsM approaches for subject 3. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 Confusion matrices obtained by the one-versus-one and OsM approaches for subject 4. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 Confusion matrices obtained by the one-versus-one and OsM approaches for subject 5. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 Confusion matrices obtained by the one-versus-one and OsM approaches for subject 6. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 Confusion matrices obtained by the one-versus-one and OsM approaches for subject 7. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 Confusion matrices obtained by the one-versus-one and OsM approaches for subject 8. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 Confusion matrices obtained by the one-versus-one and OsM approaches for subject 9. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 Confusion matrices obtained by the one-versus-one and OsM approaches for subject 10. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 Confusion matrices obtained by the one-versus-one and OsM approaches for subject 11. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . B.1 Oscillatory power for subject 1 within the alpha band. . . . . . . . . . . . . . . . B.2 Oscillatory power for subject 1 within the alpha and beta bands. . . . . . . . . . B.3 Oscillatory power for subject 2 within the alpha band. . . . . . . . . . . . . . . . B.4 Oscillatory power for subject 2 within the alpha and beta bands. . . . . . . . . . B.5 Oscillatory power for subject 3 within the alpha band. . . . . . . . . . . . . . . . B.6 Oscillatory power for subject 3 within the alpha and beta bands. . . . . . . . . . B.7 Oscillatory power for subject 4 within the alpha band. . . . . . . . . . . . . . . . B.8 Oscillatory power for subject 4 within the alpha and beta bands. . . . . . . . . . B.9 Oscillatory power for subject 5 within the alpha band. . . . . . . . . . . . . . . . B.10 Oscillatory power for subject 5 within the alpha and beta bands. . . . . . . . . . B.11 Oscillatory power for subject 6 within the alpha band. . . . . . . . . . . . . . . . B.12 Oscillatory power for subject 6 within the alpha and beta bands. . . . . . . . . . B.13 Oscillatory power for subject 7 within the alpha band. . . . . . . . . . . . . . . . B.14 Oscillatory power for subject 7 within the alpha and beta bands. . . . . . . . . . B.15 Oscillatory power for subject 8 within the alpha band. . . . . . . . . . . . . . . . B.16 Oscillatory power for subject 8 within the alpha and beta bands. . . . . . . . . . B.17 Oscillatory power for subject 9 within the alpha band. . . . . . . . . . . . . . . . B.18 Oscillatory power for subject 9 within the alpha and beta bands. . . . . . . . . . B.19 Oscillatory power for subject 10 within the alpha band. . . . . . . . . . . . . . . xiv B.20 Oscillatory power for subject 10 within the alpha and beta bands. . . . . . . . . . B.21 Oscillatory power for subject 11 within the alpha band. . . . . . . . . . . . . . . B.22 Oscillatory power for subject 11 within the alpha and beta bands. . . . . . . . . . D.1 3-axis operation mode. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . D.2 Arm control. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . D.3 Implementation example. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . D.4 Confusion matrices obtained by the OsM approach for subject 3 during an online setup. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xv

Fig. 1 . 1 :

 11 Fig.1.1: General BCI architecture. There are three main components that every BCI includes in its architecture, namely, a signal acquisition module, a signal processing stage, and the corresponding output, which can presented in many different forms according to the application of the system.

Fig. 2 . 1 :

 21 Fig. 2.1: Main EEG Rhythms. The electrical activity generated by neurons can be separated into frequency bands, which have been associated to different brain functions. In this way, by using spectral techniques, it is possible to independently analyze patterns within the different brain rhythms.

Figure 2 .

 2 Figure 2.3 shows an adaptation of Penfield's homunculus, where it is possible to see the distribution of the different body parts along their contralateral area over the primary motor cortex in accordance with electrodes C3, C z , and C4 of the international 10-20 system. Legs and feet are represented within the medial wall of the cerebrum; the trunk, upper arm, and hand are found consecutively and laterally along the hemisphere, with the face located at the most

Fig. 2 . 2 :

 22 Fig. 2.2: The standard 10-20, 10-10 and 10-5 electrode montages. The 10-20 montage is indicated by the 21 electrodes represented as black circles. The 10-10 montage, which has a total of 74 electrodes, includes the 21 electrodes of the 10-20 montage plus 53 additional ones indicated by gray circles. The 10-5 montage comprises the 74 electrodes of the 10-10 system together with 68additional ones represented by dots and white circles, which gives a total of 142 electrodes that provide a more exhaustive and uniform coverage of the head. The nomenclature for the standard 10-20 and 10-10 and 10-5 extensions is indicated. Note that electrodes on the right side have even numbers and electrodes on the left side have odd numbers, while electrodes along the midline are indicated by "z" ("z"=zero)[START_REF] Oostenveld | The 5% electrode system for high-resolution EEG and ERP measurements[END_REF].

Figure 2 .

 2 Figure 2.4 presents a series of topographic maps showing the sensorimotor rhythms over the mu band induced by a real left hand movement (fig. 2.4-A) and an imagined left hand movement (fig. 2.4-B), in comparison with the activity induced by a real right hand movement (fig. 2.4-C) and an imagined right hand movement (fig. 2.4-D). EEG signals belong to a right-handed male subject that has participated in the development of the database that will be described along the following sections. In all maps it is possible to observe the ERD phenomena characteristic of the execution of a motor task, which is represented in blue colors associated to the minimum values. For both left and right hand motor tasks, patterns are more enhanced in the case of real movements, which in comparison to the motor imageries, present darker tonalities and more defined contours around the associated electrodes over their contralateral side (C4 for left hand and C3 for right hand, see fig. 2.3 in regard of the lateralization of the somatotopically organization).

Figure 2 .

 2 Figure 2.5 presents a series of topographic maps showing the sensorimotor rhythms over the beta band induced after the termination of a real left hand movement (fig. 2.5-A) and an imagined left hand movement (fig. 2.5-B), in comparison with the activity induced after the termination of a real right hand movement (fig. 2.5-C) and an imagined right hand movement (fig. 2.5-D).EEG signals belong to the same subject considered in the previous example. In all maps it is possible to observe the ERS phenomena characteristic of the termination of a motor task (e.g., post-movement rebound), which is represented in red colors associated to the maximum values. As in the previous example, ERS patterns are observed around the associated electrodes over the contralateral side (C4 for left hand and C3 for right hand, see fig.2.3).

Fig. 2 . 4 :

 24 Fig.2.4: Topographic maps during the course of motor tasks within the mu band. A) presents the brain activity generated during a real left hand movement execution, whereas B) shows the activity generated during a left hand imagined movement. In both cases it is possible to observe the characteristic ERD patterns that are induced in the contralateral hemisphere (around C4, see fig.2.3), although these are more enhanced in the case of the real movement. C) and D) show, respectively, the real movement execution and the imagined movement of the right hand. In the same way, it is possible to observe the lateralization of the ERD patterns associated to the motor tasks (around C3, see fig.2.3). Spatial mapping was carried out by using a cubic interpolation method[START_REF] Höhne | BBCI toolbox-offline analysis[END_REF]. Note that each pair of maps was independently normalized.

Fig. 2 . 5 :

 25 Fig.2.5: Topographic maps after the termination of motor tasks within the beta band. A) presents the brain activity generated after a real left hand movement execution, whereas B) shows the activity generated after a left hand imagined movement. In both cases it is possible to observe the characteristic ERS patterns (e.g., post-movement rebound) that are induced in the contralateral hemisphere (around C4, see fig.2.3). C) and D) show, respectively, the termination of a real movement execution and the termination of an imagined movement of the right hand. In the same way, it is possible to observe the lateralization of the ERS patterns associated to the motor tasks (around C3, see fig.2.3). Spatial mapping was carried out by using a cubic interpolation method[START_REF] Höhne | BBCI toolbox-offline analysis[END_REF]. Note that each pair of maps was independently normalized.

Figure 2 .

 2 Figure 2.6 shows a detailed example of the time course of ERD/ERS% computation according to equation 2.1. For which 80 trials of right hand motor imagery recorded over electrode C3 have been considered. EEG signals belong to a right-handed female subject that has participated in the development of the database that will be described along the following sections.Figure2.6-A shows the raw EEG signals. Note that, because of space limitations, only 3 trials are displayed but that the process includes the entire set. The motor imagery starts at second "0" and lasts along 6 seconds, which is indicated by a blue segment. As mentioned, the first step consists of the band-pass filtering of the signals which, in this case is carried out throughout a 5th order Butterworth filter within the mu band in order to enhance the characteristic ERD pattern associated to this band-with (see figure 2.6-B). Subsequently, the band-pass filtered signals are independently squared across samples (figure 2.6-C), where the power decrease is already pronounced along the motor imagery period.Figure 2.6-D presents the averaging across trials, where a 2-second bar has been drawn 2 seconds before the initiation of the motor task to indicate the baseline segment that was considered in order to smooth the signal according to equation 2.1 by using a 125 ms shifting step. Finally, the resulting time course of ERD/ERS is shown in figure 2.6-E. Negative values reveal a power decrease with respect to the baseline (reference) segment, which correspond to an ERD modulation. From the figure it is possible to

Fig. 2 . 6 :

 26 Fig. 2.6: Processing method to observe the time course of ERD/ERS% according to equation 2.1. EEG segments correspond to right hand motor imageries recorded on electrode C3 (A), (B) shows the resulting band-pass filtered signals within the mu band by using a 5th order Butterworth filter, (C) illustrates the squaring of samples, (D) presents the averaging across trials, and (E) depicts the resulting time cours of ERD/ERS highlighted, respectively, in blue and red colors.

2. 6 .

 6 1 4-class database EEG signals were recorded through the openViBE 1 [85] platform from five right-handed healthy subjects at 256 Hz using a commercial REFA amplifier developed by TMS International™.Subjects were asked to perform three different motor imageries (left hand, right hand, and both hands) together with rest condition, while seated in a comfortable chair with the arms at their sides in front of a computer screen showing the task cue to be performed. The whole session consisted of 4 runs containing 20 trials per task for a total of 80 trials per class. All experiments were carried out with the consent agreement of each participant and following the statements of the WMA declaration of Helsinki on ethical principles for medical research involving human 1. http://openvibe.inria.fr subjects.

  2.7).

Fig. 2 . 7 :

 27 Fig. 2.7: Distribution of the 26 electrodes considered for the 4-class database development.The EEG cap was fitted with 26 electrodes re-referenced with respect to the common average reference across all channels and located over the extended international 10-20 system positions to cover the primary sensorimotor cortex.

Fig. 2 . 8 :

 28 Fig. 2.8: Time scheme for the 4-class database. Each trial was randomly presented and lasted for 12 seconds. During the first 6 seconds, users were asked to perform the motor imagery indicated by the task cue, which was launched throughout two simultaneously displayed panels associated, from left to right, to the left hand and right hand. The use of each body part was indicated by the presence of arrows: an arrow pointing to the left side on the left panel for left hand, an arrow pointing to the right side on the right panel for right hand, and the simultaneous combination of the two of them for the corresponding compound task. The rest condition was also considered and it was specified by the absence of arrows. After second 6, the task cue disappeared and the crosses were remaining for the next 6 seconds indicating the pause period before the next trial started.

Fig. 2 . 9 :

 29 Fig.2.9: ERD/ERS% topographic representations within the mu range during the 6-second period of the four motor tasks. The spatial mapping corresponds to the grand average across subjects (n=5) of the time course values computed in terms of equation 2.1 over all electrodes by using a cubic interpolation method[START_REF] Höhne | BBCI toolbox-offline analysis[END_REF].

Fig. 2 .

 2 Fig. 2.10: ERD/ERS% time series within the mu range for all motor tasks over electrodes C3 and C4. ERD/ERS values were estimated according to equation 2.1 over the grand average accross subjects (n=5). The oscillatory power on electrodes C3 and C4 is presented, respectively, in (A) and (B).

Fig. 2 . 11 :

 211 Fig.2.11: ERD/ERS% topographic representations within the beta range during the 6-second period after the termination of the four motor tasks. The spatial mapping corresponds to the grand average across subjects (n=5) of the time course values computed in terms of equation 2.1 over all electrodes by using a cubic interpolation method[START_REF] Höhne | BBCI toolbox-offline analysis[END_REF].

Figure 2 .

 2 Figure 2.13 shows the resulting R 2 values generated for each comparison among all classes.

Figure

  Figure 2.13-A presents the comparison between left hand motor imageries and the rest condition;

Fig. 2 . 13 :

 213 Fig. 2.13: Grand average across subjects (n=5) of correlation between the considered conditions over each electrode in terms of R 2 . The different comparisons correspond to (A) left hand MIs and rest condition; (B) left hand MIs and both hands MIs; (C) right hand MIs and rest condition; (D) right hand MIs and both hands MIs; (E) left hand MIs and right hand MIs and F) both hands MIs and rest condition. High R 2 values indicate that there is a significant difference between the corresponding MIs at the given frequency over the corresponding electrode positions.

Fig. 2 . 14 :

 214 Fig. 2.14: Topographic representations of the grand average across subjects (n=5) of the correlation values between the considered conditions in terms of R 2 .

  30 trials per class (240 trials considering the 8 classes). For the stimulation, three panels were simultaneously displayed on the screen, which were associated from left to right, to the left hand, feet and right hand. Each trial was randomly presented and lasted for 12 seconds, starting at second 0 with a cross at the center of each panel and an overlaid arrow indicating for the next 6 seconds the task to be performed: an arrow pointing to the left side on the left panel for left hand, an arrow pointing down on the central panel for feet, an arrow pointing to the right side on the right panel for right hand, and the simultaneous combination of all of them for the corresponding combined tasks (see Figure 3.1 and column 2 of Figure D.2). The rest condition was also considered and it was indicated by the absence of arrows. After second 6, the task cue disappeared and the crosses were remaining for the next 6 seconds indicating the pause period before the next trial started.

Fig. 3 . 1 :

 31 Fig. 3.1: Time scheme. Each trial was randomly presented and lasted for 12 seconds. During the first 6 seconds, users were asked to perform the motor imagery indicated by the task cue, which was launched throughout three simultaneously displayed panels associated, from left to right, to left hand, feet and right hand. The use of each body part was indicated by the presence of arrows: an arrow pointing to the left side on the left panel for left hand, an arrow pointing down on the central panel for feet, an arrow pointing to the right side on the right panel for right hand, and the simultaneous combination of all of them for the corresponding compound motor imageries. After second 6, the task cue disappeared and the crosses were remaining for the next 6 seconds indicating the pause period before the next trial started.

Fig. 3 . 2 :

 32 Fig. 3.2: Distribution of the relative oscillatory power along the sensorimotor cortex during the course of the different motor tasks for subject 3. ERD/ERS% values were obtained for all electrodes within the mu frequency range and averaged over the 4-second period starting one second after the task cue was presented. Negative values indicate ERD% modulation, which is characteristic of a motor task execution that appears over the contralateral hemisphere of the body part used in the process (broadly, electrode C4 for left hand, C z for feet, and C3 for right hand). Note that each topographic map is independently normalized to enhance the corresponding patterns.

Fig. 3 . 3 :

 33 Fig.3.3: ERD/ERS% time series within the mu range for all motor tasks over electrode C3 (i.e, the activity source associated with the use of the right hand). ERD/ERS values were estimated according to equation 2.1 over the EEG trials from subject 3. Note that all motor imageries involving the use of the right hand (i.e., right hand, right hand in combination with feet, both hands, and both hands in combination with feet) present considerably lower values that those that do not include it (i.e., rest, left hand, feet, and left hand in combination with feet).

Fig. 3 . 4 :

 34 Fig. 3.4: ERD/ERS% time series within the mu range for all motor tasks over electrode C z . ERD/ERS values were estimated according to equation 2.1 over the EEG trials from subject 3.Note that the motor imageries involving the use of the feet (i.e., feet, left hand in combination with feet, right hand in combination with feet, and both hands in combination with feet) present lower ERD/ERS% values than those that do no it include it (i.e., rest, left hand, right hand, and both hands). However, this contrast is not as significant as the one observed on electrode C3, possibly due to the fact that the regions associated to both the left and right hands are quite close to electrode C z , which might cause it to incorporate some of this activity.

Fig. 3 . 5 :

 35 Fig.3.5: ERD/ERS% time series within the mu range for all motor tasks over electrode C4. ERD/ERS values were estimated according to equation 2.1 over the EEG trials from subject 3. Note that all motor imageries involving the use of the left hand (i.e., left hand, left hand in combination with feet, both hands, and both hands in combination with feet) present considerably lower values that those that do not include it (i.e., rest, right hand, feet, and right hand in combination with feet).

Fig. 3 . 6 :

 36 Fig. 3.6: Topographic representations for one of the subjects of the correlation values between the comparison over all conditions in terms of R 2 . High R 2 values indicate that there is a significant difference between the considered motor tasks over the corresponding positions. Note that each topographic map is independently normalized to enhance the corresponding patterns.

Fig. 3 . 7 :

 37 Fig.3.7: Box plots for subject 3 comparing the power spectrum magnitude at 12 Hz within the time interval from [2-5 s] among all the considered conditions over electrode C3. The distribution of data shows that there is no significant difference between motor imageries including the use of the right hand over this region, independently from the fact of whether the motor imagery is simple or combined. In contrast there is a significant difference with all the other tasks that do not include the use of the right hand which, in turn, present a uniform distribution between them.

Fig. 3 . 8 :

 38 Fig. 3.8: Box plots for subject 3 comparing the power spectrum magnitude at 12 Hz within the time interval from [2-5 s] among all the considered conditions over electrode C z . The distribution of data reflects what was observed from the ERD/ERS% time series, and confirms that EEG activity around the central region is influenced by the surrounding sources.

Fig. 3 . 9 :

 39 Fig.3.9: Box plots for subject 3 comparing the power spectrum magnitude at 12 Hz within the time interval from [2-5 s] among all the considered conditions over electrode C4. The distribution of data shows that there is no significant difference between motor imageries including the use of the left hand over this region, independently from the fact of whether the motor imagery is simple or combined. In contrast there is a significant difference with all the other tasks that do not include the use of the left hand which, in turn, present a uniform distribution between them.

Contents 4 . 1 49 4. 1 . 1 50 4. 1 . 2 54 4. 2 58 4. 3 62 4. 5 63 4. 5 . 1 64 4. 6 66 4. 6 . 1

 414911501254258362563516466661 Feature extraction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Covariance matrix . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Normal distribution and eigenvalue decomposition . . . . . . . . . . . . 51 4.1.3 Common Spatial Patterns (CSP) . . . . . . . . . . . . . . . . . . . . . . Analytical Common Spatial Patterns (ACSP) . . . . . . . . . . . . . Feature selection . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 4.3.1 Mutual information . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 4.4 Filter Bank Common Spatial Pattern (FBCSP) . . . . . . . . . . . . 61 4.4.1 Mutual Information-based Best Individual Feature (MIBIF) algorithm . Common Spatial Pattern by Joint Approximate Diagonalization (CSP by JAD) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Information theoretic feature extraction (ITFE) . . . . . . . . . . . . . . Classification . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Discriminant functions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66 4.6.2 Distance-based classification . . . . . . . . . . . . . . . . . . . . . . . . . 74

  the following chapter the concept of feature extraction is introduced and elaborated in the context of EEG signals during motor imagery. To this end, a widely used algorithm known as Common Spatial Pattern (CSP), which has been shown to be powerful at discriminating sensorimotor rhythms and lays the basis for the methods contributed by the present work, is elucidated together with two of its variants, namely, Analytic Common Spatial patterns (ACSP) and Common Spatial Pattern by Joint Approximate Diagonalization (CSP by JAD). Additionally, we will present a strategy called Filter Bank Common Spatial Pattern (FBCSP) that optimizes the algorithm performance by decomposing the EEG signals into a series of frequency bands,

A

  fundamental machine learning task consists in searching for patterns in data; which can be achieved by tuning the parameters of an adaptive model that incorporates the characteristics of a data set that is called training set. The categories (or classes) of the elements contained in the training set are known in advance, typically as the result of conducting a supervised experiment.These categories can be expressed by a target vector t representing the identity of each element x. Then, the aforementioned model can be expressed as a function y(x); which takes a new element x as input to generate an output vector y that is encoded in the same manner as the target vectors. The specific form of the function y(x) is determined by the training set during the learning phase. Once the model is trained, it is possible to determine the identity of new elements that are said to compose the testing set, which must be processed in the same way as the training data. In most of the cases, the variability of the input vectors will be such that the training set can not represent the entire population, and thus the task to correctly categorize new elements differing from those used in the training phase becomes a challenge. For many applications, the original input vectors are transformed into a new space of variables where, hopefully, the pattern recognition problem will be easier to solve. This processing stage is called feature extraction and, basically, it consists of creating derived values (features) from the original input vectors that are meant to be informative and non-redundant, so that the subsequent generalization and interpretation are facilitated due to the reduced representation of the complete initial training set. It is worth mentioning that not all methods for feature extraction reduce dimensionality, some of them aim at improving the signal-to-noise ratio without affecting data dimension, whereas others still increase it towards a better class separability[START_REF] Bishop | Pattern Recognition and Machine Learning[END_REF].

Fig. 4 . 1 :

 41 Fig. 4.1:The distribution of data is reflected in their covariance matrix. Diagonal spread is captured by the covariance, while axis-aligned spread is captured by the variance. A) shows a white noise distribution, in this case the covariance matrix is diagonal (i.e., all non-diagonal elements are zero). B) and C) present a diagonal distribution whose correlation is reflected in the sign of the non-diagonal elements of the covariance matrix. D) and E) show different dispersions in relation to the main axis, which can also be observed from their covariance matrices.

. 8 )

 8 Substituting equation 4.8 into equation 4.3, the Mahalanobis distance quadratic form becomes

Fig. 4 . 2 :

 42 Fig. 4.2: Normal distribution. The blue curve shows the elliptical surface of constant probability density for a normal distribution within a two-dimensional space x = (x 1 , x 2 ), on which the density is exp(-1/2) of its value at x = µ. The major axes of the ellipse are defined by the eigenvectors u i of the covariance matrix, with corresponding eigenvalues λ_i.

  2} denotes the EEG single trials belonging to class k over C channels and T samples, with representing the transpose. Then, the means of the covariance matrices Σ 1 and Σ 2 are computed by averaging for each class the covariance matrices of the successive training trials. So that the spatial filters W can be computed by solving the following generalized eigenvalue decomposition problem [20]:

15 )

 15 figure 4.4-A and figure 4.4-B respectively. We can observe that the lateralization for each motor task is enhanced with respect to figures 4.3-A-B (consider the axes scale). Furthermore, if we look at the distribution of the transformed data (see 4.4-C), we observe that both groups are simultaneously uncorrelated.

Fig. 4 . 3 :

 43 Fig. 4.3: EEG trials during left hand and right hand motor imageries in electrodes C3 (A) and C4 (B).Figure (C) shows a 2D representation where it is possible to observe that both groups are normally distributed. The variance of each distribution exhibits the characteristics of the sensorimotor rhythms, and presents a smaller magnitude with respect to the dimension associated to the contralateral electrode, which reflects the desynchronization of the electrical activity during the corresponding task. However, even though both distributions present ERD lateralization, data are likely correlated.

  Fig. 4.3: EEG trials during left hand and right hand motor imageries in electrodes C3 (A) and C4 (B).Figure (C) shows a 2D representation where it is possible to observe that both groups are normally distributed. The variance of each distribution exhibits the characteristics of the sensorimotor rhythms, and presents a smaller magnitude with respect to the dimension associated to the contralateral electrode, which reflects the desynchronization of the electrical activity during the corresponding task. However, even though both distributions present ERD lateralization, data are likely correlated.

Fig. 4 . 4 :

 44 Fig. 4.4: EEG trials during left hand and right hand motor imageries in electrodes C3 and C4 after having applied the CSP algorithm, which generates uncorrelated data for both classes.

Fig. 4 . 5 :

 45 Fig. 4.5: Projected trials generated by the first two CSP pairs of spatial filters over the band-pass filtered EEG signals of two consecutive right hand and left hand motor imageries. From top to down, the first and last projections are associated to the first CSP filter, while the two in the middle are linked to the second pair. Note the simultaneous changes of the signals variance within each pair of projections.

Fig. 4 . 6 :

 46 Fig. 4.6: CSP features extracted from the changes in variance of the projected trials according to equation 4.15

Fig. 4 . 7 :

 47 Fig. 4.7: Analytic Common Spatial Pattern algorithm. The first column of images shows the band-pass filtered EEG signals over electrodes C3 and C4 during a one-second segment of left hand motor imagery (A) and right hand motor imagery (B). In the second column the corresponding projections after applying the standard CSP algorithm over the left hand (C) and right hand (D) segments are presented. The third column shows the projections by using the analytical version of the method, namely, the ACSP algorithm over the left hand (E) and right hand (F) segments.Finally, in the last column, the obtained features after applying equation 4.15 over the projected trials by using the ACSP algorithm (G) and by using the CSP method (H) are presented. Note that separation in the case of the ACSP approach is enhanced in comparison to the standard algorithm.

1 . 9 filters of 4 2 .

 1942 A filter bank comprising multiple band-pass filters, which decomposes the EEG signals into multiple frequency pass bands. Various configurations of the filter bank are effective, however, a good choice consists in covering the frequency range from 4 to 40 Hz by using Hz width, namely, 4-8, 8-12,..., 36-40 Hz. This selection allows producing stable filters' responses while covering the physiological range of frequencies that is of interest for EEG analysis. Spatial filtering in terms of equation 4.14 over the 9 frequency ranges of the EEG trials generated in the previous step. Y b,i = W b E b,i , where E b,i ∈ R c×t represents the single-trial EEG measurement from the bth band-pass filter of the ith trial; Y b,i ∈ R c×t denotes E b,i after spatial filtering, W b ∈ R c×c corresponds to the CSP filters coefficient matrix; c is the number of channels; t is the number of EEG samples per channel; and denotes the transpose operator. The m pairs of CSP features of the ith trial for the bth band-pass filtered EEG measurements are computed in terms of equation 4.15 in the following terms

2 w=1P

 2 (w|f j )log 2 P (w|f j ) j,i )log 2 p(w|f j,i ),(4.27)

2 w=1P

 2 (f j,i |w)P (w).(4.29)

  34) with I(c, f (x)) the mutual information of c and f (x), and F a function space. If we consider only linear transformations, equation 4.34 thus simplifies to w * = argmax w∈R N {I(c, w x} (4.35)

(4. 36 )

 36 If we consider p(x|c) = N (0, R(x|c)), from where it follows that p(x|c) = N (0, R(x|c)), and

Figure 4 .

 4 Figure 4.8 shows the classification space generated by the CSP by JAD algorithm for a 3-class problem involving EEG signals during left hand, right hand, and both hands motor imageries of a right-handed male subject that has participated in the development of the 4-class database

. 42 )Fig. 4 . 9 :

 4249 Fig. 4.9: Geometry of a linear discriminant function in two dimensions generated by LDA. The decision boundary is shown in red, which is perpendicular to vector w, and presents a separation from the origin that is controlled by the bias parameter w 0 . The signed orthogonal distance of a given point x from the decision boundary is obtained as y(x)/ w . Data points correspond with the features extracted by the CSP algorithm presented in section 4.1.3, which were obtained in order to discriminate between EEG signals during left hand and right hand motor imageries .

Fig. 4 . 10 :

 410 Fig. 4.10: Discriminant functions for multiple classes. Attempting to construct a K class discriminant from a set of two class discriminants leads to ambiguous regions, shown in the left side. This problem can be overcome by considering the solution expressed in equation 4.45, which leads to decision boundaries that are singly connected and convex.

(4. 47 )

 47 Let us consider the case of multivariate normal distributions, for which the term p(x|C k )p(C k ) is given by equation 4.2, and let us assume that all distributions have the same covariance matrix, then, equation 4.47 derives to g(x) = max k=1,...,K

. 53 )

 53 Equation 4.53 can be interpreted as the linear model representing the hyperplane of the decision boundary discussed in section 4.6.1, so that y(x) = w x + w 0 , with

Figure 4 .

 4 Figure 4.9 shows the decision boundary obtained through linear discriminant analysis for the features extracted by the CSP algorithm presented in section 4.1.3, which were obtained in order to discriminate between EEG signals during left hand and right hand motor imageries.

  4.11). The parameters w and b that hold this assumption are found during the training stage. The goal is

Fig. 4 . 11 :

 411 Fig. 4.11: Support Vectors Machines classification by illustrating the margin maximization..

4. 6 . 2 . 1 . 1 Algorithm 2

 62112 Classification in the Riemannian manifold Algorithm 2 (MDRM) summarizes the process to obtain the Riemannian means and the corresponding distance between them and an unknown element, which are subsequently used to assign the corresponding class label. For a K-class problem, K Riemannian means are computed (one mean per class), and the label of an unknown element is assigned according to the closest Riemannian distance between itself and the computed means. Minimum Distance to Riemannian Mean (MDRM) Input: A set of X i training trials of K different known classes. Input: A testing data set comprising EEG trials of unknown class. Input: J (k) , the set of indexes for trials belonging to the k th condition. Output: k, the estimated class of the trials comprised in the testing set. 1) Compute the covariance matrices of the training trials X i to obtain P i . 2) Compute the covariance matrices of the testing trials to obtain P. 3) for k = 1 to K do P (k) G = G(P i , i ∈ J (k) ), use equation 4.67. end for 4) k = argmin k δ R (P, P (k) G ), use equation 4.64. 5) return k.

Figure 4 .

 4 Figure 4.13 shows the features extracted by applying the MDRM algorithm over all data recorded from subject 3 of the 4-class, and by considering only motor imageries of the left and

Fig. 4 . 13 :

 413 Fig. 4.13: MDRM classification. During the training phase, Riemannian means are computed for the left hand and right hand conditions. Subsequently, during validation, class labels are assigned according to the shortest Riemannian distance separating the input covariance matrices from the aforementioned class means. Note that this example is just for illustration purposes and that the same data were used to train/validate the classification model. Data belong to subject 3 from the 4-class database. Signals were band-pass filtered in the frequency range from [8-30 Hz].

Figure 4 .

 4 Figure 4.14 shows the same data that were used in the example illustrated in Figure 4.13 of the previous section after applying the MDRM algorithm over the CSP projected trials covariance matrices. It can be observed that separation between classes has been considerably enhanced.

Fig. 4 . 14 :

 414 Fig. 4.14: MDRM classification in combination with the CSP algorithm. During the training phase, Riemannian means are computed over the CSP projected trials for the left hand and right hand conditions. Subsequently, during validation, class labels are assigned according to the shortest Riemannian distance separating the input covariance matrices computed over the CSP projected trials from the aforementioned class means. In contrast to the standard MDRM algorithm, distances are more separated. Note that this example is just for illustration purposes and that the same data were used to train/validate the classification model. Data belong to subject 3 from the 4-class database. Signals were band-pass filtered in the frequency range from [8-30 Hz].

Contents

5. 1

 1 Multiclass approaches . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80 5.1.1 One-versus-one approach . . . . . . . . . . . . . . . . . . . . . . . . . . 81 5.1.2 One-versus-all approach . . . . . . . . . . . . . . . . . . . . . . . . . . . 81 5.1.3 Hierarchical approach . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82 5.2 Multilabel approaches . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83 5.2.1 One-step Multilabel (OsM) approach . . . . . . . . . . . . . . . . . . . . 85 5.2.2 Hierarchical Multilabel (HM) approach . . . . . . . . . . . . . . . . . . 86 5.2.3 One-step Hierarchichal Multilabel (OsHM) approach . . . . . . . . . . . 89

Figure 5 .

 5 Figure 5.1 shows a representation of the one-versus-one approach for a classification problem involving 4 different classes. To solve such a problem K(K -1)/2 = 6 models are required.

Fig. 5 . 1 :

 51 Fig. 5.1: One-versus-one classification for a 4-class problem. K(K -1)/2 different models, where K represents the number of classes, are generated to separate all possible pairs of groups. Finally, labels of the validation set are assigned according to the class having the highest number of "votes".

  sets are unbalanced. For instance, if we have four classes each with equal numbers of training data elements, then the individual classifiers are trained on data sets comprising 75% of samples belonging to the group containing the K -1 classes, and only 25% samples of the group involving only the C k class, which affects the original symmetry of the task.

Figure 5 .

 5 Figure 5.2 shows a representation of the one-versus-all approach for a classification problem involving 4 different classes. To solve such a problem only K = 4 models are required.

Fig. 5 . 2 :

 52 Fig. 5.2: One-versus-all classification for a 4-class problem. Consists in constructing K separate models, where K represents the number of classes, in which the k th model is generated by grouping within one group all elements belonging to class C k and, in the second class, all elements from the remaining K -1 classes.

Figure 5 .

 5 Figure 5.3 shows a representation of the hierarchical approach for a classification problem involving 4 different classes. As it can be observed, two stages are required; the first one consists of a model that selects 2 of the classes, which are subsequently evaluated by another model that is applied in the second and last stage in order to make the final decision.

Fig. 5 . 3 :

 53 Fig. 5.3: Hierarchical approach for a 4-class problem. The first stage of the classification is shown on top, where one group corresponds to the ensemble of classes A and D, and the other group to the ensemble of classes B and C. The second stage of the process is shown on the bottom; during which, according to the predicted label assigned in the previous stage, one of two models is selected for inferring the final decision.

Fig. 5 . 5 :

 55 Fig. 5.5: Architecture of the One-step Multilabel Approach (OsM) for solving the 4-class problem. Two different binary problems are generated around each activity source, where one class corresponds to all motor imageries involving the use of the associated limb (denoted by C ERD ), and the other class to those motor imageries that do not involve it (denoted by C IDLE ). The resulting features are concatenated and used to train/validate a single classification model.

Fig. 5 . 6 :

 56 Fig. 5.6: Architecture of the One-step Multilabel Approach (OsM) for solving the 8-class problem. Three different binary problems are generated around each activity source, where one class corresponds to all motor imageries involving the use of the associated limb (denoted by C ERD ), and the other class to those motor imageries that do not involve it (denoted by C IDLE ). The resulting features are concatenated and used to train/validate a single classification model.

Fig. 5 . 7 :

 57 Fig. 5.7: Architecture of the Multilabel Hierarchical approach (MH) for solving the 4-class problem. During the training phase three feature extraction and classification models are imple-mented; one model for the first stage, which is trained to discriminate between the two motor tasks involving the use of the selected hand (denoted by C ERD ) from the two motor tasks that do not involve it (denoted by C IDLE ), and two binary models in the second stage which are trained to discriminate between the motor tasks involving the hand that has not been analyzed (denoted by C ERD ) from those that do not involve it (denoted by C IDLE ), given that the hand selected for the first stage has been already recognized as active or inactive in terms of motor imagery-related activity. During the evaluation phase only two models are needed, which are selected according to the prediction of the previous stages. Note that this scheme illustrates the case where data partition is carried out by analysing first the engagement of the left hand.

Fig. 5 . 8 :

 58 Fig. 5.8: Architecture of the Multilabel Hierarchical approach (MH) for solving the 8-class problem. During the training phase seven feature extraction and classification models are implemented;one model for the first stage, which is trained to discriminate between the four motor tasks involving the use of the selected hand (denoted by C ERD ) from the four motor tasks that do not involve it (denoted by C IDLE ), two binary models in the second stage which are trained to discriminate between the motor tasks involving the hand that has not been analyzed (C ERD ) from those that do not involve it (C IDLE ), given that the hand selected for the first stage has been already recognized as active or inactive in terms of motor imagery-related activity, and finally, four binary models in the third stage to determine whether or not the feet are involved in the motor task, given that the hands have been already identified as active or inactive. During the evaluation phase only three models are needed, which are selected according to the prediction of the previous stages. Note that this scheme illustrates the case where data partition is carried out by analysing first the engagement of the left hand.

Fig. 5 . 9 :

 59 Fig. 5.9: Architecture of the One-step Multilabel Approach (OsM) for formulated for the 8-class database. Three sets of 4 CSP filters are computed from the band-pass filtered EEG signals of the training dataset. Each one of them is obtained by solving a binary problem at each one of the three main activity sources associated to the body parts considered by the paradigm. All MIs inducing ERD modulation are gathered together in one class denoted as C ERD and, all MIs that do not show ERD activity and present similar patterns to the idle state are grouped in the so-called C ERS . The 12 resulting features for each trial are concatenated and used to train an LDA classifier to solve the 8-class problem involving all motor tasks. During the evaluation phase, the unknown elements are band-pass filtered and projected into the three feature spaces to subsequently concatenate the resulting features and classify them according to the LDA model.

Fig. 6 . 1 :

 61 Fig. 6.1: 4-fold cross-validation process for assessing performance over the 4-class database for the discussed feature extraction methods in combination with the LDA and QDA algorithms, as well as for the distance-based approaches. Throughout each fold, the testing set (represented in dark blue) was constituted by all trials belonging to one of the runs (20 trials per class), and the training set (represented in dark blue) was integrated by the remaining 3 runs (60 trials per class). After all folds, all possible choices are implemented and the performance scores for the 4 folds are averaged.

Fig. 6 . 2 :

 62 Fig. 6.2: Nested cross-validation process for assessing performance over the 4-class database for the discussed feature extraction methods in combination with the SVM algorithm. Performance is evaluated through a 4-fold cross-validation. Given that the SVM algorithm requires a series of complexity parameters, another 4-fold cross validation is performed over the training set generated at each iteration, for which a range of different values are evaluated over the classification model in order to select the optimal ones. When more than one parameter need to be estimated, then the corresponding range is evaluated for each combination of values. Once such values have been determined, the model is implemented by using the original training set and evaluated on the original testing set.

Fig. 6 . 3 :

 63 Fig. 6.3: 30-fold cross-validation process for assessing performance over the 8-class database for all discussed methods in combination with the LDA and QDA algorithms, as well as for the distance-based approaches. Throughout each fold, the testing set (represented in dark blue) was constituted by one trial consecutively selected from each class (8 trials in total), and the remaining 29 trials per class (232 trials in total) as the training set. After all folds, all possible choices are implemented and the performance scores for the 30 folds are averaged.

Fig. 6 . 5 :

 65 Fig. 6.5: Comparison of the mean accuracy results obtained by using the proposed methods in combination with the CSP algorithm for all subjects of the 8-class database. Error bars indicate the Standard Error of the Mean.

Fig. 6 . 6 :

 66 Fig. 6.6: Comparison of the mean accuracy results obtained by using the proposed methods in combination with the ACSP algorithm for all subjects of the 8-class database. Error bars indicate the Standard Error of the Mean.

Fig. 6 . 7 :

 67 Fig. 6.7: Comparison of the accuracy results obtained by using the proposed methods and the one-versus-one and one-versus-all approaches in combination with the ACSP algorithm for all subjects of the 8-class database. Note that classification was carried out by Linear Discriminant Analysis. Error bars indicate the Standard Deviation.

Fig. 6 . 8 :

 68 Fig. 6.8: Comparison of the mean accuracy results obtained by using the proposed methods in combination with the ACSP algorithm for all subjects of the 8-class database. Error bars indicate the Standard Error of the Mean.

Fig. 6 . 9 :

 69 Fig. 6.9: Comparison of the mean accuracy results obtained by using the proposed methods in combination with the ACSP algorithm for all subjects of the 8-class database. Error bars indicate the Standard Error of the Mean.

Fig. 11 :

 11 Fig.11: First two features obtained with the CSP method for each one of the binary problems formulated by the OsM method for subject 3. Each row corresponds with the CSP features generated for each one of the 3 problems associated, from top to bottom, with the left hand, feet and right hand sources. From left to right the elements belonging to each one of the 8 different classes are highlighted in red if the corresponding source is not involved in the motor task, and in blue otherwise. Rst; LH; Ft; RH and BHs correspond, respectively, with rest; left hand; feet; right hand and both hands.

Fig. 12 :

 12 Fig.12: Most discriminative CSP patterns obtained for subject 3 for each one of the binary problems formulated by the HM method starting with the left hand. From top to bottom the CSP patterns generated for the left hand (A); the right hand given that there was no ERD activity found in the source associated with the left hand (B); the right hand given that the source associated with the left hand was found to generate ERD activity (C), and the four instances formulated by the hierarchical approach involving feet: both hands versus both hands and feet (D); left hand versus left hand and feet (E); right hand versus right hand and feet (F), and rest versus feet (G).

Fig. 13 :

 13 Fig. 13: Most discriminative CSP patterns obtained for subject 3 for each one of the binary problems formulated by the MH method starting with the right hand. From top to bottom the first patterns generated for the right hand (A); the left hand given that there was no ERD activity found in the source associated with the right hand (B); the left hand given that the source associated with the right hand was found to generate ERD activity (C), and the four instances formulated by the hierarchical approach involving feet: both hands versus both hands and feet (D); right hand versus right hand and feet (E); left hand versus left hand and feet (F), and rest versus feet (G).

Fig. 14 :

 14 Fig.14: Distribution of the 2 first CSP features generated by the MH approach on data from subject 3. From top to bottom, the three stages of the hierarchical classification process are shown by starting the process considering the presence/absence of ERD activity related to the use of the right hand. (A) The eight motor tasks of the problem are considered and gathered into two groups: motor tasks that involve the use of the right hand (blue), and motor tasks that do not involve it (red). For the second stage, two models are trained to determine the presence/absence of ERD activity related to the use of the left hand, given the prediction from the previous stage: (B) Motor tasks involving the right hand that also involve the use of the left hand (blue), and motor tasks involving the use of the right hand that do not involve the left hand (red). (C) Motor tasks that do not involve the use of the right hand but that involve the use of the left hand (blue), and motor tasks that do not involve the use of any hand (red). For the third stage there are only two classes remaining and four possibilities given the previous decisions: (D) Rest (red) versus feet (blue); (E) left hand (red) versus left hand & feet (blue); (F) right hand (red) versus right hand & feet, and (G) both hands (blue) versus both hands & feet (red).

Fig. 15 :

 15 Fig. 15: Distribution of the two first CSP features generated by the OsHM approach for the four instances involving the feet motor task applied on data from subject 3. From top to bottom, the binary problems: both hands-versus-both-hands & feet (A); left hand-versus-left hand & feet (B); right hand-versus-right hand & feet (C), and rest-versus-feet (D). The classes that were used to train the corresponding model are highlighted in blue (for MIs involving the source associated with the feet) and red (for MIs that do not involve the use of the feet). Data that are evaluated by a model that did not consider the corresponding class during the training phase are coloured in gray.

Fig. 16 :

 16 Fig. 16: Confusion matrices obtained by the one-versus-one and One-step Multilabel approaches for subject 1. Note that the surface of each element represents the sum of predictions generating the same label for the trials of each class.

Fig. 17 :

 17 Fig. 17: Confusion matrices obtained by the one-versus-one and One-step Multilabel approaches for subject 2. Note that the surface of each element represents the sum of predictions generating the same label for the trials of each class.

Fig. 18 :

 18 Fig. 18: Confusion matrices obtained by the one-versus-one and One-step Multilabel approaches for subject 3. Note that the surface of each element represents the sum of predictions generating the same label for the trials of each class.

Fig. 19 :

 19 Fig. 19: Confusion matrices obtained by the one-versus-one and One-step Multilabel approaches for subject 4. Note that the surface of each element represents the sum of predictions generating the same label for the trials of each class.

Fig. 20 :

 20 Fig. 20: Confusion matrices obtained by the one-versus-one and One-step Multilabel approaches for subject 5. Note that the surface of each element represents the sum of predictions generating the same label for the trials of each class.

Fig. 21 :

 21 Fig. 21: Confusion matrices obtained by the one-versus-one and One-step Multilabel approaches for subject 6. Note that the surface of each element represents the sum of predictions generating the same label for the trials of each class.

Fig. 22 :

 22 Fig. 22: Confusion matrices obtained by the one-versus-one and One-step Multilabel approaches for subject 7. Note that the surface of each element represents the sum of predictions generating the same label for the trials of each class.

Fig. 23 :

 23 Fig. 23: Confusion matrices obtained by the one-versus-one and One-step Multilabel approaches for subject 8. Note that the surface of each element represents the sum of predictions generating the same label for the trials of each class.

Fig. 24 :

 24 Fig. 24: Confusion matrices obtained by the one-versus-one and One-step Multilabel approaches for subject 9. Note that the surface of each element represents the sum of predictions generating the same label for the trials of each class.

Fig. 25 :

 25 Fig. 25: Confusion matrices obtained by the one-versus-one and One-step Multilabel approaches for subject 10. Note that the surface of each element represents the sum of predictions generating the same label for the trials of each class.

Fig. 26 :

 26 Fig. 26: Confusion matrices obtained by the one-versus-one and One-step Multilabel approaches for subject 11. Note that the surface of each element represents the sum of predictions generating the same label for the trials of each class.

L

  'importance de cette question découle de l'intérêt porté à optimiser l'algorithme Common Spatial Pattern (CSP) pour l'appliquer à des tâches de classification multi-label. Cette méthode est un algorithme très pratique et beaucoup utilisé pour extraire des caractéristiques d'imagerie motrice à partir de signaux EEG. Il génère une série de filtres spatiaux qui sont appliqués pour décomposer des données multidimensionnelles en un ensemble de composantes non corrélés[START_REF] Blankertz | Optimizing spatial filters for robust EEG single-trial analysis [revealing tricks of the trade[END_REF], qui visent à extraire des éléments qui maximisent simultanément la variance d'une classe tout en minimisant la variance de l'autre. En raison de la manière dont l'algorithme CSP réalise un tel contraste, cette approche présente l'inconvénient de ne convenir que pour discriminer entre deux classes et, pour résoudre des problèmes impliquant de nombreuses conditions, des extensions multi-classes sont nécessaires[START_REF] Gomez-Gil | Brain computer interfaces -a review[END_REF],[START_REF] Dornhege | Boosting bit rates in noninvasive eeg single-trial classifications by feature combination and multiclass paradigms[END_REF]. Les extensions classiques telles que les approches one-versus-one et one-versus-all présentent plusieurs inconvénients. Peut-être le plus important est la quantité de modèles qui sont nécessaires à mesure que le nombre de classes augmente, ce qui a un impact négatif sur le temps de calcul et sur la fiabilité des résultats. En ce sens, si nous pouvons confirmer que l'imagerie motrice d'une combinaison de deux ou plusieurs parties du corps impliquées en même temps correspond à une simple superposition des membres individuels à l'échelle EEG, alors il sera possible de concevoir des approches plus simples pour adapter non seulement la méthode CSP mais n'importe quelle autre méthode, puisque selon ce schéma, il deviendra possible d'utiliser un ensemble de classificateurs binaires afin de déterminer, pour chaque source d'activité, si la partie corporelle qui lui est associée est ou non engagée dans la tâche et donc inférer la décision finale à partir d'une combinaison de tous les sous-problèmes[START_REF] Lindig-León | Comparison of sensorimotor rhythms in EEG signals during simple and combined motor imageries over the contra and ipsilateral hemispheres[END_REF],[START_REF] Lindig-León | A multilabel classification method for detection of combined motor imageries[END_REF].Dans la présente étude, nous fournissons quelques éléments de preuve suggérant que cette hypothèse est plausible, sur laquelle nous avons développé trois nouvelles méthodes d'extraction de caractéristiques pour discriminer les imageries motrices combinées, à savoir l'approche multilabel à une étape (OsM), l'approche hiérarchique multi-label (HM) et l'approche hiérarchique multi-calques (OsHM) en une étape, qui à notre connaissance n'ont pas été considérées auparavant.La motivation était de permettre le contrôle 3D complet d'un bras robotique commercial qui offre aux utilisateurs quatorze mouvements prédéfinis. Le paradigme proposé implique l'utilisation de trois parties du corps, (à savoir: main droite, main gauche et pieds), qui en combinaison avec la condition de repos, constitue un problème à huit classes. Pour le mettre en oeuvre, nous avons utilisé une commande qui permet de changer de mode, afin de pouvoir associer différentes actions à la même commande. Six sujets en bonne santé ont participé à une séance d'enregistrement comprenant trois séries de dix essais par classe, soit 30 essais par classe. Les résultats sont légèrement supérieurs à ceux obtenus par les approches classiques one-versus-one et one-versus-all, ce qui est particulièrement intéressant si l'on considère que, pour résoudre le problème de classe 8 susmentionné, les deux méthodes exigent respectivement 28 et 8 modèles d'extraction de caractéristiques, ce qui est beaucoup plus que les 3, 7 et 6 modèles requis par les approches OsM, HM et OsHM respectivement. Cette thèse est structurée de la manière suivante: Le premier chapitre est consacré à fournir une base générale des principes et des pratiques des systèmes d'interface cerveau-ordinateur (BCI). Ces aspects seront approfondis le long du deuxième chapitre dans le contexte de l'imagerie motrice. Au cours de ces chapitres, nous présenterons la base de données à 4 classes, qui comprend une série d'enregistrements EEG de cinq sujets en bonne santé qui ont participé à une expérience d'imagerie motrice impliquant l'utilisation de la main droite, de la main gauche et des deux mains ensemble; ce qui, en plus de la condition de repos, conduit à un paradigme de quatre classes différentes. La base de données à 4 classes a été créée pour réaliser une étude préliminaire afin d'évaluer les similitudes entre l'activité obtenue sur les électrodes associées à l'utilisation de la main gauche et de la main droite lors d'imagerie motrice simple et combinée. De cette façon, en proposant un petit nombre de classes, il est possible d'enregistrer une quantité considérable d'essais sans provoquer la fatigue des utilisateurs. Sur ces données nous présenterons une série de résultats statistiques basés sur l'analyse du spectre de puissance observé sur chaque électrode au cours des quatre tâches motrices différentes, ce qui étaye l'hypothèse que les données EEG enregistrées au cours d'une imagerie motrice combinée peuvent être caractérisées comme la superposition de l'activité générée indépendamment par chacune des sources associées aux parties du corps engagées dans la tâche. Par la suite, dans le troisième chapitre, nous allons décrire les caractéristiques du bras robotique qui a motivé les solutions pratiques développées tout au long de ce projet. Dans la même ligne, nous allons décrire le schéma d'opération que nous avons conçu, qui consiste en une implémentation intelligente basée sur un contrôle à deux modes, qui permet de gérer différentes actions en utilisant la même commande. Comme mentionné, ce dispositif offre aux utilisateurs quatorze mouvements prédéfinis, ce qui exige beaucoup plus de commandes que celles obtenues avec la base de données à 4 classes. Dans cette optique, nous avons généré une base de données à 8 classes, qui comprend des signaux EEG de 6 sujets sains produisant les huit états mentaux qu'il est possible de générer par les tâches motrices impliquant l'utilisation de la main droite, de la main gauche, des deux pieds ensemble, et le repos. De la même manière que pour la base de données à 4 classes, nous présenterons une analyse qualitative des modèles ERD qui se trouvent dans les signaux EEG au cours de huit tâches différentes, ainsi qu'une comparaison statistique du spectre de puissance sur chaque électrode pour chaque couple de conditions. Ces résultats constituent la base sur laquelle nous justifierons le regroupement des classes qui a été suggéré pour les méthodes d'extraction de caractéristiques apportées par le présent travail. Dans le quatrième chapitre, le concept d'extraction de caractéristiques est introduit et élaboré dans le contexte des signaux EEG au cours de l'imagerie motrice. A travers ces pages, l'algorithme CSP sera élucidé avec deux de ses variantes, à savoir les modèles spatiaux communs analytiques (ACSP) et les modèles spatiaux communs par Diagonisation Approximative Conjointe (CSP by JAD). De plus, nous présenterons une stratégie appelée Filter Bank Common Spatial Pattern (FBCSP) qui optimise la performance de l'algorithme en décomposant les signaux EEG en une série de bandes de fréquences, ce qui favorise la sélection d'une gamme de fréquences discriminative spécifique au sujet. Par la suite, nous présenterons la notion de classification, pour laquelle deux algorithmes couramment utilisés seront décrits, à savoir l'Analyse discriminante linéaire (LDA) et les Machines à Vecteur de Support (SVM). De plus, afin d'explorer différentessolutions, le présent travail considère aussi l'utilisation de la géométrie riemannienne, avec laquelle la classification peut être obtenue directement dans l'ensemble des matrices symétriques définies positives qui est donné par les matrices de covariance des signaux EEG. Ensuite, dans la dernière section de ce chapitre, nous introduirons une nouvelle méthode de classification des systèmes BCI basée sur l'imagerie motrice, qui utilise la combinaison de l'algorithme CSP et de la géométrie riemannienne, ce qui représente l'une des contributions du présent travail. Dans le cinquième chapitre, le concept d'approche multiclasse sera introduit en décrivant trois méthodes classiques, à savoir les algorithmes un-contre-un, un-contre-tout et hiérarchique. À ce stade, nous décrirons les trois nouvelles méthodes multilabel que nous avons mises en oeuvre, à savoir l'approche multi-label à une seule étape (OsM), l'approche hiérarchique multi-label (HM) et l'approche hiérarchique monostable (OsHM) à une étape. Chacune de ces nouvelles méthodes est fondée sur les résultats qui ont été présentés dans les chapitres 2 et 3, à partir desquels nous avons proposé une méthode pratique pour effectuer la partition des données, de sorte que tout le problème est réduit à une série de tâches de classification binaires. Enfin, au long du sixième chapitre, nous présenterons les résultats, en termes de précision, obtenus par chacune des méthodes décrites sur les deux bases de données. La discussion correspondante est présentée dans le chapitre 7 et la conclusion générale au chapitre 8.

Fig. B. 2 :

 2 Fig. B.2: Oscillatory power for subject 1 within the alpha and beta bands.

Fig. B. 4 :

 4 Fig. B.4: Oscillatory power for subject 2 within the alpha and beta bands.

Fig. B. 6 :

 6 Fig. B.6: Oscillatory power for subject 3 within the alpha and beta bands.

Fig. B. 8 :

 8 Fig. B.8: Oscillatory power for subject 4 within the alpha and beta bands.

Fig. B. 10 :

 10 Fig. B.10: Oscillatory power for subject 5 within the alpha and beta bands.

Fig. B. 12 :

 12 Fig. B.12: Oscillatory power for subject 6 within the alpha and beta bands.

Fig. B. 14 :

 14 Fig. B.14: Oscillatory power for subject 7 within the alpha and beta bands.

Fig. B. 16 :

 16 Fig. B.16: Oscillatory power for subject 8 within the alpha and beta bands.

Fig. B. 18 :

 18 Fig. B.18: Oscillatory power for subject 9 within the alpha and beta bands.

Fig. B. 20 :

 20 Fig. B.20: Oscillatory power for subject 10 within the alpha and beta bands.

Fig. B. 22 :

 22 Fig. B.22: Oscillatory power for subject 11 within the alpha and beta bands.

Tables C. 21

 21 and C.22 show the accuracy results obtained over the 4-class and the 8-class databases respectively.

Fig. D. 1 :

 1 Fig. D.1: 3-axis operation mode.The robotic arm affords a full 3D control through fourteen predefined motions; including six reaching movements along the three axis (A), six rotations around the three axis (B), and a grasping function (C). These functions are distributed in three main articulations-like that resemble the forearm, the wrist, and the fingers. Note that we have decided to use all fingers together, so that the number of required commands decreases without loss of functionality, since the grasping action by using two or three of the fingers has almost the same effect.

Fig. D. 2 :

 2 Fig. D.2: Arm control. The different MIs that users were asked to perform are shown in the first column with their associated cue stimuli in the second one. The last three columns present the different robotic movements associated to each mode. Note that the rest condition is not linked to any action, and that motor imagery including the combination of all body parts represents a switch to change from one mode to another.

Table 2 .

 2

1 shows the p-values obtained from the unpaired t-test after each comparison between each pair of conditions over electrode C3. They are marked as < 0.01 (very significant); < 0.05 (significant) or the true value for > 0.05. In table D.1 the corresponding p-values for the same comparisons over electrode C4 are presented. It can be observed that, for almost all tests and

Table 2 .

 2

		1: Unpaired T-test to evaluate significance level at 5% of the different classes at
				electrode C3
	Subject RST-RH RST-LH RST-BHs RH-LH RH-BHs LH-BHs
	1	p < 0.05	p < 0.05	p < 0.01	p = 0.75 p = 0.27	p < 0.13
	2	p = 0.20	p = 0.25	p = 0.81	p < 0.01 p = 0.15	p = 0.40
	3	p < 0.01	p = 0.80	p < 0.01	p < 0.01 p < 0.01	p < 0.01
	4	p < 0.01	p = 0.80	p < 0.01	p < 0.01 p = 0.80	p < 0.01
	5	p < 0.01	p < 0.01	p < 0.01	p = 0.26 p = 0.56	p < 0.05

Table 2 . 2 :

 22 Unpaired T-test to evaluate significance level at 5% of the different classes at electrode C4

	Subject RST-RH RST-LH RST-BHs RH-LH RH-BHs LH-BHs
	1	p = 0.17	p < 0.01	p < 0.01	p < 0.05 p < 0.01	p = 0.37
	2	p = 0.48	p = 0.33	p = 0.55	p = 0.07 p = 0.15	p = 0.65
	3	p = 0.20	p = 0.12	p < 0.01	p < 0.01 p < 0.01	p < 0.01
	4	p = 0.39	p < 0.01	p < 0.01	p < 0.01 p < 0.01	p = 0.26
	5	p = 0.10	p < 0.01	p < 0.01	p < 0.01 p < 0.05	p = 0.71
	Figure 2.15 shows a box plot comparing the power spectrum for subject 4 at 12 Hz of each
		C3			C4
	A)			B)

motor imagery over electrodes C3 (fig.

2

.15-A) and C4 (fig.

2

.15-B). It can be noticed that there is no significant difference between the rest condition and the simple motor imageries at their ipsilateral hemisphere, nor is there for the combined task and simple motor imageries at their contralateral hemisphere.

15 20 25 30 35 40 45 50 55 Rest Right hand Left hand Both hands [dB] 20 25 30 35 40 45 50 Rest Right hand Left hand Both hands [dB]

  Fig.2.15: Box plots for subject 4 comparing the power spectrum magnitude at 12 Hz among all the considered conditions over electrodes C3 (A) and C4 (B). EEG signals were taken within the time interval from [1-2 s]. The distribution of data shows that there is no significant difference between simple hand MIs at their contralateral hemisphere with the combined task in the corresponding side, nor is there between simple hand MIs in the ipsilateral hemisphere with the rest condition.

	The statistical results presented in tables 2.1 and D.1, together with figure 2.15, are consistent

  Both the signal acquisition and the stimulation process were implemented on the OpenViBE platform 1[START_REF] Renard | Open-vibe: an opensource software platform to design, test and use brain-computer inter-faces in real and virtual environments[END_REF]. The EEG cap was fitted with 26 electrodes, namely, Fp1; Fp z ; Fp2; F z ; FC5;

	FC3; FC1; FC z ; FC2; FC4; FC6; C5; C3; C1; C z ; C2; C4; C6; CP5; CP3;CP1; CP z ; CP2; CP4;
	order to justify the grouping of classes, upon which the feature extraction methods contributed
	by the present work were implemented.
	3.1 8-class database

Contents 3.1 8-class database . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 3.1.1 Paradigm and time scheme . . . . . . . . . . . . . . . . . . . . . . . . . 38 3.1.2 Oscillatory power analysis . . . . . . . . . . . . . . . . . . . . . . . . . . 38 3.1.3 Statistical analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40

As will be pointed out below, the four different brain states that are generated by a paradigm such as the one described in the previous chapter (see section 2.6.1), are insufficient to afford an appropriate control for this kind of applications. Because of this reason, throughout the first section of this chapter we will present the 8-class database; which comprises EEG signals from 11 healthy subjects inducing the eight mental states that it is possible to generate by the motor tasks involving the use of the right hand, left hand, both feet together, and rest. Similarly as with the 4-class database, we will present a qualitative analysis on the ERD patterns that are found in the EEG signals during the eight different tasks. Furthermore, a statistical comparison of the power spectrum observed at each electrode among each pair of conditions will be presented in EEG signals were registered at 256 Hz from 11 right-handed healthy subjects (7 men and 4 women, from 22 to 45 years old) using a commercial amplifier Refa developed by TMS International™.

  3, we can observe that all motor imageries involving the use of the right hand (i.e., right hand, right hand in combination with feet, both hands, and both hands in combination with feet) present considerably lower ERD/ERS% values than those that do not include it (i.e., rest, left hand, feet, and left hand in combination with feet).Similarly, Figure3.4 presents the time course of ERD/ERS% over electrode C z (i.e, the activity source associated with the use of the feet). Over this region it can also be observed that the motor imageries involving the use of the feet (i.e., feet, left hand in combination with feet, right hand in combination with feet, and both hands in combination with feet) present lower ERD/ERS% values than those that do no it include it (i.e., rest, left hand, right hand,

Table 6 . 1 :

 61 Classification results on 4-class database

	Accuracy
	SVM

Table 6 . 2 :

 62 Classification results on 8-class database

			Accuracy
	Method	LDA	QDA	MD
	CSP OVO	46.0±16.0%	45.0±15.7%	NA
	CSP OVA	42.5±15.5%	41.6±15.9%	NA
	ACSP OVO	47.0±16.4%	41.1±15.9%	NA
	ACSP OVA	43.6±16.7%	43.9±15.1%	NA
	FBCSP OVO	47.1±16.0%	45.8±16.1%	NA
	FBCSP OVA	40.4±16.3%	38.7±15.9%	NA
	FBACSP OVO	44.2±15.3%	37.6±15.5%	NA
	FBACSP OVA	36.8±15.2%	37.1±14.2%	NA
	CSPbyJAD	44.9±14.0%	38.9±16.6%	NA
	MDRM	NA	NA	40.3±13.7%
	CSP OVO /MDRM	NA	NA	41.5±15.1%%
	CSP OVA /MDRM	NA	NA	41.7±14.8%
	OsM CSP	47.7±16.2%	41.1±16.48%	NA
	OsM ACSP	49.3±16.5% 43.4±15.8%	NA
	OsM FBCSP	41.0±15.5%	37.4±15.1%	NA
	OsM FBACSP	39.8±14.9%	23.8±13.8%	NA
	OsM MDRM	NA	NA	34.1±12.9%
	OsM CSP/MDRM	NA	NA	37.5±12.9%%
	HM CSP /left	37.9±15.7%	36.8±14.7%	NA
	HM ACSP /left	43.2±16.4%	40.8±15.5%	NA
	HM FBCSP /left	42.5±16.1%	42.8±15.7%	NA
	HM FBACSP /left	42.3±15.9%	39.4±16.1%	NA
	HM MDRM /left	NA	NA	37.7±14.3%
	HM CSP/MDRM /left	NA	NA	41.5±15.3%
	HM CSP /right	43.5±16.6%	41.9±15.0%	NA
	HM ACSP /right	45.0±16.2%	41.1±16.0%	NA
	HM FBCSP /right	40.3±15.9%	37.8±15.4%	NA
	HM FBACSP /right	40.3±15.3%	37.8±15.4%	NA
	HM MDRM /right	NA	NA	35.1±14.6%
	HM CSP/MDRM /right	NA	NA	39.0±13.8%
	OsHM CSP	47.1±15.7%	43.83±17.3%	NA
	OsHM ACSP	46.7±16.3%	41.9±16.8%	NA
	OsHM FBCSP	44.1±14.4%	-	NA
	OsHM FBACSP	42.3±15.9%	-	NA
	OsHM MDRM	NA	NA	36.8±13.7%
	OsHM CSP/MDRM	NA	NA	37.0±15.1%

  classification model during the calibration phase, although results are unsatisfactory when it is applied to provide online feedback. Whereas, in other cases, it is not even possible to generate a classification model during the training stage. Besides the illiteracy problem, motor imagery based BCIs are susceptible to intra-session variability, which provokes inconstancy of data recorded from the same user throughout different dates as a result of possible changes in the electrodes location, and/or variations of the mental state of the subject. This variability can cause that the initial calibration is not longer accurate for subsequent settings, which encumbers the daily use of the system. Also in this regard, it should be noticed that performance in out-of-the-lab environments gets affected by additional noise sources, as well as by different types of distractions that may cause non-stationary features[START_REF] Brandl | Bringing bci into everyday life: motor imagery in a pseudo realistic environment[END_REF].

Finally, the most important challenge for adopting this kind of applications may lie in the establishment of effective training methods, since the existing protocols essentially provide theoretical considerations that do not stimulate the acquirement of appropriate skills to effectively modulate EEG patterns, which might be responsible for inadequate user performance and BCI illiteracy

[START_REF] Jeunet | Why standard brain-computer interface (BCI) training protocols should be changed: an experimental study[END_REF]

.

Table C .

 C 10: Accuracy results for the FBCSP method using one-versus-all classification on the 4-class database

		19±5.34%	47.19±8.06%	51.87±5.45%	52.50±3.68%	46.88±5.15%
	AVG	51.56±8.77%	49.31±7.91%	48.19±7.44%	48.00±7.28%	41.13±7.52%
			Classification method
					SVM
	Subject	LDA	QDA	Linear	Polynomial	RBF
	1	45.62±6.65%	47.19±8.00% 41.25±3.23% 41.56±3.29%	34.38±9.44%
	2	38.44±8.00%	38.75±3.54% 36.25±3.68% 34.38±3.89%	29.38±9.44%
	3	70.31±10.38% 63.44±18.99% 63.12±8.93% 62.81±11.24% 43.75±19.55%
	4	52.19±4.13%	54.37±3.89% 53.75±2.70% 52.19±3.13%	44.69±6.40%
	5	54.69±7.17%	53.75±5.10% 56.56±7.17% 48.13±6.81%	39.06±6.24%
	AVG	52.25±7.27%	51.50±7.90% 50.19±5.14% 47.81±5.67%	38.25±9.20%

4.6.2.2).

Finally, results obtained by the three contributed algorithms, namely, OsM approach, HM 52.

Table C .

 C [START_REF] Pfurtscheller | Event-related EEG/MEG synchronization and desynchronization: basic principles[END_REF]: Accuracy results for the FBCSP method using one-versus-one classification on the 8-class database

		58.75±10.75%
	4	40.62±2.92%
	5	43.75±8.35%
	AVG	45.37±6.47%
		Classification method
	Subject	LDA	QDA
	1	27.92±14.93% 29.17±16.19%
	2	41.25±15.45% 41.25±15.79%
	3	89.17±9.13%	89.17±8.42%
	4	50.42±20.36% 45.83±17.47%
	5	23.75±15.17% 22.08±14.93%
	6	46.25±16.79% 45.83±18.08%
	7	40.42±16.96% 40.00±14.46%
	8	48.33±19.07% 45.83±20.06%
	9	43.33±16.65% 38.33±18.26%
	10	49.17±15.72% 50.00±17.37%
	11	58.33±15.51% 56.25±16.33%
	AVG 47.12±15.98% 45.79±16.12%

Table C .

 C [START_REF] Blankertz | Optimizing spatial filters for robust EEG single-trial analysis [revealing tricks of the trade[END_REF]: Accuracy results for the MDRM method on the 8-class database As described in section 4.6.2.2, we have proposed a new classification approach that consists of the combination of the CSP method and the Minimum Distance to Riemannian Mean algorithm.

	Subject	Accuracy
	1	29.17±12.85%
	2	38.75±11.99%
	3	71.67±15.72%
	4	45.00±16.28%
	5	27.92±11.69%
	6	15.42±7.10%
	7	39.17±15.99%
	8	45.83±16.52%
	9	20.00±11.18%
	10	54.58±14.86%
	11	55.83±16.65%
	AVG	40.3±13.71%
	C.7 Results: CSP + MDRM algorithm

Table C .

 C 21: Accuracy results for the CSP algorithm applied in combination with the MDRM approach on the 4-class database

		19.58±9.67%
	10	47.92±17.40%
	11	45.83±15.51%
	AVG	35.15±14.58%
	Classification method

Table C .

 C , for the 8-class database, table C.57 presents the results obtained when the analysis starts on the activity generated by use of the left hand, and table C.58, when it starts with the right hand.TableC.55: Accuracy results for the HM method in combination with the CSP+MDRM approach on the 4-class database starting with the left hand 56: Accuracy results for the HM method in combination with the CSP+MDRM approach on the 4-class database starting with the right hand

	Subject	Accuracy
	1	50.00±7.36%
	2	37.19±4.13%
	3	54.69±16.94%
	4	54.06±5.98%
	5	54.38±6.17%
	AVG	50.06±8.12%
	Subject	Accuracy
	1	50.00±8.54%
	2	40.31±6.24%
	3	61.25±22.48%
	4	53.44±2.13%
	5	55.63±4.62%
	AVG	52.12±8.80%

Table C .

 C 57: Accuracy results for the HM method in combination with the CSP+MDRM approach on the 8-class database starting with the left hand

	Subject	Accuracy
	1	29.58±14.11%
	2	39.17±15.99%
	3	80.42±13.80%
	4	42.50±20.13%
	5	25.00±9.28%
	6	24.58±12.05%
	7	40.83±17.96%
	8	43.75±16.33%
	9	22.50±15.19%
	10	56.25±17.29%
	11	52.08±15.77%
	AVG	41.51±15.26%

Table C .

 C [START_REF] Pfurtscheller | Synergetics of rhythms in biological systems[END_REF]: Accuracy results for the HM method in combination with the CSP+MDRM approach on the 8-class database starting with the right hand .10.5 Results: OsHM method by using the MDRM algorithmSimilar to what we have discussed in section C.8.5, the implementation of the OsHM approach in combination with the MDRM algorithm admits different variants. However, we have decided to stay as close to the MDRM method, so that labels have been predicted from the combination of the shortest distances separating the input data from the means associated to each one of the 6 binary problems formulated by the OsHM approach. The corresponding accuracy results are shown in table C.63.Table C.63:Accuracy results for the OsHM method in combination with the MDRM approach on the 8-class database .10.6 Results: OsHM method by using the CSP + MDRM approach

		19.17±14.21%
	6	40.00±16.54%
	7	35.42±18.30%
	8	38.75±13.27%
	9	33.33±16.52%
	10	51.67±15.65%
	11	56.67±16.33%
	AVG 42.33±15.92%
	Subject	Accuracy Label
	1 Subject	32.08±13.80% Accuracy
	2 1 %	39.58±10.93%
	3 2	65.00±8.93% 37.92±10.63%
	4 3	47.50±18.10% 60.42±15.77%
	5 4	22.08±14.56% 35.42±15.07%
	6 5	25.42±14.49% 23.33±9.13%
	7 6	37.50±15.74% 24.17±12.69%
	8 7	38.33±13.51% 37.50±18.86%
	9 8	21.67±11.34% 35.42±16.11%
	10 9	54.17±18.38% 15.42±8.49%
	11 10	46.25±11.90% 48.75±14.81%
	AVG 11	39.05±13.79% 49.17±16.72%
	AVG 36.82±13.66%

CC

Table C .

 C [START_REF] Blankertz | Neurophysiological predictor of SMR-based bci performance[END_REF]: Accuracy results for the OsHM method in combination with the CSP+MDRM on the 8-class database

	Subject	Accuracy
	1	28.33±12.69%
	2	42.50±16.28%
	3	74.17±15.37%
	4	32.08±17.58%
	5	19.17±16.33%
	6	28.17±15.37%
	7	33.75±15.79%
	8	33.33±14.43%
	9	15.83±9.25%
	10	52.08±16.77%
	11	47.08±16.64%
	AVG	36.95±15.13%

Le nombre de commandes que permet un système BCI à base d'imagerie motrice dépend du nombre d'états cérébraux que le système est capable de détecter, quantité qui est limitée par le nombre de parties du corps avec lesquelles les utilisateurs sont capables de produire des modèles d'activité appropriés. De plus, la difficulté de la tâche augmente avec le nombre total de commandes disponibles. Augmenter le nombre de commandes n'est possible que si l'on peut discriminer de manière efficace un nombre correspondant d'états mentaux[START_REF] Dornhege | Increase information transfer rates in BCI by CSP extension to multi-class[END_REF]. Il existe deux types d'approches pour résoudre des problèmes multiclasses: D'une part, certaines techniques tentent de prédire le label en utilisant une méthode directe, comme les algorithmes basés sur la distance. D'autres part, certaines méthodes utilisent des extensions multiclasses afin de réduire le problème à un ensemble de décisions binaires qui peuvent être accomplies de différentes façons, par exemple, en utilisant les approches classiques one-versus-one ou one-versus-all[START_REF] Dornhege | Boosting bit rates in noninvasive eeg single-trial classifications by feature combination and multiclass paradigms[END_REF]. Une méthode multiclasse intéressante, qui permet aux utilisateurs de disposer d'un vaste répertoire de commandes différentes, est connue sous le nom d'approche multi-label. Cette approche est basée sur la détection d'imageries motrices combinées (c'est-à-dire impliquant deux ou plusieurs parties du corps à la fois)[START_REF] Weibo | Eeg feature comparison and classification of simple and compound limb motor imagery[END_REF],[START_REF] Royer | EEG) control of a virtual helicopter in 3-dimensional space using intelligent control strategies[END_REF]. En comparaison avec le schéma standard, cette approche présente l'avantage d'augmenter considérablement le nombre d'états différents du cerveau détectables, tout en utilisant le même nombre de parties du corps (de l'ordre de 2 P par rapport à P, où P est le nombre de parties du corps impliquées).Contrairement aux paradigmes multi-classes standard, auxquels laplupart des recherches ont été consacrées, moins de travail a été rapporté sur les modèles oscillatoires cérébraux induits par les imageries motrices combinées. La technique de décomposition en mode empirique combinée avec l'analyse de synchronisation de phase a été proposée et appliquée pour évaluer la synchronisation des zones cérébrales motrices supplémentaire et primaire lors d'une imagerie motrice complexe combinant l'action du corps et des membres. Les résultats de classification ont suggéré la faisabilité de cette approche pour des applications BCI de récupération motrice des membres inférieurs[START_REF] Zhou | A novel technique for phase synchrony measurement from the complex motor imaginary potential of combined body and limb action[END_REF].En outre, afin de surmonter la quantité limitée de commandes de sortie résultant du nombre limité de classes, et dans le but d'obtenir un contrôle 3D continu, l'imagerie motrice des deux mains a été adoptée afin d'augmenter le jeu d'instructions d'une interface BCI basée sur l'imagerie

Acknowledgments

6

Experimental results Along this chapter we will present the offline results generated after applying the feature extraction methods discussed in chapter 5 on the 4-class and 8-class databases. For classification, Linear Discriminant Analysis and Support Vector Machines were applied on the 4-class database in order to evaluate the performance of the proposed methods by exploring different solutions (excepting the Riemannian distance-based schema, which do not require a complementary classifier). Results obtained by the Support Vector Machines algorithm showed a considerably lower performance than the one generated by the Linear Discriminant Analysis method, presumably because the corresponding models were over-trained. Thus, for the 8-class database only results obtained by LDA are presented.

To start with, we specify the parameters that were used for signal preprocessing (i.e., filtering parameters and time window length), as well as the electrode montage that has been applied during difference between methods is significant, and that the OsM approach considerably outperforms the classical ones.

Similarly, Figure 6.7 shows the same comparison by using the analytical version of the CSP algorithm and LDA classification. Figure 6.8 shows the mean accuracy with the corresponding standard error of the mean (SEM) among all subjects for each one of the methods. And, in order to observe the performance of each method over the best subjects, Figure 6.9 present the mean accuracy together with the SEM among subject, 3, 10 and 11. Note that the difference between methods is significant, and that the OsM approach considerably outperforms the classical ones.

Discussion

Throughout this section we will analyze the characteristics of the three contributed methods together with the results presented in the previous chapter. Since all classes contained in the 4-class database are also included in the 8-class one, the following analysis will be focused only on the latter. Moreover, given that the best results were obtained by applying all methods in combination with the CSP algorithm, we will center the discussion around this solution.

To start with, we will analyze the CSP patterns [START_REF] Haufe | On the interpretation of weight vectors of linear models in multivariate neuroimaging[END_REF] and the CSP features generated by the OsM approach. Figure 10 shows a topographic representation of the first CSP patterns computed for each one of the three binary problems formulated by this method. Data were taken from subject 3 over the training set generated by one partition of the 30-fold cross validation described in section 6.2.2. It can be observed that there are strong weights over the motor cortex, more precisely over the sources that are associated to the body parts that are involved during the different motor tasks, as expected from the literature [START_REF] Pfurtscheller | Motor imagery and direct brain-computer communication[END_REF]; which confirms that such partition is neurophysiologically plausible. Fig. 10: Most discriminative CSP patterns obtained for subject 3 for each one of the binary problems formulated by the OsM approach. From left to right, the first CSP patterns generated for each binary problem considering the motor tasks inducing ERD activity from those that do not, at each one of the sources associated with the right hand (A); feet (B), and left hand (C). Note that the strongest weights are generated over the contralateral hemispheres.

For discussion purposes, Figure 11 presents the first two CSP features computed over the aforementioned training set by using the OsM approach. Note that these features correspond only to the first pair of filters and that some discriminative information may be missing. Also note that the proposed approach uses the concatenation of the three set of features by comprising them into a single vector. However, this representation is useful to observe that, in the case of the sources associated to the use of the left and right hands, features are well separated, on the contrary, features obtained from the problems involving the feet are overlapped for almost every The present work introduces a new framework for EEG motor imagery-based BCI systems that affords users with an unusual high number of commands, which facilitates a full 3D interaction that can be used in different applications. In this regard, we have implemented a system for controlling a robotic arm that provides fourteen predefined motions. Such a framework considers the use of combined motor tasks; which contrarily to standard single-label paradigms has not been extensively studied, and contributes with three different methods to solve the associated multiclass problem, namely, One-step Multilabel (OsM) approach, Hierarchical Multilabel (HM) approach, and One-step Hierarchical Multilabel (OsHM) approach, which to our best knowledge have not been considered before. These contributions are based on a particular way of grouping classes that has been inspired by some physiological considerations, such as the sensorimotor rhythms lateralization and the specific location of the activity sources along the sensorimotor cortex. This new scheme allows to transform a multilabel classification problem into a series of binary tasks, which optimizes the use of the Common Spatial Pattern (CSP) algorithm, given that it has been shown to be powerful at discriminating sensorimotor rhythms but has the inconvenience of being suitable for solving only 2-class problems.

Conclusions

We have used two databases in order to analyze the patterns of combined motor tasks, and to implement and validate the proposed methods, namely, the 4-class and the 8-class databases, which have been specifically recorded for this study. The former database involves only the use of the left and right hands; which together with the rest condition, enable four different brain states when both the single tasks and the combined one are considered (i.e., rest, left hand, right hand and both hands together). In the case of the 8-class database, we have also incorporated the use Appendices mouvements imaginés [START_REF] Avanzini | The dynamics of sensorimotor cortical oscillations during the observation of hand movements: an EEG study[END_REF]. Cette modulation prend deux formes différentes: une réduction de puissance oscillatoire au cours de la tâche motrice appelée désynchronisation événementielle (ERD), et une augmentation de puissance oscillatoire connue sous le nom de Event-Related Synchronization (ERS) qui apparaît juste après la fin de la tâche motrice, et qui dure quelques centaines de millisecondes [START_REF] Pfurtscheller | Event-related EEG/MEG synchronization and desynchronization: basic principles[END_REF]. Lors de la préparation d'une tâche motrice, certaines régions corticales et sous-corticales sont activées. A l'inverse, au cours de l'imagerie motrice, la plus grande partie de l'activité a lieu dans le cortex moteur principal, dans l'hémisphère contralatéral correspondant [START_REF] Pfurtscheller | On the existence of different types of central beta rhythms below 30 hz[END_REF]. On sait que la représentation topographique et le changement de puissance de la bande au cours de l'imagerie motrice sont localisés dans des régions spécifiques du cortex moteur primaire (en gros, les activités correspondant au côté gauche sont observées sur l'électrode C4, les activités correspondant au côté droit sur C3, les pieds sur C z , etc.) [START_REF] Pfurtscheller | Event-related EEG/MEG synchronization and desynchronization: basic principles[END_REF], [START_REF] Barachant | Multi-class brain computer interface classification by riemannian geometry[END_REF].

approach, and OsHM approach, in combination with all the aforementioned feature extraction methods are presented for both databases.

All feature extraction methods (excepting the Riemannian distance-based schema, which do not require a complementary classifier), where applied together with Linear Discriminant Analysis (LDA, see section 4.6.1.3) for classification. In the case of the 4-class database, results by using Support Vectors Machines (SVM, see section 4.6.1.4)are also presented.

C.1 Results: CSP method

In this section we present the offline results generated by applying the standard CSP algorithm discussed in section 4.1.3 over both databases and for all subjects. As mentioned, the number of the m selected pairs was set to m=2, which generates a total of four spatial filters for the corresponding feature extraction. Accuracy performance by using the discussed classification algorithms under the one-versusone approach is reported on tables C.

C.4 Results: FBACSP method

The analytical version of the FBCSP method discussed in section 4.4 has been also implemented as another contribution of the present work, to which we will refer as the Filter Bank Analytical

Common Spatial Pattern (FBACSP) algorithm. For its implementation, the same parameters of the standard FBCSP method have been selected, the main difference is that feature extraction is carried out according to the ASCP algorithm discussed in section 4.2. In this way, both feature extraction and classification, were carried out by following a one-versus-one approach (see tables C.13 and C.15 for the 4-class and 8-class databases respectively), and under a one-versus-all approach (see tables C.14 for the 4-class database and C.16 for the 8-class one). As mentioned,

given that the ACSP method generates complex features, it is necessary to use a special strategy for applying the SVM approach; therefore, we only present the results generated in combination with Linear Discriminant Analysis classification.

C.5 Results: CSP by JAD method

In this section we present the results generated after applying the CSP by JAD algorithms. As mentioned in section 4.5, feature extraction under this method is directly achieved by considering all classes at once, without having to resort to any complementary technique, such as the oneversus-one or one-versus-all approaches. Tables C.17 and C.18 show the accuracy results obtained after applying the algorithm over the the 4-class and 8-class databases respectively.

C.8 Results: OsM method

Throughout the next sections we will report the accuracy results obtained by the OsM method described in section 5.2.1, which corresponds to one of the main contribution of the present work.

As mentioned, this approach generates P binary problems, where P corresponds to the number of body parts included in the paradigm (2 problems for the 4-class database, and 3 problems for the 8-class one intended to determine whether each one of the considered limbs are engaged in the motor imagery or not).

In order to explore different solutions, we have applied this partition technique on both databases in combination with all the feature extraction algorithms that have been presented in chapter 4.

C.8.1 Results: OsM method by using the CSP algorithm

The first arrangement to be presented for the implementation of the OsM approach is in combination with the CSP algorithm. As mentioned 2 pairs of spatial filters were selected (4 in total).

The extracted features from each one of the P binary tasks are subsequently concatenated, which generates input vectors of 12 elements each.

Tables C. [START_REF] Mason | Temporal control paradigms for direct brain interfaces -rethinking the definition of asynchronous and synchronous[END_REF] and C.24 present the accuracy results that were obtained, respectively, for the 4-class and the 8-class databases after classifying the computed features. Here we present the results generated by the OsM method in combination with the ACSP algorithm. The process is similar to the one discussed in the previous section, the main difference is that in this case the analytical form of the CSP algorithm explained in section 4.2 is considered.

Given that the computed features are defined in the complex space, only results by using linear discriminant analysis classification are presented.

Tables C. [START_REF] Mak | Clinical applications of brain-computer interfaces: current state and future prospects[END_REF] and C.26 show the accuracy results after classifying the obtained features over the 4-class and the 8-class databases respectively. The results shown below were obtained by applying the same method described in the previous section, except that the CSP algorithm was considered in its analytical form. Given that the generated features are defined in the complex space, results in combination with SVM classification for the 4-class database are not reported.

Tables C. [START_REF] Graiman | Brain-computer interfaces[END_REF] and C.30 show the accuracies obtained for each subject by using linear discriminant analysis classification over the 4-class and the 8-class databases respectively. The last configuration applied under the OsM approach is the combination of the CSP algorithm with the MDRM method discussed in section 4.6.2.2. Stated briefly, it consists in the same procedure described in the preceding section with the only difference that the covariance matrices, which are used as the EEG descriptors, are computed over the trials projection obtained after applying the CSP algorithm, instead of over the band-pass filtered signal segments (as it is done by the standard MDRM). Note that the CSP filters are computed exactly as it was described in section C.8.1.

Tables C. [START_REF] Roy | On the regulation of the blood-supply of the brain[END_REF] and C.34 show the accuracy results obtained from the 4-class and the 8-class databases respectively.

C.9 Results: HM method

Throughout the next sections results obtained by the HM approach in combination with all the discussed feature extraction and classification algorithms will be presented. Note that hierarchical methods involve only binary problems that are consecutively applied (see sections 5.1.3 and 5.2.2), therefore, classification can be directly applied without using any additional grouping approach.

Also note that the HM method can start either by considering first the activity associated with the use of the left hand, or by considering the activity associated with the use of the right hand.

In order to evaluate both variants the two approaches have been implemented for all subjects.

C.9.1 Results: HM method by using the CSP algorithm

To start with, we report the performance of the HM approach in combination with the standard CSP method. Again, we have considered 2 pairs of spatial filters (4 in total) for each one of the binary problems formulated at each stage of the method (2 stages comprising 3 binary problems for the 4-class database, and 3 stages containing 7 binary problems for the 8-class database).

Tables C. [START_REF] Cohen | Magnetoencephalography: detection of the brain's electrical activity with a superconducting magnetometer[END_REF] Note that features extracted by the analytical version of the CSP method are defined in the complex space, therefore, classification with the SVM algorithm is not included. C.9.3 Results: HM method by using the FBCSP algorithm

In this section we report the performance of the HM method in combination with the FBCSP algorithm. For the implementation of all classification stages, 2 pairs of spatial filters (4 in total)

were selected from each one of the models generated over the 9 frequency bands comprised in the filter bank. Subsequently, from all features, only 4 elements were extracted by using the C.9.4 Results: HM method by using the FBACSP algorithm

In this section we present the accuracy results obtained by the HM approach in combination with the FBACSP algorithm. Its implementation was similar to the one described in the previous section with the only difference that the CSP algorithm was considered in its analytical form, therefore, given that the resulting features are complex, classification was carried out only by linear discriminant analysis. Finally, we report the performance of the OsHM approach, which corresponds to the third contributed method. As mentioned in section 5.2.3, this method is a combination of the OsM and HM approaches; it implements the same binary problems related to the sources associated to the hands, as it was done with the OsM method, and incorporates the four instances associated with the source linked to the feet that were formulated by the HM algorithm in the case of 8-class database. Each one of these problems is implemented over the whole set of classes and, subsequently, the resulting features are concatenated to create single vectors, from which classification can be directly achieved. Note the OsHM approach is equivalent to the OsM method when only the 4 classes comprised in the 4-class database are considered, therefore, we only report the results obtained over the 8-class database.

C.10.1 Results: OsHM method by using the CSP algorithm

In this section we present the performance of the OsHM method in combination with the CSP algorithm. Given that 2 pairs of spatial filters (4 in total) were considered for each one of the 6 models involved in the procedure, the resulting vectors contain 24 elements each. We have also applied the FBACSP algorithm to evaluate the performance of the OsHM approach.

Its implementation was similar to the one discussed in the previous section, the only difference is that instead of using the standard CSP algorithm for extracting features within each one of the 9 frequency bands, we have applied its analytical version. Results are shown in table C.62.

Note that only results obtained by linear discriminant analysis under the assumption of identical covariance matrices for the involved conditions are presented since, as in the previous section, in many cases their covariance matrices were also singular. In order to explore new solutions for overcoming the existing limitations of EEG-based neuroprostheses, this work has explored new alternatives to control devices requiring a high number of commands. To close such study, in the following section we present an implementation example over a commercial device that affords users with a full 3D control through fourteen predefined motions.

To start with, a description of the robotic arm and the operation scheme, which consists of a smart implementation based on a switching-mode control that has been designed to operate all its functions by using only the eight commands that are obtained trough the paradigm presented in section 3.1, are described.

D.1 Description of the robotic arm

The proposed approach has been motivated by a commercial device model JACO 3-fingers by Kinova™ 1 , which affords a full 3D control through fourteen predefined motions; including six reaching movements along the three axis, six rotations around the three axis, and a grasping function. These functions are distributed in three main articulations-like that resemble the forearm, with which the device is able of approaching objects, the wrist, with which it is possible to turn the hand towards any direction, and the fingers, which allows to open and close the hand (see [START_REF] Wolpaw | BCI: Principles and Practice[END_REF]. Regarding the use of the fingers, even when the device allows controlling each pair separately, we have decided to consider them all together. In this way the number of required commands decreases without loss of functionality, since the grasping action by using two or three of the fingers has almost the same effect.

An important aspect of this device, and probably of most of its kind, is that there is a limited range of movements once the robotic arm has reached an extreme position, so that even if the system detects a command demanding a motion that is beyond the possibilities of movement, the corresponding action will not be executed.

D.1.1 Implementation challenges

The appropriate use of this robotic arm implies several challenges that must be addressed from different perspectives. For instance, users' training is one of the most difficult aspects. Even when facilitates the interaction with the subject, but on the other hand demands a considerable number of different brain states to be properly induced and classified. In this regard there are at least two aspects to be considered; first, the number of body parts with which it is possible to offer an intuitive control and whose related sources are separated enough to distinguish the induced activity is limited, and second, the classification task becomes more difficult as the amount of classes increases, which can lead to misclassified responses that hinder the correct use of the device.

Because of these considerations we have decided to explore the use of combined motor imageries, since in contrast with the standard scheme, this approach has the advantage of considerably increasing the amount of afforded commands while using the same number of body parts (in order of 2 P compared to P , where P is the number of body parts). However, even when the use of combined movements is considered, there might be not enough commands to control devices as the one we are describing. For instance, in the present work we have decided to consider the use of the left hand, right hand and both feet together. Given that these body parts are intuitively associated with movement. It is worth mentioning that feet are considered together given the proximity of the two sources, which may be difficult to distinguish.

Under this scheme users are afforded only with eight commands if the rest condition is also considered (i.e., left hand, feet, left hand in combination with feet, right hand, both hands together; right hand in combination with feet; both hands together with feet, and rest), which are not sufficient to control the fourteen available motions. To overcome the insufficiency of commands we have designed a smart switching-mode scheme that allows controlling different actions by using the same command. Accuracy after the second run was 78.5% and 84.75% for the third run. Subsequently, an offline analysis was carried out in order to evaluate the system performance by using a 3-fold cross validation for which the test set was, at each fold, all trials belonging to one run. Table D.1

shows the obtained results.