Time-frequency analysis of electroencephalogram series. II. Gabor and wavelet transforms

In this paper we compare two methods, based on the Gabor and wavelet transforms, to quantify and visualize the time evolution of frequency contents of electroencephalogram (cid:126) EEG (cid:33) time series. We found an optimal correlation between EEG visual inspection and the proposed methods in the characterization of the frequency and energy content of characteristic activity during an epileptic seizure. The quasimonofrequency behavior observed in the epileptic EEG series, in a previous work using a Gabor analysis (cid:64) J. Inst. Electr. Eng. 93 , 429 (cid:126) 1946 (cid:33)(cid:35) , is conﬁrmed with the analysis using a wavelet. Moreover, the method based on the wavelet transform allows us to build a detector of epileptic events. Both methods are exempliﬁed with EEG series obtained with depth electrodes in refractory epileptic patients. (cid:64) S1063-651X (cid:126) 96 (cid:33) 01111-7


I. INTRODUCTION
The traditional electroencephalogram ͑EEG͒ tracing is now interpreted in much the same way it was done 50 years ago. More channels are used now and much more is known about clinical implication of the waves, but the basic EEG display and quantification of it are quite similar to its predecessors. The clinical interpretation of EEG records is made by a complex process of visual pattern recognition and the association with external and evident characteristics of the disease ͑clinical symptomatology͒. In past years the use of the Fourier transform, with the introduction of personal computers, has been generalized. The analysis of EEG signals always involves the queries of quantification, i.e., the ability to state objective data in numerical and/or graphic forms that simplify the analysis of long EEG time series. Without such measures, EEG appraisal remains subjective and can hardly lead to logical systematization ͓1-4͔.
The EEG is a complex signal whose statistical properties depend on both space and time. Regarding the temporal characteristics, no-stationarity EEG signals are everchanging; nevertheless, they can be analytically subdivided into representative epochs where the stationarity hypothesis is accomplished ͓5͔. Spectral decomposition of the EEG by computing the Fourier transform has been used since the very early days of electroencephalography. The rhythmic nature of many EEG activities leads naturally to this analysis. Fourier transform allows separation of various rhythms and estimation of their frequencies independently of each other, a difficult task to perform visually if several rhythmic activities occur simultaneously. Spectral analysis can also quantify the amount of activity in a frequency band.
The spectral analysis of EEG signals has proved to be quite useful in comparing short samples of data ͑usually 1-4 sec segment of digitized data͒ ͓6,7͔ from patients against age-matched normative values, as well as in sleep stage analysis, and quantification of drugs, metabolic effects, and various diseased states ͓6͔. However, important information about peak timing is lost. The inclusion of the time evolution in the quantification of the EEG series is an open problem. The morphology and topography of sharp transients have been correlated with seizure type and therapeutic responses to different medications and surgery. An essential component of the traditional visual interpretation of the clinical EEG is the characterization of unfrequent, morphologically variable transient events, especially those associated with the epileptic seizures ͑''spikes,'' ''spikes and waves,'' etc.͒ ͓1-4,6͔. Accordingly, a great deal of energy has been spent over the years in efforts to automatically search long recordings for these phenomena and epileptic transient detection, but with different results ͓6-9͔. The most diffused quantitative method in the clinical practice is the spectral analysis together with a visual assessment ͓1-4,6͔.
When working in the frequency domain it is useful to divide EEG activities into three different categories ͓7͔: ͑i͒ spontaneous activity or background, ͑ii͒ irregular characteristic epileptic activity ͑paroxysm͒, and ͑iii͒ activity evoked by external sensory stimulation. Consequently, it is quite obvious that in the frequency domain representation, rhythmic components are relatively enhanced at corresponding frequencies, whereas transients ͑for example, epileptic spikes, isolated paroxysm, etc.͒, are smeared over the spectrum and therefore are no longer recognizable. From this it follows that the principal field of spectral analysis is the background activity, which means the first category mentioned above, whereas in the other two categories there exist only special cases to which standard spectral analysis can be successfully applied ͓6,7͔.
In the present work we compare two different techniques that allows us to make an analysis of the EEG time series in a time-frequency space. The first is based on the Gabor transform ͓10͔ and the second on the wavelet transform ͓11-14͔. The Gabor transform is a windowed Fourier transform: a window in time is used to localize the frequencies. The window is a time function whose values are nonzero only in a finite-time interval. In particular, the Gabor transform uses a Gaussian function as the window ͓10͔. Based on the Gabor transform, recently, we introduced a procedure to character-ize and visualize the changes in the EEG signal with time ͓15͔. In particular, a good dynamical characterization of an epileptic seizure was found.
A disadvantage of the use of the Gabor transform is that the size of the window is fixed. This implies that the election of the window size is a compromise between the range of frequencies that one can expect to find and the extension in time of the phenomena to be analyzed. In fact, when looking for high frequencies a small window in time should be used and the information on low frequencies could be lost. On the other hand, when looking for slow oscillations or low frequencies a long window in time should be used, and in this case the time evolution of high frequencies events will not be accurate.
The wavelet transform was recently introduced in order to overcome the mentioned drawbacks ͓11-14͔ and its use in the EEG analysis is quite recent ͓16-19͔. In fact, this transform introduces a different type of time-frequency window: a window with variable size that allows an adaptable analysis. In particular, we can consider the Gabor transform as a special case of the wavelet transform with a fixed window ͓12͔. Thus, using the wavelet transform we can make an analysis of the EEG signal based on a window automatically adjustable to the characteristics of the events to be studied, in analogy to the analysis performed with the method based on the Gabor transform. One of the main proposals of the present work will be to check whether there is loss of relevant information using the Gabor transform instead of the wavelet transform. The importance of this is that the method based on the Gabor transform is more easily interpreted by the team of physicians by its analogies with the Fourier transform. In particular, we will confirm our previous results about the quasimonofrequency behavior during an epileptic seizure as well as paroxysmal activity. In addition, the method based on the wavelet transform lets us localize in time, in a precise way, isolated spikes and spike waves, as well as other epileptic transients. This can be used as an automatic transient detector.
This paper is organized in the following way. In Sec. II the experimental setup and the clinical data are presented. In Sec. III the time-frequency methods based on the Gabor and wavelets transforms are introduced. In Sec. IV results using both methods are compared and discussed for intracraneal epileptic EEG series. Finally, a summary is given in Sec. V.

II. EXPERIMENTAL SETUP AND CLINICAL DATA
For each patient the strategy for the use of implanted electrodes is planned in relation to the spatial and temporal organization of the epileptic discharges. This information is simultaneously correlated with clinical symptomatology. Here we present EEG signals from two patients, who were selected because they present different seizure activities. These two patients, whose EEG time series analysis are shown here, were explored with eight multilead electrodes, each one 2 mm long and 1.5 mm apart. The analysis of before and during epileptic seizure data is accomplished by visual analysis of the EEG record. Each signal was amplified and filtered using a 1-40 Hz bandpass filter. A four-pole Butterworth filter was used as a low-pass filter, serving as an antialiasing scheme. The EEG signals were digitalized at 256 Hz through a 10-bit analog-to-digital converter.
According to the visual assessment of the EEG seizure recording, patient I presented an epileptogenic area in the hippocampus with immediate propagation to the girus cingular and the supplementary motor area, on the left hemisphere. In Fig. 1 ͑signal I͒, we display the EEG signal for 16 sec, corresponding to a depth electrode in the hippocampus. In this sample, we can see an isolate paroxysm at 2 sec and a paroxysmal activity that starts around the 4 sec and finishes around 13.5 sec.
Patient II, according to the visual assessment of the EEG seizure recording, presented an epileptogenic area in the left amygdala with propagation to the contralateral amygdala. In Fig. 2 ͑signal II͒, we show the EEG signal for 40 sec, corresponding to a depth electrode in the left amygdala. In this sample, the epileptic seizure starts around the 0 sec and finishes around 34 sec.
The use of depth electrodes provides records where the noise and artifact contamination effects ͑usually present in the EEG series obtained with scalp electrodes͒ are minimized. The applicability of the proposed methods is not restricted to the use of this kind of EEG records.

A. Method based on Gabor transform
In a previous work ͓15͔ we introduce a time-frequency analysis that takes as a basic element the Gabor transform ͓10͔. We performed the Gabor transform of the EEG signal denoted by S(t) as

͑1͒
We used for g D (t) a slide Gaussian window with width D and the asterisk denotes complex conjugation. According to this algorithm, one-dimensional signals are represented in a combined time-frequency space ͓15͔. These functions are situated on a lattice in this combined space, with clearances t 0 and 0 in the time and frequency axes, respectively. In the present analysis the slide Gaussian window has a velocity displacement of 64 data points and width Dϭ2 sec. Thus the resolution in the time-frequency space was ⌬ϭ0.125 Hz and ⌬tϭ0.25 sec.
We defined the time evolution of the spectral frequency content and we denoted by B (i) this function for the frequency band i, that is, for the frequency interval ( min (i) , max (i) ), as then for the i band the power spectral intensity as a time function will be and consequently the total spectral power intensity is where the spectral intensity content is defined in the frequency interval ͑Ϫϱ,ϱ͒. We also define the power spectral intensity per band relative to the total intensity as FIG. 3. Power spectral intensity per band relative to the total intensity, as a function of time for the EEG signal shown in Fig. 1.

͑5͒
For the subsequent analysis of the EEG signal we define for the different bands the mean weight frequency value as Finally, we defined the main peak frequency in the i band at time t, M (i) , as the frequency value for which B (i) takes its maximum value in the interval ( min (i) , max (i) ): The Fourier spectrum will be represented by only one sharp peak at one frequency when there is a monofrequency signal. For this case, if we evaluate the mean weight fre-quency and the main peak frequency, they both will be the same. Therefore, when (t) is approximately equal to M (t) during an appreciable time interval in some band, we shall say that there is a quasimonofrequency engagement in that band. We stress that in our formalism a signal will be quasimonofrequent in this band if this engagement is observed during a reasonable period. Now we introduce the parameter ⌬ (i) and call it the monofrequency deviation. This parameter, as a function of time, gives us an idea about the periods in which the engagements are relevant: Moreover, in order to compare these new time series, for different bands and channels we normalized each one to its The importance of having introduced these new time series (t), (i) (t), M (i) (t), and ⌬ N (i) (t) is that they allow us to characterize the paroxysmal activity and epileptic seizure as well as its evolution with time by means of quantificable magnitudes that are independent of the signal's morphology. Also, throughout this formalism, valuable dynamical information about the epileptic seizure can be extracted ͓15͔.

B. Method based on wavelet transform
Briefly speaking, a wavelet is a rapidly decreasing oscillation function, for which we can change the scale values in order to match the frequency we are seeking ͓11-14͔. We will perform the EEG analysis using a wavelet with adequate scale values and shifting them in time in order to cover the complete EEG records ͓16-19͔.
The integral wavelet transform of a finite energy signal S(t) is defined by where a,bR are parameters; in particular, a is related to the oscillation frequency and b the localization in time of the wavelet function, respectively. We assume that the function FIG. 5. Signal (S) and signal energy distribution (E j ), as a function of time, for the EEG signal shown in Fig. 1, in the corresponding wavelet resolution levels j. S(t) is integrable square ͑Lebesgue͒, which means that S(t)L 2 ͑R͒. When a function (t)L 2 ͑R͒ is any basic wavelet, it verifies ͓11͔ ͵ Ϫϱ ϱ ͑t͒dtϭ0. ͑10͒ In the present analysis we have used a multiresolution scheme based on the spline wavelet transform, with a discretized version of the integral wavelet transform given by Eq. ͑9͒ ͓20͔. We selected this wavelet because it provides a solution to the problems considered basic in event detection: ͑i͒ optimal localization in the time-frequency domain, ͑ii͒ characterization of the different types of epileptic events, and ͑iii͒ computational efficiency of the algorithm proposed.
The wavelet function (t) ͓20͔ is where (t) is the cubic spline compactly supported ͑scaling function͒:

͑12͒
FIG. 6. Wavelet residual signal (R j ) for the EEG signal shown in Fig. 1, in the corresponding wavelet resolution levels j. ⌺ represents the reconstructed signal by summing all R j .
In order to evaluate the time-frequency localization properties of the proposed analysis, we evaluate the values of the center (m) and radius ͑⌬͒ of the time and frequency windows: where is the Fourier transform of and ʈ ʈ L 2 denotes the norm in L 2 space. From these quantities, the window area in the time-frequency plane is 1 4 ͑ area͒ϭ⌬ ⌬ ϭ0.500 67, ͑14͒ which is almost the optimal value 0.5 ͓12͔. Thus the selection of this wavelet guarantees a good localization in the time-frequency plane. For a given EEG signal S(t), initially represented by its polynomial spline coefficients at zero resolution, the wavelet decomposition can be written as where the numbers d 1 (k),d 2 (k),...,d N (k) are the wavelet coefficients and the sequence ͕c N (k)͖ represents the coarser resolution signal at resolution level N. If this decomposition is carried out over all resolutions levels, the wavelet expansion is obtained. In each level j the series expansion given by Eq. ͑16͒ has the property of complete oscillation ͓12͔, which makes this decomposition useful for event localization.
In each level j, the residual signal will be and contains the part of S(t) with the frequency range in the interval where ⌬t is the sampling time. The original signal can be recovered by S(t)ϭ ͚ j R j (t).
When the family ͕ k, j (t)ϭ(2 Ϫ j tϪk)͖ is an orthonormal basis in L 2 ͑R͒, the concept of energy is linked to the usual notions derived from the Fourier theory. Then the energy function will be the sum of the square of the coefficients of the series expansion given by Eq. ͑16͒,

͑19͒
But the wavelets that we are using in the present analysis belong to the more general class of biorthogonal wavelets ͓20͔. This means that there exists a function (t) called the dual of (t). The dual (t) is itself a wavelet ͓20͔, such that Then the family j,k ϭ(2 Ϫ j tϪk) is called the dual basis of j,k . As a consequence, every signal S(t) can be written as In this ͑biorthogonal͒ case the energy associated with the signal S(t) is given by An algorithm for detecting epileptic events can be developed based on the previously established multiresolution analysis and the energy associated with the EEG signal. In the wavelet multiresolution framework ͑described above͒, it is possible to evaluate the energy corresponding to each level, and this can be used for the detection of the characteristic epileptic events ͓19͔. Since we are using dyadic decomposition of the range of frequencies, from a signal of M samples, we have M /2 j coefficients at level j. In order to obtain an accurate event detection, we distribute uniformly the ''atoms'' of energy ͓each one of the terms of Eq. ͑23͔͒ in each level j along 2 j points. Defining for integers r in the interval (kϪ1)2 j Ͻrрk2 j , the energy in each resolution level jϭ1,...,N will be and the energy in each sampled time rϭ1,...,M will be As a consequence, different epileptic events ͑spikes and spike waves͒, isolated or a succession of them, can be characterized through the corresponding values of e j (r) in different resolution levels. Then a detection method can be implemented by looking at when the value e j (r) is greater than a threshold D j appropriately defined for each resolution level ͓19͔.

IV. RESULTS AND DISCUSSION
The results of EEG spectral analysis are often grouped into the traditional frequency bands, whose boundaries can vary a little according to the particular experiment being considered, and they can be adjusted as required ͓1-3͔. On the FIG. 8. Same as Fig. 4 for the EEG signal shown in Fig. 2. other hand, if the signal presents some remarkable frequency characteristics, these must be observed in an independent way whatever the division in bands of the total frequency range. In order to compare the two methods proposed for the analysis of EEG series in the time-frequency space, we chose six frequency bands associated with the resolution levels ͑which correspond to an octave division͒ appropriate for wavelet analysis in the scheme of multiresolution proposed. We denoted these band-resolution levels by B j ͑jϭ1, . . . ,6͒ and their frequency limits are given in Table I. Note that the maximum frequency considered in the present analysis is 32 Hz, in agreement with the bandpass filter used in the acquisition of the EEG signal.

A. Signal I
In Figs. 3 and 4 we display the power spectral intensity per band relative to the total intensity ( j) and the normalized monofrequency deviation ⌬ N ( j) for the EEG signal display in Fig. 1. In Figs. 5 and 6 we show the energy distribution for the same signal at the different resolution levels E j and the corresponding residual signal R j according with the wavelet transform analysis.
Looking at the time evolution of ( j) and E j , good agreement between the temporal changes in these quantities and the signal morphology can be established. Comparing Figs. 3 and 5, it can be seen that both methods present similar energy-band-level distributions.
The greatest energy ͑intensity͒ contribution corresponds to B 2 -B 5 . In both analysis we observe an almost constant contribution in time for B 5 . On the other hand, B 4 has important contributions between 2-9 and 12-14 sec. The most significative contribution for B 2 and B 3 appears between 4 and 13 sec and 8 and 12 sec, respectively.
As we mentioned above, for small values of ⌬ N during a reasonable time interval, a quasimonofrequency band behavior will evolve ͓15͔. In Figs. 4͑c͒ and 4͑d͒, we can see a long engagement for B 3 and B 5 between 3 and 9 sec and 2 and 10 sec, respectively. Shorter engagements can be observed in B 2 , B 4 , and B 6 ͓see Figs. 4͑b͒, 4͑d͒, and 4͑f͔͒. This quasimonofrequency behavior also can be observed by the wavelet analysis throughout the corresponding residual signals in each resolution level ͑see Fig. 6͒. Note that now the quasimonofrequency behavior is associated with almost repetitive patterns in the corresponding time intervals. In particular, this can be verified by performing a Fourier transform for the corresponding portion of the residual signal. In addition, these wavelet residual signals can show engagements that are not seen with the Gabor transform based method. This missing of engagements is associated with the fact that the Gabor transform has a fixed window and certainly can be solved if we run the described procedure for different window widths. An example of a missed engagement can be observed in B 3 between 8 and 12 sec in the corresponding wavelet residual signal ͑see Fig. 6͒, which is not observed with the running of the Gabor based method with a window width Dϭ2 sec ͓see Fig. 4͑e͔͒.

B. Signal II
The signal shown in Fig. 2 has a richer morphology than that shown in Fig. 1; in particular, bursts between 26 and 32 sec are observed. As for the previous signal, the results for the method based on the Gabor transform are displayed in Figs. 7 and 8 and the wavelet transform in Figs. 9 and 10, respectively. From Figs. 7 and 9 we can conclude that a similar global behavior is extracted by both methods. This behavior is characterized by an alternation in the dominant contribution of the different band-resolution levels. We stress that temporal resolution of the bursts ͓Figs. 9͑a͒ and 9͑b͔͒ is better in the wavelet energy diagram than in the Gabor diagram ͑Fig. 7͒. In particular, in the Gabor analysis, the bursts are observed as a global increase in the fast frequencies contribution, in B 4 between 24 and 32 sec.
From Fig. 8 it can be seen that the more important quasimonofrequency engagements are in B 2 , B 3 , and B 4 ͓Figs. 8͑a͒, 8͑b͒, and 8͑c͒, respectively͒. Short engagements can be observed in the other bands. In particular, in B 2 , these engagements correspond to bursts in the signal. Again, this monofrequency behavior can be related to a repetitive morphology of the wavelet residuals.
One of the most important characteristics of the wavelet transform is the great accuracy for the resolution of events in time; as an example of this, the isolate paroxism around 2 sec in signal I ͑Fig. 1͒ and the bursts in signal II ͑Fig. 2͒ can be detected as a high contribution at the corresponding time in the energy resolution levels 3 and 2, respectively. These characteristics can be used as a tool for automatic event detection.

V. CONCLUSION
In this work we compared two methods that provide a good global description of the EEG signals in the timefrequency space. The wavelet transform gives a more accurate temporal localization as well as a good detection of short events; in particular, this characteristic can be used as a tool for automatic event detection. The method based on the Gabor transform gives a global ''average'' description in a compact way for analyzing the EEG signal. In particular, this method is a useful tool in the analysis of long-time EEG series. When a more accurate description of the signal is required the wavelet transform must be used.
We want to stress that the time series generated by both methods analyzed here provide quantifiable and objective information about the frequency content and their relative intensities ͑energy͒ present in each interval of the EEG signal. In this way hidden information can be made evident. For example, a fast activity that is evident in a first visual inspection could be modulated by a low frequency that is not easy to detect in the EEG trace. This potential hidden information can be appreciated in a better way by analyzing the quasimonofrequency engagement parameter ⌬ N ͑the Gabor transform method͒ or wavelet residuals R. This information has a great physiological relevance for physicians because this allows them to characterize and understand in a better way the underlying dynamics of the different epileptic seizures and it can be useful to identify the possible global pacemakers. The use of these two proposed time-frequency analyses, together with the clinical patient history and the visual assessment of the EEG by an expert physician, can contribute to a better diagnosis of epilepsy.