Comp080504.qxd The Journal of Engineering Research Vol. 6, No. 2 (2009) 40-46 1. Introduction Heart failure is a common condition that usually devel- ops slowly as the heart muscle weakens and needs to work harder to keep blood flowing through the body. Heart fail- ure develops following injury to the heart such as the dam- age caused by heart attack, long-term high blood pressure, or an abnormality of one of the heart valves. Heart failure is often not recognized until a more advanced stage of heart failure, commonly referred to as congestive heart failure (CHF), in which fluid may leak into the lungs, feet, and in some cases the liver or abdominal cavity. Physicians often assess the stage of heart failure according to the New York Heart Association (NYHA) functional classification system. This system relates symptoms to everyday activities and the patient's quality of life. According to this system, heart failure is classified into 4 classes (Heart Failure Society of America, http://www.abouthf.org/questions_what_is_hf.htm): _____________________________________________ *Corresponding author’s e-mail: abhossen@squ.edu.om 1. Class I (Mild), symptoms with more than ordinary activity. 2. Class II (Mild), symptoms with ordinary activity. 3. Class III (Moderate), symptoms with minimal activity. 4. Class IV (Severe), symptoms at rest. Physicians often order a number of tests when explor- ing a possible diagnosis of heart failure. The most impor- tant of these tests is an echocardiogram. Echo cardiog- raphy is a noninvasive, entirely safe test that uses ultra- sound to image the heart as it is beating (Congestive Heart Failure, University of Maryland Medicine, http:/www.umm.edu/patiented/doc13diagnos.htm). Cardiac ultrasounds provide information about: * Accurate indications of valve function. * The amount of blood flow through the heart's chambers. * The location of the failure, whether has occurred on the left side, the right side, or both. Physicians use information from the echocardiography A Pattern Recognition Technique Based on Wavelet Decomposition for Identification of Patients With Congestive Heart Failure Abdulnasir Hossen*a and Bader Al-Ghunaimia Department of Electrical and Computer Engineering, College of Engineering, Sultan Qaboos University, P.O. Box 33, P.C 123, Al-Khoud, Muscat, Sultanate of Oman. Received 4 May 2008; accepted 15 February 2009 Abstract: A pattern recognition technique based on approximate estimation of power spectral densities (PSD) of sub-bands resulted from wavelet decomposition of R-R interval (RRI) data for identification of patients with Congestive Heart Failure (CHF) is investigated. Both trial and test data used in this work are drawn from MIT databases. Two standard patterns of the base-2 logarithmic values of the reciprocal of the probability measure of the approximated PSD of CHF patients and normal subjects are derived by averaging all corresponding val- ues of all sub-bands of 12 CHF data and 12 normal subjects in the trial set. The computed pattern of each data under test is then compared band-by-band with both standard patterns of CHF and normal subjects to find the closest pattern. The new technique resulted in an identification accuracy of about 90% by applying it on the test data. Keywords: Congestive heart failure, Pattern recognition, Wavelet decomposition, Soft-decision, Power spectral density ÊÉ≤àMC’G Ö∏≤dG õéY ≈°Vôe øY ∞°ûµ∏d áFõÛG äÉéjƒª∏d ájOOÎdG Ωõ◊G ≈∏Y áYRƒŸG IQ~≤∏d •É‰CG õ««“ á«æ≤J »ª«æ¨dG Q~H h Ú°ùM ô°UÉædG~ÑY áá°°UUÓÓÿÿGGÖ∏≤dG äÉ° Ñf ∫~©e Ò¨àe IQÉ°T’ äÉéjƒŸG áFõŒ øe áŒÉædG …OOÎdG ∞«£dG ΩõM ≈∏Y áYRƒŸG á«Ñjô≤àdG IQ~≤dG áaÉãµd •É‰C’G õ««ªàd á«æ≤J OÉéjG åëÑdG Gòg øª° àj : É°SÉe ~¡©e äÉfÉ«H I~YÉb øe ¢üëØ∏d áe~îà°ùŸG ∂∏Jh áHôéà∏d áe~îà°ùŸG äÉfÉ«ÑdG äòNCG .Úª«∏°ùdG ¢UÉî°T’G øY ºgõ««“ h ÊÉ≤àMC’G Ö∏≤dG õéY ≈°Vôe øY ∞°ûµdG πLCG øe ≈°VôŸ ɪg~MCG ,á«Ñjô≤àdG á«Ø«£dG IQ~≤dG áaÉãµd áHƒ°ùÙG á«dɪàMC’G ܃∏≤Ÿ 2-¢SÉ°SG »“QÉZƒd »æëæŸ Ú°SÉ«b Ú£‰ ¥É≤à°TG ” .I~ëàŸG äÉj’ƒdG ‘ »Lƒ∏æµàdG õà«°S ƒ°T 12 O~©d h ÊÉ≤àMC’G Ö∏≤dG õé©H ¢ jôe 12 O~©d I~FÉ©dG á«ÑjôéàdG äÉfÉ«Ñ∏d IôXÉæàŸG á«YôØdG Ωõë∏d IQ~≤dG ∫~©e OÉéjÉH ∂dPh Úª«∏°S ¢UÉî°TC’ ôNB’G ÊÉ≤àMC’G Ö∏≤dG õéY ÊÉ≤àMC’G Ö∏≤dG õéY »àdÉ◊ Ú≤à°ûŸG Ú«°SÉ«≤dG Ú£ªædG øe πc ™e áeõëH áeõM É¡àfQÉ≤e ” ~≤a ¢üëØdG â– ádÉ◊G äÉfÉ«H §‰ ÜÉ°ùM ~©H ¬fCÉa ¢üëØdG AÉæKCG ÉeCG .º«∏°S ¢üî°T .¢üëØdG äÉfÉ«H •É‰CG õ««“ ‘ %90 ¤G π°üJ ábO áÑ°ùf I~j~÷G á«æ≤àdG √òg â£YCG .É¡©e ∞æ°üàd ÜôbCG »g áÄa ájCG ¤G áaô©Ÿ º«∏°ùdG Ö∏≤dG h Gáá««MMÉÉààØØŸŸGG ääGGOOôôØØŸŸ .…OOÎdG ∞«£∏d IQ~≤dG áaÉãc ,¿ôŸG QÉ«àNC’G ,áFõÛG äÉéjƒŸG ,•É‰C’G õ««“ ,ÊÉ≤àMC’G Ö∏≤dG õéY : 41 The Journal of Engineering Research Vol. 6, No. 2 (2009) 40-46 for calculating the ejection fraction (EF), which is the per- cent of the blood pumped out during each heartbeat. EF is important for determining the severity of heart failure. People with a healthy heart usually have an EF of 50 per- cent or greater. Most people with heart failure, but not all, have an EF of 40 percent or less. The echocardiogram is the most accurate diagnostic test but an expensive one. So there is a need to a noninvasive simple test that helps in identification of patients who need an echocardiogram test. 2. Heart Rate Variability Heart rate variability (HRV) is referred to as the beat-to-beat variation in heart rate. Instantaneous heart rate is measured as the time in seconds between peaks of two consecutive R waves of the ECG signal. This time is referred to as the RRI. The beat-to-beat fluctuations in the rhythm of the heart provide us with an indirect measure of the heart health, as defined by the degree of balance in sympathetic and vagus nerve activity (Robert Nolan and Heart Rate Variability http://www.bio- medical.com/news_display.cfmnewsid=27). HRV analy- sis serves as a marker for cardiovascular disease because cardiac dysfunction is often manifested by systematic changes in the variability of the R-R interval sequence rel- ative to that of normal controls. The variation of heart rate accompanies the variation of several physiological activi- ties such as breathing, thermoregulation and blood pres- sure changes (Task Force of the European Society of Cardiology and the North American Society of pacing and Electrophysiology). Several HRV abnormalities have been described in patients with CHF. It has been shown that patients with heart failure have decreased HRV, particularly those who are at risk of cardiac death (Ponikowski et al. 1997). Analysis of the HRV using many signal processing tech- niques can help to detect abnormalities, by manipulating RRI data either in the time-domain or in the frequency domain (Teich et al. 2001). The frequency spectrum of the RRI data is divided into three main bands: 3. RRI Data Groups Two different groups of data are used as described below: 3.1 Trial Data The CHF records and Normal Sinus Rhythm (NSR) records are drawn from Massachusetts Institute of Technology (MIT) database, which includes 12 records from normal subjects (age 29-64 years) (Physionet, The MIT-BIH Normal Sinus Rhythm Database, http://www.physionet.org/physiobank/database/nsrdb) and 12 records from patients (age 22-71 years) with CHF (NYHA class 3 and 4) (Physionet, The BIDMC Congestive Heart Failure Database, http://www.phys- ionet.org/physiobank/database/chfdb). 3.2 Test Data This group contains 53 NSR (Physionet, Normal Sinus Rhythm RR Interval Database, http://www.physionet.org/physiobank/database/nsr2db) and 17 CHF (Physionet, Congestive Hert Filyre RR Interval Database, heep://www.physionet.org/phys- iobank/database/ch2db) recordings that are used to test the performance of the classification algorithm. The CHF subjects are selected from a larger set that contains 29 long-term recordings. The selected 17 records have CHF with NYHA class 3 and 4, and the remaining 12 records have NYHA class 1 and 2, and therefore have been excluded from the study. The subjects for the selected records are 8 men, aged 39 to 68, and 2 women aged 38 and 59; gender is unknown for the 7 remaining records, but aged between 35 and 64 years. The NSR data of this group contains 53 long-term (about 24 hours) RRI records. The subjects are 30 men; aged 28.5 to 76, and 24 women aged 58 to 73. The RRI data, for both trial and test groups, are gener- ated from the annotation file, which is associated with each record, using WFDB software (Physionet, an NIH- NCRR Research Resource, http://www.physionet.org/physiotools/wfdb.shtml). In the beginning, all records are truncated to a length of 75821 RR intervals, since the shortest record has such a length. Such truncation was also performed by other researchers (Teich et al. 2001). It was found in (Teich et al. 2001) that with the help of the wavelet transform standard deviation the normal sub- jects exhibited greater fluctuations at scale window between 16 and 32 samples than those afflicted with heart failure. It was possible to completely discriminate between 12 CHF and 12 NSR subjects that represent the trial data in this study. 4. Soft-Decision Technique A method was implemented for screening obstructive sleep apnea and normal controls using the estimated PSD by a soft-decision algorithm on decomposed sub-bands (Hossen, et al. 2003 and Hossen et al. 2005). The same method was also used to screen patients with CHF using power spectral ratios LF/HF and VLF/HF as screening parameters on MIT trial data only (Hossen and Al- Ghunaimi, 2004). Wavelet decompositions were used instead of subband decomposition to estimate the PSD of a signal (Hossen 2004). In this work a pattern recognition technique based on PSD estimation of wavelet decom- • The very low-frequency band (VLF): f ∈ (0.0033 – 0.04) Hz. • The low-frequency band (LF): f ∈ (0.04 – 0.15) Hz. • The high-frequency band (HF): f ∈ (0.15 – 0.4) Hz. 42 The Journal of Engineering Research Vol. 6, No. 2 (2009) 40-46 posed sub-bands of RRI data is implemented for the pur- pose of identification of patients having CHF from normal controls. A complete analysis is included on both MIT trial and test data sets. The work is examined by measur- ing the classification performances and including a con- sistency test. 4.1 Wavelet-Decomposition and Band-Selection The block-diagram of the wavelet decomposition is shown in Fig. 1. The input signal x(n) of length-N is first filtered by low- pass (LPF) and high-pass (HPF) filters and then down- sampled by a factor of 2 to produce both the "approxima- tion" a(n) and the "details" d(n). Assuming haar-filters are used, a(n) and d(n) can be obtained by: (1) If there is no information about the energy distribution of the input sequence, a band-selection algorithm (Hossen and Heute, 1993) can be used to decide (as a hard deci- sion) which band is to be computed or kept for more pro- cessing. This method depends on the energy comparison between the low- and high-frequency subsequences after the down sampling in Fig. 1. (2) According to the sign of B, the decision is taken: If B is positive, the low-frequency band is considered, and if B is negative, the high-frequency band is considered. Since we are not interested in the value of B, but only in its sign, a more-simpler equation than Eq. 2 can be obtained approximately as (Hossen and Heute, 1993): (3) 4.2 Estimation of Power Spectral Density The decomposition stages in Fig. 1 are computed with all branches up to a certain stage m to obtain 2m subbands. All estimator results up to stage m are stored, and a prob- ability measure is assigned to each path (frequency band) to bear the primary information. If J(L) is the assigned probability of the input signal being primarily low-pass, the number J(H) = 1- J(L) is the probability that the signal is primarily high-pass. One simple way to make the probability assignments is to use the ratio of the number of positive comparisons between | a(n) | and | d(n) | in Eq. (3) to the total number of comparisons for a given stage. At the following stage, the resulting estimate can be interpreted as the conditional probability of the new input sequence containing primarily low (high) frequency com- ponents, given that the previous branch was predominant- ly low (high)-pass. Using this reasoning and laws of prob- ability, the assignments for the probability measure of the resulting sub-bands should be made equal to the product of the previous branch probability and the conditional probability estimated at a given stage. Fig. 2 shows this step of probability assignment for 8 sub-bands. The probabilities P(Bi) derived from the estimator out- puts, where i is the index of the band, may be interpreted themselves as a coarse measurement of the PSD: The higher the probability value of any band, the higher is its power-spectral content. For m decomposition stages, 2m bands are resulted. Each band covers (0.5/2m) Hz of the RRI spectrum (0-0.5) Hz. 5. Pattern Recognition Technique The proposed pattern recognition technique is imple- mented using the following steps: x(n) LPF HPF v(n) u(n) 2 2 a(n) d(n) Figure 1. A single stage wevelet-decomposition 1 2 2 1 2 a( n ) [ x( n ) x( n )]= + + 1 2 2 1 2 d( n ) [ x( n ) x( n )]= − + 1 2 2 2 0 N B ( a( n )) ( d( n )) n − = −∑ = 1 2 0 N sgn( B ) sgn | a( n ) | | d( n ) | n − = −∑ = Figure 2. Probability assignments of 8 sub-bands 43 The Journal of Engineering Research Vol. 6, No. 2 (2009) 40-46 1. The probability measures are computed for each data of the trial set up to stage m to obtain 2m probabilities P(Bi), where i is the index of the band. The wavelet filters used are of the type Daubechies (db4). The selection of the wavelet filter is a matter of compro- mise between the complexity and good filter perform- ances. 2. The logarithmic based-2 values of the reciprocal of the probabilities can be obtained as: (4) A staircase approximation of the values of I(Bi) can be plotted for all bands. 3. An average plot is found by averaging all I(Bi) values of corresponding bands for all 12 CHF data and nor- mal data in the trial set to obtain two standard plots for CHF and normal as in Fig. 3 for 32 sub-bands. A clear increase in the values of I(Bi) (clear reduction in the probabilities P(Bi) ) in the LF region can be noticed in CHF plot compared to that of normal plot. 4. The I(Bi) values for each data under test are found for 32 sub-bands. A classification factor CF is then deter- mined as: (5) Where IH (Bi) and IN (Bi) are the standard patterns for CHF and normal respectively. The RRI spectrum in the frequency region (0.0625 Hz - 0.4375 Hz) corresponding to bands 5 to 28 is considered for determining the patterns. The LF and HF regions are used in the recognition, while the VLF and the upper high frequency regions are excluded. Depending on the sign of CF, the algorithm can decide whether the data is CHF or normal. The data with a negative CF is considered as nor- mal and the data with a positive CF is a CHF data. 6. Performance Evaluation The performance of a classifier is evaluated by three main metrics: sensitivity, specificity, and accuracy as defined below (Rangayyan, 2001): Sensitivity (%) = TP*100/(TP+FN) (6) Specificity (%) = TN *100/(TN +FP) (7) Accuracy (%) = (TN +TP)*100/T (8) where, TP, TN, FN, FP are defined as in Table 1, and T is the total number of data under test. Sensitivity represents the ability of a classifier to detect the positive cases, eg. CHF. Specificity indicates the abil- ity of a classifier to detect negative cases, eg. normal sub- jects. Accuracy represents the overall performance of a classifier. It indicates the percentage of correctly classi- fied positives and negative cases from the total cases (Rangayyan, 2001). Index of Band (i) 0 A 10 15 20 25 30 35 0 1 2 3 4 5 6 7 8 9 I(B i) Normal CHF Figure 3. Standard patterns of CHF and normal with 32 Sub-bands 1 2I ( B ) log ( )i P( B )i = ( ) ( )( ) ( ) ( )( ) 228 282 5 5 CF I B I B I B I Bi H i i N i i i = − − −∑ ∑ = = Table 1. The confusion matrix 44 The Journal of Engineering Research Vol. 6, No. 2 (2009) 40-46 7. Results and Discussions Figures 4 and 5 show the values of CF for all CHF and normal data using 32 Sub-bands for trial and test data respectively. The algorithm classifies CHF and normal trial data with an accuracy of 100% (Fig. 4) and classifies the test data with sensitivity of (13/17) and with a speci- ficity of (50/53), with an overall identification accuracy of (63/70) 90% (Fig. 5). 8. Consistency of Results The leave-one-out method (Marques 2002) is used to test the consistency of our results. In this method, which is also called partition method, the available set of data (94 records, 24 trial records (12 CHF, 12 normal) and 70 test records (17 CHF, 53 normal)) is divided into k=n sub- sets which rotate in their use of design (trial) and test. The classifier is designed with n-1 subsets and tested on the remaining subset as follows: 1. Divide the available data into k=n=5 subsets (Groups: G1, G2, G3, G4, G5) of randomly patterns. Each sub- set contains 19 records (6 CHF and 13 normal). G3 contains one CHF data less than other groups. 2. Design the classifier using the first n-1=4 subsets and 0 5 10 15 -15 -10 -5 0 5 10 15 Index of data CF CHF Normal Figure 4. Classification factor results of trial data with 32 sub-bands 0 10 20 30 40 50 60 -25 -20 -15 -10 -5 0 5 10 15 CF CHF Normal Index of data Figure 5. Classification factor results of test data with 32 Sub-bands 45 The Journal of Engineering Research Vol. 6, No. 2 (2009) 40-46 test it on the remaining subset. Estimate the perform- ance of the classifier. 3. Repeat the previous step rotating the position of the test subset, obtaining there by k=5 estimates. Tables 2 shows the results of applying this method on 32 sub-bands. This proves that the results obtained on test data are very consistent. This method was also tested with 16-subband wavelet decomposition. The accuracy result was almost the same with the 32 subband decomposition. For further investiga- tion of the screening capability, three different distance metrics were used: * NMSD: Negative Mean Square Distance, which meas- ures the mean square distance of the negative outputs (normal) from the threshold (zero) and is defined as: (9) * PMSD: Positive Mean Square Distance, which measures the mean square distance of the positive outputs (CHF) from the threshold (zero) and is defined as: (10) * TMSD: Total Mean Square Distance, which measures the mean square distance of all outputs (normal and CHF) from the threshold (zero) and is defined as: (11) Where, CF is the classification factor for 53 (normal) and 17 (CHF) test data. Those distance metrics are listed in Table 3. From this table, it can be concluded that increasing the sub-bands from 16 to 32, improves the results (the distances of classified points from threshold). The pattern recognition method results in identical accuracy values when it is used with 16 or 32 sub-bands, while the consistency of the method is improved as the number of sub-bands increases from 16 to 32. On the other hand the omplexity of the algorithm is increased as the number of sub-bands increases from 16 to 32 sub- bands. It has been found that the execution time of the complete program of obtaining the pattern of the signal under test and compare it with both standard patterns and finding the result of identification of the subject, is almost doubled when the number of sub-bands increases from 16 to 32, while the time needed only for comparison with both standard patterns and finding the results of identifi- cation is almost the same as the number of sub-bands increases from 16 to 32. 9. Comparison with Other Techniques An Accuracy of 100% is achieved in (Teich et al. 2001) on the same trial data at window scales of 16-samples and 32-samples using some individual wavelets parameters. The test data in our study was also investigated in a time- domain technique in (Asyali, 2003), which resulted in an accuracy of 93.2%. Both techniques in (Teich et al. 2001) and (Asyali, 2003) can not be considered as screening techniques because neither of them was implemented on both trial and test data. So their results can be considered as data-dependent. The same drawback is applied on our work in (Hossen and Al-Ghunaimi 2004), since the 100% efficiency was obtained only on trial data. On the opposite, our proposed technique in this paper is implemented on both trial data and test data with an accuracy of 100% and 90% respectively and with a total accuracy of 92.55%. Two average standard plots are obtained from trial data and then have been used to iden- tify the test data. 10. Conclusions A soft decision algorithm of PSD estimation of decom- posed wavelet sub-bands is implemented to obtain two standard plots of CHF and normal subjects for the purpose of classification between them in a new pattern recogni- tion method. The base-2 logarithmic values of the recipro- cal of the estimated PSD (with db4 wavelets) have been used to compute the classification factor. The power spec- tral density of the HRV of CHF patients is noticed to be reduced compared to that of normal subject especially at LF region. The frequency range (0.0625 to 0.4375) Hz of the PSD spectrum of the RRI data is used in the recognition. An accuracy of 100% is obtained on 24 MIT trail data and 90% on 70 MIT test data. References Asyali., M.H., 2003, " Discrimination Power of Long- Term Heart Rate Variability Measures," Proceedings of the 25th Annual International Conference of the IEEE, Vol. 1, pp. 200-203. Congestive Heart Failure, University of Maryland Medicine, available at: http://www.umm.edu/patient- ed/doc13diagnos.html. Test Group Sensitivity Specificity Accuracy G1 5/6 10/13 15/19 G2 4/6 11/13 15/19 G3 5/5 11/13 16/18 G4 5/6 13/13 18/19 G5 4/6 11/13 15/19 Table 2. Results of leave-one-out-method using 32 sub-bands 53 2 53 1 i NMSD CF / i = = ∑ = 17 2 17 1 i PMSD CF / i = = ∑ = 70 2 70 1 i TMSD CF / i = = ∑ = Sub-bands NMSD PMSD TMSD 16 14.49 15.17 12.38 32 108.85 121.91 58.65 Table 3. Distance metrics 46 The Journal of Engineering Research Vol. 6, No. 2 (2009) 40-46 Heart Failure Society of America, available at: (http://www.abouthf.org/questions_what_is_hf.htm). Hossen A., Al-Ghunaimi B. and Hassan, M.O., 2005, "Subband Decomposition Soft Decision Algorithm for Heart Rate Variability Analysis in Patients with OSA and Normal Controls," Signal Processing, Vol. 85, pp. 95-106. Hossen, A. and Al-Ghunaimi, B., 2004, " A New Method for Ccreening of Patients with Congestive Heart Failure," International Journal of Computational Intelligence, Vol. 1(3), pp. 266-270. Hossen, A., 2004, " Power Spectral Density Estimation via Wavelet Decomposition," Electronic Letters, Vol. 40(17), pp. 1055-1056. Hossen, A., Al-Ghunaimi, B. and Hassan, M. O., 2003. "A New Simple Algorithm for Heart Rate Variability Analysis in Patients with Obstructive Sleep Apnea and Normal Controls," International Journal of Bioelectromagnetism, Vol. 5(1), pp. 238-239. Hossen, A., Heute, U., 1993, "Fully Adaptive Evaluation of SB-DFT," Proceedings of IEEE Int. Symp. on Circuits and Systems, Chicago, Illinois, USA. Marques, de Sa, J.P., 2002, "Pattern Recognition, Concepts, Methods and Applications," Springer Verlag. Physionet, an NIH-NCRR Research Resource, WFDB Software package, available at: http://www.phys- ionet.org/physiotools/wfdb.shtml. Physionet, Congestive Heart Failure RR Interval Database physionet.org/physiobank/database/chf2db/. Physionet, Normal Sinus Rhythm RR Interval Database. physionet.org/physiobank/database/nsr2db/. Physionet, The BIDMC Congestive Heart Failure Database, physionet.org/physiobank/database/chfdb/. Physionet, The MIT-BIH Normal Sinus Rhythm Database, physionet.org/physiobank/database/nsrdb/. Ponikowski, P., et al. 1997, “Depressed Heart Rate Variability as an Independent Predictor of Death in Chronic Congestive Heart Failure Secondary to Ischemic or Idiopathic Dilated Cardiomyopathy,” Am J Cardiol, Vol. 79, pp. 1645-50. Rangayyan, R. M., 2001, "Biomedical Signal Analysis: A Case-Study Approach," IEEE Press, pp. 466-472. Robert Nolan and Heart Rate Variability, available at bio- medical.com/news_display.cfm?newsid=27. Task Force of the European Society of Cardiology and the North American Society of pacing and Electrophysiology, 1996 Heart Rate Variability, stan- dards of measurements, physiological interpretation, and clinical use, Circulation 93, pp. 1043-1065. Teich M.C., Lowen S.B., Jost B.M., Vibe-Rheymer K., and Heneghan C., 2001. Heart Rate Variability: Measures and Models, Nonlinear Biomedical Signal Processing,” Dynamic Analysis and Modeling, IEEE Press, New York, Chapter 6, Vo. II, pp. 159-213.