Automatic detection of epileptic slow-waves in EEG based on empirical mode decomposition and wavelet transform

2.1. Subjects and recordings

EEG data applied in this paper are collected from the subjects of the second affiliated hospital of Xi’an Jiaotong University, Xi’an, China. The subjects include 28 patients with epilepsy (12-25 year-old males and females) and two normal adults (a 22 year-old male and a 30 year-old male. They are volunteers in Xi’an Jiaotong University, Xi’an, China. They are able bodied, free from medication and any disorders of, or injuries to, the central nervous system.). 21 out of 28 patients are training set to determine the threshold of detecting slow-wave and the remaining 7 patients are test set. To obtain dataset 1 including blinking artifacts, subjects (normal subjects) are instructed to self generate blinking at ease. We also achieve multi artifacts EEG data (dataset 2) which contain blinking artifacts and Electromyography (EMG). Data acquisition system is WEE-1000A/K, Nihon Kohden Corporation, Japan. Exploring cup electrodes are fixed to the scalp at Fp1, Fp2, F7, F8, C3, C4, T5, P3, P4, T6, O1, and O2 according to the International 10-20 System. The reference electrode is located on the ipsilateral ear electrode. The epileptic EEG signals, dataset 1, and dataset 2 are recorded at a sampling frequency of 200 Hz. All the EEG data are amplified with bandpass filter of 0.5–60 Hz. 50 Hz notch filter is employed to wipe power frequency out. Block diagram of the EEG signal processing is shown in Fig. 1.

Fig. 1Block diagram of epileptic EEG signals diagnosis

2.2. Artifacts removing methods

BA may be presented throughout the recording process and the amplitude of this artifact is much greater than the EEG signals [12]. In this paper, we use EMD and regression method to eliminate it.

1) Empirical mode decomposition: empirical mode decomposition (EMD) proposed by N. E. Huang in 1998 [15] is a nonlinear signal processing method. It decomposes signals into a series of intrinsic mode functions (IMFs) indicating the bases of the signals. Each IMF satisfies two basic conditions: (1) in the whole data set, the number of extrema and the number of zero-crossing must be the same or differ at most by one; (2) the local mean value of the envelope defined by local maxima and the local minima is zero. The algorithm of EMD for extracting IMFs from $\text{x}\left(t\right)$ is been summarized as follows [16]:

(1) Compute all the local extrema including maxima and minima of $x\left(t\right)$.

(2) Obtain the upper envelope $e_n\left(t\right)$ and lower envelope $e_1\left(t\right)$ respectively by cubic spline interpolation.

(3) Compute the mean $m\left(t\right)$: $m\left(t{)=(e}_n\right(t{)+e}_1\left(t\right))/2$.

(4) Extract $h_1\left(t\right)$ from $x\left(t\right)$ as $h_1\left(t\right)=x\left(t\right)-m\left(t\right).$ Check $h_1\left(t\right)$ by the two basic conditions and repeat step 1-4 until an IMF is obtained.

(5) Compute the residual as: $r_1\left(t\right)=x\left(t\right)-c_1\left(t\right)$, when the first IMF $c_1\left(t\right)$ is extracted.

(6) Consider $r_1\left(t\right)$ as new $x\left(t\right)$ and repeat steps 1-5 until $N$ IMFs are obtained.

If residual $c_N\left(t\right)$ or $r_N\left(t\right)$ is a constant or a function with only one maxima and minima, the procedure stops.

After the procedure, $N$ IMFs and a residual $r_N\left(t\right)$ are achieved. The original signal $x\left(t\right)$ can be described as Eq. (1). The lower order IMF owns the higher frequency and vice versa.

2) Regression method: regression methods often assume that the scalp potential is a linear combination of brain and blinking potentials. By subtracting propagated blinking from EEG recordings, EEG signals can be recovered [17]. The algorithm is summarized as Eq. (2). [9, 10, 18, 19]:

where ${EEG}_{est,i}\left(t\right)$ is the reconstructed EEG. ${EEG}_{m,i}\left(t\right)$ is the EEG signal except Fp1 and Fp2. ${EEG}_{b,j}\left(t\right)$ is EEG polluted by BA in Fp1 or Fp2. $i$ is electrode index except Fp1 and Fp2. $j$ indicates Fp1 or Fp2. $\alpha$ and $\beta$ are computed in Eq. (3) and Eq. (4) respectively:

where ${\gamma }_{i,j}$ stands for the correlation coefficient of ${EEG}_{m,i}\left(t\right)$ and ${EEG}_{b,j}\left(t\right)$. $s{d}_i$ and $s{d}_j$ indicate the standard deviations of ${EEG}_{m,i}\left(t\right)$ and ${EEG}_{b,j}\left(t\right)$ respectively. $x_j$ and $y_i$ are the means of ${EEG}_{b,j}\left(t\right)$ and ${EEG}_{m,i}\left(t\right)$ respectively.

2.3. Slow-wave detecting method

Wavelet transform (WT) that is proposed in the late 1980s has the advantage of analysis signals in frequency-scale. The usage of WT is widely. A common use is wavelet filter [20, 21] besides, epilepsy detection is also an area in which the WT has been particularly successful since WT captures transient features and localizes them in both time and frequency content accurately [12]. The wavelet transform is an integral transform and the formula of wavelet transform is described in Eq. (5):

where $\phi \left(t\right)$ is the mother function wavelet. $a$ is the dilation parameter and $b$ stands for the translation parameter.

3.1. Optimal weight for combination

The combination method ($\lambda$ and $\mu$ in $\lambda \text{IMF2}+\mu \text{IMF3}$) takes the summation value between the two input IMF2 and IMF3. For searching optimal weight, we make comparison on three different ways. The first one is $\lambda =\mu =1$. The second is that $\lambda =m$ (in Eq. (6)) and $\mu =1-m$ when the energy of IMF2 is lower than IMF3, i.e., $m$ is less than 1; on the contrary, we define $\lambda =1-m$ and $\mu =m$. The last is similar with the second except that $n$ (in Eq. (7)) is employed:

where $EN\text{(IMF}i\text{)}$ and $VAR\text{(IMF}i\text{)}$ are the average energy and standard deviation of the $i$-th IMF respectively ($i=$2, 3).

In Fig. 3, box 1, box 2, and box 3 represent correlation coefficients calculated by three different combination ways mentioned above respectively. It is obviously that the first combination method ($\lambda =\mu =1$) yield higher correlation coefficients than those using other ways. This suggests that the threshold defined by the first way is higher than others and it benefits to decrease erroneous judgment. Fig. 3 also clearly shows that the first way ($\lambda =\mu =1$) holds higher robustness when another artifact is included.

Fig. 3Statistical comparison. The box 1 is significantly higher than the second (p< 0.05), this situation also appears when comparing the first with the third (p< 0.05). Figures (a) and (b) are plotted by dataset 1 and dataset 2 respectively. Each sample stands for five seconds EEG signals. The length of data should not be too short for preventing not enough maxima and minima existing or too long for avoiding reduction of the correlation coefficients between original EEG signal and IMF2 + IMF3

a)

b)

3.2. Artifact removing

Therefore, IMF2 + IMF3 can be served as useful signal in identifying BA from EEG background. It is more clearly in Fig. 2 (a)-(d). We also analyse the clinical EEGs containing BAs and the results validate this finding. Correlation coefficient is employed to quantify the similarity between original EEG signal and IMF2 + IMF3. If the correlation coefficient of original EEG signal and IMF2 + IMF3 is higher than the predefined threshold in Fp1 or Fp2, BA in the original EEG signal is detected.

The threshold discriminating blinking artifacts should be confirmed. The databases contain clinical epileptic EEG data, dataset 1, and dataset 2. Segmentation is performed for a set of thresholds between 0 and 1 with a step of 0.005. For each threshold, the number that is accordance with the prior awareness form of data style is counted. The threshold which leads to the maximum number is kept as the definitive threshold. If several threshold levels lead to the same number, the lowest one is taken. After analyzing, the threshold is set to 0.715.

For purpose of removing blinking artifacts completely, we filter the original polluted EEG signals in Fp1 or Fp2 with bandstop filter of 5-8 Hz that is the main frequency of blinking artifacts. Considering the BAs that are obtained by subtracting the filtered data from the original EEG as template, the influence of artifacts on other electrodes (except Fp1 and Fp2) is reduced by regression method. Fig. 4 shows the performance of artifact cancelling method.

Fig. 4a) Original EEG signal which contains blinking in Fp1, b) original EEG signal which contains blinking in Fp2, c) original EEG signal polluted by blinking in F3, d) the EEG signal processed by artifact cancelling method in F3. Blinkings in Fp1 and Fp2 are boxed in boxes

3.3. Slow-wave detection

To locate slow-waves, we apply the db6 wavelet function because the shape of its mother wavelet resembles slow-wave. Fig. 5(a) and (b) display the scale function and the wavelet function of db6 respectively.

The analysis of slow-wave based on discrete wavelet transform (DWT) is performed in three stages: 1) decomposition, using DWT to analysis of EEG signals in four levels; 2) feature extraction, selecting the frequency band, the fourth level, that involves the frequency of slow-wave and cubing the slow-wave wavelet coefficients; and 3) label, detecting and marking the slow-wave relying on slow-wave threshold.

Since EEG is complex in clinic and the amplitude of EEG is inter- and intra-subject variability, statistical method is the best choice. Here we select 500 hours EEG data from 21 patients which contain epileptic slow-waves. For each patient, the result of marking slow-wave is censored by two experts. The method for determining threshold is identical with the artifacts removing threshold and the optimal threshold is decided to 0.403.

Finally, the slow-wave detection method is tested on 7 clinical epileptic EEGs. The results are very promising and the average false detection rate [3] is 15.67 %. The details are presented in section 4 (Table 1).

Fig. 5The wavelet applied: a) scale function of db6, b) wavelet function of db6

Fig. 6Decomposing of epileptic EEG signal that contains slow-wave by EMD: a) the original epileptic EEG signal which contains slow-wave at about 1 s, b) IMF2, c) IMF3, d) IMF2+IMF3

3.4. Performance of artifact removal method

The artifact removing method we present should not disturb the detection of slow-wave (Fig. 6). We compare Fig. 2 with Fig. 6 and the result indicates that the similarity of original signal and IMF2 + IMF3 is greatly lower in Fig. 6, i.e., in the case of blinking artifacts contained, the signals of IMF2 + IMF3 are more similar with original signals. Consequently, the artifact removing method can’t generate misjudgment between slow-waves and artifacts. We also find that the proposed method of eliminating blinking artifacts in this paper can decrease the false detection rate of slow-waves (Fig. 7). It is more obvious in the boxes (Fig. 7(c)-(d)).

Fig. 7a) The original EEG signal which contains blinking artifacts in Fp2 (blinking which pollute P4 channel’s signals is circled), b) the original data which do not contain slow-waves in P4, c) third order of wavelet coefficients of the signal which is preprocessed by blinking cancelling method in P4, d) third order of wavelet coefficients of the original signal which do not have to be preprocessed in P4

4.1. Discussion

The advantage of third order of wavelet coefficients method comparing with traditional first order of wavelet coefficients method is depicted in Fig. 8. The original data are epileptic EEGs containing slow-waves in O2. This figure makes it clear that the difference between slow-wave and other waves expands when cubing the wavelet coefficient (it is more obvious in the boxes). Hence, the accuracy ratio of detecting slow-wave must rise up with third order of wavelet coefficients applying.

Fig. 8Comparison of third order of wavelet coefficients and first order of wavelet coefficients: a) third order of wavelet coefficients, b) first order of wavelet coefficients

Detailed results are displayed in Table 1. All the data we analyzed are reviewed by two experienced clinicians and they also determine the number and location of slow-waves in each clinical EEG segments. It is obvious that the false detection rate reduces by 53.5 % corresponding to the traditional on average.