A new algorithm for wavelet-based heart rate variability analysis

A ne w algorithm for wa velet-based heart rate v ariabil ity analysis Constantino A. Garc´ ıa a, ∗ , Abraha m Otero b , Xos ´ e V ila c , Da vid G. M ´ ar quez a a Centr o Singular de In vestiga ci ´ on en T ecnolox ´ ıas da Informaci ´ on (CITIUS), University of Santiag o de Compostela, 15782, Santia go de Compostela, Spai n. b Department of Informatio n and Communicat ions Syste ms Engin eering, Universit y San P ablo CEU, 28668, Madrid, Spain. c Department of Computer Science , Unive rsity of V igo, Campus As Lagoa s s / n, 3200 4, Ourense , Spai n. Abstract One of the most prom ising non -inv asi ve markers of the acti vity of the autonomic nervous system is Heart Rate V ariability (HR V). HR V an alysis toolkits of ten provide spectral analysis techn iques using the Fourier transfor m, which assumes that th e hea rt rate series is station ary . T o overcome this issue, the Sho rt Time Fourier T ransfor m is of ten u sed ( STFT). Howe ver , the wavelet trans- form is th ought to be a mo re suitable tool f or analyzing non- stationary signa ls than the STFT . Giv en the lack o f su pport fo r wavelet-based analysis in HR V to olkits, such analy sis mu st be implemented by the research er . This has made this techniqu e un derutilized. This p aper presents a new algorith m to per form HR V power spectrum analysis based on the Maximal Overlap Discrete W a velet P acket Transform (MO DW PT). Th e alg orithm ca lculates the power in any spectral ban d with a given tolerance for th e band’ s bou ndaries. The M ODW PT decomp osition tree is pruned to av oid calculating un necessary wavelet coe ffi cien ts, thereby op- timizing execution time. The center of en ergy shif t correction is applied to a chieve optimum alignment of the wavelet coe ffi cients. This algorithm h as been imp lemented in RHR V , an o pen- source pa ckage fo r HR V analysis. T o the best o f our knowledge, RHR V is the first HR V toolkit with suppo rt fo r wa velet-based spectral analysis. K e ywor ds: Heart Rate V ariab ility, w avelet transform, w avelet p ackets, RHR V . 1. Introduction Heart Rate V ariability (HR V) refers to th e variation over time o f the intervals between consec- utiv e heartbeats. Since the heart rhyth m is m odulated by the a utonom ic nervous system (ANS), HR V is tho ught to reflect the activity of the sympathetic a nd para sympathetic bran ches of the ANS. Th e continuou s modulation of the ANS results in continuo us v ariations in hear t rate. One of th e m ost powerful HR V analysis techn iques is based on the sp ectral an alysis o f the time se- ries o btained f rom the distances between e ach p air of con secutive heartbeats. The HR V power spectrum is a useful tool as a predictor of multiple patholog ies [1], [2 ]. ∗ Correspondi ng auth or: T el.: + 34 8818 16391 . Fax: + 34 8818 16405. Email addr esses: constantin oantonio.garcia@usc.es (Constantino A. Garc´ ıa), abraham.ote ro@usc.es (Abraham Otero), anton@uvig o.es (Xos ´ e V ila), david.gonza lez.marquez@usc.es (David G. M ´ arquez) Prep rint submitted to Biomedical Signal Proc essing and Contr ol J uly 24, 2018 Akselrod et al. [3] describ ed three compo nents in the HR V p ower spectrum with phy siolog- ical relevance: the very low frequency (VLF) compon ent (frequencies below 0.03 Hz), which is modulated by the r enin-an giotensin system; the low fr equency (LF) compo nent (0.03-0 .15 Hz), which is th ought to be of both symp athetic and para sympathetic natur e; an d th e hig h f requen cy (HF) component (0.18 -0.4 Hz), which is related to the p arasympath etic system . At p resent, ther e is no absolu te consensus o n the precise limits of the boun daries of these th ree ba nds. In the literature we can find authors who use slightly di ff erent band s’ boundaries [4]. There exist sev eral HR V sp ectral an alysis techniques. These techniq ues may be classified as n onpara metric and parametric [5]. The main advantage of the no nparame tric methods is the simplicity an d speed of the algorithm used (The Fast Fourier Transform). The main advantage of the parametric me thods is that they give smoother spectral c ompon ents. Ho wever , paramet- ric methods presen t problem s regarding to co rrect mod el order selection . Although these tech- niques are widely used, they have no temp oral resolu tion. This sev erely lim its their ability to analyze non -stationary signals and transient ph enomen a. T o alleviate this limitation tempo ral windows are often u sed, so th at small segments o f the w hole signal are analyzed . Am ong th ese technique s we may high light the Short Time Fourier T ransform (ST FT) [6]. Howe ver , time- frequen cy resolu tion o f the STFT de pends on the spread of the win dow used. Thu s, the STFT has fixed time-fr equency resolution: high fre quency resolution implies poor time resolutio n and vice versa. Con versely to Fourier , the wa velet tra nsform performs time-fre quency an alysis and it is recognize d as a powerful too l to study non-stationar y sign als [7]. HR V analy sis to olkits such as Kubio s HR V [8] o r aHR V [9] o nly enable HR V spectral analy- sis ba sed on the Fourier tran sform or parametric methods. T o the best of our kno wledge, the only option for u sing th e wa velet tran sform in HR V analysis is to manually imp lement the algorithms, probab ly with the support of some general wa velet libra ry . This is tedio us, and prone to error . Although some research ers h av e d one this [10], [11], m any mo re (especially those with a med- ical back groun d) choose to u se Fourier-based too ls, ev en when they know that the sign al being analyzed is non -stationary . A query in the PubMed database with the term s “hear t rate variabil- ity Fourier transform ” retu rns 660 articles, wh ile a q uery with the terms “h eart rate variability wa velet transfo rm” only retur ns 145 ar ticles. The lack o f to ols for car rying out HR V analysis using the wa velet tra nsform has made th is poten tially sup erior analysis tec hnique underu tilized in comp arison with t he Fourier transform. In this paper we present an algorithm to perform HR V power spe ctrum analysis based on the Maximal Overlap Discrete W avelet P acket Transform (MODWPT). The algorithm calculates the spectrogra m in any freq uency b and, allowing a certain to lerance for the p osition of the b and’ s bound aries. T he algo rithm has be en validated over simulated an d real RR series. Its cap ability fo r identifyin g fast chan ges in the RR series’ spec tral compo nents h as been compa red with the STFT and a windowed version of the Burg method, sho wing that these techniq ues miss some transient changes that are su ccessfully identified b y the wavelet transfo rm. W e have imp lemented the algorithm in RHR V , an open- source packag e f or HR V analysis publicly av ailable o n th e In ternet. A previous version of this algor ithm was published in [12]. Section 2 starts with a b rief revie w of the wav elet transfo rm, with p articular attention to the MOD WPT , and then introdu ces our algorithm to per form HR V power spec trum an alysis. Section 2.5 pr ovides a shor t d escription of the implem entation o f the algorithm in the RHR V package. Section 3 p resents a compariso n between our algo rithm, the STFT and th e windowed Burg method over simulated and real RR series. Finally , the results of this pap er are discussed and some conclusio ns are gi ven. 2 2. Material and methods A brief r evie w of some important wavelet con cepts f or our algor ithm is now given. A w avelet is a small wa ve ψ ( t ) (oscillating func tion) that is well co ncentrated in time. This functio n m ust have unitary norm k ψ k = 1 and verify the so- called adm issibility condition : R ∞ −∞ ψ ( t ) d t = 0 . ψ ( t ) can be translated and dilated in time, yielding a set of wa velet functions: ψ u , s ( t ) = 1 √ s ψ  t − u s  , (1) where s > 0 is a d ilation factor , and u is a real number representing th e translations. As ψ generates all ψ u , s function s, it is called mo ther wa velet. A contin uous wavelet transform measures the tim e-frequ ency variations of a signal f by correlating it with ψ u , s W f ( u , s ) = Z ∞ −∞ f ( t ) ψ ∗ u , s ( t ) d t . (2) In o rder to make th e wav elet tra nsform imp lementable on a c omputer, bo th dilatio n and translation factors must be discretized. This can be achie ved as follows: ( ψ j , n = 1 √ 2 j ψ t − 2 j n 2 j !) j , n ∈ Z . (3) This family is an orth onorm al basis of L 2 ( R ). Orthog onal wav elets dilated b y 2 j can be used to study sign al variations at the resolu tion 2 − j . Thu s, these families of wa velets o riginate a multire solution sign al ana lysis. M ultiresolution analy sis p rojects sign als at various resolutio n spaces V j . Each V j space contain s all possible appro ximation s at the r esolution 2 − j . Thus, each decomp osition level incr eases the spectral resolution o f the d ecompo sition, at the expense of losing temp oral r esolution. Let { V j } j ∈ Z be a mu ltiresolution appro ximation verifying V j + 1 ⊂ V j ∀ j ∈ Z and let W j be the orthogonal complement of V j in V j − 1 : V j − 1 = V j ⊕ W j . According to [13], the families ( φ j , n = 1 √ 2 j φ t − 2 j n 2 j !) n ∈ Z and ( ψ j , n = 1 √ 2 j ψ t − 2 j n 2 j !) n ∈ Z (4) are an ortho norma l basis for V j and W j , respecti vely , for all j ∈ Z . ψ j , n are th e w avelet f unctions and φ j , n are the scaling functions. Thus, we can approx imate an y function f ǫ L 2 ( R ) at resolution 2 − j as P V j f = ∞ X n = −∞ h f , φ j , n i φ j , n = ∞ X n = −∞ a j [ n ] φ j , n (5) and the orthog onal projection of f on to detail space W j is: P W j f = ∞ X n = −∞ h f , ψ j , n i ψ j , n = ∞ X n = −∞ d j [ n ] ψ j , n . (6) where a j [ n ] and d j [ n ] are called the approx imation and detail coe ffi cients, respecti vely . 3 Mallat proved [7] th at both appro ximation and detail coe ffi cients may be calc ulated using a filter bank. Let h [ n ] and g [ n ] be th e FIR filters that will b e u sed to co mpute the appr oximation and detail coe ffi cien ts, respectiv ely . It has b een proven [ 13] that the filter h [ n ] = h 1 √ 2 φ  t 2  , φ ( t − n ) i and that g [ n ] = h 1 √ 2 ψ  t 2  , φ ( t − n ) i . g [ n ] and h [ n ] can be regard ed as an appro ximation to a high-p ass filter (the wavelet filter) an d to a low-pass filter (the scaling filter), r espectively . By applying recursively over the ap prox imation coe ffi cients th e same filtering operation followed by sub-sampling by two, we obtain th e multiresolution expr ession of f . This alg orithm, known as the pyramid algorithm , is the most e ffi cient way of comp uting the Fast Orthog onal W av elet T ransfor m (FO WT) [13]. 2.1. MOD WPT Giv en that the filterin g operation is on ly a pplied over the ap proxim ation coe ffi cien ts, the FO WT only p rovides inform ation on a limited s et of frequen cy band s which need not be th e on es used in the HR V analysis. A more suitable wa velet transform is n eeded fo r ou r alg orithm: the wa velet packet de composition (WPD). Instead o f dividing only th e app roximatio n co e ffi cients a j [ n ], both detail and app roximatio n coe ffi cien ts are decomp osed successi vely by apply ing high pass and low pass filters to each set of coe ffi cien ts. Among the WPD tran sforms we h av e cho sen the MODWPT [14] because it is less sensitive to the starting point of the time series and it is applicable to non dyad ic sequences. Furthe rmore, the MODWPT av oids the su b-sampling step, and theref ore it has the same nu mber of wavelet coe ffi cients in every d ecompo sition lev el. This simplifies com putations in volving di ff e rent d e- composition levels. The j th lev el of the MODWPT decomp oses the frequen cy interval [0 , f s / 2], where f s is the sampling freq uency of the orig inal signal f , into 2 j equal width intervals (see Fig. 1). Thus, the n th node (beginning at ze ro) in the j th lev el of the decomp osition tree, th e ( j , n ) n ode, is associated with the frequ ency in terval f s 2 j + 1 [ n , n + 1]. Each node will hav e N wa velet coe ffi cients associated, N being th e length o f the sampled signal f . The N -dimensional vector W j , n will d enote the N wa velet coe ffi cients associated with the ( j , n ) n ode. W 0,0 j=0 j=1 j=2 j=3 h g W 3,4 W 3,5 W 2,2 g h W 3,6 W 3,7 W 2,3 h g W 3,0 W 3,1 W 2,0 g h W 3,2 W 3,3 W 2,1 0 1/8 1/16 3/8 1/4 3/16 7/16 5/16 1/2 W 1,0 h g W 1,1 g h h g 0 1/8 3/8 1/4 1/2 0 1/4 1/2 0 1/2 Figure 1: MODWPT decomposition tree with the nodes select ed to cov er the band [0 , 7 / 16] Hz. W 0 , 0 represent s the origina l signal, f ( t ) . 4 MODWP T coe ffi cien ts f ulfill that: k f k 2 = 2 j − 1 X n = 0 k W j , n k 2 ∀ j . (7) Therefo re, gi ven a frequ ency ban d  f 1 , f 2  = f s 2 j + 1 [ k , k ′ + 1] , f s being the samp ling freque ncy and k , k ′ and j integers, the spectra l power in  f 1 , f 2  , P (  f 1 , f 2  ), can b e calculated f rom the approp riate MOD WPT coe ffi cients. W e just n eed to find the node s ( j , k ) and ( j , k ′ ) and comp ute the spectral power in the band  f 1 , f 2  as: P (  f 1 , f 2  ) = k ′ X n = k k W j , n k 2 . (8) 2.2. F inding a pr oper cover for a s et of fr equen cy bands Equation (8) can only b e a pplied to ba nds that can be written as f s 2 j + 1 [ k , k ′ + 1], f s being the sampling f requen cy an d k , k ′ and j integers. In an HR V spectral analysis the user can be interested in bands that cannot be written this way . This for ces us to perm it a certain er ror when w e tr y to cover the bands with co e ffi cients o btained fr om the MODWPT decom position. Let [ f l , f u ] be the ban d, and let ǫ l and ǫ u be the maximu m errors allowed fo r the b eginning and the ending o f the band, respectively . W e need to find a nod e ( j , n ) whose lower frequ ency correspo nds roughly to f l with the tolerance allowed b y ǫ l , i.e.:      f l − f s 2 j + 1 n      ≤ ǫ l . (9) Analogou sly , we also need to find a n ode ( j ′ , n ′ ) wh ose uppe r frequ ency corr esponds rough ly to f u with the tolerance allowed b y ǫ u :      f u − f s 2 j ′ + 1 ( n ′ + 1)      ≤ ǫ u . (10) W e shall refer to (9) and (10) as the cover conditions . The le vel j of the d ecompo sition tree in which the no de ( j , n ) tha t fulfills (9) is fou nd needs not be th e sam e as the level j ′ in wh ich the n ode ( j ′ , n ′ ) th at f ulfills (10) is fo und. But, (8) requires that j = j ′ . T o av oid this pro blem we may think that, after the nod es ( j , n ) and ( j ′ , n ′ ) have been found , we co uld descen d the node that is at the high er level, dec omposing it to the level of the other node [12]. Howe ver , when descendin g to lo w le vels in the frequency decom position tree, frequen cy r esolution increases a nd temporal resolution d ecreases because of the Gabor- Heisenberg un certainty prin ciple fo r sign als: ∆ t ∆ f ≥ 1 / 2 [ 15]. Fur thermor e, wa velet coe ffi cients su ff er a circular time shif t as a result of the use of wa velet filters. As we descend to lo wer le vels, the circular shift will be more prono unced as we perfo rm mo re convolutions. Therefo re, to o btain a good tempor al resolution, we should a v oid descen ding to deep le vels of the tree. Let’ s suppo se that, w hen loo king for a cover for the band of inter est  f 1 , f 2  , the nod es ( j , n ) and ( j ′ , n ′ ) are selected. These node s must fu lfill th e cover condition s (9) and (1 0) with f l = f 1 and f u = f 2 . Gen erally , the nod es ( j , n ) a nd ( j ′ , n ′ ) are at d i ff erent lev els (i.e., j , j ′ ). Further- more, these nodes cover the band f s 2 j + 1 [ n , n + 1] ∪ f s 2 j ′ + 1 [ n ′ , n ′ + 1], but the band f s h n + 1 2 j + 1 , n ′ 2 j ′ + 1 i = h f ′ 1 , f ′ 2 i remains uncovered ( see nodes (1 , 0) and (3 , 6) of Fig. 1). Th e problem has been redu ced 5 to findin g a cover for the ban d h f ′ 1 , f ′ 2 i . Thu s, new n odes will be selected fulfilling th e cover condition s (9) an d (10), with f l = f ′ 1 and f u = f ′ 2 . The complete cover of  f 1 , f 2  can be achieved applying this t echn ique recursi vely . The set of nodes resulting from this cover will be referred to as Γ . There may be overlap between the band s associated with the nodes ( j , n ) and ( j ′ , n ′ ) (the original n odes fulfilling the c over co nditions ). This overlap will o ccur if a n ode is the par ent of the other no de or if b oth nodes are eq ual. If ( j , n ) = ( j ′ , n ′ ), the n ode covers all the ba nd and the search has ended . If the nod e ( j ′ , n ′ ) is a child o f the no de ( j , n ) , we sh all replace th e node ( j , n ) with its (unique) child verifying (9): the node ( j + 1 , 2 n ). Thu s, the algo rithm continues with the nodes ( j + 1 , 2 n ) and ( j ′ , n ′ ). The p rocess shall be repeated until ther e is no overlap b etween the two n odes fulfilling the cover con dition . Similarly , if the node ( j , n ) is a child of the node ( j ′ , n ′ ), we replace the nod e ( j ′ , n ′ ) with the node ( j ′ + 1 , 2 n ′ + 1) . Figure 1 illustrates the c over process for the ban d [0 , 7 / 16 ] Hz with f s = 1 Hz an d ǫ u = ǫ l = 0 . 01 u sing the “nodes at th e same le vel” criteria used in [1 2] (dash ed node s) and th e “no des at di ff erent lev els” criteria proposed here (dark nodes). There exist multiple covers of any band  f l , f u  . Our algo rithm builds a cover using the no des of the higher le vels of th e MOD WPT tree. Th is selection criteria leads to a cover that min imizes temporal shift and maximizes tem poral resolution . Furthermo re, this cover also minimizes the number of nod es. Equation (8) ca nnot be applied to this decomposition , because the decomposition uses no des from di ff erent levels. Howev er, b ecause of the wa velet coe ffi cients proper ty k W j , n k 2 = k W j + 1 , 2 n k 2 + k W j + 1 , 2 n + 1 k 2 , (11) we may write (8) as P (  f 1 , f 2  ) ≈ k ′ X n = k k W j , n k 2 = X W j , n ǫ Γ k W j , n k 2 . (12) In prac tical application s, the com plete MODWPT decomposition tree will be rarely need ed. In or der to improve the execution time of wa velet analysis we have dev elope d a prunin g algo rithm that we call Pruned MODWPT (PMODWPT). Th e PMODWPT , instead of expan ding every nod e of the decom position tree will on ly calculate th ose ch ildren requ ired to cover a frequen cy band. Figure 2 illu strates the pru ning perform ed when calculating the power of the b and [0 , 3 / 8 ] Hz being f s = 2 Hz. W 0,0 j =0 0 1 j=1 W 1,1 W 1,0 0 1/2 1 j=2 W 2,0 W 2,1 0 1/2 1/4 j=3 W 3,0 W 3 ,1 W 3,2 W 3,3 0 1/8 1/4 1/2 3/8 Figure 2: Prune procedure using the PMODWPT . The crosses indica tes which node s have been pruned. 2.3. Shift corr ection As a consequen ce of the fr equency-d epende nt phase respo nse of the wa velet filters, it is necessary to keep track of which wa velet coe ffi cients contain the energy con tribution of a g iv en 6 portion of the signa l being analy zed (see Fig. 3). Th at is, the wavelet coe ffi c ients must be corrected so that they can be accurately aligned with the original time series. Hess-Nielsen and W ickerhauser pro posed in [16] a techniqu e to comp ute shift correction s exact for linear phase filters, and showed how to estimate th e p erturbatio n that a d eviation from linear p hase prod uces. This suggests that the meth od will pr oduce better results with filters whose phase do es not de viate muc h f rom a linear phase filter . Thu s, this algorithm can b e used to compute app roximate shift corr ections for a w ide rang e of wa velets, inc luding least asymm etric and extremal phase Daubechies, Symmlet or Coiflet wa velets. T ime shifts prop osed in [16] are based on the n otion of the “center of energy ” of a filter . Let a l be a filter of length l . Its center of energy , E { a l } , is giv en by: E { a l } = l − 1 P n = 0 n · a l [ n ] 2 l − 1 P n = 0 a l [ n ] 2 (13) The time shif t need ed for the node ( j , n ) of th e MODWPT is always an advance (shift to th e left) of | p j , n | units that is given by [1 6] | p j , n | = C j , n ( E { g } − E { h } ) + (2 j − 1) E { h } , (14) where C j , n is the bit-reversal of the binary gra y code which enco des each of th e MODWPT tree nodes [16]. 0 500 1000 1500 0.00 0.04 0.08 power MOD WPT without shift correction 0 500 1000 1500 0.00 0.04 0.08 time (s) power MOD WPT with shift correction Figure 3: T ime shift using MODWPT wit hout (top) and with (bott om) “Ce nter of Ener gy” correction. The verti cal line s indica te where spectra l po wer should be. 2.4. W a velet-based HRV analysis Listing 1 shows the pseud ocode of our w av elet-based HR V analysis algorithm. Th e function waveletAnalysis takes as inpu t parameters the interpolated and filtered R R series, the boundaries 7 of the frequ ency bands, th e wa velet to b e used in the an alysis, th e sampling frequen cy o f the RR series and the tolera nce allowed when c overing the frequ ency bands to be analyzed. The algorithm starts by findin g the covers for each of th e b ands in which the user wants to calcu late the spec tral power . Th e calculate Cover fun ction implemen ts the re cursive covering alg orithm presented in sectio n 2 .2. The getNode function selects the initial nodes fu lfilling th e cover condi- tions given by (9) and (10). It expand s the nod e whose frequency interv al contains the boundary frequen cy ( f ). If none of its c hildren verify the cover co ndition f or f , the child contain ing f is selected an d the pr ocess is repeated. T he fi llGapsBetweenNod es fun ction completes the cover giv en a pair o f initial nodes: ( j , n ) and ( k , m ) th at fu lfill (9) and (10), respectively . Th is fun ction distinguishes the fou r cases discussed in section 2.2: ( j , n ) = ( k , m ); ( k , m ) is a child of the ( j , n ) node; ( j , n ) is a child of ( k , m ) or the general case where there is an uncovered band . Once the cover for the band s has been found, the wa velet coe ffi c ients are computed using the PMODW PT . Note th at we u se a single d ecompo sition tree to calculate the po wer in all the band s. The shift corr ection descr ibed in section 2 .3 is comp uted using shiftCorr ectio n . The tar get nodes are spe cified to a void app lying corre ctions to nodes th at will not be used in subseque nt calcula- tions. Since the wavelet and scaling filters are in volved in (14), the wa velet used to perfor m the analysis m ust also be specified in the shiftCorr ection f unction. Finally , the sp ectrogram o f each band is comp uted with P W j , n ǫ Γ | W j , n | 2 by the getSpectr ogram fun ction. Listing 1: Pseudocode of our wav elet -based HR V analysis algorithm. w a v e l e t A n a l y s i s ( r r , V LFmin , V LF ma x , LF min , LF ma x , H F min , H F ma x , w avelet , f s , ← ֓ bandt olerance ) { V LF nod e s = c a l c u l a t e C o v e r ( [ V L F min , V L F ma x ] , f s , band t olerance ) LF nod e s = c a l c u l a t e C o v e r ( [ L F min , L F ma x ] , f s , band t olerance ) H F node s = c a l c u l a t e C o v e r ( [ H F min , H F ma x ] , f s , band t olerance ) wrr = P M O D W P T ( r r , wav elet , V LF node s ∪ L F nod e s ∪ H F nod e s ) wrr = s h i f t C o r r e c t i o n ( w rr , wavelet , V LF nod e s ∪ L F nod e s ∪ H F nod e s ) V LF power = g e t S p e c t r o g r a m ( w rr , V LF nod e s ) LF power = g e t S p e c t r o g r a m ( wr r , LF nod e s ) H F power = g e t S p e c t r o g r a m ( wr r , H F node s ) r e t u r n ( V LF power , L F power , H F power ) } c a l c u l a t e C o v e r ( [ f l , f u ] , f s , er r or ) { i f ( f u == f l ) r e t u r n ( {} ) W j , n = g e t N o d e ( f l , f s , er r or ) W k , m = g e t N o d e ( f u , f s , er ror ) r e t u r n ( f i l l G a p s B e t w e e n N o d e s ( W j , n , W k , m , f s , er r or ) ) } g e t N o d e ( f , f s , er ror , ty pe ) { i = 1 nodeT o E x pand = 0 w h i l e ( T RU E ) { b j = nodeT oE x pand ∆ = f s 2 i + 1 f o r ( j i n (2 · b j ) : (2 · b j + 1) ) { inter val = [ j ∆ , ( j + 1) ∆ ] i f ( f ∈ inte rval ) nod eT oE x pand = j i f ( inter val v e r i f i e s b a n d c o v e r c o n d i t i o n f o r f a n d er r or ) r e t u r n ( W i , j ) } i = i + 1 } } f i l l G a p s B e t w e e n N o d e s ( W j , n , W k , m , f s , er ror ) { i f ( W j , n == W k , m ) r e t u r n ( W j , n ) i f ( W k , m i s C h i l d o f W j , n ) r e t u r n ( f i l l G a p s B e t w e e n N o d e s ( W j + 1 , 2 n , W k , m , f s , er r or ) ) 8 i f ( W j , n i s C h i l d o f W k , m ) r e t u r n ( f i l l G a p s B e t w e e n N o d e s ( W j , n , W k + 1 , 2 m + 1 , f s , er r or ) ) med iumN od e s = c a l c u l a t e C o v e r ( f s h n + 1 2 j + 1 , m 2 k + 1 i , f s , er r or ) r e t u r n ( n W j , n , W k , m o ∪ med iumN ode s ) } g e t S p e c t r o g r a m ( w rr , nod e s ) { S { x } [ m ] = 0 f o r ( W j , n i n node s ) { d j , n [ m ] = g e t W a v e l e t C o e f f i c i e n t s ( wr r , W j , n ) S { x } [ m ] = S { x } [ m ] + | d j , n [ m ] | 2 } r e t u r n ( S { x } [ m ] ) } 2.5. Implementatio n o f the algorithms in RHRV RHR V is an open -source pack age fo r the R environment for statistical comp uting that com - prise a complete set of tools for HR V a nalysis. Fu rther details about RHR V p ackage may be found in [17]. Here we sh all on ly describe the RHR V fu nctionality which is related to the algo- rithm presented in this paper . RHR V imports data files contain ing h eartbeat positions. Suppo rted form ats in clude ASCII ( LoadBea tAscii functio n), EDF ( Lo adBeatE DFPlus ), Polar ( LoadB eatP olar ), Suun to ( Loa d- BeatSuu nto ) and WFDB data files ( Loa dBeatWFDB ) . T o com pute the instantan eous heart rate series BuildNIHR can be used. A filtering operation c an be carried out in ord er to eliminate outliers or spurio us points p resent in the HR time series with unacc eptable physiolo gical v alues ( F ilterNIHR ). A un iformly sampled heart rate signal (with equ ally spaced values) is obtained using Interpola teNIHR . Spectral power HR V an alysis is perfo rmed w ith the CalculateP owerBan d fu nction. This function co mputes the spectrog ram of the heart r ate series in th e UL F , VLF , LF and HF fr e- quency ba nds. Boundar ies of th e bands may be ch osen by the user . If bo undaries are no t spec- ified, default values are u sed: U LF , [0 , 0 . 03] Hz; VLF , [0 . 03 , 0 . 05 ] Hz; LF , [0 . 0 5 , 0 . 15 ] Hz; HF , [0 . 15 , 0 . 4] Hz. Until now , Calc ulateP owerBand used the STFT to compu te the spe ctral power (window size, win dow shift and zero padd ing may be spe cified b y the u ser). The wavelet analysis algorithm p resented in this paper was included in this func tion in such a way that we maintain backward com patibility . Th us, both Fourier and wavelet a nalysis may be used with the Calcu- lateP owerBand function. T ype of an alysis c an be selected b y the user by specify ing the type parameter (“fou rier” or “w avelet”). When using wa velet analysis, in addition to the frequ ency band s, an error for the boundar ies (default value is 0.1 in absolute terms) an d a mother wav elet can be specified b y the u ser . Some of th e most used wa velets are a vailable: “h aar”, extremal phase (“ d4”, “d 6”, “d8” and “d16”) and the least asymmetric (“la8”, “la16 ” and “la20 ”) Daubechies, and th e best localized (“b l14” and “bl20”) am ong o thers. Default value is “d4”. Listing 2 sho ws all the RHR V c ode required to perfor m a ty pical wa velet-based spectral analysis. Listing 2: HR V wave let-based analysis in RHR V . m d = C r e a t e H R V D a t a ( ) m d = L o a d B e a t A s c i i ( m d , ” B e a t P o s i t i o n s . b e a t s ” ) m d = B u i l d N I H R ( m d ) m d = F i l t e r N I H R ( m d ) m d = I n t e r p o l a t e N I H R ( md , f r e q h r = 4 ) m d = C r e a t e F r e q A n a l y s i s ( m d ) 9 m d = C a l c u l a t e P o w e r B a n d ( md , i n d e x F r e q A n a l y s i s = 1 , U L F m i n = 0 , U L F m a x = 0 . 0 3 , V L F m i n ← ֓ = 0 . 0 3 , V L F m a x = 0 . 0 5 , L F m i n = 0 . 0 5 , L F m a x = 0 . 1 5 , H F m i n = 0 . 1 5 , H F m a x = 0 . 4 , ← ֓ t y p e = ” w a v e l e t ” , w a v e l e t = ” d 4 ” , b a n d t o l e r a n c e = 0 . 0 1 ) P l o t P o w e r B a n d ( m d , i n d e x F r e q A n a l y s i s = 1 ) 3. Results 3.1. T emp oral r e solution In o rder to compar e the te mporal resolution of th e HR V a nalysis technique s based o n our algorithm with those based on the Fourier transfo rm (the STFT) and par ametric estimation (th e windowed Burg metho d), simu lated RR series will b e used. Simu lated signals are used instead o f real sign als beca use when using r eal sign als, it is d i ffi cult to kn ow which is the “corr ect” result. Therefo re, we use simulated signals to know exactly what spectral componen ts they ha ve at each instant. The I ntegral Pulse Frequency Modulation (IPFM) model [18], [19] is a widely acc epted technique used to genera te RR series. The IPFM mo del simulates the sino-atr ial node (SA) modulatio n by using a mo dulating signal m ( t ) and the SA fu nction as trigger of the car diac contraction by using a threshold ˆ T . W e shall u se a signal with several fast (every 1 6 seconds) spectral changes as modulatin g sign al: m ( t ) =                0 . 3 s in (2 π · 0 . 09375 · t ) 0 ≤ t < 16 s 0 . 3 s in (2 π · 0 . 03125 · t ) 16 ≤ t < 32 s 0 . 3 s in (2 π · 0 . 09375 · t ) 32 ≤ t < 48 s 0 . 3 s in (2 π · 0 . 03125 · t ) 48 ≤ t < 64 s 0 . 3 s in (2 π · 0 . 09375 · t ) 64 ≤ t < 80 s The idea is to test the time-frequ ency tran sforms in a non-station ary scenario . T o obtain a more realistic simulation, white noise was added to the m ( t ) signal. Figure 4 shows how the wavelet transfor m co rrectly finds m ost of th e spectral p ower in th e first band between 16 s and 32 s, and b etween 48 s and 64 s; and in th e second band between 0 s and 16 s, between 32 s and 48 s and between 64 s and 80 s, while the other methods cannot track the spe ctral chang es. The wav elet analy sis was perf ormed using a lea st asy mmetric Daubech ies filter of width 8. Th e toler ance was set to 0. 01. The selection of the window parameters f or the STFT and the Burg m ethod is not trivial. Note that the min imum size of the window for the STFT should be app roximately T ≈ 1 / 0 . 031 25 = 32 s . Howe ver , there is a spectral cha nge ev ery 16 seco nds. The results shown in Fig. 4( a) were obtain ed using a 30-seco nd window with a 1-secon d shift. In the Burg method , we may use smaller windows, provided that we take enoug h points to estimate the model in eac h segment. The parametric an alysis shown in Fig. 4(b) was perf ormed using a 16-secon d window with a 1-secon d shift. The selected value for the model order was 1 6 [20]. I t can be appreciated h ow the STFT and the Burg method canno t track changes on this signal. Using smaller length analysis windows did not improve th ese results. 10 0 20 40 60 80 1000 3000 5000 7000 time (s) power [0,0.0625] Hz 0 20 40 60 80 1000 3000 5000 time (s) power [0.0625,0.125] Hz (a) Spectrogram using STFT HR V an alysis. 0 20 40 60 80 0.1 0.2 0.3 0.4 0.5 0.6 time (s) power [0,0.0625] Hz 0 20 40 60 80 0.05 0.15 0.25 time (s) power [0.0625,0.125] Hz (b) Spectrogram using parametric HR V analy- sis. 0 20 40 60 80 0 50 150 250 time (s) power [0,0.0625] Hz 0 20 40 60 80 0 50 100 200 time (s) power [0.0625,0.125] Hz (c) Spectrogram using wave let HR V analysi s. Figure 4: Spectrogra m analy sis of the simulate d RR series. T able 1: Relati ve power per frequenc y and pe r zone for ST FT , windo wed Bur g method and wa velet analysis. (a) Fourier Analysis band \ tim e(s) [0,16) [ 16,32 ) [32,48 ) [48,64 ) [64,80) VLF 0.209 0.235 0.067 0.217 0.272 LF 0.128 0.147 0.384 0.224 0.117 (b) Parametric Analysis band \ tim e(s) [0,16) [ 16,32 ) [32,48 ) [48,64 ) [64,80) VLF 0.180 0.245 0.083 0.259 0.233 LF 0.228 0.124 0.288 0.167 0.193 (c) W av elet analy sis band \ tim e(s) [0,16) [ 16,32 ) [32,48 ) [48,64 ) [64,80) VLF 0.040 0.400 0.041 0.447 0.072 LF 0.270 0.071 0.344 0.077 0.238 11 For each of the t wo band s, and fo r each of the fi ve zones with di ff erent spectral components, we calculated the ratio of the power that the band presents in ea ch zo ne divided by the overall power of the band in the fiv e zo nes. Theoretically , the power in the VLF b and sho uld be d is- tributed am ong the 2 nd and 4 th zones, wher eas the power in the LF b and should be distributed among the 1 st , 3 rd and 5 th zones. Th erefore, the ide al ratios if pe rfect time-fr equency discrimi- nation is o btained are (0 , 0 . 5 , 0 , 0 . 5 , 0) and (1 / 3 , 0 , 1 / 3 , 0 , 1 / 3). T ables 1(a), 1(b ) a nd 1(c) show the real ratios for th e STFT , the Burg method an d w avelet an alysis, respecti vely . W e can see that the wav elet ratios are closer to the theoretical values. 3.2. Computation al bur den The top of Fig. 5 compares execution time ( as a function of the signal size) of HR V analysis algorithm s b ased on MODWPT , PMODWPT (our algo rithm) an d the STFT . The input signal to these algorithm s was gene rated ran domly . In o rder to co mpare Fourier with wa velet-based analysis, two config urations of th e STFT ty pically used on HR V analysis were selected. First Fourier analysis used a window size and a displacement value of 5 minutes and 3 0 s, respecti vely (“T ypical Fourier” in Fig. 5). The second Fourier analy sis u sed a shorter win dow in ord er to achieve a highe r temp oral r esolution. W in dow size and disp lacement took a value of 30 s and 2.5 s, re spectiv ely (“High Resolution Fourier” in Fig. 5). W a velet analysis was perfo rmed using least asymmetric Daubechies of width 8 ( “la8”) an d extremal pha se Daubech ies of width 4 (“d4”) since the e ffi ciency depends on the filter length. PMODW PT and “T yp ical Fourier” STFT are muc h more e ffi cient th an MODWPT an d “High Resolution Fourier” STFT when analy zing HR V signals (see the top o f Fig. 5). Th e botto m o f Fig. 5 shows that the per forman ce of the PMODW PT is comparab le to the “T ypical Four ier” analysis. 1e+05 2e+05 3e+05 4e+05 5e+05 0 5 10 20 30 Fourier Vs. wavelet efficiency Signal length in samples ex ecution time (s) PMODWPT (la8) PMODWPT (d4) MODWPT (la8) T ypical Fourier High resolution Fourier 5e+05 6e+05 7e+05 8e+05 9e+05 1e+06 1 2 3 4 5 Signal length in samples ex ecution time (s) Figure 5: Performance of HR V analysis using PMODWPT , MODWPT and STFT . 3.3. V alida tion on r ea l data W e have tested the a lgorithms pr esented in this p aper o n th e reco rdings of the Apnea ECG database u sed in the Phy sioNet / Compu ters in Cardiology Challenge 2000 [ 21]. Obstructive Sleep 12 Apnea-Hy popne a (OSAH) Syn drome is a sleep- breathin g disord er chara cterized by the presence of total (ap neas) and / or partial (hypo pneas) cessations of respira tory airflow while the patient is asleep. There is an in terest in developing low-cost OSAH screening tests that can be carried out in the patient’ s home an d that can decrease the workload of hospitals’ Sleep Units. This is related to the go al of the Comp uters in Cardio logy Challeng e 2000: to develop a diagno stic test for OSAH f rom a single E CG lead. T he dataset for th is challenge is d ivided into a lear ning set and a test set. E ach of these sets consist of 35 reco rdings of the modified lead V2 of the patien ts ECG recorded during nocturn al rest. Both tra ining and test sets are m ade o f 20 r ecordin gs o f patients su ff ering from OSAH, 5 r ecording s of patients wh o were o n th e bord erline between normality and OSAH, and 10 recor dings of control pa tients who did not su ff er from th e disorder . Each recor ding inclu des m inute by minute annotation s indicating the presenc e o r a bsence of apneas during that minute. The goal of the secon d p art of the challeng e, the one which will be addressed here, is to detect whether or nor the p atient has su ff ered an apnea du ring each minute of noc turnal rest. T o th is end we shall use two algorithms pr eviously pu blished in the bibliography which are based on the calculation o f a ratio between HR V spectral power in two di ff erent ban ds. Specifically , we have computed the Drinnan ratio ( R d ) [22] and the Otero ratio ( R o ) [23]. Th e Dr innan ra tio and the Otero ratio are defined as R d = P ([0 . 005 , 0 . 01 H z ]) P ([0 . 01 , 0 . 05 H z ]) and R o = P ([0 . 026 , 0 . 06 H z ]) P ([0 . 06 , 0 . 25 H z ]) , respectiv ely . Both ratios were computed with RHR V using both wa velet and Fourier an alysis. The ratios obtained for each minute of the rec ording were f ed to a suppo rt vector machine (SVM) [24]. The SVM was trained using the lear ning set and validated o n th e validation set of the Apn ea ECG Database. The scores (pe rcentage o f minutes labelled co rrectly) o btained in the minute by minute apnea classification using R d , R o and ( R d , R o ) as SVM p arameters when the spectr al power in the bands was calculated using wav elets wer e: 74.9%, 6 8.4% and 78. 5%, respectively . When using Fourier, th e scores were 71.4%, 66.8% and 75.4 %, respecti vely . W a velet-based analy sis p erform s slightly better than Fourier-based analysis in all scena rios. This may be due to the n on-stationa ry natu re o f the signal bein g analyze d, an d the fact th at the higher tempo ral resolution of th e wa velet analy sis can help minimize th e spectral co ntributions of apn eas whic h h av e occurr ed outside the min ute in q uestion, but c lose to the en d o f th e p revious minute or to the beginning of the following one. 4. Discussion W e h av e presen ted an algorithm to p erform HR V power spectrum analysis b ased on the MODWP T . The computatio nal load of the our algo rithm is comp arable to the load o f widely used STFT -based algorith ms. The algorith m has been validated over simulated RR series with known spectr al compo nents. W e have shown th at the STFT and the wind owed Burg method miss some quick changes that are successfully identified by the MOD WPT . These results sug gest that wa velet-based analysis is a better tool to analy ze fast transient pheno mena in the RR time series than other technique s based on windo wing. In order to obtain optimal temporal resolution, we should av oid descend ing to d eep lev els of the PMODWPT deco mposition tree. In RHR V , a warning is gene rated if th e band cover needed for the analysis requires expanding mo re than lo g 2 ( N / ( L − 1) + 1) levels, N bein g the number of samples of the signal an d L th e filter length. A carefu l selection of the freq uency bands to be analyzed provid es some contro l over t he depth of the tre e. For example, if the RR time series is sampled at 4 Hz, an d we want to ob tain the p ower in the band [0 . 27 , 0 . 5] Hz with a tolerance in the position of the ban d’ s bo undaries of 0.01, ou r algorithm wil l need to descend s even levels on 13 the tree. Howev er, to calculate the p ower in [0 . 26 , 0 . 5] Hz with the same toleran ce o ur algorithm only needs to d escend three levels. In this way , we obtain a g ood estimate of the power in the band [0 . 27 , 0 . 5] Hz with out compromising the temporal resolution of the results. A co rollary of the phe nomeno n d escribed in th e p revious p aragraph is that, in or der to achieve optimum tempo ral reso lution (and therefor e maximize the chance o f identifyin g fast tr ansient pheno mena), the spectral bands used in HR V analysis with wa velets will prob ably ha ve to di ff er from those traditio nally associated w ith VLF , L F and HF . This open s the qu estion of what may be the pathophysiolo gical significance of spectral b ands di ff erent f rom those which h av e already been widely studied in the literature. The mother wavelet used in HR V analysis also in fluences the time- freque ncy resolution be- cause it determin es the filter shape. Fur ther study on how the m other wavelet influen ces the re- sults of the HR V analysis is requ ired. For example, our initial tests suggest that shorter wa velets have greater tim e reso lution than the larger wavelets, wh ereas the latter have a better filtering behavior than the first ones. When perf orming a wavelet analysis, a careful choice of the band s is need ed to avoid heavy computatio ns, as well as a choice of the mo ther wavelet. Wh en using STFT an alysis we also have to cho ose the f requen cy bands (altho ugh we h av e m ore flexibility in the cho ice), as well as the window size, wind ow displacemen t and z ero paddin g size. Overall, mor e p arameter choices need to b e mad e when using the STFT . Furtherm ore, the selection of the STFT p arameters is complex, not only b ecause of the computation e ffi ciency , but also because the use of a certain window size influences the temporal resolu tion and the frequen cy bands a vailable to the an alysis. T o add ress this issue, some au thors use di ff erent windows for each frequen cy ban d. T his makes the analysis process hea vy . Mo reover , th e results obtain ed in each of the b ands are not comparable since they have been obtained using di ff eren t windows. The need for tuning a lower number of kno bs, and the possibility of using the same p arameters in the analysis of all the HR V fr equency band s are two additional advantages of using wavelets in H R V studies. The algorithm described in this paper has been imp lemented in the RHR V package for the R en viron ment in its 3.0 version. T o th e b est of our knowledge, RHR V is the first HR V an alysis toolkit that supp orts wa velet-based spectral analysis of the RR time ser ies. This sof tware can be freely downloaded fro m [25]. T he availability of this packag e will enable research ers to carr y out HR V p ower spectr um an alysis based on the wavelet tr ansform in a simple manner (see Listing 2). W e hope that this will help in crease the numb er of HR V studies that u se the h igher temporal resolution wa velet-based techniques. Acknowledgements This work was supported in p art by the Eur opean Regional Development Fund (ERDF / FEDER) under the project CN2012 / 151 of the Galician Ministry of Education. NO TICE: this is the authors version of a work that was accep ted for publication in Biomedi- cal Signal Proce ssing and Control. Changes resulting fr om the p ublishing process, such as peer revie w , editing, corre ctions, structu ral fo rmatting, an d o ther qu ality co ntrol mec hanisms may not be reflected in th is docume nt. Chan ges may have be en made to this work since it was sub mitted for publication . A definitive version was subsequ ently published in Biom edical Signa l Process- ing and Contro l, V olume 8 , I ssue 6, November 2013, Pages 5 4255 0, DO I:10.10 16 / j.b spc.201 3.05.006 14 References [1] M. Malik , A. Camm, Hea rt rate varia bility , Futura Pub. Co., Armonk, NY , 1995. [2] A. Jovic, N. Boguno vic, Eva luating and compari ng performance of featur e combina tions of heart rate v ariabil ity measures for cardiac rhythm classifica tion, Biomedical Signal Processing and Control 7 (2012) 245–255. [3] S. Akselrod, D. Gordon, F . Ubel, D. Shannon, A. Ber ger , R. Cohen, Powe r spectrum analysis of heart rate fluctua- tion: a quantitat iv e probe of beat-to-b eat cardio va scular control , Science 213 (1981) 220. [4] M. J. Lewi s, M. Kingsle y , A. Short, K. Simpson, Influence of high-f requenc y bandwid th on heart rate varia bility analysi s during physical ex ercise, Biomedical Signal Processing and Control 2 (2007) 34–39. [5] Hea rt rate va riabilit y: standards of measurement, physiologi cal interpre tation and clinical use. T ask Force of the European Soci ety of Cardiol ogy and the North Americ an Society of P acing and Electrop hysiology ., Circulati on 93 (1996) 1043–1065. [6] J. V ila, F . Palaci os, J . Presedo, M. Fernandez-Del gado, P . Felix, S. Barro, Time-fre quency analysis of heart-rate v ariabili ty , E nginee ring in Medic ine and Biol ogy Magazine, IEEE 16 (1997) 119–126. [7] S. Mallat, A theory for multiresolutio n signal dec omposition: The wa vele t representati on, IE EE Transacti ons on Patt ern Anal ysis and Machine Intelligenc e 11 (1989) 674–693. [8] K ubios HR V, Heart Rate Va riabilit y Analysis Softwa re, biosi gnal Anal ysis and Medica l Imaging Group, Uni v ersity of Eastern Finland, http: // k ubios.uef.fi / , website last accessed February 2012. [9] Ne vrokard -aHR V software, http: // www .nevro kard.eu / inde x.html , website last accessed February 2012. [10] P . S. Addi son, W avele t transforms and the ecg: a re vie w , Physiological measurement 26 (2005) R155. [11] A. Hoss en, B. Al-Ghunaimi, A wa vel et-based soft decisi on techn ique for screeni ng of patients with conge sti ve heart failu re, Biomedical Signal Processing and Control 2 (2007) 135–143. [12] C. Garc´ ıa, A. Otero, X. V ila, M. Lado, An open s ource tool for heart r ate v ariabili ty wav elet-base d spectral analysi s, in: Interna tional Joint Confe rence on Bi omedical Enginee ring Systems a nd T echnologies, BIOSIGN ALS 2012 , pp. 206–211. [13] S. Mallat, A wa velet tour of signal processing. Second Edition, Academic Pr , 1999. [14] D. Perci v al, A. W alden, W avele t m ethods for time serie s ana lysis, vo lume 4, Cambridge Uni v Pr, 2006. [15] D. Gabor , T heory of communication. part 1: The analysis of informatio n, Journal of the Institution of Electr ical Engineers-P art III: Rad io and Comm unicat ion Engineeri ng 93 (1946) 429–441 . [16] N. Hess-Nielsen, M. W ickerh auser , W av elets and time-freque ncy analysis, Proceedings of the IEEE 84 (1996) 523–540. [17] L. Rodr´ ıguez-Li ˜ nar es, A. M ´ endez, M. Lado, D. Olivi eri, X. V ila, I. G ´ omez-Co nde, An open source tool for heart rate vari ability s pectra l anal ysis, Computer methods and programs in biomedicin e 103 (2011) 39–50. [18] R. Berger , S. Akselrod , D. Gordon, R. Cohen, An e ffi cient algor ithm for spectral analysis of heart rate v ariabil ity , IEEE Transa ctions on Biomedical Engineering (1986) 900–904. [19] B. Hyndman, R. Mohn, A pulse modulato r model for pac emaker acti vity , in: Digest of the 10th Int. Conf. Med. & Biol. Eng, p. 223. [20] A. Boar dman, F . S. Schlin dwein, A. P . Rocha , A. L eite, A study on the optimum order of autoregressi v e models for heart rate va riabilit y , Physiolog ical measure ment 23 (2002) 325. [21] T . Penzel, G. Moody , R. Mark, A. Goldberge r , J. Peter , The A pnea-ECG data base, in: Compute rs in Cardiol ogy 2000, IEEE, pp. 255–258. [22] M. Drinnan, J. Allen, P . Langley , A. Murray , Detectio n of sleep apnoea from frequency analysis of heart rate v ariabili ty , in: Compute rs in Ca rdiology 2000, IE EE, pp. 259–26 2. [23] A. Otero, S. Dapena, P . Felix, J. Presedo, M. T arasc ´ o, A low cost scree ning test for obstruct iv e sleep apnea that can be perfo rmed at the patient ’ s home, in: IEEE Int ernationa l Symposium on Inte lligent Signa l Processing, 2009. WISP 2009, IEEE, pp. 199–204. [24] D. Meyer , Support vector machin es, The Interfac e to libsvm in package e1071. e1071 V ignette (2012 ). [25] The Comprehensi ve R Archi ve Network (CRAN), htt p: // cran.r -project.or g / web / packages / rhrv / index.html , website last accessed February 2012. 15 Figure captions • Figure 1: MODWPT decompo sition tree with the nod es selected to cover th e band [ 0 , 7 / 16 ] Hz. W 0 , 0 represents the original signal, f ( t ) . • Figure 2: Pru ne procedur e using the PMODWPT . The cro sses indicates which nodes have been pruned. • Figure 3: Time sh ift using MODWPT withou t (to p) and with (bottom) “Cen ter of Energy” correction . The vertical lines indicate where spectral power should be. • Listing 1: Pseudocode of our wa velet-based HR V analysis algorithm. • Listing 2: HR V wa velet-based analysis in RHR V . • Figure 4: Spectrogr am an alysis of the simulated RR series. • T able 1: Relati ve power per frequen cy an d per zone for STFT , windowed Burg metho d and wa velet analysis. • Figure 5: Performa nce of HR V analy sis using PMOD WPT , MODW PT and STFT . 16