Cyril Touze

Antoine Chaigne

Lyapunov Exponents from Experimental Time Series: Application to Cymbal Vibrations

Lyapunov exponents are among the most relevant and most informative invariants for detecting and quantifying chaos in a dynamical system. This method is applied here to the analysis of cymbal vibrations. The advantage of using a quadratic fit for determining the Jacobian of the dynamics is presented. In addition, the interest of using a time step for the evolution of the neighbourhood not equal to the timelag used for the reconstruction of the phase space is underlined. The robustness of the algorithm used yields a high degree of confidence in the characterization and in the quantification of the chaotic state.

To illustrate these features in the case of cymbal vibrations, transitions from quasiperiodicity to chaos are exhibited. The quasiperiodic state of the system is characterized together by the power spectrum of the experimental signal and by calculation of the Lyapunov spectrum.

Introduction

Cymbals belong to the category of non linear percussion instruments. In their normal use, these instruments are struck by a mallet. Pairs of cymbals can also be struck against each other, especially in bands or in orchestral use. Due to the rather heavy strokes, the vibrations of cymbals generally exhibit large amplitude, compared to their thickness, which leads to nonlinear effects. As a consequence, the pitch of such instruments is not clearly defined. Only at the end of the decay, as the magnitude of the motion becomes sufficiently small, one can clearly hear the most salient eigenfrequencies of the normal modes predicted by the linear theory. Previous investigations on cymbals are mostly concerned with the linear regime [ 1], and only a few references address the question of nonlinearity which is of prime importance for these instruments [START_REF] Wilbur | Subharmonic generation in cymbals at large amplitude[END_REF][START_REF] Fletcher | Nonlinear dynamics and chaos in musical instruments[END_REF]. In order to tackle this problem, one strategy consists in investigating the underlying physics of the vibrating structure. This method has been applied to gongs in the past [START_REF] Fletcher | Nonlinear frequency shifts in quasisphericalcap shells: pitch glide in chinese gongs[END_REF][START_REF] Legge | Nonlinearity, chaos, and the sound of shallow gongs[END_REF]. One complementary approach is to investigate the dynamical behaviour of these instruments with the help of nonlinear signal processing tools. This is the purpose of the present paper.

The main difficulty for the analysis of struck cymbals follows from the broadband excitation of the initial pulse. This pulse excites simultaneously a very large number of frequencies, in a very short amount of time, which makes it almost impossible to analyse the successive bifurcation mechanisms and the routes to chaos. In addition, the magnitude and location of the pulse is generally hard to reproduce. Therefore, it seems more appealing to drive the cymbal with slowly and continuously increasing magnitude, so that transitions from the linear to the chaotic state can be observed, and to limit the spectrum of the excitation. In practice, the cymbal is driven here sinusoidally at frequencies close to the most salient linear eigenfrequencies of the instrument, or intentionally apart from these eigenfrequencies [START_REF] Touze | Analysis of cymbal vibration using nonlinear signal processing tools[END_REF]. The velocity of one selected point of the instrument is recorded during about three minutes (see Figure 5). At large amplitude, it is interesting to notice that the sound obtained by driving the cymbal with a sinusoid instead of a pulse is clearly comparable to characteristic cymbal sounds, although the tone quality may be somewhat poorer. This is a strong argument in favor of considering the typical sounds of the instrument in ordinary performances as a combination of the nonlinear effects due to a limited number of active modes.

In this paper, particular attention is paid to the calculation of the Lyapunov exponents from experimental time series. The method used here for computing the Lyapunov spectrum relies on an idea conjointly developped by Eckmann and Ruelle [START_REF] Eckmann | Ergodic theory of chaos and strange attractors[END_REF][START_REF] Eckmann | Liapunov exponents from time series[END_REF] and Sano and Sawada [START_REF] Sano | Measurements of the Lyapunov spectrum from chaotic time series[END_REF]. This method consists of approximating the matrix of the linearized flow in the reconstructed tangent space. One major advantage of this method is that, in theory, the complete spectrum of the exponents can be obtained. Other methods, in which the divergence of nearby trajectories is used, are also possible candidates, but their principal.limitation follows from the fact that, in practice, only the maximal exponent can be obtained [START_REF] Wolf | Determining Lyapunov exponents from a time series[END_REF][START_REF] Kantz | A robust method to estimate the maximal Lyapunov exponents of a time series[END_REF]. Other methods, based on a global reconstruction of the vector field and especially designed for short or noisy data sets, are also available [START_REF] Kadtke | Global dynamical equations and Lyapunov exponents from noisy chaotic time series[END_REF][START_REF] Brown | Calculating Lyapunov exponents for short and/or noisy data sets[END_REF]. These methods will not be considered here.

The main features of the algorithm presented in this paper are the following: first, the local Jacobian of the evolution function in the reconstructed tangent space is approximated by a higher-order local polynomial [START_REF] Bryant | Lyapunov exponents from observed time series[END_REF][START_REF] Briggs | An improved method for estimating Liapunov exponents of chaotic time series[END_REF][START_REF] Brown | Computing the Lyapunov spectrum of a dynamical system from an observed time series[END_REF]. As a consequence, it is expected that the accuracy in the determination of the negative exponents will be improved, and that the presence of spurious exponents, multiples of the true exponents, will be avoided. One example of the existence of such spurious exponents was first noted in [START_REF] Eckmann | Liapunov exponents from time series[END_REF]. More recently, an analytic proof has been given by Sauer et al. [START_REF] Sauer | Spurious Lyapunov exponents in attractor reconstruction[END_REF]. Second, a time step for the evolution of the neighbourhood is selected, which is different from the timelag used for the reconstruction of the attractor. This improves significantly the detennination of the exponents. These properties of the algorithm will be illustrated numerically on standard computer-generated time series originated from chaotic discrete and continuous dynamical systems.

Interpretation of the Lyapunov exponents related to real world data is a rather difficult task which should be done with great care. However, a number of successful applications can be mentioned, which cover various fields of physical systems: hydrodynamics [START_REF] Eckmann | Liapunov exponents from time series[END_REF], c hemistry [START_REF] Wolf | Determining Lyapunov exponents from a time series[END_REF], fundamental physics [START_REF] Abarbanel | Analysis of observed chaotic data[END_REF][START_REF] Kantz | Nonlinear time series analysis[END_REF] and acoustics [START_REF] Holzfuss | Liapunov exponents from a time series of acoustic chaos[END_REF][START_REF] Frison | Chaos in ocean ambient noise[END_REF]. including musical acoustics [START_REF] Wilson | Characterizing the clarinet tone: Measurements of Lyapunov exponents, correlation dimension, and unsteadiness[END_REF][START_REF] Lee | Chaos in segments from Korean traditional singing and Western singing[END_REF], and speech [START_REF] Banbrook | How to extract Lyapunov exponents from short and noisy time series[END_REF]. In section 3 of the paper, the computation of the Lyapunov exponents is appbed t.o experimental time series obtained from measurements performed on a thin crash cymbal, which is driven sinusoidally at its center by means of a shaker. Similar experiments were conducted in the past by Fletcher. Rossing et al. [START_REF] Fletcher | Nonlinear dynamics and chaos in musical instruments[END_REF][START_REF] Legge | Nonlinearity, chaos, and the sound of shallow gongs[END_REF][START_REF] Wilbur | Subharmonic generation in cymbals at large amplitude[END_REF]. By increasing slowly the amplitude of the shaker motion, transiti ons from linear to chaotic motion are observed. A complete description of the experimental setup can be found in [6}, together with the results of the algorithms used for estimating the dimensions of the system (correlation dimension, and false nearest neighbours). The calculation of the Lyapunov exponents makes it possible to detect the presence of chaos, and to quantify the "degree of disorder" in the chaotic regime. In addition. this analysis gives us a better understanding on the route to chaos observed in the system, through testing of the quasi periodicity of the measured signal. This point will be developed in section 4.

Calculation of the Lyapunov spectrum

Algorithm

The methods of current use in nonlinear signal processing are based on the Takens reconstruction theorem [START_REF] Takens | Detecting strange attractors in turbulence[END_REF]. Recent surveys present both the theoretical background and applications of these techniques for chaotic signals [START_REF] Abarbanel | Analysis of observed chaotic data[END_REF][START_REF] Kantz | Nonlinear time series analysis[END_REF][START_REF] Ott | Coping with chaos: Analysis of chaotic data and the exploitation of chaotic systems[END_REF]. Therefore, the fundations of this theory will not be developed further in this paper.

In order to compute the Lyapunov spectrum of a given time series, it is necessary to obtain the local Jacobian of the underlying now from the trajectory in the reconstructed phase space. The principle is the following: Given a time series {s(n)}n= l. .. N. a multivariate trajectory is formed with a timelag T and a so-called embedding dimension d5. The value of the timelag is usually taken as the first minimum of the Average Mutual Information function [START_REF] Fraser | lndependant coordinates for strange attractors from mutual information[END_REF]. A lower bound for the embedding dimension can be given by using the test of the false nearest neighbours [START_REF] Kennel | Determining embedding dimension for phase-space reconstruction using a geometrical construction[END_REF]. This yields the vector:

y(k) = [s(k), s(k + T), ... , s (k +(dE -l)T)] .
Denoting F the evolution function which maps the vector y (k) onto the vector y(k + Tp), one can write:

y(k + TF) = F(y(k)), {1)
where TF is the time step for this local neighbourhood-to- neighbourhood mapping. The main task of the algorithm is to compute the local Jacobian matrix DF(k) in the tangent space, since the Lyapunov exponents are computed from the eigenvalues of the Jacobian. Defining further:

• { y{k;) } i=I. .. Nv : the Nv nearest neighbours ofy(k),

• TJ(k, i) = y (k;) -y (k), i ;:; 1 .. . N u: the small dis- placements in the neighbourhood ofy(k),

• TJ(k+Tp,i) = y(k; +Te)y (k +Tp),i = 1 ... N 11 :

the small displacements in the evolved neighbourhood, after T F time steps, the local Jacobian is then calculated via a Taylor expansion in equation (1). The standard method consists in Iinearizing F , but it has been found that higher-order Taylor expansion give better resu lts [START_REF] Brown | Computing the Lyapunov spectrum of a dynamical system from an observed time series[END_REF].

A second-order ex pansion, for example, with i = 1 ... N 11 , yields:

TJ(k + Tp , i) = (2) DF(k)ry(k, i) + ~ry(k, i)tH(k)ry(k, i) where ( 8F•) DF (k);,; = B y; (y(k))
is the Jacobian of the evolution function F , and ( 82F; ) H (k);,;,t = fJyjl)y, (y(k))

is the Hessian of F . Equation ( 2) is then projected onto a subspace of dimension dL 5: de. de-is the local dimension corresponding to the dimension of the manifold on which the dynamics takes place, which should not be mistaken for the embedding dimension dE . The dimension dE should be large enough so as to avoid the presence of false neighbours [START_REF] Abarbanel | Analysis of observed chaotic data[END_REF][START_REF] Kantz | Nonlinear time series analysis[END_REF]. All subsequent calculations are made in a dL-dimension subspace, which provides us with d1, Lyapunov exponents. The next step of th e algorithm consists in calculating DF and H through matrix inversion corresponding to equa- tion (2) (the Hessian part can be written in a matrix form. See, for example. [START_REF] Brown | Computing the Lyapunov spectrum of a dynamical system from an observed time series[END_REF] for more details). The coefficients of H are useful only for assessing the trajectory, since the Lya- punov spectrum is computed from the eigenvalues of DF.

Thus, only the Jacobian part is finally retained.

The last step of the algorithm is the classical QR decomposition of the Jacobian, which is evaluated along the phase space trajectory, in order to calculate the dL Lyapunov ex- ponents [START_REF] Eckmann | Ergodic theory of chaos and strange attractors[END_REF]. With K the number of points on the attractor for which the Jacobian has been estimated, R(k) the R matrix from the QR decomposition at y (k ), and T the sampling time step, the itlt Lyapunov exponent is given" by: This summary of the main steps of the algorithm shows that many parameters are involved: T and dE for the reconstruction, Tp , dr" K and Nv for the fit. For most of them, simple rules are clearly established to avoid pitfalls. A large number of references exist in the literature devoted to non linear signal processing where the rules for selecting appropriately these parameters are specified (see, for example, [START_REF] Abarbanel | Analysis of observed chaotic data[END_REF][START_REF] Kantz | Nonlinear time series analysis[END_REF][START_REF] Grassberger | Nonlinear time sequence analysis[END_REF]).

Higher-order Taylor expansion

Most of the algorithms used for estimating the local Jacobian matri x restrict the Taylor expansion of equation (l) to the fi rst-order (linear) term. In the present study, the use of higher-order expansions has been investigated. The major expected improvement is to increase significanrly the accuracy of the result~, especially for the negative ex ponents for which the linear algorithms present a number of difficulties [START_REF] Briggs | An improved method for estimating Liapunov exponents of chaotic time series[END_REF][START_REF] Brown | Computing the Lyapunov spectrum of a dynamical system from an observed time series[END_REF]. Two simultaneous burdens result from a linear fit on the Jacobian since, first, the trajectory must be approximated and, second, the Lyapunov exponents must be given. Usi ng a higher-order expansion implies that all terms in the matrices of equation ( 2) are involved in the approximation of the local neighbourhood-to-neighbourhood mapping.

The accuracy of the results does not increase monotonically with the degree dp of the expansion. This feature can be explained as follows: with increasing dp. it is necessary to take a larger number of neighbours into account in order to perform the matrix inversion in good conditions. As a consequence, the assumption of small neighbourhood for the Taylor expansion is not verified any more. In other words, the number of selected neighbours Nv for computing the matrix must be large enough to overcome numerical instability, but, at the same ti me, this number is limited in order to keep the neighbourhood as small as possible. A good compromise consists in selecting N v is given by [START_REF] Brown | Computing the Lyapunov spectrum of a dynamical system from an observed time series[END_REF]:

Observations made by computing a third-order expansion of equation (1) show that the accuracy of the results is not significantly improved, compared to the second-order one. Therefore, it has been decided to select a second-order expansion.

Another important fact that j ustifies the use of a secondorder expansion is linked to the problem of spurious exponents. One major problem, when dealing with experimental time series, is that the local dimension of the attractor, and thus the number of Lyapunov exponents, is not known. As a consequence, if the dimension is too large, then the algorithm presented in section 2. 1 generates spurious exponents. With a first-order expansion , for a noise-free case, it is foun d that some of the spurious exponents arc integer multiples of the true ones, as noted in [START_REF] Eckmann | Liapunov exponents from time series[END_REF] and demonstrated in [START_REF] Sauer | Spurious Lyapunov exponents in attractor reconstruction[END_REF].

With a second-order expansion, still for a noise-free case, the presence of tJ1ose spurious exponents which are multiples of the true ones, is avoided. With a linear fi t, a spurious ex ponent equal to twice the positive one appears as the local Table I. Calculated Lyapu nov exponents for the Henon map, with linear and quadratic fits. The length of the sequence is equal to 40000 points, with T = TFTT = 1 for a discrete map, dE dimension dL increases. This phenomenon may lead to the incorrect conclusion that two positive exponents exist in the system. By using a quadratic fit, those false multiple positive exponents are avoided. Despite our efforts, no analytical proof of this property could be given. Indeed, the mathematical derivations leading to the analytic result of [START_REF] Sauer | Spurious Lyapunov exponents in attractor reconstruction[END_REF], when extended to the quadratic case, give a system which bears an implicit solution [START_REF] Touze | Techniques d'ordre superieur pour I' elimination d' exposants de Lyapunov fallacieux[END_REF]. However, a number of numerical calculations performed on different systems (Henon map, logistic map, Lorenz system) show that the quadratic fit eliminates the multiple spurious exponents. This feature is illustrated for the Henon map in Table I.

. Selection of the evolution time Tp

Selecting an appropriate value for Tp is a particularly hard task in the case of continuous-time dynamical systems. Previous algorithms use Tp = T as a general feature [START_REF] Eckmann | Liapunov exponents from time series[END_REF], because it substantially simplifies the implementation. However, there are no other justifications for this choice, and thus it is not recommended to select Tp = T, since there are no relations between these two q uanti ties, as mcntionned in [START_REF] Brown | Computing the Lyapunov spectrum of a dynamical system from an observed time series[END_REF]. The choice of Tp is directly connected to the choice of the discrete map which best fits the continuous How, and no general rules are available for this problem.

However, for too small evolution time Tp, the local contraction or expansion properties of the How arc not well taken into account, thus leadi ng to a bad estimate of the whole spectrum. On the contrary, when Tp approaches T, there is some probability for the trajectory to be too smoothed, hence leading to underestimate the negative exponents.

By studying the effect of Tp on a number of continuous maps, it has been observed that using a too small Tp leads to poor results, and that the precision of the estimated Lyapunov spectrum becomes acceptable as TF 2: T / 2 . This feature is illustrated on a time series generated from the forced Duffing x

+ 8x + x + x 3 = a cos(wt)
with 8 = 0.2, a = 40, w = 1. For this choice of parameters, the Lyapunov exponents are .A1 = 0.11, .A2 = 0, A3 = -0.31 (in s-1 ) [START_REF] Parlitz | Identification of true and spurious Lyapunov exponents from time series[END_REF]. The time series was calculated with initial conditions: x = 1 and ;i; = 0.7. The length of the series is equal to 40000 samples. The first 5000 samples were cut in order to avoid the influence of transients. A fourthorder Runge-Kutta method is applied to the system, with a sampling time step of 0.04 s. For this value, the first minimum of the Average Mutual Information function is obtained for T = 10 samples. In this example, Figure I shows that the best spectra are obtained for Tp 2: T /2. It has been found that this empirical rule remains valid for higher embedding dimensions.

In a second set of experiments, a time series obtained from the first coordinate of the Lorenz system, sampled by a fourth-order Runge-Kutta scheme, has been used. The sampling time step is equal to T = 0.01 s. In this case, the first minimum of the Averaged Mutual Information function is equal to 0.1 s, which corresponds to T = 10 samples. Figure 2 shows that in this case, the best spectrum is obtained

forTp = T/2.
In conclusion, it turns out that even if no general rule can be asserted, all particular cases investigated show that selecting Tp :::; T /2 leads to bad results.

Results on computer-generated time series

In order to check its validity and limits, the algorithm has been tested on computer-generated time series derived from well-known dynamical systems (logistic, Henon and Ikeda maps, Lorenz and Rossler systems, Duffing equation, superposition oftwo incommensurate sinusoids). In each situation, will be presented below. The first example is aimed at illustrating the performance of the method on a chaotic system, and the second example shows that the method used also yields excellent results with a non-chaotic signal.

Figure 3 shows the computed Lyapunov exponents for increasing values of dL, with the time series used in the previous section. The evolution time is here Tp= 6 samples. The convergence ofthe true exponents (marked by dotted lines on Figure 3) with increasing dL is clear, whereas the spurious ones vary with d£. This can help in identifying the spuri-ous exponents. This figure shows that the absolute value of the false exponents is systematically larger than the absolute value of the largest true one. However, this property cannot be generalized. Tests performed on the Lorenz system, for example, show that the spurious exponents obtained in this case are also all negatives, but that they are located between the vanishing exponent and the most negative one. These results are in accordance with those previously obtained by other authors [START_REF] Eckmann | Liapunov exponents from time series[END_REF][START_REF] Brown | Computing the Lyapunov spectrum of a dynamical system from an observed time series[END_REF][START_REF] Holzfuss | Liapunov exponents from a time series of acoustic chaos[END_REF][START_REF] Stoop | Calculation of Lyapunov exponents avoiding spurious elements[END_REF]. From these observations, one is led to conclude that it is not possible to exhibit a general rule for the location of the spurious numerical exponents in the Lyapunov spectrum.

A signal composed of two incommensurate audible frequencies has been also tested. The results show that the algorithm is robust for finding the two-torus on which the trajectory takes place, by exhibiting two clearly vanishing exponents whatever the local dimension d£. The signal is made of the following frequencies: fi = 440Hz and h = 2v'2JI = 1244.5 Hz, with equal amplitudes. A trajectory of 40000 points has been used, with a sampling frequency Fe = 48000 Hz, a standard value for audio signals. The delay T for the reconstruction is equal to 16 samples. Figure 4 shows the convergence of the exponents, for a local dimension dL equal to 3. The signature of a two-torus is assessed by the presence of two zero exponents, while the remaining ones are all negatives.

A commonly accepted method to discriminate the spurious exponents consists in computing the Lyapunov spectrum for a given time series, read first forward and then backward in time [START_REF] Parlitz | Identification of true and spurious Lyapunov exponents from time series[END_REF]. The expected behaviour is that the signs of the true exponents change with time reversal, whereas the spurious ones do not. This method has been found to give excellent results with time series issued from discrete maps (Henon, logistic and Ikeda map has been tested in this manner), thus leading to an undoubted identification of the spurious exponents. Unfortunately, in the case of continuous flows, the expected behaviour is not so clearly visible, and the identification of true and spurious exponents becomes slightly impossible. Thus, this method was not used in the present study.

Application to the cymbal

Experimental set-up

The previously described method for determining Lyapunov exponents has been applied to experimental time series on a Zildjian thin crash cymbal (diameter: 41 cm), in order to assess and quantify the chaotic regime. In this experiment, the cymbal is clamped at its center to an impedance head (B&K 8001) mounted on a LDS vibration exciter. The amplitude of the driving force is slowly and linearly increased by modulating the sinusoidal driving signal by a triangle waveform of very low frequency. The magnitude of the driving force lies within the range 1 to 20 N, depending on the driving frequency. The magnitude of the displacement at the edge is generally between 0.1 and 1 mm. However, this magnitude 0. 02,-------~-----~------- may become very large for some values of the excitation frequency. In ce1tain cases, edge displacements up to 1 cm have been observed. The analyzed signal is delivered by an accelerometer (B&K 4393) glued either at the edge or at a distance of r /2 from the center, where r is the radius of the cymbal, and recorded on a DAT-tape with sampling frequency Fe = 48000 Hz. The duration of the recordings lies between 2 and 5 minutes. One example of recorded signal is shown in Figure 5. Sharp transitions can be seen on this figure, as well as the relaxation effect after suppression of the excitation, around t = 2 min 40 s. In order to observe various types of transitions, the forcing frequency is either close to one linear eigenfrequency of the cymbal, or significantly apart from eigenfrequencies and from low-order linear combinations of these frequencies. In a previous work, calculations of the correlation dimension and of the embedding dimension, by means of the false nearest neighbours method, have been performed on these signals. The results show, among other things, that the dimension of the underlying dynamics is less than or equal to 7 [START_REF] Touze | Analysis of cymbal vibration using nonlinear signal processing tools[END_REF].

As it is generally the case for experimental signals, when compared to computer-generated series, for example, the present signals exhibit a finite signal-to-noise ratio (typically 50 dB). This signal-to-noise ratio is essentially due to the accelerometer. Moreover, the reconstruction does not display a clear fractal form. Since the algorithm makes extensive use of the geometrical properties of the phase portrait (neighbouring trajectories), some diffic ulties may be encountered in the detennination of the exponents. To illustrate this, a two-dimensional reconstruction is shown in Figure 6, for a driving frequency Fexc = 526Hz, after the onset of chaos.

Results

The results presented below were obtained in the large amplitude vibration regime, characterized by a dense and broadband Fourier spectrum. In fact, one primary goal of this study was to assess that the cymbal vibration was really chaotic. Thus, the convergence of a positive Lyapunov exponent was expected. Calculations of the Lyapunov exponents were made for four different driving frequencies and for two different locations of the accelerometer. In each situation, a portion of 40000 points of the signal is selected (which corresponds to a duration of 0.83 s), once the chaotic regime is well established. In Figure 5, for example, a portion of 0.83 s was selected around 2mn30s. At that time, the driving force is typically equal to 10 N.

A typical curve of convergence is shown in Figure 7. This curve has been o btained for particular values of the involved set of parameters. It has been checked that the best convergence is obtained for Tp around T /2 as shown previously on calculated time series. It is observed that the convergence is less pronounced as Tp is equal to T. One can notice, in particular, the convergence of the largest exponent, for increasing values of dL, which clearly indicates the presence of a positive exponent in the system. Thus, it. can be concluded that the vibration is chaotic. One problem arising from the blurred trajectory of the experimental signal is that spurious exponents appear at positive, though relatively small, values.

In order to validate the analysis, the results obtained are checked against theoretical rules given by different authors:

• The Kaplan-Yorke conjecture [START_REF] Kaplan | Chaotic behaviour in multidimensional difference equations[END_REF] yields a dimension dKY from the spectrum of the Lyapunov exponents defined by:

d N L:~1 Ai I<Y = + IAN+ll
where N is such as : L:~1 Ai > 0 and L:~t 1 Ai < 0.

This dimension dJ<y can be used in comparison to the correlation dimension d 2 [START_REF] Grassberger | Measuring the strangeness of strange attractors[END_REF][START_REF] Touze | Analysis of cymbal vibration using nonlinear signal processing tools[END_REF] in order to get confidence in the calculated exponents. This also yields a guideline for the• selection of the local dimension dL .

• Theoretical results show that, when computing the Lyapunov spectrum from a time series derived from a continuous dynamical system, at least one vanishing exponent must be present which represents the continuity of the trajectory. ... g However, applications of theses rules show that one cannot conclude with absolute certainty on the validity of all Lyapunov exponents. Table 11 shows the results obtained for an acceleration signal recorded at the edge of the cymbal, with Fexc = 526 Hz. These results are also shown in Figure 7 for dL E [START_REF] Fletcher | Nonlinear dynamics and chaos in musical instruments[END_REF][START_REF] Eckmann | Liapunov exponents from time series[END_REF]. The correlation dimension for the same signal has already been calculated, which gives d2 = 3.6 [START_REF] Touze | Analysis of cymbal vibration using nonlinear signal processing tools[END_REF]. Comparing this value with dKY calculated with the Lyapunov exponents shows that, for dL 2:. 6, dKY becomes too large compared to d 2 = 3.6. Thus, the best candidates for the spectrum are obtained for dL equal to 4 or 5. At this stage, it cannot be concluded with certainty which of these two spec-Table 11. Calculated Lyapunov exponents for three different values of the local dimension dL. for the cymbal driven at 526 Hz, and with the acceleration measured at the edge. The length of the experimental sequence is equal to 40000 points. The other fixed parameters are: ds = 8. T = 16, TF = 10, K = 2000. The exponents are given here in normalized units with respect to the sampling time, i.e. in sample- Ill shows the results for dL = 4 and 5 , which are here also the best possible candidates. Application of the rule related to the p resence or not of a vanishing exponent in the spectrum could lead to the conclusion that the exponents obtained with dL = 4 are right. However, the correlation dimension of this signal is equal to 3.4, and thus the Lyapunov spectrum calculated with dL = 5 cannot be completely rejected.

••••••••••••e••••• •••••~ 0 ............ •-•••o•••••
In view of the results obtained in the two typical cases presented in Table ll and in Table Ill, one can be confiden t in the estimation of the largest exponent but not in the relevance of the other exponents in the spectrum. For these reasons, it has been decided to limit the present analysis to the convergence of the largest exponent.

T he results obtained for four different driving frequencies are summarized in Table IV, together with the correlation dimension d 2 and the minimum embedding di mension dE to use, given by the false nearest neighbours. This Table shows that a positive Lyapunov exponent has been found in each case, which is certainly the strongest indication that the system really exhibits chaos, since the calculation of Lyapunov exponent is more reliable than the correlation dimension.

-0.08425 3.06 -0.02460 -0.09147 4.02

The precision given in this Table for the largest Lyapunov exponent is just the stan dard deviation calculated on the last I 000 iterations. The calculations are made typically with K E [1 500, 2500] iterations, and the convergence is obtained in less than 1000 iterations. Figure 4 shows an example for a theoretical case (two incommensurate frequencies), where it turns out that 1500 iterations are largely sufficient for the convergence. Finally, it can be noticed in Table IV that the exponent<: are of the same order of magnitude for all different driving frequencies.

Route to chaos

Ruelle-Takens theory

In this section, attention is paid to the characterization of experimental transition scenarii from quasiperiodicity to chaos observed in the vibrations of the cymbal submitted to sinusoidal driving force with increasing amplitude. Calculation of the Lyapunov spectrum during the transition, as well as power spectrum analysis, are used for identifying the route to chaos observed. T his section starts with a brief survey on the Ruelle-Takens theory. The interested reader is referred to dedicated references for a complete description of the most classical routes to chaos and of their general properties (see, for example, [START_REF] Schuster | Deterministic chaos. 3rd augmented edition[END_REF][START_REF] Manneville | Structures dissipatives, chaos et turbulence[END_REF]). Specific routes to chaos (such as period doubling, quasiperiodiciry, in termittencies, crises, .. . ) have been extensively studied in the past, on both theoretical and experimental viewpoints. T hese routes show, in particular, that the underlying mechanisms are governed b y some universal laws. In this context, the Ruelle-Takens theory postulates that a limited number of Hopf bifurcations are involved in the transition from a stationary state to a chaotic state, for any physical system [START_REF] Ruelle | On the nature of turbulence[END_REF][START_REF] Newhouse | Occurrence of axiom-A attractors near quasi-periodic flows on Tm, m 2: 3. Communications in[END_REF]. In other words, this theory establishes that, if a system exhibits three Hopf bifurcations, then there is high probability for this system to be chaotic with a strange attractor. In physical systems, the following characteristic features are thus expected, as the control parameter of the system increases (in our case, the magnitude of the driving force): Initially, The power spectrum of the observed variable (here, the acceleration of a given point on the cymbal) shows one peak at a given frequency !I (plus, eventually, some harmonics). The phase portrait will be a simple closed curve, theLyapunov spectrum will display a single zero exponent, the other ones being negatives. Then, after a first Hopf bifurcati on, the Fourier spectrum is made of two frequencies fi and fz, and, in addition, of all combinations f = md 1 + m2!z. with (m 1 , mz) E Z 2 . The phase space trajectory is located on a two-torus. At this stage, one should make a clear distinction between two cases. The first case corresponds to mode locking which occurs when the ratio between the two frequencies ft and fz is rational. As a consequence, the trajectory on the torus will close itself after a finite number of cycles. In the second situation, the ratio between these two frequencies is irrational: the frequencies are incommensurate, and thus the torus is densely covered by the trajectory which never closes itself. In this case, the Lyapunov spectrum exhibits clearly two zero exponents. When the third Hopf bifurcation is about to appear, the phase portrait shows typical foldings, indicating that the two-torus will be destroyed, and hence that the system will gain an additional dimension [START_REF] Berge | L'ordre dans le chaos[END_REF][START_REF] Curry | A transition from Hopfbifuraction to chaos[END_REF]. If another Hopf bifurcation arises (i.e. if a third frequency fa is to appear in the power spectrum), then a "broad-band noise", which is characteristic of chaos, is observed, and the torus breaks up in favor of a strange attractor. It will be shown in the next section that the experimental signals recorded on the cymbal show, to variable extent, all the abovementioned characteristics.

A clear survey of the transition scenarii can be found in [START_REF] Schuster | Deterministic chaos. 3rd augmented edition[END_REF]. Details on the quasiperiodicity routes to chaos are studied in [START_REF] Jensen | Transition to chaos by interaction of resonances in dissipative systems. I. circle maps[END_REF][START_REF] Bohr | Transition to chaos by interaction of resonances in dissipative systems. II. Josephson junctions, charge-density waves, and standard maps[END_REF].

Mode-locking and quasiperiodicity

Figure 8 shows the phase portrait and the power spectrum of the acceleration signal recorded on the cymbal at r /2, at three successive instants of time. The driving frequency is F exc = 440 Hz. In the first raw of the figure, the vibratory motion of sponding to the continuity of the flow), the other ones being negatives. By increasing the amplitude of the driving force, a Hopf bifurcation occurs. Therefore, a second frequency and combinations between Fexc and this new frequency become visible. The motion is quasi periodic. The phase portrait displays typical foldings indicating that the trajectory is located on a two-torus which is about to be destroyed. With increasing amplitude, a second bifurcation is observed, and a broadband spectrum indicates the onset of chaos, as it can be seen on the phase portrait also.

A zoom on the low-frequency part of the power spectrum in the quasiperiodic state is shown in Figure 9. The identification of a quasiperiodic motion can be made here by noticing that all frequency peaks, in this part of the spectrum, can be written as the combination f = m 1 h + m2h, with h = F exc :::: 440 Hz, and another fixed frequency h :::: 698 Hz. This frequency h corresponds here to one particular eigenfrequency of the cymbal. This result is in accordance with one possible interpretation of the Ruelle-Takens scenario presented in [START_REF] Manneville | Structures dissipatives, chaos et turbulence[END_REF]: the Hopf bifurcation occurs when one mode of the system is destabilized and becomes active (ie the real part of the eigenvalue of this eigenmode becomes positive). The vibratory motion is then governed by two incommensurate frequencies. Calculations of the Lyapunov exponents were made j ust before the onset of chaos. The Lyapunov spectrum clearly displays two vanishing exponents (see Table V). This indicates that the trajectory is located on a manifold topologically equivalent to a two-torus and confirms that the transition to chaos is obtained by the break-up of a two-torus. If another eigenfrequency becomes active with increasing amplitude, then the motion becomes chaotic, as shown in the last raw of Figure 8, and as has been shown by the presence of one positive Lyapunov exponent (see section 3).

Similar observations and calculations were made for other driving frequencies, for which there is no evident relationship with the eigenfrequencies of the cymbal (Fexc = 412, 526 and 385 Hz have been tested). Moreover, experiments at a fixed level of excitation (which is sufficient to reach the quasi periodic state), and with slight variations of the excitation frequency, have been conducted. Hysteresis cycles were clearly observed, wi th different resonance curves, depending on F exc is increased or decreased.

Finally, some experiments were conducted with a driving frequency equal to twice one particular eigenfrequency of the cymbal. In this case, as expected, the following behaviour is observed: a mode-locking occurs between Fe:xc and the eigenfrequency defined by fa. = Fexc/2. This result can be seen in Figure I 0, where the driving frequency is Fe:xc ::::248Hz, and the eigenfrequency corresponding to the (4,0) mode of the cymbal is 124 Hz. An energy transfer from the low-frequency to the high-frequency range is clearly visible in the power spectrum shown in figure I 0. Here, the motion is phase-locked in a ratio 2: I. As the amplitude ofthe forcing frequency increases, the system suddenly changes to chaotic behaviour which is characterized by a positive Lyapunov exponent, as it has been shown in Section 3 (,\ = 0.3 ± 0 .1 (ms)-1 . This value was not reported in table IV because the recorded signal was to short to yield sufficient accuracy to be extremely confident in the result). This result suggests that the system is strucLUrally unstable after the second Hopf Bifurcation, and that no quasiperiodicity involving three different frequencies occurs. With increasing amplitude of the driving force, the two-torus directly breaks up in favor of a strange attractor.

In conclusion, the whole set of analysis tools used for characterizing the transition scenario (Fourier and Lyapunov spectrum, phase portrait) indicates that the main features of a -4

~-60 • ~ -80 IL.I!U.J....J...I_J.....J _ _ _ _ _ .....:..J --4 -2 0 2 500 1 000 1500 2000

S{l)

x to'

Frequency {Hz)

Ruelle-Takens scenario are present. However, at this stage of the work, the number of analyzed signals is limited, and thus no g(o6a( (aws sucn as tne Farey tree at tne onset of cnaos [START_REF] Jensen | Transition to chaos by interaction of resonances in dissipative systems. I. circle maps[END_REF]421. can be clearly exhibited.

Conclusion

In this paper, a study on the vibratory motion of a cymbal has been presented, using nonlinear signal processing tools. Special emphasis has been put on the computation of the Lyapunov exponents from experimental time series, since it is one of the most re liable and most relevant invariant in order to detect and quantify chaos. In this context, an algorithm for estimating the Lyapunov spectrum has been developed. Specific features were used for improving the quality of the estimation. First, a second-order Taylor expansion has been used in order to improve the accuracy with which the exponents are given. It has been shown numerically that this quadratic fit avoids the presence of spurious exponents which are multiples of the true ones. Second, a time step for the evolution of the neighbourhood, which is different from the time step used for the reconstruction, has been selected. This also improves the accuracy of the results. T he ability of the method has been clearly shown on theoretical time series. In the case of experimental time series, only the largest Lyapunov exponent has been clearly exhibited, like it is most often the case in adverse experimental context.

The cakulations of the Lyapunov exponents were made before the chaotic regime. This yields additional arguments for the identification of the observed transition scenario. The exponent~. together with the observation of spectral peaks combinations, allows c lear confidence in the fact that the system is subjected to a transition from quasiperiodicity to chaos. Depending o n the value of the driving frequency, mo<fe-(ocKing and quasiperiodic states have been observed, as well as the succession of two Hopf bifurcations from the Linear to the chaotic vibration.

The following general features of the phenomena have been underlined. First, a mode-locking is observed in the case of a driving frequency equal to twice the value of one particular cigenfrequency of the cymbal. Second, a quasiperiodicity is observed for other drivi ng frequencies. This quasiperiodicity is associated to the combination of spectral peaks and to typical foldings of the trajectory in the phase space. In this case, the presence of two vanishing exponents in the Lyapunov spectrum indicates the presence of a two-torus. When the ampli tude of the forcing frequency increases, the two-torus breaks in favor of a strange attractor, whose presence is clearly identified by a positive Lyapunov exponent. This completes the description of the transition scenario for the cymbal, and quantifies the chaotic regime in terms of magnitude of the largest Lyapunov exponent.

= 8 and

 8 K=l200. The apparition of a spurious exponent at twice the true positive one >.1 = 0.418, with the linear fit, is clear. The second theoretical exponent is >.2 = -1.621. With the quadratic fit, the spurious exponents are all negatives. 41866 -0.3 1485 -0.44684 -0.8641 1 -1.61900

Figure 1 .

 1 Figure 1. Lyapunov equation for the Duffing equation, as a function of the evolution time TF. The parameters are K = 2000, T = 10, dE= 6, and dL = 3. 40000 points have been used for the calculation, with a quadratic fit. It can be seen that the acceptable values for the evolution time are obtained for TF :::: T /2. The true exponent are marked by dotted lines.

0

  -~: :::::::::~: :::::::::~: :::::::::~:.::•::::::~: :.::.::•:~: ::::::•::~: :::::::::~: :::::::::~: :::::::::~• --~3~~4--~5--~6--~7--~8--~9--~1~0 Evolution time Tp

Figure 2 .Figure 3 .

 23 Figure2. Lyapunov exponents for the Lorenz system, as a function of the evolution time TF. The optimal value for TF is T /2. The calculation has been made with 40000 points from the first coordinate of the Lorenz system, with K = 1500, dE = 4, dL = 3, and T =10. The true exponents are marked by dotted lines. A quadratic fit has been used. The influence of TF is particularly visible on the negative exponent (marked by a(+)), whereas the positive one and the vanishing one((*) and (o)) do not vary significantly with TF.

Figure 4 .Figure 5 .

 45 Figure 4. Convergence of the Lyapunov exponents with the number K of iterations for a signal composed by two incommensurate frequencies. The two first exponents are equal to zero, the third one (spurious) is negative, which is consistent with the theoretical results. The fixed parameters are: dE= 6, dL = 3, T = 16, TF= 10.

Figure 6 .

 6 Figure 6. Trajectory in the reconstructed phase space, with T = 16, for nn acceleration signal of the cymbal d riven at 526Hz, once the chaotic regime is established.

7 Figure 7 .

 77 Figure 7. Lyapunov exponents in (sample)-1 , as a function of the local d imension dL, for the acceleration signal recorded at the edge of the cymbal driven with Fexc = 526Hz. 30000 points are used for the analysis. The parameters have the following values: dE = 8, K = 2000, T = 16, TF = 10. A clear convergence to 0.020 is observed for the largest exponent.

Table V . Figure 8 .

 .8 Figure 8. Reconstructed traj ectory and power spectrum at three successive instants of time, showing three typical phases of the transition to chaos. The cymbal is driven at Fexc = 440Hz and the acceleration is recorded at r /2.

  weakly nonlinear. The power spectrum shows a peak at f = Fexc. and some harmonics. The reconstructed trajectory is a simple closed loop. The Lyapunov spectrum exhibits one vanishing exponent (,\ = 0.0 ± w-4 , corre-

Figure 9 .

 9 Figure 9. Power spectrum of the previous signal below 1200 Hz during the quasiperiodic state, with Fexc =440 Hz. All frequency peaks can bewrittenas linear combination of h = Fe:r.c and h =698Hz, where /2 corresponds to one eigenfrequency of the cymbal.

Figure 10 . 2 ~

 102 Figure 10. Reconstructed traj ectory and power spectrum of the acceleration signal at three successive instants of time, with Fe:r;c =248 Hz. The eigenfrcquency of mode( 4,0) is ut 124 liz. It can be seen that the moti on is phase-locked afrer the first Hopf bifurcation . By increasing the amplitude of the driving frequency, the motion becomes chaotic.

Table Ill .

 Ill 1 . Calculated Lyapunov exponents in sample-1 for the cymbal driven at 412Hz, with a time series of 40000 points. The fixed parameters are: ds = 8. T = I 0, TF = 8, K = 2000.

	dL			Lyapunov exponents (in sample -I)			dKY
	4	0.01936	0.00157	-0.00130	-0.0657 1			3.12
	5	0.01979	000349	-0.00334	-0.01989	-0.0927 1		4.02
	6	0.02060	0.00625	-0.00088	-0.00938	-0.03393	-0.10969	4.55
	dL		Lyapunov exponents (in sample -I)		dKY
	4	0.02325	0.00001	-0.01779
	5	0.02457	0.00636	-0.00582

Table IV .

 IV Correlation dimension d 2 , embedding dimension ds and greatest Lyapunov exponents .X. in (ms)-1 , for four different driving frequencies. Each Lyapunov exponent calculation has been made with a time series of 40000 points, using a quadratic fit.

	F eo;c (Hz)	Recording	d2	dE	.X (ms)-1
		Position
	412	r/2	3,4	5	1,15 ± 0.01
	412	edge	2,9	5	0.82 ± 0.02
	440	r/2	3,8	5	I, 38 ± 0.02
	526	edge	3,6	5	0, 94 ± 0.01
	385	edge	4,6	7	0, 75 ± 0.01
	385	r/2	4,5	6	1.42 ± 0.06
	tra is the most valid. A similar discussion can be made for

Fexc = 412 Hz. Table

Acknowledgement

The authors would like to thank Paul Manneville and Olivier Michel for fruitful discussions related to the physics of the system, and to non linear signal processing problems. We also wish to thank David Ruelle for giving favourable consideration to the problem of the spurious exponents which are multiples of the true ones, and Den is Matignon for his help in some tedious calculations. Many thanks also to Staffan Schedin for his contribution in the experimental part of the work.


Eduardo Cardoso

Arnaldo Batista

Rui Rodrigues

Manuel Duarte Ortigueira

Cristina Bárbara

Cristina Martinho

Raul Rato

Manuel Ortigueira

A Contribution for the Automatic Sleep Classification Based on the Itakura-Saito Spectral Distance

Keywords: Sleep Staging, EEG signal processing, AR modelling, Itakura-Saito Distance

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

According to the new criteria based on the Rechtschaffen and Kales(R&K) rules, the human sleep is divided in two phases, the rapid eye movement (REM) and the nonrapid eye movement (NREM).sleep [START_REF] Akaike | Power spectrum estimation through autoregressive model fitting[END_REF]. The NREM is divided into three stages, N1 N2 and N3, being N3 the deepest state of sleep. The various states are distinguished on the basis of pattern changes in the electroencephalogram (EEG), electromyogram (EMG) and electro-oculogram (EOG). The depth of non-REM sleep stages are arranged from the lighter to the deepest, as follows: N1, N2 and N3. The awake state is characterized by alpha frequencies , to score an epoch as awake is just necessary that over 50% of the epoch has alpha rhythm. The REM state has low chin EMG tone, saw-tooth waves, trains of sharply contoured or triangular, often serrated, 2-6 Hz waves maximal in amplitude and symmetrical EOG deflections. The N1 is characterized by mixed frequency activity (4-7 Hz) and low amplitude waves and slow eye movement. Following the N1 is N2 when spindles or K-complexes disassociated with arousals occur. A spindle is a train of distinct waves with frequency 11-16 Hz (most commonly 12-14 Hz) with duration 0.5 seconds, usually maximal in amplitude in central derivations. K-complex is a well-delineated negative sharp wave immediately followed by a positive component standing out from the background EEG, with total duration around 0.5 seconds, usually maximal in amplitude when recorded from frontal derivations. For an arousal to be associated with a K complex, it must begin no more than 1 second after the termination of the K complex. The N3 stage shows waves of frequency 0.5-2Hz and peak-to-peak amplitude higher than 75 µV, when these waves occupy more than 20 to 50% of the 30 seconds epoch, this is scored as N3.

The five stages of sleep, including their repetition, occur cyclically in a normal subject. The first sleep cycle which ends after the conclusion of the first REM stage, usually lasts for about 100 minutes. Each subsequent cycle lasts longer, as its respective closing REM stage extends in duration. So a subject with a regular sleep can complete five cycles in a typical sleep nigh. This sleep cycle is variable, influenced by several agents such as sleep deprivation, frequently changing sleep schedules, stress and environment effects. Another factor that affects the sleep cycle is age: during infancy and childhood the percentage of REM sleep is highest, decreasing with age [START_REF] Ebrahimi | Assessment of Itakura Distance as a valuable feature for computeraided classification of sleep stages[END_REF]. Figure 1 shows the sleep cycles in childhood, early adulthood and old age [START_REF] Estrada | Itakura Distance: A useful similarity measure between EEG and EOG signals in computer-aided classification of sleep stages[END_REF]. Figure 2 shows the relative distribution of the wake, REM and Slow wave sleep (N3) along age [START_REF] Gray | Distortion Measures for Speech Processing[END_REF]. REM period decreases and stabilizes at around 16 years of age, while the awaking period steadily increases with age. REM sleep is associated with the development of the brain in the early life.

The aim of our work is to contribute to the automatic sleep classification using the EEG signal (C3-A2) signal and the Itakura-Saito [START_REF] Iber | The American Academy of Sleep Medicine Manual for the Scoring of Sleep and Associated Events[END_REF] spectral distance estimator. A spectral template for the wake state is firstly obtained. This template is then compared with all the spectral representations of all night epochs. Without any further processing a hypnogram is obtained where N2, N3 and wake states are distinguishable. However N1 and REM stages are often indistinguishable, requiring other detection strategies, namely using the EOG. This scoring falls in the category of a continuous unsupervised sleep staging for an all night session.

-Contribution to Technological Innovation

This paper presents the following main contributions to innovation:

1. An automated process to the calculation of the awake spectrum template based on the spectrogram and the AR modeling. This procedure is adequate for a unsupervised sleep classification scheme. 2. Validation of the ISD as a tool for the Sleep Classification Implementation of a software tool that incorporates the developed algorithms in a user friendly interface (SleepLab).

Materials and Methods

Pre-Processing

We have dealt with signal saturation commonly occurring in continuous EEG recordings by replacing each of these epochs by the previous one, a procedure also used in the visual scoring. Baseline removal was performed through wavelet filtering. Movement artefacts that lead to abnormal signal amplitudes were dealt in the same way as in saturation cases. In all cases the affected epochs were replaced by the previous ones, a procedure also used in the visual scoring.

Awake Template Estimation

An awake template is required for the purpose of the spectral similarity test. Awake epochs were obtained via the power spectrum (Welch method) of each epoch, and then selecting the epochs whose power spectrum peak lie between 8 and 11 Hz, and is above a selected threshold. Due to spindle activity and other EEG components, this power peak must be predominant in order for the epoch classified as awake. An autoregressive (AR) spectral estimation of all the awake classified epochs are calculated, and their mean represents the awake template. Figure 3 represents the power spectrum of the template (solid) and of a sleep epoch (solid) relative to a study patient. The power spectrum of the EEG epochs typically exhibit two main peaks.

Itakura-Saito Distance

This tool [START_REF] Iber | The American Academy of Sleep Medicine Manual for the Scoring of Sleep and Associated Events[END_REF] [8] is used to compare each EEG epoch with the mentioned template, using their spectra, since the Euclidian distance is inadequate for this case. In this study we have used AR modelling for obtaining the spectra [START_REF] Rechtschaffen | Sleep and Dreaming[END_REF] [START_REF] Soong | A Frequency-Weighted Itakura Spectral Distortion Measure and It's Application to Speech Recognition in Noise[END_REF]. The Itakura-Saito spectral distance (ISD) d IS is defined by the general expression:

2 2 2 2 ( ) ( ) log 1 2 ( ) ( ) IS X S d d S X π π ω ω ω π ω ω -       = +   -         ∫ (1)
where is the spectrum of the epoch under study and is the template spectrum. This distance is sensible to the relative spectra amplitudes which can be a benefit since some sleep staging can be performed based o the EEG amplitude, namely the higher amplitudes normally present in the N3. However this becomes drawback if the EEG also has amplitude fluctuations due, for instance, to electrode displacement. To overcome this, we have normalized all the spectra to have unity gain, at zero frequency, being the price paid for this to reduce sensibility to the signal amplitude and therefore decreasing N3 detection rate. After some experimentation we have found that this option produced better results in our cases. Itakura-Saito distance has been used in the context of speech recognition and measures the dissimilarity of two spectra and , having the functions peaks a predominant weight in the d IS value, which makes sense in our case where peaks in certain frequency bands define the respective sleep stage.

Results

Figure 4 shows a patient hypnogram expert scored (blue) and the Itakura-Saito Spectral (ISD) distances obtained as described above (green). The ISD has not been processed, besides a 4 coefficient moving average filter. We confirm that the awake states follow the hypnogram therefore validating our template.

Stages N2 and N3 are also detected. REM is also detected although with lesser accuracy. This patient, as most of the others, exhibit only a few N1 epochs being therefore not possible to fully assess the ISD to detect this state. It is also patent that the ISD plot allows for the evaluation of the way the sleep state continuously progress from stage to stage, unlike the visual scoring that produces discrete stage transitions. Figure 5 shows the mean and standard deviation of the ISD values for each visual scored stage, for the same patient. As expected ISD is lower in the awake state and higher in the N3 state. N1 and REM stages are similar spectral wise, which can be seen in this figure. A software tool named SleepLab was developed in Matlab [11] to streamline the data pre-processing, template estimation and Itakura-Saito distance. The user interface is shown in Figure 6. Other versions of the Itakura distance are also obtained for comparison purposes. This platform is opened to future developments as ICA (Independent Component Analysis) to remove the EOG artifacts, and sleep spindle detection. This is an open platform that will include newer developments. Fig. 6. A screen-shot of the user interface. The Save button exports all the results and models for comparison and statistical studies.

Conclusions and Future Work

Continuous EEG signal analysis for automatic sleep classification is a difficult task. This kind of signal is typically affected by the usual artifacts related to the patient unpredictable movements and possible electrode displacements during the sleep period. We have introduced an automatic procedure to obtain the awake template using the spectrogram of all the night´s epochs, which is computationally very efficient and essential for the ISD calculation. We have used the ISD formulation (distortion measure) which we have found to produce better results than the gainoptimized Itakura-Saito distortion measure generally used [START_REF] Gray | Distortion Measures for Speech Processing[END_REF]. It is still unclear which of the distortion measures better suits the sleep scoring via EEG signal analysis, and we intend to continue the research in order to differentiate and assess the different spectral distances measures for this task. In this work the ISD produced automatic hypnograms that could detect the awake, N2 and N3 states. The assessment of the N1 is compromised by the fact that this state is relatively rare among the tested patients. The REM is identical to the N1, spectral wise, and therefore requires other features to be automatically detected, namely the EOG.

We also intend to submit the ISD obtained results to a Neuro-Fuzzy classifier, as this seems to be an obvious step to follow. Furthermore, we plan to use the EOG as a classifier of the REM stage and apply further artifact removing techniques such as ICA. EMG chin activity is planned to be introduced as an important indicator of sleep deepness.

  Sleep EEG data from Hospital Pulido Valente was used in this work. EEG signals were recorded at a rate of 100Hz. Visual sleep scoring of sleep stages were scored by Hospital experts, according to Rechtschaffen & Kales based on 30-second epochs. These hypnograms stand as our reference.

Fig. 1 .

 1 Fig. 1. Reference hypnograms for childhood, early adulthood and old age [3].

Fig. 2 .

 2 Fig. 2. Wake, REM and deep sleep time along age. [4].

Fig. 3 .

 3 Fig. 3. Spectrum of the obtained EEG awake template (solid) and a sleep epoch for the same patient (dash). See text. Its relative location and value varies according to the sleep/awake state. The order of the AR model has been estimated based on three factors: the Minimum Description

Fig. 4 .Fig. 5 .

 45 Fig. 4. Visually scored hypnogram (blue) and ISD automatic scoring (green). See text.