Abstract

A new method for analysis of electroencephalogram (EEG) signals using empirical mode decomposition (EMD) and Fourier-Bessel (FB) expansion has been presented in this paper. The EMD decomposes an EEG signal into a finite set of band-limited signals termed intrinsic mode functions (IMFs). The mean frequency (MF) for each IMF has been computed using FB expansion. The MF measure of the IMFs has been used as a feature in order to identify the difference between ictal and seizure-free intracranial EEG signals. It has been shown that the MF feature of the IMFs has provided statistically significant difference between ictal and seizure-free EEG signals. Simulation results are included to illustrate the effectiveness of the proposed method.

1. Introduction

Epileptic seizures are the outcome of the transient and unexpected electrical disturbance of the brain. The electroencephalogram (EEG) signal has been the most utilized signal to clinically assess brain activities. The detection of epileptic seizures in the EEG signals is an important part in the diagnosis of epilepsy [1].

Parameters extracted from EEG signals are highly useful in diagnostics. Spectral analysis is a common technique used for the analysis of EEG signals, as it reveals the frequencies present in the signal. An underlying assumption of the Fourier transform, however, is that the signal being analyzed is stationary (i.e., the mean value, variance, and frequency content of the signal do not change over time).

Recently, nonlinear methods have been proposed to extract new parameters linked to the electrical activity of the brain. Among these parameters, the Lyapunov exponent provides clinically useful information about the signal [2]; the correlation dimension techniques can contain information about the different neurological states of the brain [3]; the fractal dimension (FD) and entropy measure the complexity or the degree of disorder of the EEG signal [4, 5], while correlation integral, the measure sensitive to wide variety of nonlinearities, used in [6], could be used to characterize the epileptogenic regions of the brain during the interictal period. However, recent work shows that the EEG signals exhibit nonstationary behavior [7, 8].

In this paper, a new technique of EEG signal analysis is presented, which is based on the empirical mode decomposition (EMD) developed specially for nonlinear and nonstationary time-series analysis [9] and the Fourier-Bessel (FB) expansion suitable for nonstationary signal representation [10]. The EMD extracts the local oscillations composing the signal, referred to as the intrinsic mode functions (IMFs), as well as the residual representing the local trends. The IMFs can be considered as a set of narrow-band nonstationary signals. The coefficients of the FB expansion have been used to compute the mean frequency of the IMFs. The FB coefficients are unique for a given signal in the same way that Fourier series coefficients are unique for a given signal. However, unlike the sinusoidal basis functions in the Fourier transform, the Bessel functions are aperiodic, and decay over time. These features of the Bessel functions make the FB series expansion suitable for analysis of nonstationary signals when compared to simple Fourier transform [11, 12]. The MF measure of the IMFs has been used as a feature in order to discriminate seizures from seizure-free intervals in intracranial EEG data recordings.

2. Empirical Mode Decomposition

Empirical mode decomposition (EMD) represents any temporal signal into a finite set of amplitude and frequency modulated (AM-FM) oscillating components which are bases of the decomposition. The decomposition is an intuitive and adaptive signal-dependent decomposition. Moreover, the decomposition does not require any conditions about the stationarity and linearity of the signal. The principle of the EMD technique is to decompose a signal π‘₯(𝑑) automatically into a set of the band-limited functions 𝐷𝑝(𝑑) named intrinsic mode functions (IMFs) [9]. Each IMF satisfies two basic conditions: (i) in the complete dataset, the number of extrema and the number of zero-crossings must be the same or differ at most by one, (ii) at any point, the mean value of the envelope defined by the local maxima and the envelope defined by the local minima is zero. The first condition is similar to the narrow-band requirement for a stationary Gaussian process and the second condition is a local requirement induced from the global one, and necessary to ensure that the instantaneous frequency will not have redundant fluctuations as induced by asymmetric waveforms. The EMD algorithm [13] for the signal π‘₯(𝑑) can be summarized as follows. (1)Set 𝑔1(𝑑)=π‘₯(𝑑).(2)Detect the extrema (both maxima and minima) of 𝑔1(𝑑).(3)Generate the upper and lower envelopes π‘’π‘š(𝑑) and 𝑒𝑙(𝑑), respectively, by connecting the maxima and minima separately with cubic spline interpolation.(4)Determine the local mean as π‘š(𝑑)=(π‘’π‘š(𝑑)+𝑒𝑙(𝑑))/2.(5)IMF should have zero local mean; subtract π‘š(𝑑) from the original signal as 𝑔1(𝑑)=𝑔1(𝑑)βˆ’π‘š(𝑑).(6)Decide whether 𝑔1(𝑑) is an IMF or not by checking the two basic conditions as described above.(7)Repeat steps (2) to (6) and end when an IMF 𝑔1(𝑑) is obtained. Once the first IMF is derived, define 𝐷1(𝑑)=𝑔1(𝑑), which is the smallest temporal scale in π‘₯(𝑑). To find the rest of the IMF components, generate the residue π‘Ÿ1(𝑑) of the data by subtracting 𝐷1(𝑑) from the signal as π‘Ÿ1(𝑑)=π‘₯(𝑑)βˆ’π·1(𝑑). The sifting process will be continued until the final residue is a constant, a monotonic function, or a function with only maxima and one minima from which no more IMFs can be derived. The subsequent basis functions and the residues are computed as π‘Ÿ1(𝑑)βˆ’π·2(𝑑)=π‘Ÿ2(𝑑),…,π‘Ÿπ‘€βˆ’1(𝑑)βˆ’π·π‘€(𝑑)=π‘Ÿπ‘€(𝑑),(1) where π‘Ÿπ‘€(𝑑) is the final residue. At the end of the decomposition, the signal π‘₯(𝑑) is represented as π‘₯(𝑑)=𝑀𝑝=1𝐷𝑝(𝑑)+π‘Ÿπ‘€(𝑑),(2) where 𝑀 is the number of IMFs and π‘Ÿπ‘€(𝑑) is the final residue.

Matlab codes are available at http://perso.ens-lyon.fr/patrick.flandrin/emd.html. An example of the application of EMD on the 23.6 seconds EEG time series is shown in Figure 1.

3. Mean Frequency Computation Using Fourier-Bessel Expansion

The zero-order Fourier-Bessel series expansion of a signal π‘₯(𝑑) considered over some arbitrary interval (0,π‘Ž) is expressed as in [10] π‘₯(𝑑)=π‘„ξ“π‘š=1πΆπ‘šπ½0ξ‚€πœ†π‘šπ‘Žπ‘‘ξ‚,(3) where πœ†π‘š,π‘š=1,2,3,… are the ascending-order positive roots of 𝐽0(πœ†)=0, and 𝐽0((πœ†π‘š/π‘Ž)𝑑) are the zero-order Bessel functions. The sequence of Bessel functions {𝐽0((πœ†π‘š/π‘Ž)𝑑)} forms an orthogonal set on the interval 0β‰€π‘‘β‰€π‘Ž with respect to the weight 𝑑, that is, ξ€œπ‘Ž0𝑑𝐽0ξ‚€πœ†π‘šπ‘Žπ‘‘ξ‚π½0ξ‚€πœ†π‘›π‘Žπ‘‘ξ‚π‘‘π‘‘=0,forπ‘šβ‰ π‘›.(4) Using the orthogonality of the set {𝐽0((πœ†π‘š/π‘Ž)𝑑)}, the FB coefficients πΆπ‘š are computed by using the following equation πΆπ‘š=2βˆ«π‘Ž0𝑑π‘₯(𝑑)𝐽0πœ†ξ€·ξ€·π‘šξ€Έπ‘‘ξ€Έ/π‘Žπ‘‘π‘‘π‘Ž2𝐽1ξ€·πœ†π‘šξ€Έξ€»2,(5) with 1β‰€π‘šβ‰€π‘„, where 𝑄 is the order of the FB expansion and, 𝐽1(πœ†π‘š) are the first-order Bessel functions. The FB expansion order 𝑄 must be known a priori. The interval between successive zero-crossings of the Bessel function 𝐽0(πœ†) increases slowly with time and approaches πœ‹ in the limit. If order 𝑄 is unknown, then in order to cover full signal bandwidth, the half of the sampling frequency, 𝑄, must be equal to the length of the signal.

It should be noted that the FB series coefficients πΆπ‘š are unique for a given signal, similarly as the Fourier series coefficients are unique for a given signal. However, unlike the sinusoidal basis functions in the Fourier series, the Bessel functions decay over time. This feature of the Bessel functions makes the FB series expansion suitable for nonstationary signals.

The mean frequency is calculated as in [14] 𝑓mean=βˆ‘π‘„π‘š=1π‘“π‘šπΈπ‘šβˆ‘π‘„π‘š=1πΈπ‘š,(6) where πΈπ‘š=𝐢2π‘šπ‘Ž22𝐽1ξ€·πœ†π‘šξ€Έξ€»2𝑓=(energyatorderπ‘š),π‘š=πœ†π‘š2πœ‹π‘Ž=(frequencyatorderπ‘š).(7) The selection of the optimum window size π‘Ž is required for a good resolution. A larger window provides a finer resolution in frequency, which also means that a greater number of FB coefficients will be needed to cover the same signal bandwidth. The mean frequency of the IMFs was computed using FB expansion. Mean frequency represents the centroid of the spectrum, and thus characterizes the frequency components of the intrinsic mode functions of the EEG signal.

4. Results and Discussion

The use of EMD before calculating mean frequency was necessary owing to the nonstationary and nonlinear nature of the EEG signals. Mean frequency (MF) estimation was performed using the Fourier-Bessel expansion method.

In this section, the MF estimate of the IMFs has been considered as a feature in discriminating EEG patterns in intracranial EEG data recordings. For this purpose, EEG recordings having seizure-free intervals and seizures are considered. The EEG dataset which is available online in [15] is used. In this section, a short description is given and please refer to [15] for further detail. The dataset includes single channel EEG data from healthy and epileptic subjects. The data has five subsets denoted as A, B, C, D, and E, each containing 100 single-channel recordings, each recording of 23.6 seconds in duration. The subsets A and B have been acquired using surface EEG recordings of five healthy volunteers with eyes open and closed, respectively. The signals in two sets have been measured in seizure-free intervals from five patients in epileptic zone (subset D) and from hyppocampal formation of opposite hemisphere of the brain (subset C). Finally, the subset E contains seizure activity selected from all recording sites exhibiting ictal activity. The subsets A and B have been recorded extracranially. The sampling frequency of the data is 173.61 Hz. From the dataset, the subsets C, D, and E have been selected because these have been acquired intracranially. The subsets C and D are combined to form one class and subset E forms the other class.

The MF values have been estimated for both the classes using intrinsic mode functions. The value of MF is small in seizure intervals when compared with that for seizure-free intervals. The class discriminating ability of MF feature is quantified using Kruskal-Wallis statistical test. The MF values are significantly different among the two classes of EEG signals (𝑝<0.01). The results are shown in Figure 2 for the first four intrinsic mode functions. The result suggests that MFs are effective for discriminating seizure from seizure-free intervals for intracranial EEG recordings.

The use of empirical mode decomposition in the present study was justified on the basis of the lack of stationarity in EEG signals. In this way, the entire signal could be analysed simultaneously, at least for a given frequency band. However, although EMD does decompose a signal into different frequency bands (IMFs), the interpretation of the results derived from these IMFs is problematic. Unlike other methods such as wavelets, the number of bands (IMF) is dependent on the frequency content of the signal being analysed, with the number of IMF varying for each signal. Such a property is problematic, as the IMFs used for comparison might relate to different frequency bands. For instance, IMF4 for two subjects might be the fourth of four IMFs for one subject, but the fourth of eight IMFs for another. In such an example, the former would represent the low-frequency components in the signal, while the latter would represent frequency components from the middle of the spectrum.

5. Conclusion

The use of EMD to decompose signals into IMFs is a promising method. However, caution should be exercised when interpreting results from individual IMFs. It would be of interest to develop a method to standardise the comparison of individual or summed IMFs, in order to better compare seizure and seizure-free intervals in the EEG signals.

To establish the clinical use for the proposed technique, it is necessary to test on out-of-sample datasets. It may require the collection of a very large database of recordings of sufficient duration (many hours).