Automatic detection of epileptiform activity by single-level wavelet analysis

Clinical Neurophysiology 110 (1999) 239–249

Automatic detection of epileptiform activity by single-level wavelet analysis Flavio Sartoretto a ,*, Mario Ermani b a

Dipartimento di Matematica Applicata e Informatica, Universita ‘Ca’ Foscari’ di Venezia, Via Torino 155, 30171 Mestre VE, Italy b Dipartimento di Scienze Neurologiche e Psichiatriche, Universita di Padova, Via Giustiniani 5, 35128 Padova PD, Italy Accepted for publication: 10 July 1998

Abstract We describe a new strategy to automatically identify epileptiform activity in EEG. Our scheme is based upon detecting epileptic spikes, via multiresolution analysis, a relatively new tool in signal processing, which allows for dramatic improvements in the efficiency of basic wavelet analysis. We perform a single-level analysis, which is fast and delivers satisfactory results, provided a wise strategy is adopted. Key points are: the identification of suitable wavelets, in order to gain high computational efficiency; the recognition of a proper resolution level; the computation of an appropriate dynamic threshold, in order to pick out the pathological events. Using a suitable wavelet as the model of a threshold-event proved to be a good choice for devising an algorithm which efficiently performs automatic analysis at high-sensitivity levels. The proposed algorithm was implemented into a C++ multiplatform code having an user-friendly interface, which runs on generalpurpose PCs. Results obtained on a set of test tracings, show that the sensitivity of the automatic analysis can be as high as 96%, while less than 5% of the overall recording time is marked. The computational complexity of our algorithm is O (N). Its highly efficient implementation allows for the analysis of up to 310 s of 8 channel EEG, by spending one mere CPU second on a standard PC.  1999 Elsevier Science Ireland Ltd. All rights reserved. Keywords: Automatic EEG analysis; Wavelet analysis; Seizure detection; Spike detection; Computational efficiency

1. Introduction Performing EEG analysis by human readers is a timeconsuming, high-cost activity, thus devising programs for semi-automatic analysis started a number of years ago. When treating epileptic pathologies, the development of dynamical, mobile long-term EEG recording techniques has allowed to detect sporadic events and to evaluate the percentage of critical episodes during all-day activity. On the other hand, such long EEG recordings have enlarged the amount of data to be analyzed, thus increasing the demand for automatic analysis devices (Burr and Stefan, 1987). Unfortunately, dynamical recordings report everyday activity undertaken by the patient, and noise, as well as nonpathological events which resemble pathological ones, are recorded. Physiological and instrumental artifacts, changes * Corresponding author. Fax: +39-049-827-5995; e-mail: [email protected]

in the physiological state of the subject, produce non-stationary EEG events which are easily mistaken by an automatic device for clinically relevant events. Many approaches were adopted for automatic EEG analysis. The recognition algorithms documented in the literature can be roughly divided into three classes: (i) statistical analysis; (ii) mathematical analysis; and (iii) neural network approach. Statistical analysis assumes that the EEG is a basic stationary stochastic process, summed up with nonstationary ones. Unpredictable events, like spikes, are detected by parameter analysis. Some results obtained in analyzing epileptic patterns have been reported (Lopes da Silva et al., 1976; Faure et al., 1980). The classification of non-stationary events is problematic. Unfortunately, its high computational cost discourages use of this approach. Mathematical analysis relies upon derivatives, integrals, the Fourier transform and other mathematical tools to identify peculiar portions of the EEG signal, viewed as a mathematical function of time. This is the so-called mimetic

0013-4694/99/$ - see front matter  1999 Elsevier Science Ireland Ltd. All rights reserved. PII: S00 13-4694(98)001 16-3

CLINPH 97803

240

F. Sartoretto, M. Ermani / Clinical Neurophysiology 110 (1999) 239–249

approach, because it mimics the human approach to the detection of peculiar structures. A well-known automatic procedure was proposed in Gotman and Gloor (1976) where a number of mathematical features is studied. The method was evaluated and it was shown that automatic and visual analysis are comparable, but automatic methods allow for more consistent interpretations (Gotman et al., 1978; Gotman, 1980; Gotman and Wang, 1991; Gotman and Wang, 1992). Alternative mimetic approaches are presented in Ehrenberg and Penry (1976), where a good agreement between computer analysis and human inspection is recorded. A low efficiency was reported using general purpose microcomputers when analyzing long-term recordings. The neural network approach is a powerful tool used for attacking a large variety of problems, such as pattern recognition and graph matching. EEG analysis has recently started to be considered. Preliminary results are presented in Gabor and Seyal (1992), and a deeper insight is given in Webber et al. (1994) where the problem of wisely selecting the input data to the network is analyzed. The results seem to be encouraging, but the efficiency of the method is really poor on standard PCs, demanding implementation on adhoc hardware to provide the fast analysis required by practical applications. The strategy adopted in the present paper is of the second type. It can be viewed as a frequency analysis approach. It goes beyond the use of Fourier transforms, exploiting the concepts of wavelets and multiresolution analysis (Chui, 1992; Daubechies, 1992; Kaiser, 1994). These are relatively new tools for signal processing, which provide time-frequency resolution by decomposing the signal at different scales. While Fourier analysis is not apt for non-repetitive body signals, wavelet transforms allow for processing such signals. In Schiff et al. (1994), applications of these concepts to the analysis of surrogate EEG data are reported, exploiting a cluster-detection technique used in astronomy. The potential capabilities of multiresolution in analyzing a variety of EEG transients are shown in Clark et al. (1995). The wavelet transform has also been recently proposed for discrimination in spike superposition (Zouridakis and Tam, 1997), and its applicability for real-time applications has been suggested. These results show that specific problems in EEG analysis are worth attacking by multiresolution. In D’Attellis et al. (1997) such procedure is applied to the detection of epileptic events. The approach is the same as that exploited in the present paper, relying upon the concept of ‘energy’ of the signal, as provided by multiresolution analysis. Their approach is computationally more expensive than the one adopted here. Moreover, no detailed analysis of the ‘threshold’ computation, which is a key point in the detection algorithm, is provided. The real life performances of the detection procedure when analyzing clinical records are not recorded. In this paper, we report the real-life applicability of single-level analysis, in a multiresolution context, for reducing the amount of EEG data to be analyzed by

human readers when dealing with epileptiform activity. Single level analysis will be shown as an extremely efficient tool in identifying epileptiform events by pointing out epileptic spikes, both in paroxysmal and detached intercritical activities. Detailed results on a set of test recordings show that nearly all the epileptiform activity is pointed out. Misinterpretations due to artifacts are observed, nevertheless a small percentage of the overall recording is picked out by the automatic analysis. The performance of the codes for EEG analysis is not exhaustively measured in the available literature. The implementation of our algorithm has been carefully performed in order to attain high efficiency. We show that running our procedure on a general-purpose computer allows the analysis of an 8 channel EEG tracing, each channel being 128 Hz-sampled, at up to 310 times the real-time.

2. Materials and methods 2.1. Wavelets and multiresolution analysis In principle, spikes and spike-and-wave detection in EEG can be attacked by time-frequency analysis. The Fourier transform (FT) is a commonly used tool to perform frequency analysis, but it is not suitable to provide resolution in time, so the Windowed Fourier transform (WFT) is preferred (Daubechies, 1992). Unfortunately, FT and WFT are not suitable to analyze non-repetitive signals, like EEG ones. The continuous wavelet transform can be thought of as a generalization of the WFT, which eliminates many drawbacks of the latter (Daubechies, 1992), and allows for the analysis of EEG signals (Schiff et al., 1994; Clark et al., 1995). The Fourier transform is a projection of the signal on the dilations of exp(it), while the CWT is obtained by projecting over the dilations and translations, w((t − b)/a), of a function w(x), called mother wavelet. In time-frequency analysis only positive scaling is considered, so a . 0 is assumed in the sequel. In order to devise efficient algorithms to compute wavelet decompositions and/or reconstructions of a discrete signal, like a sampled EEG, the parameters a and b must belong to a discrete set of values. Usually a = am0 , m [ Z, a0 . 1, is chosen. The resulting wavelets are defined as − m=2

wm, n (x) = a0

w(a0− m (x − nb0 am 0 ))

m,n [ Z. Note that in the sequel, a0 = 2, b0 = 1 is set, being the most usual choice. The corresponding discrete wavelet transform (DWT) is fˆ w (m, n) =

…

+∞ −∞

f (t)wm, n (t)dt

Orthonormal wavelets are computationally the most efficient wavelets for signal processing. Multiresolution analysis is a natural framework for studying and constructing orthonormal wavelet bases. The idea underlying multireso-

F. Sartoretto, M. Ermani / Clinical Neurophysiology 110 (1999) 239–249

lution analysis (Mallat, 1989) is to decompose L2(R) into nested subspaces …V2 ⊂ V1 ⊂ V0 ⊂ V − 1 ⊂ V − 2 ⊂ … which satisfy

[ V =L

2

j

j[Z

\ V = {0}

…R†,

j

j[Z

Let Pj be the orthogonal projection onto Vj. For a given f [ L2(R), Pjf is called the approximation at resolution 2 − j or at resolution level j. Let Wj be the orthogonal complement of Vj in Vj − 1. Let Qj be the orthogonal projector onto Wj. Given a space Vj − 1 in a multiresolution analysis n = {Vk, k [ Z}, any f [ Vj − 1 admits an orthogonal decomposition into an approximated signal Pjf and a detail signal Qjf. This is the key point of multiresolution analysis. A signal f = P0f can be analyzed, taking advantage of a given multiresolution n, f and w being its scaling function and a mother wavelet, respectively. One can obtain (Daubechies, 1992) (j − 1) c(j) k = 〈f , fj, k 〉 = ∑ hn − 2k 〈f , fj − 1, n 〉 = ∑ hn − 2k cn n

n

(1)

where hn = , f, f − 1,n .. Assume that the sequence h = (hk)k [ Z is available and we want to perform the wavelet decomposition of the input signal P0f. Using Eq. (1), one can compute the coefficients ck(j) at any level of analysis, j. Thus, the approximation at the resolution 2 − j is obtained by plain sums and products over the coefficients in the approximation at level 2 − j + 1. When h and P0f are finite n-tuples, this strategy is called pyramidal algorithm, because the number of non-zero coefficients decreases going up from one level to the other. This is the strategy which allows us to reach a very high computational efficiency in EEG processing. The same approach is worth using to compute the details Qjf of the signal. It can be proved that (Daubechies, 1992) dk(j) = 〈f , wj, k 〉 = ∑ gn − 2k 〈f , wj − 1, n 〉 = ∑ gn − 2k dn(j − 1) n

n

(2)

where gn = (−1)nh − n + 1. Using this relation, once h and the coefficients ck(j) for appropriate levels are available, the coefficients dk(j) can be efficiently computed. The DWT of the signal is given by the set of coefficients ck(j) and dk(j), k = 0, 1, −1, 2, −2,..., j = 0, 1,.... The wavelet decomposition given by Eqs. (1) and (2) can be viewed as a sequence of filtering operations on the signal (Mallat, 1989; Vetterli and Herley, 1992). From the point of view of signal processing, a scaling function is a low-pass filter, while a wavelet is a band-pass filter (Mallat, 1989; Vetterli and Herley, 1992). High computational efficiency is easily attained when scaling functions and wavelets with compact support are chosen, since in this case convolutions (1) and (2) apply on a finite number of terms. In D’Attellis et al. (1997),

241

biorthogonal wavelets are exploited for EEG analysis, which belong to a more general class than orthogonal wavelets. Note that orthogonal wavelets provide faster analysis, thus we focused on finding orthonormal compactly supported wavelets, as biorthogonal ones are not proven as having any greater ability in detecting epileptic signals. The orthogonal wavelets that we considered have no analytical representation. They are characterized by their own set of coefficients, hn. The smaller is the number of such coefficients, the faster is the computation of the DWT via Eq. (1) and Eq. (2). Like in D’Attellis et al. (1997), we first considered one typical family of approximation functions and attached wavelets, introduced by I. Daubechies, which will be here denoted by Nf and Nw, respectively (they are the functions Nf and Nw plotted in Fig. 6.3 in Daubechies (1992), par. 6.4). Recall that a wavelet can be considered as a band-pass filter. A good wavelet should approximately fit an ideal band-pass filter. The wavelets of the previous family does not match this requirement well, unless N is large, which implies a high computational cost for the analysis. Daubechies proposed families of wavelets which better comply an ideal filter, also when they are characterized by a small number of coefficients. Fig. 1 shows a scaling function gf3 , and wavelet gw3 (after Daubechies (1992), par. 6.4), which were obtained in order to chase this feature. The wavelet gf3 has 8 coefficients hn, n = 0,…, 7. which are reported in Table 1. We also considered other orthogonal wavelets, proposed in the literature. Vetterli and Herley (1992) built a scaling function Vf, characterized by 22 hn coefficients, and attached wavelet Vw, which provide good time and frequency localization. We performed EEG analysis exploiting several types of wavelets. Extensive testing lead us to focus on 5w, gw3 , and Vw wavelets. The wavelet gw3 has nearly the same frequency localization features as 5w, (see Fig. 2), the latter having 10 coefficients hn thus being computationally less profitable. Vw gives better frequency localization, but the ensuing decomposition is not efficient, being characterized by far more coefficients. In accordance with these considerations, gw3 turned out to provide the best trade-off between accurate time-frequency localization and profitable computational efficiency. Further reasons, which are explained in Section 2.2, reinforce choosing this wavelet for spike analysis. 2.2. EEG analysis EEG recordings include a number of channels. Our approach is based upon a channel-wise analysis. Let (i) (i) (i) f(i) s = (fs,1,fs,2...fs,Ns) be the sampled, finite-length signal on the i-th channel. Assume f(i)(t) = Skf(i) s,kfr,k(t), for a suitable dilation level r, is the theoretical continuous signal. We rescale the computations, so that level r is the ground level 0. Our recording device had a sampling rate of 128 Hz on each

242

F. Sartoretto, M. Ermani / Clinical Neurophysiology 110 (1999) 239–249

Fig. 1. Scaling function gf3 (left) and wavelet gw3 (right).

channel. Assuming optimal sampling at Nyquist frequency was performed, f(i) is supposed to have frequencies in the interval (0–64) Hz. Multiresolution analysis provides its time-frequency decomposition. Computing Pjf, j = 1,2,…, (i.e. going up through the approximation levels) smoothes out some frequencies, while the detail signal, Qjf, retains the dropped components (Mallat, 1989; Vetterli and Herley, 1992). Note that the detail signal represents the difference in neighboring levels of approximations (Mallat, 1989), thus it is well suited to detect non-stationary components in the signal, such as spikes. Starting from this observation, we focused on detecting epileptic spikes, leaving to further studies the analysis of other phenomena such as waves in a spike-and-wave assembly. A level of analysis is expected to be reached where the detail signal includes the frequency components characterizing spikes. For a given function g, ˆ ∩ R≥ 0 be its (positive) frequency band. let B(g) = supp(g) Going up a level by Multiresolution Analysis with dyadic wavelets roughly halves the bandwidth of the signal, while the detail retains the components smoothed out (Mallat, 1989). Thus. starting from B(f(i)) = (0.64). wt ideally obtain B(Pjf(i)) = 2-j(0,64), and B(Qjf(i)) = B(Pj-1f(i))\B(Pjf(i)), j = 1,2,... Table 2 shows the ideal frequency intervals for the first four resolution levels. At level j = 2, we obtain B(Q2f (i)) . (16, 32), which our experimental results suggest to encompass most frequency components of epileptiform spikes. Thus, the detail signal at level j = 2 can be expected to be large when spikes are displayed. Fig. 3 shows an 8 channel EEG lay-out. Seconds 10 and 11 are influenced by multiple-spikes affecting all channels. Let f(2) s be obtained by Table 1

Fig. 2. Fourier spectrum of the wavelet gw3 (solid line) and those of 5w (dotted line) and Vw (dashed-dotted line).

sampling channel 2 at a 128 Hz rate. Fig. 4 reports the absolute values of the detail signals, Qjf (2), for j = 1–5, when gf3 and gw3 are exploited to perform signal processing. Peaks can be seen at the resolution level 2 in the time slice where multiple spikes were recorded, confirming that a close connection between spike presence and a large detail signal at level 2 is to be expected. These considerations suggest that a threshold is to be computed, so that when the detail signal at a suitable level j exceeds this threshold, an epileptic spike is likely to have occurred. A single-level analysis is thus chased. Indeed, different types of epileptic events can be detected by considering the DWT at different levels of analysis. We focused on a single level for detecting only spikes. Using a single-level analysis seems to reduce the need for multiresolution analysis. Though considering more levels is actually exploited in D’Attellis et al. (1997), the one-level approach allows for a very fast analysis yet providing a good skimming of the tracings (see Section 2.3), provided a suitable filter is chosen, which amounts to picking up a suitable wavelet for the analysis. The latter proved to be the winning strategy in the present analysis. The computation of the threshold remains the most delicate step in the procedure. In D’Attellis et al. (1997) a sketch of such a step is given. Due to its importance in the whole detecting algorithm, our approach in the sequel is described in detail.

The coefficients hn for gf3 n

hn

Table 2

0 1 2 3 4 5 6 7

−0.1135418 −0.0343775 0.51308565 0.77863292 0.32390632 −0.0911275 −0.0163434 0.05397888

Ideal frequency intervals (Hz), analyzed at resolution levels j = 1,2,3,4, when a 128 Hz sampled signal is considered Level j

P jf

Qjf

1 2 3 4

0–32 0–16 0–8 0–4

32–64 16–32 8–16 4–8

F. Sartoretto, M. Ermani / Clinical Neurophysiology 110 (1999) 239–249

Fig. 3. Eight-channel EEG signal. Sampling rate was 128 Hz.

Assume we have a suitable multiresolution analysis with attached scaling function f and an appropriate model event wavelet w. We exploit wj,k(t), as the sample model threshold event at level j, centered at t = kaj0 = k2j. Define Dwj = maxtwj,k(t) − mintwj,k(t). Recall that w(t) is assumed to be (continuous) compactly supported, so Dwj is welldefined. Moreover, Dwj is independent of k, since translations do not affect the values of the maxima and minima. The sample threshold-event wj,k is standardized by taking Sj, k =

wj, k Dwj

(3)

which can be thought of as a normalized threshold-event. The normalized threshold event is independent of the EEG lay-out under analysis, while the characteristics of the EEG signals dramatically change when different subjects are considered. To tune the model event to the EEG under analysis, we proceed in the following manner. Define the mean stationary jump for a given EEG signal as the length of the usual 95% probability interval DE = 2.0 × 1.96j, j being the standard deviation of the raw data. Tuning the threshold-event for a given EEG recording, amounts to scale the normalized signal given by Eq. (3) by multiplying to the jump DE. Note that we normalize the threshold by considering the signal by itself, rather than its DWT, as in D’Attellis et al. (1997). This approach can allow an easier calibration. The adapted threshold-event is then defined by Sj,a k =

243

(1995) note that one reason which makes automatic EEG analysis a difficult task to perform is a lack of ‘... features that reflects the relevant information’. Our results suggest that, for spike detection, a good feature can be an appropriate wavelet, which in some sense resembles an epileptic spike, and moreover provides good time-frequency localization, yet supplies a computationally efficient decomposition of the signal. The wavelet gw3 has these characteristics. Combining its satisfactory features in frequency decomposition shown in Section 2.2, with the fact that it can mildly resemble an epileptic spike, this latter wavelet seems a suitable choice. The results shown in the sequel were obtained by exploiting gw3 . The threshold Ta is likely to allow the identification of all events modeled by a suitable wavelet w, but recall that in our case we are interested in epileptic events. They are characterized by high energies, so a scaling parameter ws ≥ 1 must be used, in order to improve the selectivity of the automatic analysis. The static threshold T = ws ⋅ T a

(4)

is thus obtained. Fig. 5 shows the static threshold, obtained at level j = 2, by setting ws = 2. Note that not all spikes have details which are larger than the threshold, only epileptic spikes are likely to have enough energy to be marked out. To analyze long-lasting recordings, the static threshold must be periodically recomputed, because the signal behavior changes with time. The values m and j can be affected by a 100% variation during a 5 s time interval, and even by a 400% variation in a 2 min interval (Gotman and Gloor, 1976; p. 516). By experiments, we found that as previously suggested in Gotman and Gloor (1976), computing the threshold each Dt = 5 s, a tuning to brain activity can be obtained, yet avoiding too much recalculation. Every Dt seconds a new dynamic threshold T¯n is computed using Eq. (4). The evaluation is based upon the raw data contained in the interval I, encompassing the last Dt analyzed seconds. To avoid the overestimation of erratic sharp variations, the

wj, k DE: Dwj

…

Its DWT, at level j, is a Sˆ W (j, n) =

+∞ −∞

Sj,a k (t)wj, n (t)dt =

DE d Dwj k, n

The adapted threshold Ta = DE/Dwj is thus defined. Peculiar events can be identified in time by choosing the appropriate resolution level, j, and picking out those time values ¯t = n¯ aj0 = n¯ 2j when lfˆ W (j, n¯ )l ≥ T a holds true for the DWT of the signal, ˆfw(m,n). Summarizing, our main idea is to compute the threshold by choosing a model threshold event, whose DWT will provide the threshold value. Clark et al.

Fig. 4. Multiresolution analysis of channel 2 in the layout reported in Fig. 3. The absolute value of the detail signal is shown. Five resolution levels were computed.

244

F. Sartoretto, M. Ermani / Clinical Neurophysiology 110 (1999) 239–249

Fig. 5. (a) EEG layout showing epileptic spikes. (b) Absolute value of the detail signal at resolution level 2. The dashed line represents the threshold computed by setting ws = 2.

old threshold, To, is reconsidered, the final dynamic threshold, Tf, being computed by the formula Tf = (1 − pm )To + pm T¯ n The memory parameter pm , 1, is intended to properly weight the importance of To. By experiments, we found that a suitable setting is pm = 0.75. On the other hand, we discovered that when spikes have been detected in I, Tf is likely to be overestimated. In this case, we increase the weight given to To, by lowering pm, to avoid overestimation of spikes on the oncoming analysis. By extensive tests, we found that when spikes have been detected in I, a suitable (0) setting is p(1) m = pm = 0.2, while pm = pm = 0.75 holds elsewhere, as previously suggested. Recall that at level j = 2, the detail signal Q2f(i) has frequencies in the interval (16, 32) Hz. Dyadic wavelets provide frequency analysis whose bandwidth is split by a factor of two at each level, which can be too coarse, as observed in Clark et al. (1995), p. 381. A splitting trick (Daubechies, 1992; par. 10.5) proves useful to improve the detection quality. This device is a special case of wavelet packets (Chui, 1992; par. 7.4). The idea is to perform a further wavelet decomposition on the just decomposed signal. By applying the splitting trick, the frequency interval (16, 32) Hz, is ideally halved into the bands (16, 24) and (24, 32) Hz (Daubechies, 1992). This corresponds to adjunctive lowpass and high-pass filtering, respectively. Using low-pass filtering, high frequency components can be filtered out,

usually reducing the number of non-epileptic spikes incorrectly detected as epileptic ones. In practice, the frequency separation is not so precise as in theory, so our code allows the user to choose either low-pass or high-pass adjunctive filtering, by the same wavelet which produced the initial decomposition. Our experiments suggest that low-pass filtering is to be preferred, so this choice was done in our tests. EEG recordings list several channels. As previously quoted, wavelet decomposition and spike identification is performed channel-wise. All computations are performed separately on each channel. In some cases, however, high energy epileptic events are chased, which surely affect more than a single channel. The user-setting of the least number of channels, Cn, was enabled, thus an event affecting less than Cn channels is not taken into consideration. The results shown in the following section were obtained setting Cn = 1, because we are interested in detecting both generalized and focal epileptiform activity. The main task of this work is to provide an automatic skimming of EEG tracings, thus the adopted testing procedure is as follows. A time slice starting at a given recording second, whose length is a whole number of seconds, is called a time interval (TI). Each epileptic event or set of epileptic events which are separated by less than 1 s gaps, are identified by the smallest encompassing TI, which is called an epileptiform interval (EI). At the same time, when an event is pointed out by the automatic analysis, we say that the smallest TI including that event was marked. An epileptiform interval, I, is considered to be detected (by automatic analysis) when at least one TI, say J, was marked by our code, so that either the set difference I \ J is non-void or I is contiguous to J. In other words, the marked TI either overlaps or touches the EI. The idea behind this criterion is that, given either a portion of an EI, or a contiguous TI, a human reader can easily re-collect the entire interval. The performance of our automatic procedure has been tested essentially on the ground of the number of detected epileptiform intervals. From the computational point of view, recall that a discrete WFT, using fast Fourier transform techniques has a computational cost which is O(NslogNs), Ns being the number of samples in the finite signal (Wickerhauser, 1994). A key point is that to detect epileptic spikes, a given number of multiresolution analysis levels, which depends only on the sampling rate of the digital signal has to be performed. The computational cost of such a procedure is O(Ns), as can be inferred from Eqs. (1) and (2). This result suggests that high efficiency can be attained by exploiting our approach. We shall see in Section 2.3 that wavelet analysis is practically very efficient, requiring a small number of levels to be computed. 2.3. Materials An appropriate set of long-term EEG layouts, recorded by an 8-channel MEDILOG 9000 system have been processed.

F. Sartoretto, M. Ermani / Clinical Neurophysiology 110 (1999) 239–249

Available EEG devices can count 16 or more channels, but in this paper we analyze the problem of channel-wise event identification, not the localization of foci, thus the number of available channels is not crucial. Each tracing was obtained by a variable length continuous recording. The montage consisted of 8 bipolar channels, placed according to the 10/20 international system: F2–F8, F1–F7, F8–T4, F7–T3, T4–T6, T3–T5, T6–O2, T5–O1. Analog signals were converted into digital ones by an A/D converter device with a 128 Hz sampling rate on each channel. The digital signals were stored on magnetic devices and later analyzed by a physician who marked the EI values. Highly expensive human analysis is needed as a reference basis, thus we selected a reasonable amount of representative test tracings. Moreover, the automatic analyses were performed after downloading analog MEDILOG 9000 tracings into PC files. This operation is highly time-consuming, practically preventing the automatic analysis of entire 24-h tracings. Portions with different durations were selected, in order to accomplish with the different durations of the embedded relevant events inside them (i.e. epileptic events, as well as artifacts and other events which can be mistaken for epileptic activity). EEG tracings from 10 epileptic patients have been considered. They are representative of the clinically relevant activity in terms of morphology (spikes, polyspikes, classical 3 Hz spike-andwave complexes, slow spike-and-wave complexes, polyspike-and-wave complexes), spatial distribution (focal and generalized), and discharge duration (ranging from less than 1 to 12 s). Several types of artifacts are also included. Note that a number of tracings were selected, where a portion of epileptic events is embedded into artifacts. In the sequel, each fragment pertaining to the i-th subject is labeled by EEGiX, where ‘X’ is an uppercase letter. Multiresolution analysis was applied to our set of EEG tracings, by exploiting the wavelet gw3 , the level of analysis being j = 2. A summary of the values finally chosen for the computing parameters follows: recalculation time interval Dt = 5 s, scaling factor qs = 2, memory parameters p(0) m = 0.75, p(1) m = 0.2. Least number of channels was Cn = 1. An additional splitting trick were performed, using low-pass filtering, via gw3 . Level of analysis, recalculation time interval and low-pass filtering were chosen according to the theoretical considerations and suggestions from literature, which were reported in the previous sections. The scaling factor and memory parameter were set by the following procedure. We considered the first EEG fragment coming from patient 1. We identified the parameter values which allow satisfying sensitivity and selectivity levels. Then we analyzed the second fragment from patient 1, changing the parameter values, if necessary, so that the automatic analysis yields valid sensitivity and selectivity levels for both tracings. This procedure was carried out for all tracings coming from patients 1, 4, 7 and 10. These tracings are representative of all relevant types of events and physiological conditions. The final parameter values were used in the

245

analysis of all our tracings. An average performance was obtained, which is not turned on a particular patient. We developed a complete multiplatform code which implements the algorithm described in Section 2.2. Our code allows for easy window-guided analysis of EEG tracings. The tracing can be inspected at any time on a window and after the analysis was performed, stripes are superimposed to the signal, showing the events detected by the program. The width of the stripes depends upon the computed level of DWT analysis. A palette is given for setting the user parameters, such as the level of DWT analysis, threshold recomputing time, etc.

3. Results Tables 3– 5 summarize the agreement between our code and the human reader. The first column reports the label of each fragment, the second column shows the recording daytime, the third column its length, in seconds. The EIs marked by the human reader, are reported in the fourth column, the TIs pointed out by the program can be found in the fifth one. An integer number of seconds ‘m’ reported in columns 4 and 5, denotes the time interval m − 1 , t ≤ m while ‘m − n’ denotes the interval m − 1 , t ≤ n. Inside column 5, the detected (in the sense defined in the previous section) epileptiform intervals are underlined. The last column shows the classification of the phenomena which are peculiar to the EEG fragment.

4. Discussion Concerning the ability of our code in detecting epileptiform spikes, let us see some sample layout. Fig. 6 shows the output of an automatic analysis performed on a portion of tracing EEG4B. The user-friendly interface supplied into our code easily allows for the visual inspection of the results. Dark stripes show the time slices pointed out by the program. An interictal SW activity, mainly located in the frontotemporal regions, can be observed. Our automatic analysis properly detects the epileptiform activity. Fig. 7 shows a clearer activity, starting with multiple spikes. Also, in this case our automatic analysis was able to point out the epileptiform activity. Note that, unlike other approaches, for example that exploited in Ermani et al. (1990), our code selects both detached inter-critical activities, such as single spikes and multiple spikes. Fig. 8 shows the performance of our code when epileptic activity is mixed with noise. Epileptic activity mainly affecting channels 2, 3, and 4 is correctly detected. Let us analyze the test results which were summarized in the previous section. Inspecting Tables 3 and 4, and Table 5, one observes that practically all EIs were detected. The overall sensitivity of the automatic analysis, i.e. the percen-

246

F. Sartoretto, M. Ermani / Clinical Neurophysiology 110 (1999) 239–249

Table 3 Human-code accordance. The first column shows the label given to the tracing. The second column reports the day-time of the recording. The length of the fragment is in column three. The EIs marked by the human reader are shown in the fourth column, the TI pointed out by the automatic analysis are in the fifth one. Epileptiform TIs have been underscored. ‘–’ indicates that no TI was marked. The last column shows a classification of some peculiar phenomena. Abbreviations: SW, spike-and-waves; SD, spike discharges; SWD, spike-and-wave discharges; PSWD, polyspike-and-wave discharges File

Time

Length

Human

Program

Classification

EEG1A EEG1B EEG1C EEG1D EEG2A EEG2B EEG2C EEG2D EEG2E EEG2F EEG2G EEG2H EEG3A EEG3B EEG3C EEG3D EEG3E

16:30 17:16 17:22 17:15 18:01 19:15 20:01 00:20 01:40 01:49 02:09 02:21 17:15 17:19 17:37 18:20 20:25

280 72 248 96 56 72 72 72 56 144 160 56 96 64 88 96 160

44–45, 210–211 22–28 25–29,33–38,71–85, 89–95, 217–222 32–34, 81–82 28 36–38 – 41, 43 32–34 39, 112 31, 66, 143 21 23–28 33–36 – – –

34, 44–45, 210–211 21–25, 45, 48 27–30, 34–36, 74–85, 90–95, 218-219 6–7, 9, 33, 50, 82, 96 28 37–38 19, 27, 38, 39–40 41, 43 32–33 39, 112 31, 66, 143 20, 47, 66 23–24, 96 33–35 20, 47, 66 24, 25, 26, 49, 60 2, 84–89, 92, 93

Short SD Long PSWD SD and PSWD Short SD, noise Short SD Short SD Noise 2 SW PSWD Short SD Short SD Short SD SWD PSWD Noise Noise Chew artifacts

tage of EIs marked out of their total number, was 46/ 48 = 96%. This seems a satisfactory percentage, which compares well with those of the methods reported in the literature. In a number of cases, non-epileptic spikes were pointed out. Inspecting Table 3, we see that tracings containing chew artifacts (e.g. EEG3E) together with tracings displaying noise (e.g. EEG3D) can be misinterpreted by automatic analysis. The marked events have some affinity with epileptiform activity, in some cases being confusing also for a human reader. The total selectivity, i.e. the percentage of detected EIs out of the total number of marked TIs is 46/ 123 = 37%, which can not be rated very good. Selectivity can be improved by raising the threshold, but at the expense of losing some sensitivity, and/or by integrating more levels

of DWT analysis. On the other hand, sensitivity is more important than selectivity in clinical applications, as the recordings must eventually be analyzed by a human reader. Note that finding the parameter values (threshold, recalculation interval, least number of channels, etc.); which provide the best performance or whichever tracing one could encounter, is not feasible. A large number of non-epileptic events was marked by the program, but this drawback proved to be not so serious in practice, as our code marked only 216 s over 4760 of total recording, i.e. less than 5%, thus dramatically reducing the amount of data. Table 6 shows the high efficiency of our code on some computers running under two operative systems (OS). The processing speed, SEEG, is reported, being the number of 8channel input tracing seconds, analyzed by spending a sin-

Table 4 Human-code accordance (see Table 3) File

Time

Length

Human

Program

Classification

EEG4A EEG4B EEG4C EEG4D EEG5A

13:48 01:45 18:56 04:25 08:15

232 168 176 80 616

00:44 23:59 13:40 14:11 09:55 10:05 10:30

416 168 72 48 272 64 176

112–120, 172, 174, 203–208 18, 52–53, 90 6, 44, 82–84, 125, 152 42 70, 85, 156–162, 295, 358, 367, 488, 540, 556, 575, 579, 593 – – 4–7, 9–10, 14, 19, 48–50, 60 32, 46 9, 12–13, 64, 91, 143–145, 183, 203–205, 245, 247 22–23, 46 11, 21–23, 57–58, 82–83, 93, 105, 126, 139, 141, 158–161

SWD (2–3 Hz) SWD in sleep Eye artifacts Sleep SD and PSWD

EEG6A EEG6B EEG7A EEG7B EEG8A EEG8B EEG8C

113–124, 203–209 52–53, 90–92 – – 70–71, 156–164, 295–296, 358–359, 467 – – 49–52 – 12, 61, 64 22–23, 45 22–23, 58, 158–162

Sleep Sleep, noise SD, noise Noise PSWD, eye artifacts PSWD PSWD

247

F. Sartoretto, M. Ermani / Clinical Neurophysiology 110 (1999) 239–249 Table 5 Human-code accordance (see Table 3). Asterisks in the third column mean that the recording day-time is not available File

Time

Length

Human

Program

Classification

EEG9A EEG9B EEG9C EEG10A EEG10B EEG10C EEG10D

16:57 17:48 18:02 * * * *

152 80 80 16 20 16 20

29–35 36–42 29–36 6–7 – 10, 12–13 11–12

29–34, 54, 67, 112, 113, 140, 144 36–41 11, 14, 30 6–7 3, 8, 11, 12 10, 12 11–12, 14

SWD SWD SWD PSWD Chew artifacts Short PSWD PSWD

gle CPU second. The speed ranges from 95 to 310, i.e. it is higher than the standard value SEEG = 60 estimated for a trained physician. This is perhaps the most interesting result of our test, which encourages further research in this direction. Such a high performance suggests that to improve sensitivity and/or selectivity, further post-processing multiresolution stages can be added to our automatic procedure, yet in principle maintaining its efficiency comparable with that of a human reader. Although multilevel analysis is expected to improve the identification of complex events, a simpler single-level analysis come forward as primary test for the efficiency of identification by wavelet, and it proves to be very efficient indeed. Another important research goal is the automatic localization of epileptic foci, which demands simultaneous 16 (or more) channel analysis. This is beyond the scope of the present work, being scheduled for future research.

5. Summary A detailed strategy for detecting epileptiform activity in EEG tracings was studied, developed and implemented. Our strategy relies upon the concept of multiresolution analysis proposed by Mallat (1989), which allows for the efficient detection of suspected epileptic spikes. An appropriate choice of the wavelet space and resolution level was performed. The exploitation of a suitable wavelet as a model for an epileptic spike allowed for an effective time resolution of this event. A satisfactory detection level was attained, by carefully computing a dynamic threshold. We implemented our strategy into an efficient multiplatform C++ code, equipped with a user friendly interface. The analysis of a test set of tracings, coming from several patients suffering from severe forms of epilepsy, has been presented. The aptitude of the program in detecting epilepti-

Fig. 6. Detail of the layout EEG4B. Dark stripes point out the time slices marked by the automatic analysis.

248

F. Sartoretto, M. Ermani / Clinical Neurophysiology 110 (1999) 239–249

Fig. 7. Fragment EEG10D. Dark stripes show the time slices pointed out by the program.

form discharges, independently marked out by a physician, was measured. The overall sensitivity of the code turned out to be 96%, which seems a satisfactory result. Although the measured selectivity is not high, our test results suggest that

a good skimming of the tracings is performed by our code, as it marks less than 5% of the global amount of recorded seconds. Running our code on standard PCs, up to 310 s of 8-

Fig. 8. Fragment EEG3B. Epileptic activity mixed with noise.

F. Sartoretto, M. Ermani / Clinical Neurophysiology 110 (1999) 239–249 Table 6 The number of 8-channel, 128-Hz/channel sampled, EEG seconds, SEEG, processed by spending a mere CPU second, is shown for several standard machines, running under two operative systems (OS) Computer

Processor

OS

SEEG

Macintosh 6100 Macintosh 5400 Olidata Green 740 PC Winner

PowerPC 601/66 PowerPC 603/180 486 DX2/66 Pentium 166

Mac OS Mac OS Windows 95 Windows 95

177 280 95 310

channel EEG can be analyzed by a mere CPU second, i.e. far more than an average-skilled human reader. This encouraging result confirms the real-life applicability of the method. Moreover, such high efficiency suggests that sensitivity and selectivity can be increased, yet in principle maintaining a good processing performance. Further work along this research direction is in progress.

Acknowledgements Gianfranco Boccalon efficiently implemented the algorithm. Roberto Rinaldo and Giancarlo Calvagno provided helpful hints. References Burr, W. and Stefan, H. Automatic analysis of spike-wave and background activity in mobile long-term EEG recording. Adv. Epileptol., 1987, 16: 265–272. Chui, C.K. An Introduction to Wavelets. Wavelet Analysis and Its Applications. Academic Press, Boston, MA, 1992. Clark, I., Biscay, R., Echeverria, M. and Virues, T. Multiresolution decomposition of non-stationary EEG signals: a preliminary study. Comput. Biol. Med., 1995, 25 (4): 373–382. D’Attellis, C.E., Isaacson, S.I. and Sirne, R.O. Detection of epileptic events in electroencephalograms using wavelet analysis. Ann. Biomed. Eng., 1997, 25: 286–293.

249

Daubechies, I. Ten Lectures on Wavelets. Society for Industrial and Applied Mathematics, Philadelphia, PA, 1992. Ehrenberg, B.L. and Penry, K. Computer recognition of generalized spikewave discharges. Electroenceph. clin. Neurophysiol., 1976, 41: 25–36. Ermani, M., Pellegrini, A., Tagliapietra, F. and Azzalini, A. Use of discriminant analysis, for automatic spike and wave detection. Electroenceph. clin. Neurophysiol., 1990, 75: S41. Faure, C., Dubuisson, B. and Samson-Dollfus, D. Aide au diagnostic de l’epilepsie. Int. J Bio-Med. Comput., 1980, 11: 21–32. Gabor, A.J. and Seyal, M. Automated interictal EEG spike detection using artificial neural networks. Electroenceph. clin. Neurophysiol., 1992, 83: 271–280. Gotman, J. Quantitative measurements of epileptic spike morphology in the human EEG. Electronceph. clin. Neurophysiol., 1980, 48: 551–557. Gotman, J. and Gloor, P. Automatic recognition and quantification of interictal epileptic activity in the human scalp EEG. Electroenceph. clin. Neurophysiol., 1976, 41: 513–529. Gotman, J. and Wang, L.-Y. State-dependent spike detection: concepts and preliminary results. Electroenceph. clin. Neurophysiol., 1991, 79: 11– 19. Gotman, J. and Wang, L.-Y. State dependent spike detection: validation. Electroenceph clin. Neurophysiol., 1992, 83: 12–18. Gotman, J., Gloor, P. and Schaul, N. Comparison of traditional reading of the EEG and automatic recognition of interictal epileptic activity. Electroenceph. clin. Neurophysiol., 1978, 44: 48–60. Kaiser, G. A Friendly Guide to Wavelets. Birkhauser, Berlin, 1994. Lopes da Silva, F.H., Van Hulten, K., Lommen, J.G., van Leewen, S., Van Veelan, C.W.M. and Vliegenthart, W. Automatic detection and localization of epileptic foci. Electroenceph. clin. Neurophysiol., 1976, 43: 1– 13. Mallat, S.G. A theory for multiresolution signal decomposition: the wavelet representation. IEEE Trans. Pattern Anal. Mach. Intell., 1989, 11 (7): 674–693. Schiff, S.J., Aldroubi, A., Unser, M. and Sato, S. Fast wavelet transformation of EEG. Electroenceph. clin. Neurophysiol., 1994, 91: 442–455. Vetterli, M. and Herley, C. Wavelets and filter banks: Theory and design. IEEE Trans Signal Proc., 1992, 40 (9): 2207–2232. Webber, W.R.S., Litt, B., Wilson, K. and Lesser, R.P. Practical detection of epileptiform discharges (EDs) in the EEG using an artificial neural network: a comparison of raw and parametrized data. Electroenceph. clin. Neurophysiol., 1994, 91: 194–204. Wickerhauser, M.V. Adapted Wavelet Analysis from Theory to Software. IEEE Press, Piscataway, NJ, 1994. Zouridakis, G. and Tam, D.C. Multi-unit spike discrimination using wavelet transforms. Comput. Biol. Med., 1997, 27 (1): 9–18.