Patent Publication Number: US-2005124901-A1

Title: Method and apparatus for electrophysiological and hemodynamic real-time assessment of cardiovascular fitness of a user

Description:
CROSS REFERENCE TO RELATED APPLICATIONS  
      This patent application is based upon Provisional Patent Application Ser. No. 60/527,179 filed on Dec. 05, 2003 
    
    
     BACKGROUND OF THE INVENTION  
      This invention relates to the field of heart rate monitors, and in particular to heart rate monitors used for electrophysiological and hemodynamic assessment of the cardiovascular fitness of athletes, sport enthusiasts and health conscience people.  
      Cardiovascular fitness refers to the quantity of work that can be performed by the muscles, which is critically dependent on the volume of blood that can be delivered by the heart.  
      A large number of products are available commercially from manufacturers like Polar, Timex, Acumen, Cardiosport, Performance, Sensor Dynamics, Sports Instruments and Vetta, such as the S810i by Polar, for the evaluation of cardiovascular fitness based on the user&#39;s heart rate. These devices, which are called heart rate monitors (HRM), detect the electrical activity of the heart through electrocardiogram type electrodes mounted in a chest strap. The filtered and amplified electrical signals are then transmitted to a wristwatch-type device or a device mounted on exercise equipment for further processing and display of calculated heart rate and other relevant information.  
      Use of heart rate monitors for the evaluation of cardiovascular fitness is based on the assumption that the volume of blood, and hence the volume of oxygen delivered to the muscles is directly proportional to the heart rate. For example, a higher heart rate results in more blood being supplied to the muscles. Unfortunately, this assumption is only correct for the ideal situation and it&#39;s not true in many occasions such as endurance training, maximum workloads, cardiovascular abnormalities, extreme or unusual environments, effect of caffeine, and people taking medications.  
      While it is useful to know your heart rate when exercising, it is not the only parameter that should be monitored for assessment of cardiovascular fitness. Heart rate does not necessarily give an accurate picture of cardiovascular fitness since the main factor in such a determination, how much blood the heart is actually able to supply to your body per minute (cardiac output), is not taken into account. Cardiac Output is calculated using the following equation: 
 
 CO=SV*HR    (1) 
 
 Where: 
          CO—Cardiac Output (liter/min)     SV—Stroke Volume, the volume of blood ejected from the heart due to contraction of the left ventricular,     HR—Heart Rate, number of heart beats per minute.        

      Obtaining the actual stroke volume in a clinical setting is complicated. Consequently, cardiac output is traditionally derived in a clinical setting use any of the following four methods.  
      Fick Method  
      Cardiac Output=[oxygen absorbed by the lungs (ml/min)]/[arteriovenous oxygen difference (ml/liter of blood)].  
      Oxygen consumption is derived by measuring the expired gas volume over a known period of time. The arteriovenous oxygen difference is obtained by taking blood samples or by examining the oxygen context of the user&#39;s expired and inspired gas. These methods are difficult to implement. For example, unless the patient has an endotracheal tube, measurements may be faulty because of leakage around the facemask or mouthpiece. Those of ordinary skill in the art will fully appreciate cardiac output measurement by the Fick method. Further details may be found in  The Textbook of Medical Physiology,  7 th  Ed., p. 284 by Arthur C. Guyton, which is hereby incorporated by reference.  
      Thermodilution Method  
      This technique for measuring the cardiac output requires the monitoring of temperature changes after bolus injection of a cold liquid into the blood stream. The injection is made via a catheter that contains a thermistor mounted at its tip. The thermistor measures the sequential changes in temperature which are then plotted over time. The cardiac output is inversely related to the area under the thermodilution curve. This is the standard method for monitoring cardiac output in an intensive care unit. Those of ordinary skill in the art will fully appreciate this method. Similar methods include the Indicator Dilution Method as described in  The Textbook of Medical Physiology,  7 th  Ed., p. 284-285 by Arthur C. Guyton, which is hereby incorporated by reference.  
      Echocardiography  
      This method can be used to derive cardiac output from the measurement of blood flow velocity by recording the Doppler shift of ultrasound reflected from blood cells as they pass through a vessel. The time/flow integral, which is the integral of instantaneous blood flow velocities during one cardiac cycle, is obtained for the blood flow in the left ventricular outflow tract (other sites can be used). This is multiplied by the cross-sectional area of the tract and the heart rate to give cardiac output. The main disadvantages of this method are that a skilled operator is needed, the probe is large and therefore heavy sedation or anesthesia is needed, the equipment is very expensive and the probe cannot be fixed so as to give continuous cardiac output readings without an expert user being present. Those of ordinary skill in the art will fully appreciate this method.  The Textbook of Medical Physiology,  7 th  Ed., p. 284 by Arthur C. Guyton, which is hereby incorporated by reference.  
      Thoracic Bioimpedance Technology  
      This method has the advantages of providing continuous cardiac output measurement at limited risk to the patient. A small, high frequency current is passed through the thorax via electrodes placed on the skin. Contraction of the heart produces a cyclical change in thoracic impedance. Sensing electrodes are used to measure the changes in impedance within the thorax. A constant current generator establishes a fixed level for I(o). The resulting voltage change, V(t), is used to calculate impedance. Because the impedance is assumed to be purely resistive, the total impedance, Z, is calculated by Ohm&#39;s Law. The normal impedance value for an adult is 20-48 ohms with a current frequency of 50-100 Hz. The total impedance is derived from a constant base impedance, Z o , and time-varying impedance, Z(t), as shown in equation (2) below and  FIG. 1 :  
             Z   =         V   ⁡     (   t   )         I   ⁡     (   o   )         =       Z   o     +     Z   ⁡     (   t   )                   (   2   )             
 
      Z o    140  reflects constant resistivity of tissue and bones. Time-varying impedance Z(t)  130  reflects changes in resistivity of portions of the arterial system as blood flows through the aorta.  
      The aforementioned impedance values may then be used to calculate stroke volume using the Kubicek equation (3) or any of its modifications like, for example, the Gundarov equation (4):  
             SV   =     ρ   ⁢           ⁢     (       L   2       Z   o       )     ⁢           ⁢   LVET   ⁢           ⁢       ⅆ     Z   ⁡     (   t   )           ⅆ     t   max                   (   3   )               SV   =     0.9   ⁢           ⁢   ρ   ⁢           ⁢   κ   ⁢           ⁢         Q   2     ⁢   L       Z   o   2       ⁢   LVET   ⁢           ⁢       ⅆ     Z   ⁡     (   t   )           ⅆ     t   max                   (   4   )             
 
 Where: 
          L=distance between sensing electrodes (cm),     LVET=left ventricular ejection time,     Z o =base impedance in (ohms),  
                 ⅆ     Z   ⁡     (   t   )           ⅆ     t   max         =       ⁢     magnitude   ⁢           ⁢   of   ⁢           ⁢   the   ⁢           ⁢   largest   ⁢           ⁢   negative   ⁢           ⁢   derivative   ⁢           ⁢   of                     ⁢       the   ⁢           ⁢   impedance   ⁢           ⁢   change     ,               
    Z(t) occurring during systole (ohms/s)     ρ=resistivity of blood (ohms/cm),     δ=weight correction factor,     H=height of the user,     κ=torso correction factor     Q=circumference of torso.        

      LVET is the time between the opening and closing of the aortic valve. These times can be identified using two characteristic points of the ICG. Point B corresponds with the aortic valve opening. Point X corresponds with the aortic valve closing ( FIG. 4 ).  
         ⅆ     Z   ⁡     (   t   )           ⅆ     t   max           
 
 corresponds to characteristic point C of the ICG ( FIG. 4 ), which is point of maximum amplitude in the cardiac cycle. 
 
      As indicated in equation 1, once stroke volume is known, it may be multiplied by the heart rate to determine cardiac output and consequently, cardiovascular fitness.  
      In contrast to the Fick method and the other cardiac output methods previously described, only the thoracic bioimpedance technology, also called impedance cardiography, is practically useful for real-time assessment of cardiovascular fitness during workload. Still, there are drawbacks associated with this technology: 
          Unfortunately, the cyclical change in impedance cardiography is about 0.5% to 1%. This yields a low signal to noise ratio. Consequently, the impedance signal is vulnerable to noise and artifacts.     The user should lie still in order to limit noise and artifacts associated with muscle contractions.     The impedance cardiography procedure requires positioning of at least 4 pairs of electrodes in the neck and thoracic areas. These electrodes must then be connected to a signal processing unit, thereby restricting the user&#39;s freedom of movement.     The calculation of stroke volume, and consequently cardiac output, using impedance cardiography requires ECG derived data such as RR interval and Q wave onset for detection of left ventricular ejection time and maximum impedance. ECG-based technologies incorporated in current heart rate monitors only detect electrical spikes associated with the heartbeat and not values such as RR interval and Q wave onset therefore, such monitors cannot be used as an input for impedance cardiography.     Current signal processing methods used in impedance cardiography and electrocardiography employ complex calculations and filtering techniques that necessarily impose significant hardware and software requirements. This restricts the ability to design a lightweight, portable, low cost monitoring device.        

      The present invention overcomes the aforementioned problems by offering an effective method and apparatus for simultaneous real-time electrophysiological and hemodynamic evaluation of cardiovascular fitness.  
    
    
     BRIEF DESCRIPTION OF THE DRAWINGS  
       FIG. 1  illustrates basic concepts of thoracic impedance cardiography.  
       FIG. 2  shows a heart rate monitor with hemodynamic and electrophysiological assessment features.  
       FIG. 3  is a flow chart for calculation of heart rate, heart rate variability, stroke volume and cardiac output.  
       FIG. 4  illustrates characteristic points on ECG and ICG waveforms.  
       FIGS. 5   a,    5   b,    5   c  and  5   d  show detection and refining of point R.  
       FIG. 6  illustrates point and intervals of two successive RR intervals used in the calculation of noise level and point T.  
       FIGS. 7   a  and  7   b  illustrate saw-type and spike-type noise.  
       FIGS. 7   c  and  7   d  show ECG fragment before and after cubic spline smoothing.  
       FIGS. 8   a,    8   b  and  8   c  illustrate detection of point Q.  
       FIGS. 9   a,    9   b  and  9   c  illustrate detection of point S.  
       FIGS. 10   a,    10   b,    10   c,    10   d,    10   e,    10   f  and  10   g  illustrate detection of points T and T e .  
       FIG. 11  is a flow chart for obtaining (i) point C on an ICG and (ii)  
           ⅆ     Z   ⁡     (   t   )           ⅆ     t   max         .         
       FIG. 12  is a flow chart describing an iterative algorithm for filtering and smoothing ICG signals.  
    
    
     DETAILED DESCRIPTION OF THE DRAWINGS  
      In the following description, numerous specific details are set forth to provide a thorough understanding of the present invention. However, it will be obvious to those skilled in the art that the present invention may be practiced without such specific details. For the most part, details concerning specific non-essential materials and the like have been omitted inasmuch as such details are not necessary to obtain a complete understanding of the present invention and are within the skills of persons of ordinary skill in the relevant art.  
      The fundamental aspect of the invention is a real-time simultaneous measurement of electrophysiological and hemodynamic parameters incorporated in a device, similar to the popular heart rate monitor, which can be used for comprehensive assessment of cardiovascular fitness.  
      These objectives are reached by using ECG signals, rather than ICG signals, for the measurement of key parameters needed for the calculation of electrophysiological and hemodynamic data, thus mitigating the influence of noise and artifacts on ICG signal.  
      As shown in  FIG. 1 , AC current source electrodes  120  and sensing electrodes  110  are positioned across the torso substantially parallel to subclavian arteries  160 . This configuration is in contrast to other thoracic bioimpedance set-ups whereby electrodes are placed along the torso in parallel to carotid  150  and renal arteries  170 . According to equation (5), this positioning of electrodes provides lower resistivity Z o  because, for the average person, the width of one&#39;s torso is less than the length of one&#39;s torso. therefore, the distance between electrodes (L) will be less  
             R   =       ρ   ⁢           ⁢   L     A             (   5   )             
 
 Where: 
          R=resistivity of torso     L=the distance between sensing electrodes placed at opposite sides of the torso     A=cross section area of the torso (cm 2 )     ρ=resistivity of blood (ohms/cm)        

      Z o  is considered equal to R because thoracic impedance technology is based on the assumption that the total thoracic impedance is totally resistive (i.e., the reactance component is equal zero). Consequently, Z o  will be reduced providing better sensitivity to the time-varying component Z(t) of equation (2).  
      The heart rate monitor with hemodynamic and electrophysiological assessment capabilities is shown in  FIG. 2  wherein a belt  200 , with electrodes  210 ,  220 ,  230 , and  240  and signal processing unit (“SPU”)  250  are shown. Electrodes  210  and  240  provide a path for constant alternating current (AC) application and electrodes  220  and  230  provide for voltage sensing. The electrodes  220  and  230  sense both the alternating current from electrodes  210  and  240  as well as the user&#39;s own ECG. Both of said signals are conveyed to the SPU  250 , where they are filtered, amplified, digitized and transmitted to a wristwatch like device  260  or any other kind of processing unit for calculation and visualization of electrophysiological and hemodynamic parameters. Commonly used hardware filters, utilizing the difference of predefined frequencies for the ECG and ICG signals, work to separate the ECG and ICG signals. Electrodes  210 ,  220 ,  230 , and  240  are placed in a flexible, stretchable belt that provides placement of the electrodes at specified locations regardless of the torso size.  
       FIG. 3  illustrates how heart rate (HR), heart rate variability (HRV), stroke volume (SV), and cardiac output (CO) are calculated. After activating the heart rate monitor  310  and acquiring ECG  320  and ICG  330  signals, the characteristic points (a.k.a. fiduciary points) Q, R, S, J, T and T e  of ECG and fiduciary points B and X of ICG are ascertained in step  340  ( FIGS. 3 and 4 ). Point C is ascertained on ICG in step  350 . Point B corresponds to the time of the aortic valve opening while point X corresponds to the time of aortic valve closing. Point C corresponds to the maximum of  
           ⅆ     Z   ⁡     (   t   )           ⅆ     t   max         .         
 The time interval between point B and point X defines LVET. Because point B and X are hard to identify due to the low signal to noise ratio, the present method correlates point B with the time at which point J occurs. Furthermore, point X is correlated with the points in time when T e  occurs. Point C is much less vulnerable to noise and artifacts due to its distinct location. Consequently, point C may be measured directly on the ICG instead of being derived from points on the ECG. 
 
      The detection of characteristic points Q, R, S, J, T and T e  is illustrated in  FIG. 4  to  FIG. 10   a - 10   e.    
      The detection of characteristic points starts from extraction of point R as the most distinctive point of ECG.  
      When progressing along axis t ( FIG. 4 ) and comparing amplitude, V i , of a point at current time, t i , and amplitudes, V 1  at time t i −d 1  and V 2  at time t i −d 2  ( FIG. 5   a - 5   d ), the approximate location of point R a  is found, when the following is true: 
 
( V   i   −V   1 )&gt; A   1  OR ( V   i   −V   2 )&gt; A   2  
 
 Where: 
          V i =amplitude of the current point at time t i ;     V 1 =amplitude at time t i −d 1 ;     V 2 =amplitude at time t i −d 2 ;        

      A 1 =0.25 mV and d 1 =75 ms may be used for strongly expressed (high) R waves ( 5   b ).  
      A 2 =0.15 mV and d 2 =40 ms may be applied for weakly expressed (short) R waves ( 5   a ).  
      Theses values are commonly selected by those of skill in the art because the amplitude of point R normally increases 0.25 mV within a period of 75 ms for high R waves and it increases 0.15 mV within a period of 40 ms for short R waves. However, persons of ordinary skill in the art realize certain physiological conditions may require the alteration of values A 1 , A 2 , d 1  and d 2 .  
      After a point R a  is found, the location R t  of point R is ascertained analyzing a time interval [t i , t 1 +d r ] ( FIG. 5   d ). As soon as the amplitude of point R is decreases by more than A r , the location of point R is considered found and further analysis of interval [t i , t i +d r ] is stopped. The point R r  is one step back from the point of decrease. The empiric values of d r  and A r  are 200 ms and 0.05 mV respectively, because the maximum width of R wave is not more than 200 ms and, within this period, R wave descends by not less than 0.05 mV. However, persons of ordinary skill in the art may successfully use other values.  
      RR interval is the time between two successive R points ( FIG. 6 ). QRS fragment is defined as two successive RR intervals (RR) i−1  and (RR) i  ( FIG. 6 ). Each current (RR) i  interval is tested for the level of noise.  
      First the noise level, N 1 , of saw-type noise ( FIG. 7   a ) is calculated within time interval [R i−1 +e 1 , R i−e   1 ], ( FIG. 6 ), where e 1  may typically be 75 ms.  
      Starting N 1 =0, for each point j of interval [R 1−1 +e 1 , R i −e 1 ] and each m, if |V j −V j−1 |&gt;2 m  AND |V j −V j+1 |&gt;2 m , then N 1 =N 1 +2 m , where m=3,2,1,0. At the next step the level N 2  of spike-type noise ( FIG. 7   b ) is calculated within time interval [R i−1 +e 1 , R i −e 1 ], ( FIG. 6   b ), where e 1  may typically be 115 ms.  
      Starting N 2 =0, for each point j of interval [R i−1 +e 1 , R i −e 1 ] and each m, if |V j −V j−1 &gt;m AND |V j −V j+1  |&gt;m, then N 1 =N 1 +m, where m=30, 20.  
      The total noise level of current interval (RR) i , N i =N 1 +N 2 . If the noise level N i &gt;N limit  where N limit  may be 20, then current interval (RR) i  is considered unreliable and excluded from further calculations.  
      The values e 1 =75 ms and e 2 =115 ms are empirically derived and commonly selected by those of skill in the art because indentations 75 ms and 115 ms from R point exclude Q, S and T waves from mistakenly considering these points as saw-type noise as well as point R as a spike-type noise. However, persons of ordinary skill in the art may successfully use other values.  
      N limit =20 provides a sufficient noise filtering for disclosed application however, persons of ordinary skill in the art may successfully use other values.  
      After noise filtering of RR interval, RR interval is smoothed using cubic spline interpolation algorithm included in Matlab Version 3.2 spline toolbox. However, persons of ordinary skill in the art may successfully use other smoothing techniques.  FIG. 7   c  illustrates ECG fragment before smoothing, while  FIG. 7   d  shows ECG fragment after cubic spline smoothing was applied.  
      QRS fragment ( FIG. 6 ) is defined as two successive reliable RR intervals, (RR) i−1  and (RR) i .  
      Point Q, S, J, T, and T e  are ascertained within QRS fragment.  
      Point Q is calculated from the graphs in  FIGS. 8   a - 8   c  in the following manner:  
      Referring to  FIGS. 8   a  and recalling that R has already been located as shown above, point Q may be obtained as follows:  
      A time interval to be analyzed is defined as [t DQ , t R ] where, typically, t DQ =t R −75 ms. t R  corresponds with point R as was derived above. A 75 ms period is commonly selected by those of skill in the art because the onset of Q wave is normally within 0 to 75 ms of R. However, persons of ordinary skill in the art may select other value, which may be successfully applied.  
      When sampling this period, starting at t R  and progressing towards t DQ , Q is found when the following is true: 
 
 A   i−1&gt; and (   A   R   −A   i )&gt; A   RQ  
 
 Where: 
          A denotes amplitude     A RQ =0.1 mV     A R =the amplitude at point R        

      A RQ  =0.1 mV is commonly selected by those of skill in the art because a typical R peak is at least 0.1 mV “above” the Q. However, persons of ordinary skill in the art may select other value, which may be successfully applied. This approach is valid for normal Q waves as shown on  FIG. 8   a.    
      Next, if the above conditions are not met and Q is not ascertained, the following conditions are evaluated to determine Q (See  FIG. 8   b ):  
      A time interval to be analyzed is defined as [t DQ , t R ] where, again, typically, t DQ =t R −75 ms. A 75 ms period is commonly selected by those of skill in the art because the onset of Q wave is normally within 0 to 75 ms of R. However, persons of ordinary skill in the art may select other value, which may be successfully applied.  
      When sampling this period, starting at t R  and progressing towards t DQ , Q is found when the following is true: 
 
( A   −A   i−3 )&gt; A   d  and ( A   R   −A   i )&gt;A RQ  
 
 Where: 
          A denotes amplitude     A d =0.025 mV     A RQ =0.1 mV        

      A RQ =0.1 mV is commonly selected by those of skill in the art because a typical R peak is at least 0.1 mV “above” the Q. However, persons of ordinary skill in the art may select other value, which may be successfully applied.  
      A d =0.025 mV is commonly selected by those of skill in the art because this amplitude difference is typical for abnormal Q wave shown on  FIG. 8   b.  However, persons of ordinary skill in the art may select other value, which may be successfully applied.  
      Next, if the above conditions are not met and Q is not ascertained, the following conditions are evaluated to determine Q (See  FIG. 8   c ):  
      A time interval to be analyzed is defined as [t DQ , t R ] where, again, typically, t DQ =t R −75 ms. A 75 ms period is commonly selected by those of skill in the art because the onset of Q wave is normally within 0 to 75 ms of R. However, persons of ordinary skill in the art may select other value, which may be successfully applied.  
      When sampling this period, starting at t R  and progressing towards t DQ , Q is found when the following is true:  
             A   i     -     A     i   -   3             A     i   +   3       -     A   i         ≥       Q   r     ⁢           ⁢   and   ⁢           ⁢     (       A   R     -     A   i       )       &gt;     A   RQ         
 
 Where: 
          A denotes amplitude     A RQ =0.1 mV     Q r =0.45        

      A RQ =0.1 mV is commonly selected by those of skill in the art because a typical R peak is at least 0.1 mV “above” the Q. However, persons of ordinary skill in the art may select other value, which may be successfully applied.  
      Q r =0.45 is commonly selected by those of skill in the art because it&#39;s a typical ratio for abnormal Q wave related to group of premature beats ( FIG. 8   c ). However, persons of ordinary skill in the art may select other value, which may be successfully applied.  
      Referring to  FIG. 9   a  and recalling that R has already been located as shown above, point S may be obtained as follows:  
      A time interval to be analyzed is defined as [t R , t DS ] where, typically, t DS =t R +75 ms. A 75 ms period is commonly selected by those of skill in the art because the onset of S wave is normally within 0 to 75 ms of R. However, persons of ordinary skill in the art may select other value, which may be successfully applied.  
      When sampling this period, starting at t R  and progressing towards t DS , S is found when the following is true: 
 
 A   1+1   &gt;A   i  and ( A   R   −A   i )&gt;A RS  
 
 Where: 
          A denotes amplitude     A RS =0.1 mV        

      A RS =0.1 mV is commonly selected by those of skill in the art because a typical R peak is at least 0.1 mV “above” the S. However, persons of ordinary skill in the art may select other value, which may be successfully applied. This approach is valid for normal S waves as shown on  FIG. 9   a.    
      Next, if the above conditions are not met and S is not ascertained, the following conditions are evaluated to determine S (See  FIG. 9   b ):  
      A time interval to be analyzed is defined as [t R , t DS ] where, again, typically, t DS =t R +75 ms. A 75 ms period is commonly selected by those of skill in the art because the onset of S wave is normally within 0 to 75 ms of R. However, persons of ordinary skill in the art may select other value, which may be successfully applied.  
      When sampling this period, starting at t R  and progressing towards t DS , S is found when the following is true: 
 
( A   i   −A   i+3 )&lt; A   d  and ( A   R   −A   i )&gt; A   RS  
 
 Where: 
          A denotes amplitude     A d =0.025 mV     A RS =0.1 mV        

      A RS =0.1 mV is commonly selected by those of skill in the art because a typical R peak is at least 0.1 mV “above” the S. However, persons of ordinary skill in the art may select other value, which may be successfully applied.  
      A d =0.025 mV is commonly selected by those of skill in the art because a this amplitude difference is typical for abnormal S wave shown on  FIG. 9   b.  However, persons of ordinary skill in the art may select other value, which may be successfully applied.  
      Next, if the above conditions are not met and S is not ascertained, the following conditions are evaluated to determine S (See  FIG. 9   c ):  
      A time interval to be analyzed is defined as [t R , t DS ] where, again, typically, t DS =t R +75 ms. A 75 ms period is commonly selected by those of skill in the art because the onset of S wave is normally within 0 to 75 ms of R. However, persons of ordinary skill in the art may select other value, which may be successfully applied.  
      When sampling this period, starting at t R  and progressing towards t DS , S is found when the following is true:  
             A   i     -     A     i   +   3             A     i   -   3       -     A   i         ≥       S   r     ⁢           ⁢   and   ⁢           ⁢     (       A   R     -     A   i       )       &gt;     A   RS         
 
 Where: 
          A denotes amplitude     A RS =0.1 mV     S r =0.3        

      A RS =0.1 mV is commonly selected by those of skill in the art because a typical R peak is at least 0.1 mV “above” the S. However, persons of ordinary skill in the art may select other value, which may be successfully applied.  
      S r =0.3 is commonly selected by those of skill in the art because it&#39;s a typical ratio for abnormal S wave related to group of premature beats ( FIG. 9   c ). However, persons of ordinary skill in the art may select other value, which may be successfully applied.  
      Recalling that point S has already been located as shown above, point J ( FIG. 4 ) may be obtained as follows:  
      A time interval to be analyzed is defined as [t S , t DJ ] where, typically, t DJ =t S +75 ms. A 75 ms period is commonly selected by those of skill in the art because the point J is normally within 0 to 75 ms of S. However, persons of ordinary skill in the art may select other value, which may be successfully applied.  
      When sampling this period, starting at t R  and progressing towards t DJ , J is found when the following is true: 
 
( A   i+3   −A   i )&lt;A d  
 
 Where: 
          A denotes amplitude     A d =0.025 mV        

      If point J is not ascertained, point J may be considered to have time coordinate equal t Q +50 ms. t Q +50 ms is selected because the point J normally coincides with the onset of Pre-ejection Period, which normally starts within 50 ms after point Q.  
      Recalling that point J ( FIG. 4 ) has been located and two successive RR intervals, (RR) i−1  and (RR) i  ( FIG. 6 ) have been identified as shown above, point T ( FIG. 4 ,  FIGS. 10   a - 10   e ) may be obtained as follows:  
      A time interval to be analyzed is defined as [t J , t N ], where t N =60%*[R i−1 , R i ] ( FIG. 6 ). 60% is commonly selected by those of skill in the art, because the point T is normally located within interval [t J , t N ]. However, persons of ordinary skill in the art may select other value, which may be successfully applied.  
      When sampling this period, starting at t J  and progressing towards t N , point T is found if the distance from moving point (t i , A i ) to straight line (J i , J i−1 ), which exists between point J i  of interval [R i−2 , R i ] and point J i−1  of interval [R i−1 , R i−1 ] ( FIG. 6 ), is more than T Amax . Where T Amax  is a maximum point and=2 mm. Those of skill of the art commonly select this value, however, persons of ordinary skill in the art may select other value, which may be successfully applied. This approach is valid for normal T wave as it is shown on  FIG. 10   a.    
      If T wave is inverted ( FIG. 10   b ), then, when sampling this period, starting at t J  and progressing towards t N , point T is found if the distance from moving point (t i , A i ) to straight line, (J i , J i−1 ), which exists between point J i  of interval [R i−2 , R i ] and point J i−1  of interval [R i−1 , R i−1 ] ( FIG. 6 ), is more than T Amin . Where T Amin  is a minimum point and=8 mm. This value is commonly selected by those of skill of the art, however, persons of ordinary skill in the art may select other value, which may be successfully applied. This approach is valid for a normal T wave as it is shown on  FIG. 10   b.    
      For flat T waves ( FIG. 10   c ), sampling interval [t j , t N ] as defined above and progressing toward t N , point T is found when the following is true: 
 
( A   i   −A   i+5 )&gt; A   d  
 
 Where: 
          A denotes amplitude     A d =0.025 mV        

      A d =0.025 mV is commonly selected by those of skill in the art because experimental data well correlated with this value. However, persons of ordinary skill in the art may select other value, which may be successfully applied.  
      If point T is not identified, the point is considered undetectable. Time coordinates of point T and point T e  are then considered equal to t N .  
      When point T has been located, then point T e  ( FIGS. 10   d ,  10   e ) is ascertained as follows:  
      A time interval to be analyzed starting from t T  and progressing to direction of time gain.  
      Point T e  is found when one of the following is true: 
          A i &gt;A i−1 , if T wave has shape as shown on  FIG. 10   d       or     A i &lt;A i−1 , if T wave has shape as shown on  FIG. 10   f       or     angle α i , between straight line (T, T i ) and axis t, becomes less than α i−1  and this condition remains for the period of time not less then d=40 ms for T waves shaped as shown on  FIG. 10   e  and  FIG. 10   g.          

      d=40 ms is commonly selected by those of skill in the art because experimental data well correlates with this value. However, persons of ordinary skill in the art may select other value, which may be successfully applied.  
      Other methods such as spectral analysis, signal averaging, wavelet transform, and Fourier transform may be successfully used for detection of characteristic points of ECG.  
       FIG. 11  illustrates detection of point C, which defines parameter  
         ⅆ     Z   ⁡     (   t   )           ⅆ     t   max             
 of equations (10). First, noise level q is identified  1110 . The noise level q is defined as the average change of ICG signals between two adjacent points within n time intervals [B,X]/3 using formula (6):  
             q   =         ∑   1     n   -   1       ⁢           ⁢            Y   j     -     Y     j   +   1                  2   ⁢   n               (   6   )               
 Where: 
          q=noise level of ICG signal;     j=1 to n−1, where n is total number of points (samples) in interval [B,X]/3;     Y j =j th  point (sample) of ICG signal.     n=the number of points within time interval [B,X1/3; 
 
 n=F*[B,X]/ 3 
       
      Where F is the sampling frequency (number of samples per second) of digitization performed in the processing unit  260 . Typical sample frequencies are from 200 Hz to 1,000 Hz.  
      Next, to further ensure a signal with a viable point C has been obtained, impedance change, ΔZ  1120  is deemed undetectable if:  
                 Δ   ⁢           ⁢   Z         2   ⁢           ⁢     Z   o       +     Δ   ⁢           ⁢   Z         &lt;     q     2   ⁢     Z   o                 (   7   )             
 
 Where: 
          ΔZ=impedance change;     Z o =base impedance;     q=noise level 
 
Δ Z=k*Z   o  
       

      In the present embodiment, k=0.01, although other values may be successfully used. k is typically in a range of 0.005 to 0.02.  
      In subsequent steps, acceptable ICG signals are qualified using ensemble averaging  1130  and then subjected to smoothing processes  1140 . Point C is then calculated as the peak height  
         ⅆ     Z   ⁡     (   t   )           ⅆ   t         
 
 on the ICG curve  1150  within time interval [B, X], which is equal to interval [J, T e ]. The number of consecutive [B,X] i  intervals n included in the ensemble averaging process  1130  depends on a predefined threshold signal-to-noise ratio (SNR), N, and on a permissible number, M, of consecutive [B,X] intervals used for ensemble averaging. In one embodiment, N=1,000 and M=10, although other values may be successfully used. 
 
       FIG. 12  illustrates iterative algorithm of filtering and smoothing of ICG signals.  
      Each point S j  is calculated as the averaged value of correspondent points of n consecutive [B,X] i  intervals (8):  
               S   j     =         ∑   1   n     ⁢           ⁢     Y     j   i         n             (   8   )             
 
      The expected improvement in signal-to-noise ratio after ensemble averaging process is expressed by {square root}{square root over (n)}.  
      The smoothing of the ICG signal  1140  in interval [B,X] is performed using triangular 5-points smooth method. The signal value, S j , is calculated by equation (9):  
               S   j     =         Y     j   -   2       +     2   ⁢     Y     j   -   1         +     3   ⁢     Y   j       +     2   ⁢     Y     j   +   1         +     Y     j   +   2         9             (   9   )             
 
 for j=3 to n−2, where S j  is the j th  point in the smoothed signal, Y j  is the j th  point in the original signal, and is the total number of samples (points) in the interval [B,X]. The process increases the signal-to-noise ratio by {square root}{square root over (5)}. 
 
      The total improvement of signal-to-noise ratio after ensemble averaging and smoothing is {square root}{square root over (5)}*n.  
      Returning to  FIG. 3 , and aside from calculating point C, Left Ventricular Ejection Time (LVET) 355 must be calculated. Instead of using point B and X of ICG, the method uses point J, the point of junction between the S wave and T wave, and point T e , which is at the end of the T wave of the ECG. Point J corresponds to aortic valve opening and point T e  corresponds aortic valve closing. LVET is equal to the time interval between point J and point T e  of the ECG waveform. Choosing points J and T e  provides more reliable and reproducible results, because the exact time domain location of points J and T e  are much easier to identify than the transient location of points B and X of an ICG. Also, calculation of LVET as the interval between points J and T e , automatically corrects LVET which may be significantly affected by increased or decreased heart rate.  
      Heart rate is calculated  360  ( FIG. 3 ) as a moving average of n consecutive RR intervals, where n is derived from iterative filtering and smoothing process  1130  and  1140  described in  FIG. 12 .  
      Anthropometrical data  370  comprise of torso height and torso perimeter.  
      Calculation of stroke volume  380  may employ any of the equations (3) and (4) or their modifications. In the present invention, stroke volume is calculated using Gundarov modified equation (10).  
             SV   =     0.9   ⁢           ⁢   ρ   ⁢           ⁢   κ   ⁢         H   2     ⁢   L       Z   o   2       ⁢   LVET   ⁢       ⅆ     Z   ⁡     (   t   )           ⅆ     t   max                   (   10   )             
 
 Where: 
          L=distance between sensing electrodes (cm),     LVET=left ventricular ejection time,     Z o =base impedance (ohms),  
             ⅆ     Z   ⁡     (   t   )           ⅆ     t   max         =     magnitude   ⁢           ⁢   of   ⁢           ⁢   the   ⁢           ⁢   largest   ⁢           ⁢   negative   ⁢           ⁢   derivative   ⁢           ⁢   of   ⁢           ⁢   the   ⁢           ⁢   impedance   ⁢           ⁢   change       ,       
    Z(t) occurring during systole (ohms/s),     ρ=resistivity of blood (ohms/cm),     H=height of the user&#39;s torso,     κ=torso correction factor.        

      Torso correction factors (κ) are shown in the table (2):  
                           TABLE 2                                   Torso size in cm   k                                                    80   0.0036           82   0.0035           84   0.0034           86   0.0033           88   0.0032           90   0.0032           92   0.0031           94   0.0030           96   0.0029           98   0.0028           100   0.0027           102   0.0027           104   0.0026           106   0.0025           108   0.0025           110   0.0024                      
 
      Once stroke volume  380  has been calculated, cardiac output  390  is calculated using equation (1). 
 
 CO=SV*HR  
 
 Where: 
          CO=cardiac output     SV=stroke volume     HR=heart rate        

      Although the present invention and its advantages have been described in detail, it should be understood that various changes, substitutions and alterations can be made herein without departing from the spirit and scope of the invention.