EEG signal improvement with cascaded filter based on OWA operator

Noise reduction methods have a great impact on the performance of all EEG signal processing systems. The work involves the reduction in impulsive interferences. It is impossible to suppress effectively such a noise using linear filtering approach. The work presents the properties of a cascaded filter based on the ordered weighted aggregation (OWA) operator and its application to improve the electroencephalogram (EEG) signal. The OWA operator is a class of mean-like aggregation operators. By introducing a nonlinear sorting process and by assigning appropriate values of weights to the sorted samples, the improvement in filtering in the impulsive environment is achieved. The structure of the proposed cascade consists of two layers. The first layer contains two OWA filters that share a certain number of signal samples. The second layer averages the outputs of the first layer. The algorithm for the rise of filtering efficiency has been introduced. The performance of the new method has been experimentally compared with the traditional methods of impulsive noise reduction, using synthetic as well as real signals from the CHB-MIT database. The obtained results demonstrate that in the field of impulsive noise suppression and preprocessing, the proposed cascaded OWA filter brings about a significant improvement in the signal-to-noise ratio.


Introduction
Smoothing of biomedical signals should be carried with particular attention. The requirement for efficient suppression of noise must be reconciled with the least possible effect of filtration on the signal. Special attention should be given to the filtering of such waveforms as the electroencephalogram (EEG) signal. EEG is a recording of the electrical activity of the brain from the scalp. The recorded waveforms reflect the cortical electrical activity. It is an invaluable measurement aimed to assess brain activities in terms of information relating to different physiological states of the brain. EEG waveforms allow us to explore non-invasively the electrical activity of the brain. It has been used as a valuable source of information for clinical diagnoses and monitoring, cognitive science as well as the development of new brain-computer B Tomasz Pander tpander@polsl.pl 1 Institute of Electronics, Faculty of Automatic Control, Electronics and Computer Science, Silesian University of Technology, Akademicka St. 16, 44-100 Gliwice, Poland interfaces (BCI). However, it is very important to process reliable signals. The sharpness of the EEG signal must be retained efficiently during the preprocessing stage, but any artefacts or outliers that appear in the signal should be effectively suppressed [1]. EEG signals are contaminated with various artefacts, such as signals associated with muscle activity, eye movements and blinks, body motions or power grid noise, which are not of cerebral origin [1,2]. These types of disturbances appear as high-amplitude or high-frequency bursts of activity and look like impulsive types of noise, whereas brain activities have generally lower energy than artefacts.
The first step of the EEG signal analysis is a preprocessing stage which consist of two processes such as artefact removal and data filtering. Artefacts may be triggered by various factors such as head motion during signal acquisition, physical problems in the electrode/channel/lead or connectivity problems between the head and the device. These noises create also abnormal frequencies and shapes. Numerous signal preprocessing methods have been proposed for the suppression of noise affected by eye movements or muscle activity. Lowand high-pass filters are commonly applied to remove muscle artefacts and drift, respectively. However, the spectral pattern of noises usually overlaps brain signals [3]. The use of frequency methods brings about the removal the noises but also the removal of informative parts of the EEG signal. Recently, advanced machine learning methods have been developed to deal with the artefacts. They include such methods as independent component analysis (ICA), canonical correlation analysis (CCA) [3], discrete wavelet transform (DWT), empirical mode decomposition (EMD) or cascaded Savitzky-Golay filter [1,2,4]. Neurologists have been working on the interpretation of EEG signals, but there are a few cases when real-time monitoring or automatic interpretation is particularly important. For that reason, the acquisition of a high-quality EEG signal from background interferences is an incentive to look for clear, reliable and effective methods of noise suppression. An important requirement underlying the interpretation of the registered brain activity involves proper signal quality. A common way to reach this objective consists in recording the EEG in a shielded laboratory and in preparing the subjects' skin to reduce the electrodeskin impedance before the electrodes are placed. However, these standard procedures limit the accessibility of an EEG device and narrow down its application outside the laboratory [5]. The challenge is to design a filtering method which can ensure the highest quality and low computational effort for possible wearable EEG devices.
The main objective of the work is to demonstrate the structure and properties of the cascaded filter based on the ordered weighted aggregation (OWA) operator and the applicability of the cascade OWA-based filtering in enhancing the SNR during the preprocessing of the EEG signal. All the considered filters are based on the moving data window. The outline of the paper is as follows. The fundamentals of OWA operators are described in Sect. 2.1. The proposed cascaded OWA (COWA) filter and its properties are presented in Sect. 2.2. Numerical experiments, results and comparison are presented in Sects. 3 and 4. Section 5 provides conclusions.

OWA filters
OWA operators were introduced by Yager [6], and they have been subjected to numerous research studies [7][8][9][10]. The forecasting model involving the distribution of electric power load is investigated in [7]. Different aggregation of OWA operators is presented in [8]. The methods of data smoothing are described in [9], whereas the application of OWA operator for noise suppression in ECG signal processing is presented in [10]. Smoothing is an operation which removes high-frequency fluctuations from a signal. It can be expressed as where x i denotes a discrete input signal and y i is a signal after low-pass filtering. The parameter N specifies the observation window length equal to M = 2N + 1 as well as the filter low cut-off frequency. Equation (1) defines the well-known moving average filter. In spite of its simplicity, the moving average filter is optimal for common problems: reducing random noise while maintaining a sharp step response. This makes it the premier filter for time-domain signals.
As to the aggregation aspect of data smoothing, the OWA operator describes a family of averaging operators parametrized by its weighting vector [7]. For a fixed normed weighting vector w = [w 1 , w 2 , . . . , w M ] ∈ [0, 1] M and M i=1 w i = 1, the related OWA operator is given by [6,8,9] OWA (x 1 , where (M) and M = 2N + 1. Though this permutation need not be unique (in the case of some ties), formula (2) yields the same result for any acceptable permutation [8]. The OWA operator is parametrized by the selection of weights.
if constant weights are used, then the OWA filter is a simple average [7]. In [9], several parametrized families of OWA operators were suggested for aggregation of time series data. It is important to emphasize that the weights are associated with a particular ordered position of sample rather than a specified element [7]. The necessity of signal smoothing, enhancing or preserving the shape of a signal in the cases when impulsive noises appear involves the application of a robust method as the only alternative to linear filtering. The considered problem of the OWA filter output is related to the filtering process of a time series x i using a sliding window. The output y i can be interpreted as an OWA based on the input samples [11] where w is the vector of weights and M = 2N + 1 stands for filter window length. The operation of the filter given by Eq. (3) is as follows. First, the samples in the observation window are sorted with respect to their values (rank-ordered). Then, the outliers are located at the ends of the observation window. Next, the weighted averaging operator is applied. This operation may be interpreted as an aggregation of knowledge about noise included in the samples [10].
The filtering requirements enforce that the selection of weights should be done with a special care. In the case of a simple average, the weights w i are equal to 1/M for i = 1, . . . , M. Such a window function has crisp boundaries. After nonlinear operation, i.e. the ordering of samples x i , the resulting value of the OWA filter output is affected by the outlier(s). It seems to be more reasonable to use the window function with a soft transition between the inclusion and exclusion of the input signal samples [10]. Many window functions fulfil such requirements. It means that the impact of an outlier on the output value will be small, and the impact of the remaining samples will increase. The overview of the weighted window functions for OWA filters is presented in [11]. For aims of this work, the Gaussian weighted function is used according to the following formula w i = exp − i 2 2σ 2 , where σ is the standard deviation. The tail of the Gaussian window decays rapidly, and the corresponding weight values are small at the ends of the window. The outlier samples after nonlinear operation of sorting are placed at the ends of the data window, and it results in a greater impact of the samples in the middle of the data window on the output value of the OWA robust filter [11].

Cascaded OWA filter
The idea of a cascaded OWA filter comes from the cascaded moving average (CMA) filter presented in [12]. The CMA filter is a hierarchical model that combines two moving average filters placed one after the other. The first layer of the cascaded filter contains an M-order and an N -order moving average (MA) filters. However, in this structure K -samples are common for these two moving average filters. There are samples that overlap at the end of the M-order MA filter as well as at the beginning of the N -order MA filter. The objective of the second layer of the CMA filter is to smooth the piecewise linear trends obtained from the outputs of the two MA filters in the first layer. Finally, the output of the CMA filter is the average value of both MA filters from the first layer [12].
The cascaded OWA (COWA) filter is based on the idea of the CMA filter, and the structure of the COWA filter is  Fig. 1. The first layer consists of two OWA filters. (However, it can also consists of more than two OWA filters; for simplicity, only the case of two OWA filters is considered.) The symbols M and N denote two OWA filters of the order (length), respectively, M and N defined by Eq. (2). The weight vectors are determined with the help of the window weighted functions presented in [11]. The two input windows have the common part of K samples satisfying the following condition 0 K M and 0 K N . There are overlapping samples at the end of the M-order OWA filter as well as in the beginning of the N -order OWA filter. The objective of the second layer of the COWA filter is to smooth the outputs of two OWA filters in the first layer. Finally, the output of the COWA filter is given as where w 21 , w 22 are the weights in the second layer of COWA filter, and they satisfy the following condition: w 21 + w 22 = 1, and y i,M , y i,N are the outputs of the OWA filters in the first layer. For simplicity, we assume in the work that w 21 = w 22 = 1 2 . However, the operation presented by Eq. (4) can be regarded as a special case of the OWA operator (a simple average value), and different OWA operators can be also applied. The second assumption is that M = N . As a result, the symmetry between the outputs of the first layer of the COWA filter is achieved.

Selection of number of overlapping samples
The number of the overlapping samples K satisfies the following condition: 0 K M and 0 K N . When there are no overlapping samples (K = 0), the first layer of the OWA filters produces two shapes that are shifted in time domain with respect to each other by M samples. As a result, the output of the COWA filter is the average signal of these two components. The case of full overlapping (K = M or K = N ) results in two identical signal components (for M = N ), and then, the output of the proposed COWA filter is the same as that in the single OWA filter. However, it is necessary to select K between these two limit cases.
To verify the influence of K value selection, an artificial EEG signal was used. It was generated according to the mathematical description given in [1]. The obtained composite waveform contains typical variabilities of the EEG signal [1].
Additionally, the artificial input signal was disturbed by the Gaussian noise. The signal-to-noise ratio (SNR) was , where σ 2 d is the variance of the uncorrupted signal and σ 2 ν is the variance of the observation noise. The impulsive noise ζ i is usually modelled as a Bernoulli-Gauss (BG) process, i.e. ζ i = ω i · ε i , where ω i is a Bernoulli process with the probability density function described by P(ω i = 1) = P and P(ω i = 0) = 1 − P, and ε i is a Gaussian process with zero mean and variance σ 2 ε = κ · σ 2 n , κ 1 and σ 2 ν = σ 2 n + Pσ 2 ε = (1 + Pκ)σ 2 n , where σ 2 n = 1 is the variance of zero-mean Gaussian noise and κ = 10 [13]. The results of filtering operation are evaluated by the geometric noise reduction factor (GNRF) where s i is the desired signal, x i − s i is the noise added, y i is the result of filtering and y i − s i is the residual noise. The testing procedure was the following. For each value of P (P = {0.1%, 0.5%, 1%}) in the Bernoulli-Gauss model of impulsive noise and for each value of the SNR (SNR= {0 dB, 5 dB, 10 dB, 20 dB}), 100 corrupted EEG signals were generated. The mean square error (MSE) between the deterministic and the filtered signals as well as the GNRF factor was calculated for each result of filtering with the proposed COWA filter. Finally, the mean value of MSE and GNRF was calculated. The parameter K , which determines the overlapping, was searched in the set {0, 1, . . . , N }. In the simulation N = M = 11. The width of the Gaussian weighted function (in the MATLAB environment) is controlled by the coefficient υ which is inversely proportional to the standard deviation σ and σ = (M − 1)/(2υ) where M is the length of the Gaussian weighted function. The COWA filter was used with υ = 4.5. The obtained results are presented in Figs. 2 and 3.
The case P = 0.1% represents the highest level of impulsiveness, i.e. a relatively high probability of outliers to occur, whereas for P = 1% the level of impulsiveness is relatively low. The COWA filter leads to almost the same values, respectively, for MSE and GNRF, for all values of P. The largest differences are in the case of the lowest value of SNR (0 dB) and for larger values of the number of overlapping samples K for MSE factor. For higher SNR, differences between MSE values are smaller. The impact of the selection of K is particularly clear for the GNRF factor (see Fig. 3) for the highest SNR (SNR = 20 dB), especially for odd values of K . This results in the acquisition of a symmetrical number of samples around the i-th sample. For that reason, K should be an odd number. It can be noticed that the smallest values of MSE as well as the highest values of GNRF are achieved by the COWA filter, but for a different number of the overlapping samples K . For lower value of SNR, the smallest values of MSE (see Fig. 2) are obtained for K = 1, whereas for a higher value of SNR K tends to M (or N ), but K should not reach the case of full overlapping (K = M or K = N ).
Similar observations can be made for GNRF. Its highest value is achieved in the case of lower SNR and for a small number of overlapping samples. If SNR is higher, then GNRF is higher for larger K (see Fig. 3). When K is even, the COWA filter response is shifted in the time domain by half of a sampled period with respect to the desired signal. This causes the GNRF to oscillate. In the case of odd K , the filtered signal is not affected with this problem. The above observations can be useful when selecting the COWA filter parameters.

Tuning of the width of the weighted function
The results presented in Figs. 2 and 3 were obtained for a relatively large value of the υ parameter of the Gaussian window. In effect, not only the outliers are "rejected", but also "ordinary" samples can be suppressed. To avoid such situations and to allow for the decrease in υ value when the impulsiveness level is smaller, the following algorithm is proposed.
In this case, the applied noise estimator is the median absolute deviation (MAD), which defines the standard deviation of noise as where L is the signal length and the factor 1.4826 makes the MAD scale an unbiased estimate of the standard deviation for the Gaussian data. The value τ A can be compared with the locally estimated MAD value τ i for M +N −K samples in the filter moving window according to the following condition where A, B are constants (A < B) and L is the signal length. The locally calculated τ i is a simple outlier detector. If there are no outliers, τ i is smaller than τ A , but when the outliers are present, τ i exceeds τ A . In such a case, the υ value should be larger to suppress the outliers more effectively. Figure 4 presents an example of the artificial EEG signal (SNR=10 dB) disturbed by an impulsive noise (upper plot) and the changes in τ i in the bottom plot. The surpassing of the threshold τ A by τ i leads to the change in the width of the Gaussian weighted function according to the condition (7). In this work, A = 4.5 and B = 6.

The real EEG signals
To test the proposed COWA filter, 40 signals from 5 young subjects (age ranging from 1.5 to 22 years) from the CHB-MIT database were used. This database consists of EEG recordings from paediatric subjects with intractable seizures [1,14]. To evaluate the quality of filtering, the SNR improvement factor was introduced where x raw is the raw EEG signal and x f denotes EEG waveform after filtering. The results are presented in Table 1. The reference filters are: the median filter (MedFilt) and the cascaded Savitzky-Golay filter (casSGol). The parameters for the evaluated filters are the following. The length of median filter is 7; for the COWA filter, M = N = 4 and K = 3; for the OWA filter, N = 4; and for the cascaded Savitzky-Golay filter, the polynomial order is 8 and the length is 21.
The obtained results are presented in Table 1.

Results
The considered signals contain seizure episodes that make the filtering process more demanding for the proposed COWA filter because it is designed to suppress impulses, and at the same time, the certain kind of peaks are important and should remain unchanged. The obtained results revealed the usefulness especially of the COWA filter for a given value of parameter K . (Statistically, this filter yields the highest value of SNR imp in 70% of the investigated cases of real waveforms.) When considering other cases, the highest SNR imp is achieved by the cascaded Savitzky-Golay filter. This superiority is especially evident for the subject chb03 and the fourth recording. For other subjects and waveforms, the COWA filter outperforms the reference filters. The main advantage of the proposed COWA filter is that it requires relatively small numbers of samples to work properly. However, the COWA filter performs two nonlinear operations. (The first one involves sorting of samples.) The typical sorting complexity varies from O(N log N ) to O(N 2 ), depending on the used sorting algorithm. These operations can be carried out using specialized devices like microcontrollers or the fieldprogrammable gate arrays (FPGA). An exemplary filtering of the signal from the subject chb08_02 (channel F3-C3) with the evaluated filters is presented in Fig. 5. In this case, the highest value of SNR imp is obtained with the COWA filter.
In the proposed filter, it is also necessary to determine the values of N , M, K and υ for Gaussian window weights. The experiments and evaluation of filters were performed using the MATLAB environment.

Conclusions
In the paper, the COWA filter for EEG signal improvement is proposed. It is based on a two-layer structure with two OWA filters combined in a cascade in the first layer. The robust COWA filter uses the ordered weighted aggregation operator which is used on signal samples after nonlinear sorting operation. In this method, the signal samples are not rejected but multiplied by appropriate weights calculated with the help of the Gaussian weighted function. The numerical experiments have proved that the amount of overlapping samples in the filter windows has an evident influence on the effectiveness of noise suppression. The tests of the proposed method were carried out using a synthetic signal and real EEG signals with seizure episodes. The results show the usefulness of the impulsive noise filtering method with a smoothing phase. The proposed COWA filter introduces the smallest distortions to the filtered signal and guarantees the highest values of noise reduction factor. For short-length windows in the An example of filtering the EEG signal for subject chb08_02, channel F3-C3. a Real EEG signal, b signal filtered with the COWA filter, c signal filtered with the cascaded Savitzky-Golay, d signal filtered with the median first layer, the COWA filter is much more useful for online filtering than, for example, the cascaded Savitzky-Golay filter. The proposed filter can be adapted for wearable EEG devices as the preprocessing stage. In future works this filter can be extended to system of more complicated structure, e.g. more than two OWA filters in the first layer.