We propose a new method for estimating the change-points of heart rate in the
orthosympathetic and parasympathetic bands, based on the wavelet transform in
the complex domain and the study of the change-points in the moments of the
modulus of these wavelet transforms. We observe change-points in the
distribution for both bands.Comment: proceeding of the workshop 'Fouille de donn\'ees temporelles et
  analyse de flux de donn\'ees' EGC'2009, january 27, Strasbourg, Franc

Bertrand, Pierre R.

Boudet, Gil

Chamoux, Alain

Teyssière, Gilles

English

arXiv

proceeding of the workshop 'Fouille de données temporelles et analyse de flux de données' EGC'2009, january 27, Strasbourg, FranceWe propose a new method for estimating the change-points of heart rate in the orthosympathetic and parasympathetic bands, based on the wavelet transform in the complex domain and the study of the change-points in the moments of the modulus of these wavelet transforms. We observe change-points in the distribution for both bands

Bertrand, Pierre, Raphael

HAL Clermont Université

Detection of Change–Points in the Spectral Density.With Applications to ECG DataPierre, Raphael Bertrand, Gilles Teyssie`re, Gil Boudet, Alain ChamouxTo cite this version:Pierre, Raphael Bertrand, Gilles Teyssie`re, Gil Boudet, Alain Chamoux. Detection of Change–Points in the Spectral Density. With Applications to ECG Data. proceeding of the workshop’Fouille de donne´es temporelles et analyse de flux de donne´es’ EGC’200.. 2008. <hal-00354416>HAL Id: hal-00354416https://hal.archives-ouvertes.fr/hal-00354416Submitted on 19 Jan 2009HAL is a multi-disciplinary open accessarchive for the deposit and dissemination of sci-entific research documents, whether they are pub-lished or not. The documents may come fromteaching and research institutions in France orabroad, or from public or private research centers.L’archive ouverte pluridisciplinaire HAL, estdestine´e au de´poˆt et a` la diffusion de documentsscientifiques de niveau recherche, publie´s ou non,e´manant des e´tablissements d’enseignement et derecherche franc¸ais ou e´trangers, des laboratoirespublics ou prive´s.Detection of Change–Points in the Spectral Density.With Applications to ECG DataPierre R. Bertrand, Gilles Teyssière, Gil Boudet, Alain ChamouxINRIA Saclay, APIS TeamUMR CNRS 6620, Clermont-Ferrand University, France.E-mail: Pierre.Bertrand@inria.fr,,CREATES, Aarhus University, DenmarkE-mail: stats@gillesteyssiere.net,Institut de Médecine du travail, UFR Médecine, Clermont-Ferrand UniversityGil.Boudet@wanadoo.fr,Institut de Médecine du travail, UFR Médecine, Clermont-Ferrand UniversityCHU Clermont-FerrandE-mail: alain.chamoux@u-clermont1.frRésumé. Nous proposons une nouvelle méthode pour estimer les ruptures du rythme cardiaque dans lesbandes parasympathiques et orthosympathiques basée sur la transformée en ondelettes dans le domainecomplexe et l’étude des ruptures des moments des modules de ces transformées en ondelettes. Nous obser-vons des ruptures dans la distribution pour les deux bandes de fréquence.1 IntroductionECG signal analysis has a long story after the implementation of the ambulatory monitoring by Holter at the beginningof the fifties. Recent measurement methods, due to the size reduction of the measurement devices, see Chamoux (1984),allow us to record ECG data for healthy people over a long period of time : long distance (marathon) runners, individualsdaily (24 hours) records, etc. We then obtain large datasets that allow us to characterize the variations of the heartbeat ratein the two components of the nervous autonomous system : the parasympathetic and the orthosympathetic ones.5h50 9h43’ 13h38’ 16h32’ 20h07’ 0h14’ 4h52’050010001500FIG. 1 – RR interval for a healthy subject during a period of 24 hoursDetection of Change–Points in the Spectral DensityThe data studied in this paper have been recorded the 29th of January 2008, at the Clermont-Ferrand UniversityHospital (CHU) by G. Boudet. Unlike Boudet et al. (2004) and Cottin et al (2006), these data are not recorded in theframework of a laboratory experiment, but during the real life of the individual under investigation.We have daily datasets, over 100,000 observations, and in the near future we will have data over several days, i.e., overthat sample size. We are measuring the interval between two RR peaks, i.e., if we denote by (ti)i=1,...,N , the sequence ofpeak times measured with a precision of 1.E-03 second, we consider the seriesX(ti) =(ti− ti−1)measured in seconds.Several indicators are related to heartbeat data. The most popular is the instantaneous average frequency, i.e., X(ti)−1,which is displayed by runners watches. This quantity is informative on daily activity of observed individual : sleeping andwaking up times, physical activity, e.g., sport, manual work, but does not summarize all relevant information. Note thatheartbeat data display large variations, clustering, etc, as only individuals with serious disease display a regular heart rate.Cardiologists are interested in the study of this signal in two frequency bands : the orthosympathetic and parasympa-thetic bands, i.e., the frequency bands (0.04Hz, 0.15Hz) and (0.15Hz, 0.5Hz) respectively. The definition of thesebands is the outcome of research works, see e.g., Task force of the European Society of Cardiology and the North Ameri-can Society of Pacing and Electrophysiology (1996), and is based on the fact that the energy contained inside these bandswould be a relevant indicator on the level of the stress of an individual.Indeed, for the heart rate, the parasympathetic system is often compared to the brake while the orthosympatheticsystem would be a nice accelerator ; see e.g. Goldberger (2001). At rest there is a permanent braking effect on the heartrate. Any solicitation of the cardiovascular system, any activity initially produces a reduction of parasympathetic brakefollowed by a gradual involvement of the sympathetic system. These mechanisms are very interesting to watch in manydiseases including heart failure, but also rhythm disorders that may fall under one or other of these two effects, monitoringthe therapeutic effect of several Medicines including some psychotropic. In the field of physiology such data are crucialfor measuring the level of vigilance and particularly the risk of falling asleep while driving a vehicle, the level of stressinduced by physical activity or level of perceived stress, which can be considered as a criterion of overtraining in sport.Fractals models have been used in cardiology after the works by Ivanov et al. (1999), who applied the multifractalspectrum analysis advocated by Frisch (1996), for modeling RR series and classifying individuals according to this mul-tifractal spectrum, as this spectrum discriminates between individuals who experienced hearth trouble, and those who didnot. However, this tool has some shortcomings as it requires huge samples common in turbulence analysis, and is theninappropriate for studying phenomena occurring at a resolution lower than the daily time interval, such as the variationsof the parasympathetic and orthosympathetic systems inside the day.Wavelets based methods have been used in biostatistics by Diab et al. (2007) for uterine EMG signal analysis. Ho-wever, they consider that the process studied is homogeneous, and used these methods as a classification tool. Anothersignificant difference is the fact that they use discrete wavelet decomposition, i.e., a frequency decomposition on a dyadicwavelet basis, the choice for the frequency bands is made without reference to a biological phenomenon. In our case, thechoice for the frequency band is justified by biological considerations, and we fit the wavelets inside these bands.This is why the continuous analysis of both systems and their quantification is a particularly promising research area.So, the example in Figures 4, 5, 6 shows at the observation 28,220 i.e., 13h40’50”, a simultaneous variation of bothsystems.2 Mathematical formulationFor biological reasons, the heart rate is within the interval [20, 250] bpm (beats per minute), i.e., X(t) belongs to theRR interval 60/250s. < X(t) < 60/20s.. This leads to modeling by stationary or locally stationary processes. We wishto decompose this signal in a sequence of homogeneous intervals having the same mean and/or variance. We make theassumption that the series is Gaussian, and then find the optimal segmentation of the process using the same approach asin Lavielle and Teyssière (2006) or Lavielle and Moulines (2000).Bertrand, Teyssière, Boudet, ChamouxWe assume that the signal is the sum of a piecewise constant function and a Gaussian process, centered and locallystationary. We then have the following representation :X(t) = µ(t) +∫IReitξf1/2(t, ξ) dW (ξ), for all t ∈ R, (1)where– ξ 7→ f(t, ξ) is an even and positive function, called spectral density piecewise constant, i.e., there exists a partitionτ1, . . . , τK such that f(t, ξ) = fk(ξ) for t ∈ [τi, τi+1[– the function t 7→ µ(t) is also piecewise constant for another partition τ˜1, . . . , τ˜L with µ(t) = µℓ if t ∈ [τ˜ℓ, τ˜ℓ+1[The wavelet coefficient associated with ψ isWψ(b) =∫IRψ(t− b)X(t) dt unit in second2. (2)In Bardet and Bertrand (2007) or Bardet et al. (2008), one can find a theoretical study of the wavelet coefficient forstationary (or with stationary increment) centered Gaussian processes, i.e., for X given by (1) with µ(t) = 0. Thisapproach can be generalized to locally stationary Gaussian processes, this will be detailed in a forthcoming paper.According to recommendations of Task force of the European Soc. Cardiology and the North American Soc. of Pacingand Electrophysiology (1996), we use the following notations :– [ω1, ω2] = (0.04Hz, 0.15Hz) denotes the orthosympathetic frequency band,– [ω2, ω3] = (0.15Hz, 0.5Hz) denotes the parasympathetic bandThe energy associated with each of these frequency bands and localized around the time b, is measured by the modulus ofthe complex wavelet coefficients |Wi(b)|2 for i = 1, 2, withb 7−→Wi(b) =∫IRψi(t− b)X(t) dt, i = 1, 2.In this work, these wavelets coefficients are computed at each second, i.e., the difference between two consecutive valuesfor b is equal to 1 second.How to choose the wavelets ψ1 and ψ2 ?In the idealistic case, one would use two filters ψ1 and ψ2 with compact support, the Fourier transforms of which havesupport inside the orthosympathetic and parasympathetic bands.Unfortunately, it does not exist a non null function ψ with compact time domain support and compact frequencysupport, see for instance Mallat (1998) Th 2.6 p.34. Therefore, the best we can do is to choose between a filter with acompact frequency support and a filter with a compact time domain support. The first choice is well suited for stationarymodels, see Bardet and Bertrand (2007). But, in this work, we are interested by locally stationary models, thus ourspecifications are a filter with compact time domain support as a Daubechies wavelet. The price to pay for the compactnessof the time domain support is the loss in the compactness of the frequency support. However, the frequency support is"almost compact" in the following sense :Definition 2.1 (ρ pseudo support) Let 0 < ρ < 1, be a map g ∈ L2(IR) that admits the compact interval I as a ρ pseudosupport if∫I|g(t)|2 dt∫IR|g(t)|2 dt= ρ.Fourier transform of Daubechies wavelet have a reasonably small ρ pseudo support with a ratio ρ close to 1. Moreover, thelarger the number of the Daubechies wavelet is, the closer to 1 the ratio ρ is ; see the example below with the Daubechieswavelet D6.By scaling and modulation, one can adjust the pseudo support inside a specified frequency band as stated by thefollowing propositionDetection of Change–Points in the Spectral DensityProposition 2.1 Let ψ be a filter with compact support [L1, L2] and a frequency ρ pseudo support [Λ1,Λ2], for anyfrequencies band [ω1, ω2], the map ψ1(t) = µ× eiηt ψ(λt) withµ > 0, λ =ω2 − ω1Λ2 − Λ1and η =ω1 + ω22− (ω2 − ω1)Λ2 + Λ1Λ2 − Λ1has a ρ pseudo support [ω1, ω2] and a time domain support[Λ2 − Λ1ω2 − ω1L1,Λ2 − Λ1ω2 − ω1L2].Proof. Since ψ̂1(ξ) = µ × ψ̂(ξ − ηλ), one can deduce ρ pseudo suppψ1 = η + λ × ρ pseudo suppψ and then theproposition.Different choices for the wavelets ψ1 and ψ2From Proposition 2.1, one can deduce the different possible choices of the filters ψ1 and ψ2– Daubechies wavelet D6 : In this case, we have Λ1 = 0.08, Λ2 = 1.75, ρ = 0.9999 and one can setψ1(t) = µ× eiη1tD6(λ1t) and ψ2(t) = µ× eiη2tD6(λ2t)with η1 = −0.0255, λ1 = 0.0659, η2 = −0.0585, and λ2 = 0.2096.The length of the time support are |Suppψ1| and |Suppψ2|. One can see on Fig. 2 (a) that the Fourier transformsψ̂1(x) and ψ̂2(x) have almost disjoint supports.λ1 =ω2 − ω1Λ2 − Λ1and η1 =ω1 + ω22− (ω2 − ω1)Λ2 + Λ1Λ2 − Λ1λ2 =ω3 − ω2Λ2 − Λ1and η2 =ω2 + ω32− (ω3 − ω2)Λ2 + Λ1Λ2 − Λ1|Suppψ1| =Λ2 − Λ1ω2 − ω1× |SuppD6| and |Suppψ2| =Λ2 − Λ1ω3 − ω2× |SuppD6|– Gabor wavelet : For computational reasons, we will use the Gabor wavelet defined asψ(t) = eiηtg(t), g(t) =1(σ2pi)1/4e−t22σ2 (3)see, e.g., Mallat (1998). This wavelet has the same time and frequency ρ pseudo support [−L,L] = [−3.5, 3.5] withρ = 0.9995. In the spectral domain, we haveψ̂(t) = gˆ(ξ − η), gˆ(ξ) = (4piσ2)1/4e−σ2ξ22 (4)We fit the Gabor wavelet inside the orthosympathetic or the parasympathetic frequency bands by using Prop. 2.1 ordirect calculations. One obtains the two Gabor wavelets defined by (3) with the following choice for the parameters :η1 =ω1 + ω22and σ1 =2Lω2 − ω1η2 =ω2 + ω32and σ2 =2Lω3 − ω2moreover |ρ pseudoSuppψ1| =4L2ω2 − ω1and |ρ pseudoSuppψ2| =4L2ω3 − ω2with ρ = 0.9995One can see on Fig. 2 that the Fourier transforms ψ̂1(x) and ψ̂2(x) still have almost disjoint supports.Using the Gabor wavelet is more efficient in terms of computing time, as it is at least 8 times faster. Figure 3 belowdisplays the Gabor wavelets coefficients in the orthosympathetic and parasympathetic bands respectively for the sampleplotted in Figure 1.Bertrand, Teyssière, Boudet, Chamoux−0,2 −0,1 −0 0,1 0,2 0,3 0,4 0,5 0,6 0,70,04 0.15024681012  ψ̂1ψ̂2FIG. 2 – The Fourier transforms ψ̂1(x) (left) and ψ̂2(x) (right)0246810x 104−1000−50005001000−1000−50005001000timeW1 real partW1 complex part0246810x 104−5000500−400−2000200400600timeW2 real partW2 complex partFIG. 3 – The wavelet coefficients in the orthosympathetic band (left) and in the parasympathetic band (right) of the samehealthy subject during a period of 24 hours3 Segmentation analysisWe assume that the process {Xt} is abruptly changing and is characterized by a parameter θ ∈ Θ that remains constantbetween two changes. We use this assumption to define our contrast function J(τ ,X). Let K be some integer and letτ = {τ1, τ2, . . . , τK−1} be an ordered sequence of integers satisfying 0 < τ1 < τ2 < . . . < τK−1 < n.For the detection of changes in the mean and variance of a sequence of random variables, i.e., a change in distribution,Detection of Change–Points in the Spectral Densitywe use the following contrast function, based on a Gaussian log–likelihood function :Jn(τ ,X) =1nK∑k=1‖Xτk − X¯τk‖2σˆ2k+ nk log(σˆ2k). (5)where nk = τk − τk−1 is the length of segment k, σˆ2k is the empirical variance computed on that segment k, σˆ2k =n−1k∑τki=τk−1+1(Xi − X¯)2, and X¯ is the empirical mean of X1, . . . , Xn.In fact, we minimize a penalized version of this contrast function, the penalty parameter being selected in an adaptiveway so that the obtained segmentation does not depend too much on the penalty parameter ; see also Birgé and Massart(2007) for another choice for the penalty parameter. Although this method has been devised for Gaussian processes, itempirically provides relevant results for non-Gaussian processes, e.g., financial time series, see Lavielle and Teyssière(2006) for further details.3.1 Orthosympathetic bandFor that frequency band, we obtain the following segmentation :τ = {28220, 33366, 71048},i.e., 13h40’50”, 15h06’36”, 1h34’38”.5h50’30" 13h40’50" 15h06’36" 1h34’38"050100150200250300350400450500FIG. 4 – Segmentation in the mean and variance (change in the distribution) of the modulus of the wavelet coefficients inthe orthosympathetic band3.2 Parasympathetic bandFor that frequency band, we obtain the following segmentation :τ = {11620, 21912, 28054, 31540, 36022, 40172, 52622, 70550},i.e., 9h04’10”, 11h45’42”, 13h38’04”, 14h36’10”, 15h50’22”, 17h00’02”, 20h27’32”, 1h26’40”.4 ConclusionThe example illustrated by Figures 4, 5 and 6 shows at around t = 28, 220 a simultaneous variation of both systems.But, one can also observe change-points in the orthosympathetic and parasympathetic bands occurring at different times.Bertrand, Teyssière, Boudet, Chamoux5h50’30" 9h04’10" 11h45’42" 13h38’ 14h36’ 15h50’ 17h00’02" 20h27’32" 1h26’40"050100150200250300350400FIG. 5 – Segmentation in the mean and variance (change in the distribution) of the modulus of the wavelet coefficients inthe parasympathetic bandIn the future, we will study the existence of a possible causality or sequentiality between these change–points in differentbands.5h50’30" 9h04’10" 11h45’42" 13h38’ 14h36’ 15h50’ 17h00’02" 20h27’32" 1h26’40"0100200300400500600700800900  Orthosympathetic bandParasympathetic bandHearthbeat processFIG. 6 – Mean of the modulus if the wavelets coefficients for the orthosympathetic and parasympathetic bands, and of theRR process X(t)Références[1] Bardet, J.M. and Bertrand, P.R. (2007). Identification of the multiscale fractional Brownian motion with biomecha-nical applications. J. Time Ser. Anal. 28, 1–52.[2] Bardet, J.M., Bertrand, P.R. and Billat, V. (2008). Estimation non-paramétrique de la densité spectrale d’un processusgaussien observé à des instants aléatoires. Pub. ISUP 52, 123–138.[3] Birgé, L. and Massart, P. (2007). Minimal penalties for Gaussian model selection. Probab. Theory Related Fields138, 33–73.Detection of Change–Points in the Spectral Density[4] Boudet, G., Albuisson, E., Bedu, M. and Chamoux, A. (2004). Heart rate running speed relationships during exhaus-tive bouts in the laboratory. Can. J. Appl. Physiol. 29(6), 731–742.[5] Chamoux, A. (1984). Le système Holter en pratique. Médecine du Sport 58, 43-273, 54–284.[6] Cottin, F., Leprêtre, PM., Lopes, P., Papelier, Y., Médigue, C. and Billat V. (2006). Assessment of ventilatory thre-sholds from heart rate variability in well-trained subjects during cycling, Int. J. Sports Med. 27, 959–967.[7] Diab, M.O., Marque, C. and Khalil, M.A. (2007). Classification for uterine EMG Signals : Comparison between ARmodel and statistical classification method. Int. J. Computational Cognition 5, 8–14,[8] Frisch, U. (1996). Turbulence. Cambridge University Press.[9] Goldberger, A.L. (2001). Heartbeats, hormones and health : is variability the spice of life ? Am. J. Crit. Care Med.163, 1289–1290.[10] Ivanov, P.C., Amaral, L.A.N., Goldberger, A.L., Havlin, S., Rosenblum, M., Struzik, Z, and Stanley, H.E. (1999).Multifractality in human heartbeat dynamics. Nature 399, 461–465.[11] Lavielle, M. and Moulines, E. (2000). Least-squares estimation of an unknown number of shifts in a time series. J.of Time Series Anal. 21, 33–59.[12] Lavielle, M. and Teyssière, G. (2006). Detection of multiple change-points in multivariate time series. LithuanianMath. J. 46, 287–306. [Détection de ruptures multiples dans des séries temporelles multivariées.] (Translation inFrench) Liet. Mat. Rink. 46, 351–376.[13] Mallat, S. (1998). A wavelet tour of signal processing. Academic Press.[14] Task force of the European Soc. Cardiology and the North American Society of Pacing and Electrophysiology(1996). Heart rate variability. Standards of measurement, physiological interpretation, and clinical use. Circulation93, 1043–1065.SummaryWe propose a new method for estimating the change-points of heart rate in the orthosympathetic and parasympatheticbands, based on the wavelet transform in the complex domain and the study of the change-points in the moments of themodulus of these wavelet transforms. We observe change-points in the distribution for both bands.

Detection of Change--Points in the Spectral Density. With Applications to ECG Data

HAL Descartes

Hal-Diderot

INRIA a CCSD electronic archive server

HAL Id: hal-00354416https://hal.archives-ouvertes.fr/hal-00354416Preprint submitted on 19 Jan 2009HAL is a multi-disciplinary open accessarchive for the deposit and dissemination of sci-entific research documents, whether they are pub-lished or not. The documents may come fromteaching and research institutions in France orabroad, or from public or private research centers.L’archive ouverte pluridisciplinaire HAL, estdestinée au dépôt et à la diffusion de documentsscientifiques de niveau recherche, publiés ou non,émanant des établissements d’enseignement et derecherche français ou étrangers, des laboratoirespublics ou privés.Detection of Change–Points in the Spectral Density.With Applications to ECG DataPierre, Raphael Bertrand, Gilles Teyssière, Gil Boudet, Alain ChamouxTo cite this version:Pierre, Raphael Bertrand, Gilles Teyssière, Gil Boudet, Alain Chamoux. Detection of Change–Pointsin the Spectral Density. With Applications to ECG Data. 2008. ￿hal-00354416￿Detection of Change–Points in the Spectral Density.With Applications to ECG DataPierre R. Bertrand, Gilles Teyssière, Gil Boudet, Alain ChamouxINRIA Saclay, APIS TeamUMR CNRS 6620, Clermont-Ferrand University, France.E-mail: Pierre.Bertrand@inria.fr,,CREATES, Aarhus University, DenmarkE-mail: stats@gillesteyssiere.net,Institut de Médecine du travail, UFR Médecine, Clermont-Ferrand UniversityGil.Boudet@wanadoo.fr,Institut de Médecine du travail, UFR Médecine, Clermont-Ferrand UniversityCHU Clermont-FerrandE-mail: alain.chamoux@u-clermont1.frRésumé. Nous proposons une nouvelle méthode pour estimer les ruptures du rythme cardiaque dans lesbandes parasympathiques et orthosympathiques basée sur latransformée en ondelettes dans le domainecomplexe et l’étude des ruptures des moments des modules de ces transformées en ondelettes. Nous obser-vons des ruptures dans la distribution pour les deux bandes de fréquence.1 IntroductionECG signal analysis has a long story after the implementation of the ambulatory monitoring by Holter at the beginningof the fifties. Recent measurement methods, due to the size reduction of the measurement devices, see Chamoux (1984),allow us to record ECG data for healthy people over a long period of time : long distance (marathon) runners, individualsdaily (24 hours) records, etc. We then obtain large datasetstha allow us to characterize the variations of the heartbeat ratein the two components of the nervous autonomous system : the parasympathetic and the orthosympathetic ones.5h50 9h43’ 13h38’ 16h32’ 20h07’ 0h14’ 4h52’050010001500FIG. 1 – RR interval for a healthy subject during a period of 24 hoursDetection of Change–Points in the Spectral DensityThe data studied in this paper have been recorded the 29th of January 2008, at the Clermont-Ferrand UniversityHospital (CHU) by G. Boudet. Unlike Boudett al. (2004) and Cottinet al (2006), these data are not recorded in theframework of a laboratory experiment, but during the real life of the individual under investigation.We have daily datasets, over 100,000 observations, and in the ear future we will have data over several days, i.e., overthat sample size. We are measuring the interval between two RR peaks, i.e., if we denote by(ti)i=1,...,N , the sequence ofpeak times measured with a precision of 1.E-03 second, we consider the seriesX(ti) =(ti− ti−1)measured in seconds.Several indicators are related to heartbeat data. The most popular is the instantaneous average frequency, i.e.,X(ti)−1,which is displayed by runners watches. This quantity is informative on daily activity of observed individual : sleepingandwaking up times, physical activity, e.g., sport, manual work, but does not summarize all relevant information. Note thatheartbeat data display large variations, clustering, etc,as only individuals with serious disease display a regular he t rate.Cardiologists are interested in the study of this signal in two frequency bands : the orthosympathetic and parasympa-thetic bands, i.e., the frequency bands(0.04Hz, 0.15Hz) and(0.15Hz, 0.5Hz) respectively. The definition of thesebands is the outcome of research works, see e.g., Task force of th European Society of Cardiology and the North Ameri-can Society of Pacing and Electrophysiology (1996), and is based on the fact that the energy contained inside these bandswould be a relevant indicator on the level of the stress of an indiv dual.Indeed, for the heart rate, the parasympathetic system is often compared to the brake while the orthosympatheticsystem would be a nice accelerator ; see e.g. Goldberger (2001). At rest there is a permanent braking effect on the heartrate. Any solicitation of the cardiovascular system, any activity initially produces a reduction of parasympathetic brakefollowed by a gradual involvement of the sympathetic system. These mechanisms are very interesting to watch in manydiseases including heart failure, but also rhythm disorders that may fall under one or other of these two effects, monitori gthe therapeutic effect of several Medicines including somepsychotropic. In the field of physiology such data are crucialfor measuring the level of vigilance and particularly the risk of falling asleep while driving a vehicle, the level of stressinduced by physical activity or level of perceived stress, which can be considered as a criterion of overtraining in sport.Fractals models have been used in cardiology after the worksby Ivanovet al. (1999), who applied the multifractalspectrum analysis advocated by Frisch (1996), for modelingRR series and classifying individuals according to this mul-tifractal spectrum, as this spectrum discriminates between individuals who experienced hearth trouble, and those whodidnot. However, this tool has some shortcomings as it requireshuge samples common in turbulence analysis, and is theninappropriate for studying phenomena occurring at a resolution lower than the daily time interval, such as the variationsof the parasympathetic and orthosympathetic systems inside the day.Wavelets based methods have been used in biostatistics by Dia et al. (2007) for uterine EMG signal analysis. Ho-wever, they consider that the process studied is homogeneous, and used these methods as a classification tool. Anothersignificant difference is the fact that they use discrete wavelet decomposition, i.e., a frequency decomposition on a dyadicwavelet basis, the choice for the frequency bands is made without reference to a biological phenomenon. In our case, thechoice for the frequency band is justified by biological considerations, and we fit the wavelets inside these bands.This is why the continuous analysis of both systems and theirquantification is a particularly promising research area.So, the example in Figures 4, 5, 6 shows at the observation 28,220 i.e., 13h40’50”, a simultaneous variation of bothsystems.2 Mathematical formulationFor biological reasons, the heart rate is within the interval [20, 250] bpm (beats per minute), i.e.,X(t) belongs to theRR interval60/250s. < X(t) < 60/20s.. This leads to modeling by stationary or locally stationaryp ocesses. We wishto decompose this signal in a sequence of homogeneous intervals having the same mean and/or variance. We make theassumption that the series is Gaussian, and then find the optimal segmentation of the process using the same approach asin Lavielle and Teyssière (2006) or Lavielle and Moulines (2000).Bertrand, Teyssière, Boudet, ChamouxWe assume that the signal is the sum of a piecewise constant function and a Gaussian process, centered and locallystationary. We then have the following representation :X(t) = µ(t) +∫IReitξf1/2(t, ξ) dW (ξ), for all t ∈ R, (1)where– ξ 7→ f(t, ξ) is an even and positive function, called spectral density piecewise constant, i.e., there exists a partitionτ1, . . . , τK such thatf(t, ξ) = fk(ξ) for t ∈ [τi, τi+1[– the functiont 7→ µ(t) is also piecewise constant for another partitionτ̃1, . . . , τ̃L with µ(t) = µℓ if t ∈ [τ̃ℓ, τ̃ℓ+1[The wavelet coefficient associated withψ isWψ(b) =∫IRψ(t− b)X(t) dt unit in second2. (2)In Bardet and Bertrand (2007) or Bardetet al. (2008), one can find a theoretical study of the wavelet coeffici nt forstationary (or with stationary increment) centered Gaussian processes, i.e., forX given by (1) withµ(t) = 0. Thisapproach can be generalized to locally stationary Gaussianprocesses, this will be detailed in a forthcoming paper.According to recommendations of Task force of the European Soc. Cardiology and the North American Soc. of Pacingand Electrophysiology (1996), we use the following notations :– [ω1, ω2] = (0.04Hz, 0.15Hz) denotes the orthosympathetic frequency band,– [ω2, ω3] = (0.15Hz, 0.5Hz) denotes the parasympathetic bandThe energy associated with each of these frequency bands andlocalized around the timeb, is measured by the modulus ofthe complex wavelet coefficients|Wi(b)|2 for i = 1, 2, withb 7−→Wi(b) =∫IRψi(t− b)X(t) dt, i = 1, 2.In this work, these wavelets coefficients are computed at each se ond, i.e., the difference between two consecutive valuesfor b is equal to 1 second.How to choose the waveletsψ1 andψ2 ?In the idealistic case, one would use two filtersψ1 andψ2 with compact support, the Fourier transforms of which havesupport inside the orthosympathetic and parasympathetic bands.Unfortunately, it does not exist a non null functionψ with compact time domain support and compact frequencysupport, see for instance Mallat (1998) Th 2.6 p.34. Therefore, the best we can do is to choose between a filter with acompact frequency support and a filter with a compact time domain support. The first choice is well suited for stationarymodels, see Bardet and Bertrand (2007). But, in this work, weare interested by locally stationary models, thus ourspecifications are a filter with compact time domain support as a Daubechies wavelet. The price to pay for the compactnessof the time domain support is the loss in the compactness of the requency support. However, the frequency support is"almost compact"in the following sense :Definition 2.1 (ρ pseudo support) Let0 < ρ < 1, be a mapg ∈ L2(IR) that admits the compact intervalI as aρ pseudosupport if∫I|g(t)|2 dt∫IR|g(t)|2 dt= ρ.Fourier transform of Daubechies wavelet have a reasonably small ρ pseudo support with a ratioρ close to1. Moreover, thelarger the number of the Daubechies wavelet is, the closer to1 the ratioρ is ; see the example below with the Daubechieswavelet D6.By scaling and modulation, one can adjust the pseudo supportinside a specified frequency band as stated by thefollowing propositionDetection of Change–Points in the Spectral DensityProposition 2.1 Let ψ be a filter with compact support[L1, L2] and a frequencyρ pseudo support[Λ1,Λ2], for anyfrequencies band[ω1, ω2], the mapψ1(t) = µ× eiηt ψ(λt) withµ > 0, λ =ω2 − ω1Λ2 − Λ1and η =ω1 + ω22− (ω2 − ω1)Λ2 + Λ1Λ2 − Λ1has aρ pseudo support[ω1, ω2] and a time domain support[Λ2 − Λ1ω2 − ω1L1,Λ2 − Λ1ω2 − ω1L2].Proof. Sinceψ̂1(ξ) = µ × ψ̂(ξ − ηλ), one can deduceρ pseudo suppψ1 = η + λ × ρ pseudo suppψ and then theproposition.Different choices for the waveletsψ1 andψ2From Proposition 2.1, one can deduce the different possiblechoices of the filtersψ1 andψ2– Daubechies wavelet D6 :In this case, we haveΛ1 = 0.08, Λ2 = 1.75, ρ = 0.9999 and one can setψ1(t) = µ× eiη1tD6(λ1t) and ψ2(t) = µ× eiη2tD6(λ2t)with η1 = −0.0255, λ1 = 0.0659, η2 = −0.0585, and λ2 = 0.2096.The length of the time support are|Suppψ1| and|Suppψ2|. One can see on Fig. 2 (a) that the Fourier transformsψ̂1(x) andψ̂2(x) have almost disjoint supports.λ1 =ω2 − ω1Λ2 − Λ1and η1 =ω1 + ω22− (ω2 − ω1)Λ2 + Λ1Λ2 − Λ1λ2 =ω3 − ω2Λ2 − Λ1and η2 =ω2 + ω32− (ω3 − ω2)Λ2 + Λ1Λ2 − Λ1|Suppψ1| =Λ2 − Λ1ω2 − ω1× |SuppD6| and |Suppψ2| =Λ2 − Λ1ω3 − ω2× |SuppD6|– Gabor wavelet : For computational reasons, we will use the Gabor wavelet defined asψ(t) = eiηtg(t), g(t) =1(σ2π)1/4e−t22σ2 (3)see, e.g., Mallat (1998). This wavelet has the same time and frequencyρ pseudo support[−L,L] = [−3.5, 3.5] withρ = 0.9995. In the spectral domain, we haveψ̂(t) = ĝ(ξ − η), ĝ(ξ) = (4πσ2)1/4e−σ2ξ22 (4)We fit the Gabor wavelet inside the orthosympathetic or the parasympathetic frequency bands by using Prop. 2.1 ordirect calculations. One obtains the two Gabor wavelets defined by (3) with the following choice for the parameters :η1 =ω1 + ω22and σ1 =2Lω2 − ω1η2 =ω2 + ω32and σ2 =2Lω3 − ω2moreover |ρ pseudoSuppψ1| =4L2ω2 − ω1and |ρ pseudoSuppψ2| =4L2ω3 − ω2with ρ = 0.9995One can see on Fig. 2 that the Fourier transformsψ̂1(x) andψ̂2(x) still have almost disjoint supports.Using the Gabor wavelet is more efficient in terms of computing me, as it is at least 8 times faster. Figure 3 belowdisplays the Gabor wavelets coefficients in the orthosympathetic and parasympathetic bands respectively for the sampleplotted in Figure 1.Bertrand, Teyssière, Boudet, Chamoux−0,2 −0,1 −0 0,1 0,2 0,3 0,4 0,5 0,6 0,70,04 0.15024681012  ψ̂1ψ̂2FIG. 2 – The Fourier transformŝψ1(x) (left) andψ̂2(x) (right)0246810x 104−1000−50005001000−1000−50005001000timeW1 real partW1 complex part0246810x 104−5000500−400−2000200400600timeW2 real partW2 complex partFIG. 3 – The wavelet coefficients in the orthosympathetic band (left) and in the parasympathetic band (right) of the samehealthy subject during a period of 24 hours3 Segmentation analysisWe assume that the process{Xt} is abruptly changing and is characterized by a parameterθ ∈ Θ that remains constantbetween two changes. We use this assumption to define our contrast functionJ(τ ,X). LetK be some integer and letτ = {τ1, τ2, . . . , τK−1} be an ordered sequence of integers satisfying0 < τ1 < τ2 < . . . < τK−1 < n.For the detection of changes in the mean and variance of a sequence of random variables, i.e., a change in distribution,Detection of Change–Points in the Spectral Densitywe use the following contrast function, based on a Gaussian log–likelihood function :Jn(τ ,X) =1nK∑k=1‖Xτk − X̄τk‖2σ̂2k+ nk log(σ̂2k). (5)wherenk = τk − τk−1 is the length of segmentk, σ̂2k is the empirical variance computed on that segmentk, σ̂2k =n−1k∑τki=τk−1+1(Xi − X̄)2, andX̄ is the empirical mean ofX1, . . . , Xn.In fact, we minimize a penalized version of this contrast function, the penalty parameter being selected in an adaptiveway so that the obtained segmentation does not depend too much on the penalty parameter ; see also Birgé and Massart(2007) for another choice for the penalty parameter. Although this method has been devised for Gaussian processes, itempirically provides relevant results for non-Gaussian processes, e.g., financial time series, see Lavielle and Teyssière(2006) for further details.3.1 Orthosympathetic bandFor that frequency band, we obtain the following segmentation :τ = {28220, 33366, 71048},i.e., 13h40’50”, 15h06’36”, 1h34’38”.5h50’30" 13h40’50" 15h06’36" 1h34’38"050100150200250300350400450500FIG. 4 – Segmentation in the mean and variance (change in the distribution) of the modulus of the wavelet coefficients inthe orthosympathetic band3.2 Parasympathetic bandFor that frequency band, we obtain the following segmentation :τ = {11620, 21912, 28054, 31540, 36022, 40172, 52622, 70550},i.e., 9h04’10”, 11h45’42”, 13h38’04”, 14h36’10”, 15h50’22”, 17h00’02”, 20h27’32”, 1h26’40”.4 ConclusionThe example illustrated by Figures 4, 5 and 6 shows at aroundt = 28, 220 a simultaneous variation of both systems.But, one can also observe change-points in the orthosympathetic and parasympathetic bands occurring at different times.Bertrand, Teyssière, Boudet, Chamoux5h50’30" 9h04’10" 11h45’42" 13h38’ 14h36’ 15h50’ 17h00’02" 20h27’32" 1h26’40"050100150200250300350400FIG. 5 – Segmentation in the mean and variance (change in the distribution) of the modulus of the wavelet coefficients inthe parasympathetic bandIn the future, we will study the existence of a possible causality or sequentiality between these change–points in differentbands.5h50’30" 9h04’10" 11h45’42" 13h38’ 14h36’ 15h50’ 17h00’02" 20h27’32" 1h26’40"0100200300400500600700800900  Orthosympathetic bandParasympathetic bandHearthbeat processFIG. 6 – Mean of the modulus if the wavelets coefficients for the orthosympathetic and parasympathetic bands, and of theRR processX(t)Références[1] Bardet, J.M. and Bertrand, P.R. (2007). Identification of the multiscale fractional Brownian motion with biomecha-nical applications.J. Time Ser. Anal.28, 1–52.[2] Bardet, J.M., Bertrand, P.R. and Billat, V. (2008). Estimation non-paramétrique de la densité spectrale d’un processusgaussien observé à des instants aléatoires.Pub ISUP52, 123–138.[3] Birgé, L. and Massart, P. (2007). Minimal penalties for Gaussian model selection.Probab. Theory Related Fields138, 33–73.Detection of Change–Points in the Spectral Density[4] Boudet, G., Albuisson, E., Bedu, M. and Chamoux, A. (2004). Heart rate running speed relationships during exhaus-tive bouts in the laboratory.Can. J. Appl. Physiol.29(6), 731–742.[5] Chamoux, A. (1984). Le système Holter en pratique.Médecine du Sport58, 43-273, 54–284.[6] Cottin, F., Leprêtre, PM., Lopes, P., Papelier, Y., Médigue, C. and Billat V. (2006). Assessment of ventilatory thre-sholds from heart rate variability in well-trained subjects during cycling,Int. J. Sports Med.27, 959–967.[7] Diab, M.O., Marque, C. and Khalil, M.A. (2007). Classification for uterine EMG Signals : Comparison between ARmodel and statistical classification method.Int J. Computational Cognition5, 8–14,[8] Frisch, U. (1996).Turbulence. Cambridge University Press.[9] Goldberger, A.L. (2001). Heartbeats, hormones and healt : is variability the spice of life ?Am. J. Crit. Care Med.163, 1289–1290.[10] Ivanov, P.C., Amaral, L.A.N., Goldberger, A.L., Havlin, S., Rosenblum, M., Struzik, Z, and Stanley, H.E. (1999).Multifractality in human heartbeat dynamics.Nature399, 461–465.[11] Lavielle, M. and Moulines, E. (2000). Least-squares estimation of an unknown number of shifts in a time series.J.of Time Series Anal.21, 33–59.[12] Lavielle, M. and Teyssière, G. (2006). Detection of multiple change-points in multivariate time series.LithuanianMath. J.46, 287–306. [Détection de ruptures multiples dans des sériestemporelles multivariées.] (Translation inFrench) Liet. Mat. Rink.46, 351–376.[13] Mallat, S. (1998).A wavelet tour of signal processing. Academic Press.[14] Task force of the European Soc. Cardiology and the NorthAmerican Society of Pacing and Electrophysiology(1996). Heart rate variability. Standards of measurement,physiological interpretation, and clinical use.Circulation93, 1043–1065.SummaryWe propose a new method for estimating the change-points of heart rate in the orthosympathetic and parasympatheticbands, based on the wavelet transform in the complex domain and the study of the change-points in the moments of themodulus of these wavelet transforms. We observe change-points in the distribution for both bands.

Detection of Change--Points in the Spectral Density. With Applications to ECG Data

Abstract

Similar works

Full text

Available Versions

HAL Clermont Université

HAL Descartes

Hal-Diderot

INRIA a CCSD electronic archive server