A New PLV-Spatial Filtering to Improve the Classification Performance in BCI Systems

Objective: The performance of an EEG-based brain-computer interface (BCI) system is highly dependent on signal preprocessing. This manuscript presents a filtering method to improve the feature classification algorithms typically used in BCI. Methods: A graph Laplacian quadratic form using the Phase Locking Value (PLV) is applied to generate a new filtered signal in the preprocessing stage. Results: The accuracy of the classification algorithms improved significantly (up to 27.18% in the BCI Competition IV dataset, and up to 42.56% with records made with an Emotiv EPOC+). In addition, the proposed filtering algorithm has similar or better results when compared with the Filter Bank Common Spatial Pattern (FBCSP), which has disadvantages in a multiclass classification. Conclusion: This paper shows how our PLV-based filtering between EEG channels could improve the performance of a BCI.


I. INTRODUCTION
B RAIN-COMPUTER interfaces are of special interest due to their potential in a wide range of situations [1]. The most practical contexts are video games, in which the user can control a game character or the environment, or control systems for peripherals that help people, such as wheelchairs or prostheses for people with motor impairments [2], [3]. BCI systems currently suffer from a lack of standardization in terms of measurement devices-their configuration, preprocessing methods and signal conditioning-which has a decisive influence on their result and accuracy. There is also no standardization in terms of the tools to use to extract and classify information. All this means that this technology is not yet an everyday tool that can be applied by the end user without expert supervision [4]. A BCI user also requires their EEG device to have the same features as those provided by clinical measurement devices, but at a low cost and with ease of use. However, these requirements are detrimental to the quality of the final product. With today's technology, easy-to-use wireless EEG measurement devices do not offer the same performance as more sophisticated systems, such as those used in clinical applications. Low-cost measurement equipment is usually limited in structure (fewer sensors, lower sampling rates, etc.) and has a poor signal-to-noise ratio. This makes the extraction of information more complex and causes problems obtaining a good final classification.
Extensive techniques have been applied to alleviate the effects of this noise in motor imaging task classification and improve the BCI performance. For example, Rahman et al. [5] with their approach, Space-Frequency Spatial Localized Filtering (SFLSF), improved their classification results. Their method begins by dividing the scalp into local overlapping spatial windows to apply local frequency bands through a filter bank, and later a spatial filter obtains the features to be classified. Yongkoo Park et al. applied local Common Spatial Pattern (CSP) to the best channel regions that are defined using several methods: an interquartile range (IQR) or an "above the mean" rule based on a variance ratio dispersion score (VRDS) and an inter-class feature distance (ICFD) [6]. Jin et al. [7] proposed a novel time filter that introduces local time weighting into the objective function of the taskrelated component analysis (TRCA)-based method and uses singular value decomposition. Also, they proposed another method based on the Dempster-Shafer theory to take into consideration the distribution of features of motor imagery trials [8]. Other examples are the novel Correlation-Based Temporal Window Selection (CTWS) algorithm for MI-based BCIs proposed by Feng et al. [9], or the Discriminative Canonical Pattern Matching (DCPM) algorithm proposed by Xiao et al. [10].
This manuscript presents our PLV-Spatial Filtering (PLV-SF) method, which is a novel process based on a graph model to capture the spatial dependencies between EEG device systems and reconstruct the signals by exploiting their spatio-temporal relationships and removing the existing noise.
The graph model is based on the phase synchronization metric between channels; specifically, the phase locking value (PLV), which is widely employed to measure synchronization between EEG acquisition sensors in BCI, has been used. The PLV represents the interaction between two processes or systems with dynamic states. Moreover, this metric is statistically indicated by a non-random phase or phase difference distribution, and is therefore less affected by noise than amplitude synchronization [11]. Moreover, as Kolev et al. [12] demonstrated in their study to differentiate the brain responses of two age groups, the PLV works much more accurately to assess dynamic brain processes. In the field of neuroscience, the PLV is typically extracted between several sensors that record brain signals. Benzy and Vinod et al. [13] select the most significant pairs of electrodes based on the PLV in different frequency bands to differentiate the movement of both hands. Wang et al. fuse the PLV with the one-versusthe-rest filter-bank common spatial pattern (OVR-FBCSP) to improve the robustness of motor imagery (MI) classification. Tosun et al. [14] use sliding windows to obtain the eigenvalues of the PLV in the intrinsic mode functions (IMFs) of each channel, classifying four motor image tasks and demonstrating the high relevance of the PLV in these types of tasks.
Graph signal processing has become increasingly interesting in recent years [15], [16], [17]. There are works in which signals are sampled and reconstructed [18], and it is important to consider that, in a real context, signals vary over time, but they evolve smoothly and do not exhibit abrupt changes.
For this reason, a typical measure of the smoothness of graphs is the graph Laplacian quadratic form. In this study, we propose our approach, which includes a new signal preprocessing step. This method combines the graph Laplacian quadratic form and the PLV between sensors. The signal at each electrode location is determined as the average of the remaining electrode signals weighted according to their phase synchronization levels in the original signal. To verify that this preprocessing method improves the classification results, well-known classification algorithms were applied to the extracted feature. In addition, the PLV-SF was compared with the Common Spatial Pattern method. Since the trend in recent years has been to process signals in certain specific bands, both the application of the PLV-SF and the calculation of the CSP for the different bands studied have been carried out, the latter being known as Filter Bank CSP [19], [20], [21], [22], [23].
This manuscript describes the steps involved in this process realization, starting with the Materials and Methods section, which describes the data and the software used in this experiment, the concrete description of the proposed method and the classification algorithms that will validate its effectiveness. Next, in the Results and Discussion section, the results obtained are described, as well as the results of other authors with similar methodologies. Finally, the Conclusions section highlights the key points found in this study.

II. MATERIALS AND METHODS
The materials and procedures involved in the proposed filtering model methodology are listed and described below.

A. Data and Analysis Software
Two datasets were used. The first one was the BCI Competition IV 2a [24], which was provided by the Institute for Knowledge Discovery (Laboratory of Brain-Computer Interfaces), Graz University of Technology. The device used has a sampling rate of 250 Hz, and its sensors follow the 10-20 distribution, as shown in Fig. 1. Additionally, these data were previously filtered with a band-pass filter between 0.5 and 100 Hz, with an additional 50 Hz notch filter to remove line noise.
The second dataset was recorded using a wireless low-cost device, Emotiv EPOC+, and the Emotiv Xavier Testbench v3.1.21 software. The device consists of 16 sensors and two reference electrodes at P3 (CMS, Common Mode Sense) and P4 (DRL, Driven Right Leg)), as per the 10-20 system as well, at an acquiring frequency of 128 Hz, Fig 1. The data were analyzed using Matlab® (The MathWorks, Inc.) and Fieldtrip [25]. Fieldtrip is a toolkit for preprocessing and advanced analysis methods for MEG, EEG, iEEG and NIRS recordings, which was developed by the Donders Institute for Brain, Cognition and Behavior in Nijmegen (The Netherlands), in collaboration with other institutes. Fast-FC was also used, which is an open-source toolbox to efficiently compute connectivity indices [26].

B. Participants
The BCI Competition IV 2a dataset consists of EEG recordings from 9 people, and contains data on the baseline state, motor imagery tasks for both hands, the feet and the tongue [24].
The Emotiv EPOC+ dataset was recorded from 13 people between 18 and 51 years old who do not have any pathologies. The dataset consists of three classes: the baseline state and the motor imagery of both hands. This study was approved by the ethics committee of the University of La Laguna (register number: CEIBA2020-0405).

C. Experimental Protocol
The BCI Competition IV 2a dataset consists of four classes of motor imagination tasks: imagination of the left hand, right hand, both feet and tongue. Each subject performed two sessions (on different days) in which the following protocol was applied: a cross appears on a black screen at startup, in addition to a short warning tone. After two seconds, a signal appears indicating which of the four tasks to perform (up, down, right, and left arrow, associated with the motor imagery task for the tongue, feet, right and left hand, respectively). This signal is held for 1.25 seconds, but the user must hold the task for up to 6 seconds. This task was followed by a 1.5-second rest, after which the protocol was rerun. For more information on the dataset and protocol used, see the detailed description provided in [24]. Because only basal states and motor imagination of both hands were considered, the tongue and foot tasks have been omitted from this study.
The experimental paradigm shown in Fig. 2 was designed for the Emotiv EPOC+ device. Similar to the BCI Competition dataset, in the first part a cross is displayed, and the user does not have to do anything. This is recorded and corresponds to the baseline state (this section is 20-seconds long). Then the cross disappears, and a 10-second rest period starts. While the screen is black, the user does not have to do anything. Finally, a right-or left-arrow is displayed for 20 seconds, corresponding to the motor imagery of the hand (right or left). While the arrow appears on the screen, the user has to perform the task, which is recorded and ends when the arrow disappears. The complete experimental protocol was repeated four times for each motor imagery task, generating the final Emotiv EPOC+ Dataset.

D. Data Preprocessing
Both datasets undergo the same signal preprocessing. First, a band pass filter between 1-40 Hz is applied. Second, the REBLINCA procedure defined by Di Flumeri et al. [27] is carried out to reduce the blinking effect. This procedure is a modified version of the Gratton algorithm [28], where an EOG channel is not needed because the front sensor placed on the midline sagittal plane of the skull is used.
In this paper, the Fz channel is used for the BCI Competition IV, while in the case of the Emotiv EPOC+ dataset, as it did not have a central sensor, a reference signal (AFz) is generated as the average of AF3 and AF4.
These signals (Fz or AFz) are used to calculate the REBLINCA regression signal (RegrSignal) using a Butterworth filter between 1-7 Hz for each dataset (frequency range used in the Gratton algorithm with the V-EOG channel). Then, RegrSignal is derived and normalized to have a mean of zero and a variance of one. Finally, a moving average is applied to the square of this signal. The outcome is a new signal called the threshold signal (ThresSignal), which is used to intensify and locate the blink disturbance. Then, in the blinking areas, the original signal is corrected by subtracting the weighted regression signal (RegrSignal) for each channel: where we previously estimated the weight of the proportion of the main reference signal present in each channel as: Finally, an automatic artefact rejection as defined in the Fieldtrip tutorial [29] is applied: the signal of each channel (x ch (t)) is standardized (z ch (t)) with respect to its standard deviation (σ ch ) and its average (μ ch ) over time (t). After that, a new signal (z sum (t)) is created as the average of all the standardized data of each channel, where N is the number of channels: A threshold (th rej ) of 3σ sum around the mean was chosen. This threshold is applied to remove the outliers, which were considered artefacts and noise.

E. Frequency Bands
Motor imagery is described as the cognitive process of a person thinking that they are performing a task when they are not actually performing any kind of movement. As described above, our case is focused on both hands. Therefore, the alpha (7.5-12.5 Hz) and beta (12.5-30 Hz) bands were used for the study. As can be seen in the literature, these bands have been extensively studied to identify the tasks defined (examples include McFarland et al. [30], Bao-Guo Xu and Song [31] and Hermosilla et al. [32]). Due to the above, it is necessary to pass our EEG signal to a frequency-time context, so the Morlet Wavelet Transform is applied [33] to yield the power for both bands studied (alpha and beta).

F. PLV-SF Method
Many EEG devices have a low signal-to-noise ratio and signal coupling between sensors. This all leads to undesired behavior in a BCI classification system, which could produce errors in a real use case. The method proposed here reduces the noise and diminishes artifacts present in the signal.
The main idea of this method is to apply the synchronization information to exclusively highlight the signal patterns in each task studied. Note that each of the tasks will have a unique synchronization pattern associated with it. Each single sensor signal is filtered considering the functional influence of the other sensors that affect it to a greater or lesser extent, depending on their phase synchronization level.
The procedure has the following phases: 1) Phase Lock Value: The signal obtained from the EEG source has to be transformed into a complex value to compute the synchronization metric for the alpha and beta band in order to apply the filter to each of them.
The phase definition used in this paper is based on the Hilbert Transformation (HT) [26]. The first part of this transformation filters the signal in the frequency band of interest by using a Finite Impulse Response (FIR) filter (once forward and once backward). In the second part, a circular fill is performed to ensure that the filter produces zero phase distortion and minimal edge effect. The HT of the signal is defined as: where p.v. denotes the principal value of Cauchy and the analytical representation of x k (t) is the complex-valued signal defined as: x a k (t) = x k (t) + i x k (5) from which the phase is defined as: Formally, the phase synchronization (PS) of two signals (x k and x l ) is stated as: On a practical level, it is useful to define a time-independent value that characterizes the entire interval of interest. One of the most used methods in BCI is the mean phase coherence, also known as the phase locking value (PLV): where stands for its average value, yielding outcomes in the range from 0 (no PS) to 1 (perfect PS).
2) Graph Laplacian Matrix: From the previous PLVs, the following adjacency matrix is generated: where N is the number of channels. Note that the adjacency matrix is a symmetric matrix.
Once calculated, the adjacency matrix is used to determine the diagonal matrix (D) as: These two matrices are computed to obtain the graph Laplacian matrix as: The graph Laplacian matrix can be interpreted as a difference operator for signals defined on the PLV graph. The nodes of the PLV graph are the sensors (x k ), and the edges are the connections between them, with each one weighted with its PLV (P LV kl for nodes x k and x l ).
3) Graph Laplacian Quadratic Form: Alpha and beta frequency bands were calculated using the Morlet Wave Transform [33]. As a result, each node has a temporal power signal associated with M samples for each band, where M corresponds to 500 samples for the BCI Competition IV dataset and 256 samples for the Emotiv EPOC+ dataset (corresponding to two seconds in each system). Therefore, the matrix time-series signal can be modeled as: where x is a column vector of X. From this, the graph Laplacian quadratic form of the graph can be expressed as: This equation represents the weighted sum of the neighborhood variations across all nodes (see reference [15], on graphical signal processing, for a detailed definition). From this definition, it is possible to calculate the graph Laplacian quadratic form of the trial over time as: 4) Convex Problem: Finally, the signals with the channel to be filtered can be modeled as follows: where B corresponds to a binary matrix in which the channel to be filtered is eliminated, and • denotes the Hadamard product operator. Therefore, the X signal is generated by solving the following convex problem: where F denotes the Frobenius (or L2) norm for matrices. This process of eliminating and generating each channel's signal with the information from the remaining signals is carried out individually. Therefore, applying this convex problem to all the channels of a trial yields a new matrix with the same dimensions. Its content is the filtered signal based on the actual values of the remaining channels. Furthermore, the application of the synchronization-based adjacency matrix places greater importance on sensors that have a higher PS. In addition, we studied whether the filtered signal (plvSig) at each location lost information contained in the original real signal (realSig) that could improve the classification. Therefore, to see the effect of PLV-SF, the classification accuracy is calculated with a new rebuilt signal (fSig) that is generated with different forgetting factors (FF). By increasing this factor, the influence of the filtering increases at the cost of decreasing the contribution of the original signal:

G. Classification
The classification procedure carried out seeks to contrast the effectiveness of the proposed filtering method previously defined by using various classification algorithms that are well known in the literature, Table I.
The classification algorithms are applied to each user's data, using each temporal sample of all the channels. On the one hand, unfiltered data and, on the other hand, data with the application of PLV-SF are used in order to compare this proposed methodology. Therefore, the classification algorithm will have as input a feature for each channel and frequency band, yielding a dataset with the total size of all the registered samples.
Because the split of each user's original sets is randomly assigned to the cross-validation and testing processes, the procedure is repeated 20 times. To compare the algorithms, the mean population accuracy is estimated, where the accuracy of each individual has been calculated as the average of the 20 test procedures.

III. RESULTS AND DISCUSSION
The results obtained in this research show how the proposed method, based on PLV graph Laplacian quadratic form, improves certain classification algorithms. The new values for the power are influenced by the channels with which they have the highest synchronization. If the power of these channels is high, the filtered channel has a high value; if not, it will have a low value. Generically, the filtering decreases the power where the original signal is extremely high whereas in cases where this signal is low, its value is increased.
The results obtained with different FF values are shown in Fig. 3 and Fig. 4, where the original data has a 0-FF value, and the complete filtered signal has a 1-FF value. The methodology applied was repeated several times, as defined in the classification in the Materials and Methods section. These figures show how the accuracy value is highest when the proportion of the original signal used is the lowest (FF value equal to 1). Specifically, the accuracies of the algorithms at these extreme FF values, where the original signal and the fully filtered signal are used, show increases between −1.62% and 27.18% in the BCI Competition IV, and between −2.37% and 42.56% in the Emotiv EPOC+ dataset.
Since not every algorithm showed a significant improvement, the Wilcoxon test (a non-parametric test) was performed to confirm any gain, Table II. The p-values obtained by comparing the results of the raw and filtered data show whether there is a significant improvement for each algorithm. In both datasets, the algorithms of the Discriminant Linear Analysis, SVM with a Linear or Quadratic kernel, and the KNN with a Fine, Medium, Cubic, Weighted or Cosine kernel stand out. Although in the case of the Emotiv EPOC+ dataset, all the algorithms show significant improvements, except the Quadratic Discriminant Analysis.
These results obtained by the PLV-SF method are comparable with recent studies on synchronization metrics, including that of Rodrigues et al. [34]. They used the BCI Competition IV dataset and classified its four main classes (the motor imagery task of the right/left hand, tongue and feet) by using  Table I shows which algorithm corresponds to each number (#). different features, such as Pearson's correlation, Spearman's correlation, phase coherence, STR, etc. Finally, they obtained mean accuracy values in the range of 37% to 57% using a least squares classifier. We also find the paper by Ai et al. [35], in which their proposed method, the combination of features of functional brain networks with the common spatial pattern (CSP) algorithm and local feature scale decomposition (LCD), yields an accuracy of 79.7% with this same dataset. The values obtained in our research with the three classes selected (basal and motor imagery of the right-and  left-hand state) show similar accuracy results with two algorithms (SVM with Quadratic kernel and KNN with Cosine kernel).
Elsewhere, the research by Khan et al. [36] aims to improve multiclass classification accuracy for motor image movement using common sub-band spatial patterns with the sequential feature selection method (SBCSP-SBFS). They used two datasets, one from Emotiv EPOC+ and the other from wet gel electrodes, both to classify three classes (motor imagery of the left and right hand and the rest state). Their results yielded a maximum accuracy of 60.61% for the Emotiv EPOC headset and 86.50% for wet gel electrodes. Another example is the research of Lin et al. [37] using a BiLSTM neural network as a binary classifier (motor imagery right and left hand), yielding an accuracy of 87.14%. Working in similar contexts, the paper by Mwata-Velu et al. [38] presents a hybrid architecture of convolutional neural networks (CNN) and long short-term memory networks (LSTM), which achieve an accuracy of 84.69% and 79.2%. They applied their own data and a state-of-the-art dataset of motor imagery EEG in an effort to create their own BCI system with an Emotiv system. These accuracies are comparable to the PLV-SF method using the Linear Discriminant Analysis, SVM (Linear and Quadratic), and KNN (with any kernel except the Coarse one), whose results exceed 80% accuracy.
Furthermore, PLV-SF method was contrasted with a wellknown spatial filter method, the FBCSP. In this case, the MI of both hands of the Emotiv EPOC+ dataset is classified. Only these two classes were used in order to apply the original CSP algorithm, which is based on a binary comparison. Four different cases were studied, in which the classification algorithms were used with: the raw data, the data filtered by CSP, the data filtered by the PLV-SF and the combination of both filters (PLV-SF-CSP). Fig. 5 shows the accuracies obtained from this analysis. If we apply the non-parametric Wilcoxon test, the p-values show significant differences in five algorithms: Discriminant Analysis with a Linear (p < 0.01) and a Quadratic (p < 0.01) kernel, and the SVMs, which use a Linear (p < 0.01), Quadratic (p < 0.01) and Medium Gaussian (p < 0.01) kernel. We can see from Fig. 5 that the significance in the SVM algorithm with Quadratic and Gaussian kernel algorithms shows that FBCSP is better, while in the case of Linear Discriminant Analysis and the SVM with the Linear and Quadratic kernel, the PLV-SF obtains significantly better values. If the same non-parametric test is applied with the FBCSP method versus the combination of both, the PLV-SF-CSP method, there are no significant differences in the Decision Trees with a Medium (p = 0.649) and Coarse (p = 0.2381) kernel and with a Quadratic Discriminant Analysis algorithm (p = 0.6081). Therefore, by observing Fig. 5 we can affirm that the combination of both methods is better in all the algorithms except in these three, which yield similar results. Similarly, it is possible to apply the Wilcoxon test between the values obtained with the PLV-SF and the PLV-SF-CSP method, which shows significant differences in all the algorithms (p < 0.05), except with a Fine Decision Tree (p = 0.1008). Looking again at Fig. 5, all the algorithms show improvements except two, the Linear Discriminant Analysis and the Linear SVM, in which simply applying the PLV-SF method is better. Yu et al. [39] show that our results are comparable to theirs. They perform a binary comparison with users selected from several datasets, obtaining average accuracies ranging from 73.75 to 86.11. Elsewhere, the binary classification carried out by Ghanbar et al. [40] obtains accuracy ranges between 58.61% and 83.48%. In our case, the Linear Discriminant Analysis, Linear SVM and Quadratic SVM algorithms stand out, with values exceeding 90% accuracy.
We must emphasize the importance of our method in terms of data processing, since in the case of CSP, it is necessary to work not only with the class that we want to filter, but also with the remaining classes, increasing its complexity in a multiclass classification context. The PLV-SF method works individually with each class.

IV. CONCLUSION
There is an extensive body of literature on the development of BCI systems that require high accuracy to ensure proper results. Specifically, in this manuscript, a well-known metric, namely PLV, which is commonly used as a feature in classifiers, has been shown to be an effective metric to use in the graph Laplacian quadratic form to filter the signal, improving the performance of the BCI.
The PLV-SF method significantly facilitates the task of a wide set of classifiers commonly used in this type of application, which is an essential aspect for good performance in a real context. Moreover, similar or better results (depending on the algorithm used) were obtained compared to the FBCSP method.
In summary, the application of PLV as a filter method reduces the amount of noise and artifacts in the signal, while highlighting the characteristics of the EEG signal for each task. This makes the data preprocessing more accurate, and the classifier performs better.