POWER SPECTRUM ANALYSIS OF HUMAN
HEART RATE VARIABILITY:
A QUANTITATIVE MEASURE OF CARDIOVASCULAR
FUNCTION AND CONTROL
by
Kenneth Mark Riff, M.D.
B.S., University of Michigan, 1973
M.D., University of California at San Francisco, 1977
A thesis submitted to the
Faculty of the Graduate School of the
University of Colorado in partial fulfillment
of the requirements for the degree of
Master of Science
Department of Electrical Engineering and Computer Science
1988
This thesis for the Master of Science degree by
Kenneth Mark Riff, M.D.
has been approved for the
Department of
Electrical Engineering and Computer Science
by
Date.
Riff, Kenneth Mark (Master of Engineering, Electrical Engineering)
Power Spectrum Analysis of Human Heart Rate Variability: A
Quantitative Measure of Cardiovascular Function and Control
Thesis directed by Professor Jochen Edrich
Analysis of heart rate variability (HRV) can yield
significant information about the activity of fundamental
cardiovascular control mechanisms modulated by the autonomic
nervous system (ANS). HRV is best evaluated by analysis of its
power spectrum, as this method fractionates total HRV into its
separate frequency components and permits discrimination between
the different physiologic sources of HRV.
Various techniques of deriving the power spectrum of HRV are
reviewed. Stochastic point process theory leads to the derivation
of two nonequivalent spectra, the spectrum of counts and the
interval spectrum, and methods of estimating these spectra are
reviewed.
A review of cellular electrophysiology, which emphasizes the
role of integrative control mechanisms, leads to a model of heart
beat generation termed the integral pulse frequency modulator
(IPFM). The properties of the IPFM are derived and examined in
detail, including its frequency response, inputoutput relation
ships, and analysis of its nonlinear behavior. It is shown that a
signal derived from the IPFM's output pulse train, the
instantaneous frequency, retains a linear relationship to the input
modulating signal, while the stochastic spectral estimators contain
significant frequency distortion due to the nonlinear behavior of
iv
the IPFM. Two methods of converting the instantaneous frequency
signal into a stationary continuous signal suitable for spectral
analysis, linear interpolation and lowpass filtering, are derived
and analyzed. The sources of error in the use of the instantaneous
frequency, nonlinearities and aliasing, are examined in detail.
The behavior of an IPFM simulated on a digital computer is
examined. It is shown that the frequency range for valid spectral
analysis is approximately half the Nyquist frequency. An unusual
type of aliasing, termed signal aliasing, is derived and examined.
A complete system, including hardware and software, for
performing clinical HRVSA is described. The role of certain signal
processing algorithms, including digital filtering, autocorrelation
techniques, and time domain lag windowing, are examined.
HRVSA is performed and analyzed in a variety of clinical
situations, including changes in posture, aerobic exercise, and
treatment of abnormal autonomic function. It is shown that HRVSA
provides a powerful noninvasive tool to quantitate ANS activity
and function.
This thesis is dedicated to my wife, Julie, and my daughter, Emily,
whose encouragment, patience, and love made it all possible.
It is also dedicated to Dr. Henry Roth, whose generous financial
support allowed this concept to become a reality.
CONTENTS
CHAPTER
I. INTRODUCTION: THE ANALYSIS OF HEART RATE VARIABILITY......1
1.1 Methods of Analysis of Heart Rate Variability......... .2
1.1.1 Statistical Methods of Analysis..................3
1.1.2 Power Spectrum Analysis.........................5
II. THE PHYSIOLOGY OF HEART RATE VARIABILITY...................9
2.1 Basic Electrophysiology................................9
2.1.1 Membrane Passive Electrical Properties..........11
2.1.2 Membrane Active Electrical Properties..........13
2.1.3 The Role of Spatial and Temporal Integration....14
2.1.4 Pacemaker Cells................................15
2.1.5 Cardiac Electrophysiology......................16
2.2 The Autonomic Nervous System.........................17
2.2.1 ANS Control of the Cardiovascular System.......18
2.2.2 The Origin of Heart Rate Variability...........23
2.2.3 Clinical Implications..........................24
III. METHODS OF DETERMINING THE SPECTRUM OF
HEART RATE VARIABILITY..................................25
3.1 Spectral Analysis Using Stochastic
Point Process Theory................................27
3.1.1 Stationary Point Processes.....................27
3.1.2 Second Order Properties of Intervals...........28
3.1.3 Second Order Properties of Counts..............29
vii
3.1.4 Estimation of the Interval Spectrum.............31
3.1.5 Estimation of the Spectrum of Counts............32
3.1.6 Comparison of the Interval Spectrum
with the Spectrum of Counts....................34
3.2 Spectral Analysis Using Time Series Methods...........34
3.3 Spectral Analysis Using a Communications Model...36
IV. THE INTEGRAL PULSE FREQUENCY MODULATOR....................39
4.1 Description of the Integral Pulse
Frequency Modulator.................................40
4.2 Analysis of InputOutput Relationships................42
4.3 Physiologic Justification.............................46
4.4 System Properties.....................................50
4.5 The Spectrum of a SinusoidallyModulated Pulse Train..55
V. SPECTRAL ESTIMATION UTILIZING
THE INSTANTANEOUS FREQUENCY.............................59
5.1 LowPass Filtering....................................60
5.2 Polynomial Interpolation..............................63
5.3 Sources of Error: Aliasing and Nonlinear Effects.....68
5.3.1 Aliasing........................................68
5.3.2 Nonlinear Effects..............................70
5.4 Spectral Analysis of the Instantaneous Frequency......76
VI. SPECTRAL ANALYSIS OF SIMULATED DATA.......................83
6.1 Simulations of Linear and Nonlinear
Frequency Response Effects..........................85
6.1.1 A Single Modulating Sinusoid....................85
6.1.2 Two Modulating Sinusoids........................87
6.1.3 Three Modulating Sinusoids of Equal Phases...89
6.1.4 Three Sinusoids of Differing Phases.............91
viii
6.1.5 Two Sinusoids in the Second Order
Nonlinear Range..............................93
6.1.6 Two Sinusoids in Higher Order Nonlinear Range..95
6.1.7 Frequency Response Simulation Summary...........97
6.2 Simulation of Aliasing Effects.......................98
6.2.1 Classical Aliasing Effects......................98
6.2.2 Nonclassical Aliasing Effects.................100
VII. TECHNIQUES OF CLINICAL DATA ACQUISITION AND ANALYSIS... .'.102
7.1 Data Acquisition System..............................102
7.1.1 Patient Interface..............................102
7.1.2 Analog Signal PreProcessing...................103
7.1.3 AnalogtoDigital Conversion...................104
7.1.4 Preliminary Data Processing and Storage........105
7.2 Data Analysis System.................................107
7.2.1 Hardware Specifications....................... 107
7.2.2 Software: General Considerations...............107
7.2.3 Preliminary Data Processing....................107
7.2.4 Digital Filtering..............................108
7.2.5 Record Segmentation and Autocorrelation........109
7.2.6 Spectral Estimation............................110
7.3 Spectral Display and Analysis.......................Ill
7.3.1 The Compressed Spectral Array.................Ill
7.3.2 Statistical Description of the Spectrum........Ill
ix
VIII. CLINICAL HEART RATE VARIABILITY SPECTRAL ANALYSIS.........114
8.1 The Normal Resting Heart Rate Variability Spectrum...115
8.2 The Heart Rate Variability Spectrum on
Attaining the Upright Posture......................118
8.2.1 The Normal Response to Standing................118
8.2.2 Response to Standing with Abnormal
Autonomic Function...........................126
8.3 The Heart Rate Variability Response to Exercise.....133
8.4 Normalization of Autonomic Activity.................140
8.4.1 Sympathicotonic Orthostatic Hypotension........140
8.4.2 Sympathicotonic Orthostatic Hypotension
After Drug Therapy...........................146
IX. DISCUSSION, SUMMARY, AND DIRECTIONS FOR FUTURE RESEARCH..151
9.1 Discussion...........................................151
9.2 Summary.............................................155
9.3 Directions for Future Research..................... 159
REFERENCES......................................................161
FIGURES
Figure
21 Electrical model of the passive
cell membrane properties...............................12
22 The cell pacemaker potential...............................16
23 Block diagram of the autonomic nervous system's
influences on heart rate variability...................18
24 Control model of the baroregulatory system.................21
41 The integral pulse frequency modulator.....................41
42 The inputoutput relationship of the integral pulse
frequency modulator and the derivation of the
instantaneous frequency........................................44
4 3 Spectrum of a sinusoidallymodulated pulse train
produced by the integral pulse frequency modulator.....56
5 1 Derivation of the data points used in linear
interpolation of the instantaneous frequency...........66
5 2 Construction of the linearly interpolated estimate
of the instantaneous frequency.........................67
6 1 IPFM simulation: a single modulating sinusoid...........86
62 IPFM simulation: two modulating sinusoids...............88
63 IPFM simulation: three modulating sinusoids of
equal phases...........................................90
64 IPFM simulation: three modulating sinusoids of
different phases.......................................92
65 IPFM simulation: two modulating sinusoids, one in the
second order nonlinear range of the IPFM.............94
66 IPFM simulation: two modulating sinusoids, one in the
higher order nonlinear range of the IPFM.............96
67 IPFM simulation: classical aliasing effect..............99
xi
68 IPEM simulation: nonclassical aliasing effect.............101
81 Instantaneous heart rate of a healthy subject
resting in the supine position.........................117
82 Filtered instantaneous heart rate of
supine resting heart rate data.........................117
83 Compressed spectral array of supine resting
heart rate data........................................118
84 Instantaneous heart rate of a healthy subject
attaining the upright posture..........................121
85 Filtered instantaneous heart rate of healthy
subject attaining the upright posture.................121
86 Compressed spectral array of healthy subject
attaining the upright posture........................122
87 Total AC power response of healthy subject
attaining the upright posture........................124
88 Mean power frequency response of healthy subject
attaining the upright posture........................124
89 Standard error of the mean power frequency response of
healthy subject attaining the upright posture.........125
810 Midband power percentage response of healthy subject
attaining the upright posture........................125
811 Maximum power frequency response of healthy subject
attaining the upright posture........................125
812 Instantaneous heart rate of subject with autonomic
dysfunction attaining the upright posture.............128
813 Filtered instantaneous.heart rate of subject with
autonomic dysfunction attaining the upright posture....128
814 Compressed spectral array of subject with autonomic
dysfunction attaining the upright posture..............129
815 Total AC power response of subject with autonomic
dysfunction attaining the upright posture..............131
816 Mean power frequency response of subject with autonomic
dysfunction attaining the upright posture..............131
xii
817 Standard error of the mean power frequency response of
subject with autonomic dysfunction attaining
the upright posture...........................................132
818 Midband power percentage response of subject with
autonomic dysfunction attaining the upright posture...132
819 Maximum power frequency response of subject with
autonomic dysfunction attaining the upright posture...133
820 Instantaneous heart rate response to exercise...........135
821 Filtered instantaneous heart rate
response to exercise................................ 135
822 Compressed spectral array to exercise...................136
823 Total AC power response to exercise.....................138
824 Mean power frequency response to exercise...............138
825 Standard error of the mean power frequency
response to exercise..................................139
826 Midband power percentage response to exercise...........139
827 Instantaneous heart rate response to standing in
untreated sympathicotonic orthostatic hypotension.....142
828 Filtered instantaneous heart rate response
to standing in untreated sympathicotonic
orthostatic hypotension.......................................142
829 Compressed spectral array on standing in
untreated sympathicotonic orthostatic hypotension.....143
830 Total AC power response to standing in
untreated sympathicotonic orthostatic hypotension.....144
831 Mean power frequency response to standing in
untreated sympathicotonic orthostatic hypotension.....145
832 Midband power percentage response to standing in
untreated sympathicotonic orthostatic hypotension.....145
833 Instantaneous heart rate response to standing in
treated sympathicotonic orthostatic hypotension.......147
834 Filtered instantaneous heart rate response
to standing in treated sympathicotonic
orthostatic hypotension
147
Xlll
835 Compressed spectral array on standing in treated
sympathicotonic orthostatic hypotension..............148
836 Total AC power response on standing in treated
sympathicotonic orthostatic hypotension..............149
837 Mean power frequency response to standing in treated
sympathicotonic orthostatic hypotension..............150
838 Midband power percentage response to standing in
treated sympathicotonic orthostatic hypotension.......150
CHAPTER I
INTRODUCTION:
THE ANALYSIS OF HEART RATE VARIABILITY
When the normal human heart rate is accurately measured, it
is found to be not quite regular, but rather undergoes spontaneous
beattobeat variations of varying frequency and magnitude. The
most convenient method of accurately detecting and timing the
occurrence of a heart beat is through the use of the surface
electrocardiogram (EOG). The heart beat generates a large positive
deflection in the EOG termed the "R wave," and the study of heart
rate variability (HRV) has come to be known in the medical
literature as the study of "RR variability."
The most obvious variation in heart rate occurs at the
frequency of respiration, with the heart fate varying cyclically at
the respiratory frequency (approximately 0.20 to 0.25 Hz,
corresponding to 1215 breaths/minute), with a magnitude and phase
relationship which is a function of respiratory frequency (1). This
variation is termed "respiratory sinus arrhythmia". While the
mechanisms responsible for respiratory sinus arrhythmia are not
precisely known, it is felt to be related to neural reflexes
generated by changing pressure gradients inside the chest during
respiration (24). Corresponding oscillations in continuously
recorded blood pressure measurements have been termed "Traube
2
waves" (5).
Slower heart rate oscillations at a frequency of 0.1 Hz are
also found. A corresponding oscillation in the continuously
recorded blood pressure has been termed "Mayer waves," and appears
to be related to blood pressure regulatory mechanisms (6).
Finally, very low frequency variations in heart rate, on the
order of 0.05 Hz and less, can also be detected, and have been
ascribed to mechanisms involved in temperature regulation (7).
The variability found in the heart rate is precisely mirrored
in the similar variability of continuouslyrecorded blood pressure
measurements (8,9), suggesting that this variability may carry
significant information about the fundamental status and regulatory
activity of the entire cardiovascular system. As it is much easier
methodologically to measure HRV than to obtain the corresponding
blood pressure measurements (which require the use of an indwelling
intraarterial catheter), and, as the variability in both
parameters appears to be virtually identical, investigation has
focussed on the analysis of HRV as a probe of the status of the
cardiovascular system.
1.1 Methods of Analysis of Heart Rate Variability
Measurement of HRV is based on measuring sequential RR
intervals, or the instantaneous heart rate (IHR) derived from
these, where the IHR is defined as
3
IHR (beats/minute) = 60/RR interval (seconds)
11)
Factors affecting the decision to use measurements of period (i.e.,
RR intervals) or of frequency (i.e., IHR) are discussed more fully
in Chapter IV.
1.1.1 Statistical Methods of Analysis
Traditional methods of quantitating HRV have relied upon
statistical indices of the variability (10,11). Some of the more
widely used indices of HRV, and their respective advantages and
limitations, are summarized below.
a) Variance or standard deviation (SD) of the RR intervals,
where the SD is defined in the usual fashion:
While the SD yields a measure of the spread of the RR
intervals around the mean RR interval, to be meaningful it
requires that the mean heart rate remain constant over the sampling
period, which is rarely true in practice. The results are also
highly dependent on the mean heart rate, which is undesirable (10).
Finally, the results may be significantly affected by a small
number of erroneous data points.
b) Mean square successive difference (MSSD) of the RR
intervals, where
12)
[RR^RR,]2
N 1
13)
MSSD = Y.
4
While this method is relatively insensitive to slow trends in
the mean heart rate, it has been found to have poor correlation
with clinical results (10).
c) Difference between the maximum and minimum heart rate
(MaxMin HR) during the sampling period, where
MaxMin HR IHRmax IHRmin 14)
This statistic measures the total amplitude of heart rate
variation and is particularly applicable where a stimulus to
increase heart rate variability is applied. It is, however,
critically dependent on maintaining a constant mean heart rate
during the sampling interval and, in addition, a single erroneous
data point may completely invalidate the statistic.
d) Ratio of the longest to the shortest RR interval (Max/Min
RR) during the sampling period, where
Max/Min RR RRmax/RRmin 15)
This statistic is also critically dependent upon the mean
heart rate remaining constant during the sampling interval and may
also be completely dominated and invalidated by a single erroneous
data point.
This discussion serves to illustrate the problems in
adequately quantitating HRV by simple statistical parameters. In
addition to the limitations discussed for each specific statistic,
all suffer from a conmon shortcoming: only total HRV is
quantitated, allowing no discrimination between the various
components that comprise the total variability. Power spectrum
analysis of HRV overcomes many of the limitations associated with
these statistical parameters.
5
1.1.2 Power Spectrum Analysis of Heart Rate Variability
Analysis of the power spectrum of a timevarying quantity is
employed in diverse engineering problems. In this application,
power spectrum analysis allows quantitation of total HRV as well as
fractionating HRV into different spectral regions, allowing
discrimination of each of the components contributing to the total
variability. Filtering algorithms may be used to isolate important
components or to remove undesirable variability, such as slow
trends in the mean heart rate. Overall, it allows a much more
detailed analysis and quantitation of HRV than any of the
statistical parameters discussed above.
Power spectrum analysis of HRV was first described in 1973
(1214), but received little attention for another decade, perhaps
due to the sophisticated data processing resources required for its
implementation. Since 1981, there has been a resurgence of interest
in the technique, and a number of papers have been published in the
medical and engineering literature examining methods of performing
HRV spectral analysis and describing typical results. Since 1983,
papers have begun to appear in the medical literature describing
its applicability to actual clinical situations, with virtually all
investigators concluding that the technique may represent a major
advance in noninvasive cardiovascular monitoring in a variety of
clinical situations. Unfortunately, careful examination of these
papers reveals that most authors have not considered some basic
6
problems in performing heart rate variability spectral analysis
(HRVSA).
HRVSA is complicated by a fundamental problem: the signal to
be analyzed is a point process. That is, the data consists of a
discrete ordered set of points on the time axis representing tiroes
of occurrence of heart beats, and the power spectrum of a point
process is not uniquely defined (15). While various investigators
have used different methods to attempt to circumvent this problem,
a comprehensive analysis of the mathematical and physiologic bases
of HRV, and their implications in designing an appropriate spectral
analysis algorithm, has not been performed. The goals of this
thesis are to examine the physiologic basis of HRV, to explore the
mathematical models of the generation of HRV, and to incorporate
these fundamental considerations into the design of a technique to
perform HRVSA. Once a physiologically and mathematically sound
technique for performing HRVSA is derived, and a thorough
understanding gained of its capabilities and limitations, the
technique is then used to perform HRVSA in a variety of clinical
situations.
As an overview of this thesis, Chapters II, III, and IV are a
wideranging review and compilation of material from a broad range
of sources that are important in the understanding of HRV. These
chapters are a synthesis of material from widely disparate
disciplines, including cellular physiology, neurophysiology,
7
stochastic process theory, and communications theory, and establish
a firm theoretical framework for the material that follows. While
original criticism and conments are included throughout this
discussion, these chapters primarily review known material.
Chapters V through IX sill present original work on the theory and
practice of HRV analysis.
Chapter II discusses the physiology of HRV, including cardiac
electrophysiology and the role of the autonomic nervous system in
the regulation of the cardiovascular system and in the generation
of HRV.
Chapter III discusses the mathematical basis of methods of
determining the power spectrum of a point process such as a heart
rate event series, and describes the three fundamentally different
approaches that have been used to derive the power spectrum of HRV.
Chapter IV describes a particular model of HRV generation,
the integral pulse frequency modulator (IPFM), and reviews some of
its system properties.
Chapter V introduces new material to the literature on the
IPFM, including an analysis of nonlinear effects, unique aliasing
effects, and an analysis of methods of demodulating a signal
produced by an IPFM.
Chapter VI presents spectral analysis of data from an IPFM
simulated on a digital computer, which validates the mathematical
analyses presented in the previous chapters.
8
Chapter VII describes the techniques developed to perform
HRVSA, based on the mathematical model of the previous chapters,
and describes a complete system implemented on a general purpose
digital microcomputer.
Chapter VIII demonstrates experimental results obtained with
the system described in Chapter VII in a variety of clinical
situations.
Chapter IX discusses the results demonstrated in Chapter VIII
and suggests clinical implications and some potential uses for the
technique.
9
CHAPTER II
THE PHYSIOLOGY OF HEART RATE VARIABILITY
2.1 Basic Klectrophysiology
Living cells are surrounded by a cell membrane: a complex
layer of lipids and proteins which act as a semipermeable membrane
to ions found on both sides of the membrane. Specialized proteins
in the cell membrane act as ionic pumps, which, together with other
proteins which control the membrane's permeability to different
ions, generate an ionic concentration gradient across the membrane.
The major ions involved in this process include sodium, potassium,
and chloride. The ionic gradients produce an electrical potential
difference across the membrane in accordance with the Nemst
equation (16):
V = RT/ZF In [Ci/Co] 21)
where R = Universal gas constant
T = Absolute temperature
Z = Valence of the ion
F = Faraday's constant
Ci = Ionic concentration inside the membrane
Co = Ionic concentration outside the membrane
10
The relationship between the total membrane potential Em and
the concentrations of Na+ K+, and Cl is given by the
GoldmanHodgkinKatz equation (17)
ftT P,[jr]. + P,.[A;q). + fc.[Ct], *
where:
Pk, PNa, and Pci are the permeabilities of K+, Na+, and Cl
[K+]o, [Na+]o, [Cl]o are their extracellular concentrations
[K+]i, [Na+]i, [Cl~]i are their intracellular concentrations
It is clear from equation 22 that changes in the
permeability of the cell membrane to the different ions will lead
to changes in the membrane potential, and it apparently by this
method that cells control transmembrane potential.
Due to the ionic concentration gradients maintained by the
ionic pumps in the cell membrane and the permeabilities maintained
by the proteins in the cell membrane, the interior of a normal
human cell is approximately 70 millivolts moire negative than the
exterior of the cell, which is known as the resting potention Er.
In addition to these properties which are characteristic of
all human cells, certain cells (notably nerve, muscle, and cardiac
conduction tissue, which resembles nerve cells) also demonstrate
special properties of electrical excitation. These properties are
divided into passive and active processes.
11
At the junctions and interconnections of electrically
excitable cells, there are specialized connections called synapses.
These synapses release specific chemicals, known as neurotransmit
ters, in response to electrical stimulation. These neurotransmit
ters diffuse across the cell junction and alter the permeability of
the cell membrane on the other side which, in turn, alters the
transmembrane potential of the cell. Some neurotransmitters
increase (or hyperpolarize) the resting transroembrane potential,
which others decrease (or depolarize) this potential. A typical
nerve cell maintains connections with numerous other cells, each
joined through synapses dispersed over the cell membrane.
2.1.1 Membrane Passive Electrical Properties
At rest, the cell membrane can be modeled as a distributed
resistorcapacitor network similar to that used in cable theory and
depicted in Figure 21, and the passive electrical properties of
the cell are known as its cable properties.
12
Figure 21
Electrical model of a cable applicable to the resting cell
membrane. See text for physiologic equivalents of electrical
components.
The behavior of a cable such as that in Figure 21 is
described by the following differential equation (17):
+ EEr0
23)
where:
A.
^ m
f ifi^ m
is known as the membrane length constant, and
is known as the membrane time constant.
13
If current is injected into the cable at point X=0 and time
t=0 in Figure 21, then the voltage Em at point X=0 at time t is
given by
E n = E r + AE 2~4>
where:
Er is the resting membrane potential
a Eo is the steady state depolarization at t = 0
erf is the error function
and thus the cell is said to demonstrate temporal summation or
temporal integration.
After a steadystate distribution of potential has been
reached, the potential at any point x along the cable is given by
Em = Er + AEexp ( x/k) 25)
and thus the cell is said to demonstrate spatial summation or
spatial integration. The importance of these concepts will become
clear in the model of heart beat generation presented in Chapter
IV.
2.1.2 Membrane Active Electrical Properties
In addition to the passive electrical properties described
above, electically excitable cells demonstrate a special property
of active electrical excitation. When any point of the membrane is
depolarized to approximately 40 mv (the threshold) by the passive
measures described above, the cell membrane abruptly becomes highly
14
permeable to Na+ ions, which then enter the cell and change the
transmembrane potential to +40 mv, generating a spike known as the
action potential. The action potential propagates along the
membrane initially through its passive cable properties, which
causes other parts of the membrane to depolarize actively, thus
propagating an active spike throughout the membrane. The action
potential can also pass across synapses to affect other cells by
the mechanisms previously described. It is this active electrical
conduction that permits long distance conmunication (i.e., greater
than the few millimeters permitted by passive cable conduction) in
the nervous system. The action potential is an allornone
phenomenon; that is, there is no gradation in action potentials,
and they thus resemble digital signals in that information is
carried in the time sequence of the spikes rather than in the
structure of the spikes themselves.
2.L3 The Role of Spatial and Temporal Integration
The generation of an action potential in an electrically
excitable cell is controlled by the synaptic integration mechanisms
previously described. The membrane receives both depolarizing (or
excitatory) and hyperpolarizing (or inhibitory) stimuli separated
in time and space. The generation of an action potential depends on
the integration of influences combining at a given point on the
membrane in sufficient strength to generate an action potential.
This is thought to be the fundamental method of control and
15
integration of nervous system function (18). As will be discussed
in Section 2.1.5, this is also the model used to describe the
control of heart beat generation.
2.1.4 Pacemaker Cells
Certain cells exhibit spontaneous action potential generation
in the absence of stimulation. These cells are known as pacemaker
cells and are widely distributed throughout the nervous system and
heart. Pacemaker activity results from a special property of the
membrane: the resting membrane potential has no stable value, due
to an increased membrane permeability to sodium. The resting cell
membrane potential falls slowly toward zero rather than remaining
constant as it does in nonpacemaker cells; this slow depolarization
is termed the pacemaker potential.
When the pacemaker potential reaches the theshold voltage, an
impulse is generated, The rate of initiation of impulses depends on
the slope of the pacemaker potential, as shown in Figure 22. This
slope is highly dependent on the temperature, the ion concentra
tions, and the presence or absence of small concentrations of the
neurotransmitters acetylcholine and epinephrine (19).
16
B
Figure 22
Demonstration of how varying the slope of the pacemaker potential
varies the rate of spontaneous depolarization of a cell. Er is the
resting membrane potential and Eth is the membrane threshold for
action potential generation. A) is the unstimulated rate of
depolarization. B) represents the effect of adding a depolarizing
agent such as epinephine, while C) represents the effect of adding
a hyperpolarizing agent such as acetylcholine.
2.1.5 Cardiac Electrophvsiology
Heart muscle is electrically excitable tissue; depolarization
of the cardiac muscle cell membrane causes the cell to contract
through a process known as excitationcontraction coupling.
Normally, each heart beat is initiated by the depolarization of a
specialized group of pacemaker cells known as the sinoatrial (SA)
node. The wave of depolarization then spreads throughout the heart
17
through a specialized conduction system, which ensures that each
part of the heart depolarises, and thereby contracts, at the
appropriate time for maximum mechanical efficiency.
Consisting of pacemaker cells, the SA node is inherently
rhythmic, and the isolated heart will contract rhythmically (and
regularly) under the sole influence of the SA node. However, the SA
node receives synaptic input from nerve fibers, some of which
release acetylcholine, which flattens the slope of the pacemaker
potential and thereby slows heart rate, while others release
epinephrine and norepinephrine, which increases the slope of the
pacemaker potential and thereby increases heart rate. It is by this
mechanism that heart rate is controlled by the nervous sytem.
2.2 The Autonomic Nervous System
The autonomic nervous system (ANS) is the part of the nervous
system which functions unconsciously to control vital bodily
functions such as heart rate, blood pressure, bowel and bladder
function, temperature regulation, sexual function, and regulation
of blood flow to the various parts of the body.
The ANS is divided into two functional components: the
parasympathetic nervous system and the sympathetic nervous system.
The parasympathetic nervous system releases the neurotransmitter
acetylcholine at synapses, and thus tends to be an inhibitory
influence. The sympathetic nervous system releases the neurotrans
mitters epinephrine and norepinephrine, and thus tends to be an
excitatory influence. Most organs receiving input from the ANS
18
receive ennervation from both branches; thus an increase in heart
rate, for example, can result from either decreased parasympathetic
stimulation or from additional sympathetic stimulation.
2.2.1 ANS Control of the Cardiovascular System
The ANS is a nonlinear, negativefeedback control system. A
simplified block diagram of the physiological interactions which
affect heart rate variability is demonstrated in Figure 23.
Figure 23
A schematic block diagram of the structure of the autonomic nervous
system's interactions and its influences on heart rate variability.
The sinus node (center), which initiates each heart beat, receives
input from both the sympathetic and parasympathetic branches of the
system, both of which receive multiple influences from multiple
physiologic parameters. (Adopted from (20))
19
Nonlinear control systems exhibit characteristics which are
quite different from that of linear systems, such as frequency
entrainment and limit cycles, both of which are exhibited in the
behavior of the ANS (12).
Blood pressure regulation is a relatively simple and
wellunderstood example of the negativefeedback control maintained
by the ANS. Blood pressure is generally modelled try a simple Ohm's
Law relationship:
BP = CO x TPR
26)
where
BP = Blood pressure
CO = Cardiac output
TPR = Total peripheral resistance
Blood pressure, the controlled variable, is sensed by
receptors in the large blood vessels leaving the heart, and a
control signal is transmitted by special nerves to integrating
centers in the brain. A fall in blood pressure is compensated by an
increase in heart rate (thus raising cardiac output) and
contraction of the smooth muscles of the blood vessels, narrowing
their diameter and thus raising TPR. A rise in blood pressure is
compensated in inverse fashion. This nonlinear control system is
believed to oscillate in a limit cycle at a frequency of 0.1 Hz
(12), which is felt to be the origin of the Mayer waves described
in Chapter I, and is later examined in more detail.
20
As an example of the complexity of the ANS, other controlled
variables can also influence heart rate and blood pressure, as seen
in Figure 23. Arterial oxygen and carbon dioxide tensions (p02 and
pC02), as well as arterial pH, are also sensed by special
chemoreceptors, and changes in any of these quantities will lead
the ANS to change the heart rate, blood pressure, or TPR in an
attempt to return these quantities to normal.
Because of the importance of the oscillations in the blood
pressure regulatory system (also known as the baroreceptor system),
which, as will be seen, contribute a prominent and important
component in the spectrum of HRV, it is now examined in more
detail.
It is known experimentally that these oscillations exhibit
entrainment with oscillations at other frequencies (e.g.,.
thermoregulatory oscillations at 0.05 Hz and respiratory oscilla
tions at 0.150.25 Hz), so that any model of the baroreceptor
system must account for frequency entrainment. A general model
which has been developed to simulate physiologic oscillations is
shown in Figure 24. In the case of the baroreceptor system, there
are three basic sections: the linear characteristics of the
peripheral vascular bed, a pure time delay, and a nonlinear switch
element.
21
Linear filter
Srltdi eleiBit
Itfamce 7 'N ^ . M %
 + M
Paata*,
Incorporating
cure time delay
G(juj)
Output jf(t)
Effector
diaractert sties
of uiuulli nuscle
embroiling psipheral
vascular resistance
or afferent
pattwy
Stlnulus
2(t)
Figure 24
Control model of the baroreceptor system which displays nonlinear
characteristics such as spontaneous oscillations and frequency
entrainment. The pire time delay in the linear filter gives rise to
spontaneous oscillations, while the nonlinear switch element
permits frequency entrainment.
The characteristics of the peripheral vascular bed have been
studied by a number of workers, and the frequency response, using a
least squares fitting technique (20), has been found to be
22.0(1 + 12 you )
(1 +92/cu)(l + 2ju>)
G(joj)
27)
22
Such a frequency response does not exhibit oscillatory activity.
However, a modification which incorporates a pure time delay
changes the situation dramatically. The open loop transfer function
of the model including a pure time delay is
GO'ou)
22.0(1 + 12 j(v)e~,a>T
(1 +92;o;)(l + 2y'u>)
28)
where T is the length of the pure time delay. The point of
oscillation occurs when the locus crosses the negative real axis of
the Nyquist plane, and can be shown (20) to be the solution to
944cu 2tan (cu 7) + 82cu + 2208cu3 + tan (o> 7) = 0
To produce an oscillation of 0.1 Hz, equation 29 requires a time
delay T of 3.3 seconds. Interestingly, studies have shown that the
latent period between electrical stimulation of baroreceptor nerves
and a change in the heart beat period in anesthetized dogs is on
the order of 2 to 3 seconds (21), lending support to the validity
of the model.
Frequency entrainment can be shown to occur by describing
function analysis (20) when a nonlinear switch element is added to
the control loop. Therefore, the model of Figure 24 gives
excellent agreement with experimental results.
Thus the ANS is seen to be a complex nonlinear control
system regulating fundamental bodily functions. It has been
extremely difficult to obtain adequate methods of evaluating ANS
function and activity, as the nerves involved are deep inside the
23
body and not accessible to examing techniques. However, over the
past decade, investigation has focussed on HRV as a probe of ANS
function and activity.
2.2.2 The Origin of Heart Rate Variability
The SA node is ennervated by both parasympathetic and
sympathetic fibers from the ANS, as diagraimoed in Figure 23. As
described above, in the absence of external stimulation, the heart
beats regularly under the influence of the spontaneously
depolarizing pacemaker cells of the SA node. However, variability
in the heart rate arises from the influence of external influences
from the ANS, and evaluation of HRV is a simple, yet powerful,
method to obtain information on the activity and function of the
ANS.
At any given instant in time, the SA node is receiving
inhibitory influences from the parasympathetic nervous system and
excitatory influences from the sympathetic nervous system, which in
turn are being controlled by integration and control mechanisms in
the brain and spinal cord. These stimuli are integrated both
spatially and temporally, as described in Section 2.1.3, at the SA
node cell membranes, sunroing in a complex fashion to form a total
instantaneous ANS influence on the SA node, which will be denoted
the autonomic activity signal (AAS). The goal, therefore, of
measuring HRV is to discern the characteristics of the AAS, and
investigation of HRV has become an accepted tool in the measurement
of ANS function.
24
2.2.3 Clinical Implications
Knowledge of the characteristics of the MS is useful in at
least two clinical situations. First, the ANS itself can be tested
by detecting abnormal HRV responses to stimuli known to affect the
ANS, thus permitting the diagnosis of ANS dysfunction. Many
diseases (including diabetes, kidney failure, and alcoholism) can
damage the ANS with wideranging clinical consequences, and
measurement of HRV and comparison with normal values allows the
detection of ANS damage (22,23). Secondly, if the ANS is
functioning normally, then shortterm changes in the AAS may be an
indication that the ANS is compensating for some unknown
physiologic stimulus such as changes in the arterial p02, pC02, or
pH (24,25), falling blood pressure, abnormal cardiac function
(2628), etc, all of which may be useful in situations such as
anesthesia (29,30) or intensive care unit monitoring (31,32).
Finally, detailed knowledge of the characteristics of the AAS may
serve as a tool for investigating fundamental cardiovascular
regulatory mechanisms (33,34).
The goal, therefore, of the remainder of this thesis, is to
integrate the knowledge of the physiology of HRV, from the
biochemical level up to the level of the intact organism, in the
design of a physiologically plausible scheme to analyze the MS.
The next chapter discusses methods that have been used to derive
the power spectrum of HRV.
25
CHAPTER III
METHODS OF DETERMINING THE POWER SPECTRUM
OF HEART RATE VARIABILITY
Heart rate data is a point process; i.e., the relevant data
consists of times of occurrences of heart beats (the point events)
on the time axis. As a point process, the individual
characteristics (e.g., shape and amplitude) of each beat are
neglected, and its time of occurrence or, equivalently, the
interval between successive beats, is taken as the only relevant
variable. Unlike that of continuous signals, spectral analysis of
point processes is complicated by the fact that the notion of a
spectrum of a point process is not intuitively obvious, and
different methods are employed to define and thereby estimate its
spectrum. This chapter reviews the three fundamentally different
approaches that have been used in the literature to estimate the
spectrum of HRV data. These methods fall into three basic
categories:
i) Stochastic point process theory, which uses statistical
methods to estimate the spectrum, and ignores the
specific dynamics and models relevant to the process of
heart beat generation. Stochastic theory defines two
26
relevant but nonequivalent spectra of random point
processes: the interval spectrum and the spectrum of
counts.
ii) Time series methods, which view the successive intervals
between beats as values in a time series, and uses
classical methods of deriving spectral extimates of a
time series (i.e., parameter estimation and estimation
of the periodogram) in deriving the spectrum of HRV.
iii) Comnunications methods, which view the series of heart
beats as a modulated spike train and attempt to recover
the original modulating signal from the event series.
This method models a modulator based on the specific
physiology of heart beat generation.
Section 3.1 reviews in some detail the theory of stochastic
point processes, particularly with regard to the problem of
spectral estimation, as this method has been widely utilized in the
literature of HKVSA, as well as in the related and mathematically
equivalent problem of nerve fiber spike trains.
Section 3.2 briefly covers the use of time series methods.
This approach has not been widely utilized and is included here
primarily for completeness, although it does have certain
advantages.
Section 3.3 introduces the concept of demodulating the heart
beat event series in order to recover the original modulating
signal, the AAS. This section serves primarily to set the stage for
27
Chapters IV and V, which cover in detail a specific model used in
the communications approach to this analysis: the integral pulse
frequency modulator.
3.1 Spectral Analysis Using Stochastic Point Process Theory
3.1.1 Stationary Point Processes
A point process is stationary if the joint distribution
function of the number of events in k fixed intervals
(ti', ti"], (t2*, t2"].......(tk', tk"]
is invariant under translation for all k = 1, 2, ... .
There are two complementary aspects of a stationary point
process. The first is the counting process Nt, the cumulative
number of events in an interval t following the arbitrarily
selected event where the observation of the process begins. The
second aspect is the sequence of random variables {Xi, X2, ....}
where the Xi are the times between successive events. Note that Nt
is a continuous function of time that can take on only integer
values, while X takes on a continuous range of values, but occurs
as a sequence of intervals {Xi} in time. These two aspects are
connected by the following fundamental relationship:
Prob(/V, t) 31)
r = 1, 2, ...
Consequently, given the distribution of counts, it is
theoretically possible to obtain the joint or marginal distribu
tions of {Xi}, and vice versa. However, the two aspects of the
28
process are equivalent only through their complete distributions
(equation 31). If a statistical analysis is based only on first
and second order properties, the implication of equation 31 is
that the first and second order properties of both the counts and
the intervals are informative. In other words, they are not
equivalent analyses.
3.1.2 Second Order Properties of Intervals
The series of times between events {Xi} is a stationary
sequence of random variables, the subscript i denoting the order in
which the random variables appeared in time. The autocorrelation
function or sequence is
It can be shown (15) that a sequence (p is the autocorrelation
sequence for a stationary random sequence {Xi} if and only if each
p, can.be represented in the form
where F(o) is a symmetrical distribution function called the
spectral distribution function with derivative f(o) called the
spectral density function. Equation 33 then becomes
cov(;rt,AVt)
var (X)
k = ..., 1, 0, 1, ... 32)
k = ..., 1, 0, 1, ... 33)
k = ..., 1, 0, 1, ... 34)
29
and it is clear that the p*'s are the Fourier coefficients of f(<).
The inverse relationship is
2 it
1 + 2 ^ ptcos(fcco)
ti
n < uo < n
35)
Since p = p for real sequences {Xi}, and, from equation 35, f(to) =
f(to), a spectral density for positive to can be defined:
/.() =2/(co)
= ^/ 1 + 2 ^ ptcos(fcu;)
0 < to Â£ n
36)
3.1.3 Second Order Properties of Counts
The second order properties of the counting process Nt are
more complex than those of the intervals. An alternative and
mathematically equivalent way to examine these properties is to
consider the differential process Nt}, where & Nt is the number
of events in (t, t + rt]; i.e., and taken in the limit as rt
approaches 0. This generates a function p(t) which is a series of
delta functions at the event tiroes, so that the second order
properties of counts are equivalent to the second order properties
of a random train of impulses (3537).
Define
30
m,(t) = A/^r)
= lim
4*0
prob{event in t + x ,t+ x + At  event at t}
77
37)
and defining m as the mean rate of the process, then
m = 1 / E {X} = lim
AT* o
e[nm)
AT
38)
Using the definitions in equations 37 and 38, it can be
shown (15) that the covariance density of the differential process
is
y + (T) = lim
Al*0
cov{AN ttAN t.t)
MO*
= m[m m + 6(z)}
39)
A complete spectral density function for the counting process, the
Fourier transform of ?.(r), may be defined as
g(cu ) f m6{t)e~,a>xdt + f m[m/[t) m}e~,a>xdr
2 7t Jm
+ f [m.{r)m)e~/a>x dx
2n 2nJy n ' 1
Again, it is convenient to define a spectral density for
nonnegative to by:
31
g + (co) = 2 g(o>)
m
e',wtdr; cu Â£ 0
311)
n
Since the stationary processes are real, mf(*) = mf(r) and mf(t)
tends to m as t*, An alternative form of equation 311, therefore,
is
where mj(s) is the onesided Laplace transform of rnf (t).
3.1.4 Estimation of the Interval Spectrum
The interval spectrum is typically estimated by the
periodogram, which is defined through the following derivation.
Define:
312)
313)
and
314)
Also define:
32
Anoj(V') + jB^joo)
2{n
i a.
'Â£xaeJmt
61
jis;
which is the finite Fourier transform of the Xs's. Since the phase
between Xs and '** is not of interest, then define:
!*>(.&) =\Hw(cv) 2
_ Ata(cu) + 3~1
4 it
The plot of i(oo) versus
While the periodogram can be shown to be an asymptotically
unbiased estimate of f(w), it is not a consistent estimate. The
properties of i(ou) are summarized below:
Â£{/ (co)/(co);
var{/no(cu)} = /(ou)2;
cov{/no(co1),/no(cu2)} = 0;
It is possible to obtain consistent estimates of f (<*>) by
using spectral smoothing procedures (15), at the cost of
introducing bias and nonzero covariance into the estimate.
3.1.5 Estimation of the Spectrum of Counts
0 < ou < it
0 < ao < n
co i co 2 318)
33
The spectrum of counts can be estimated by a straightforward
calculation based on the event series. Each event is replaced by a
delta function, so the process can be represented by a signal p(t),
where
where tk is the time of the kth occurrence of an event. For equal
intervals, (i.e., I = tk tki), the spectrum of counts consists
of an infinite series of delta functions spaced at distance 1/1 on
the frequency axis. Usually, one is interested only in frequencies
much lower than the mean repetition frequency of the events, so
only the low frequency part of the spectrum needs to be considered.
One method of estimating the low frequency portion of the
spectrum is to pass p(t) through an ideal lowpass filter with
cutoff frequency fmax; this is equivalent to convolving p(t) with
the function einc(n/.t) and amounts to replacing each delta function
at time tk by the function . The result is a continuous
signal termed the lowpass filtered event series (38). An efficient
algorithm for performing this convolution was described by Holden
and French for neuronal spike train data (3942). The socalled
Holden French algorithm takes advantage of the synmetry and
periodicity of the sine function to reduce the number of
computations by a factor of 2no from that required by direct
319)
34
convolution. The continuous signal derived from this convolution is
then regularly sampled, and the spectrum is calculated by a
discrete Fourier transform.
3.1.6 Comparison of the Interval Spectrum
with the Spectrum of Counts
That these two spectra are fundamentally different can be
seen in the fact that f+(co) is a periodic function, while g + (a>) is
not periodic. The lack of periodicity in g+(>) arises because the
corresponding function in the time domain, y.U) is a function in
continuous time. Nevertheless, it can be shown (43) that if the
deviation in the intervals is small compared to the mean interval
length, the two spectra are equal to first order over their common
frequency ranges.
3.2 Spectral Analysis Using Time Series Methods
Heart rate data may be viewed as a time series, in which the
RR intervals are considered to be a set of data points obtained in
discrete time. One complication with this method is that the usual
methods of time series analysis require a uniform sampling
interval, while hart rate data occurs at irregular times.
Therefore, the sampling interval is directly related to the value
of the sampled function (a situation known as intrinsic sampling).
It can be shown (4447) that if the mean sampling frequency exceeds
the Nyquist rate, then the signal can, in theory, be completely
35
reconstructed. Interpolation schemes have been described, using
simple polynomial interpolators, that provide better performance
than a WienerKolmogorov filter (48,49).
A simple way of circumventing the requirement of regular
sampling is to consider each sample as a function of interval
number, rather than of time; i.e, RRi is associated not with a
time, but rather with interval number i. Spectral techniques then
yield data with units of cycles/beat rather than the usual
cycles/second (Hz), however the data can be converted to a function
of time by normalizing the abcissa by the mean interval length.
HRV has been analyzed in this way (50,51), using a general
autoregressive (AR) model, in which the generating process is given
by:
where p is the order of the autoregressive process and e(k) is a
zero mean white noise process with variance i.e.,
320)
E{e(fc)> = 0
E = o26,k
321)
The system function H(z) is the Z transform of 320:
36
H (z)
Y(Â£)
E(z)
1
1 +
Â£
(tnZ
322)
The power density spectrum S(f) is then determined as
S(/) = G7At  1 +
n1
2/n/iJ
2
Various methods may be used to estimate the set of
autoregressive coefficients {ai} and the white noise variance a*,
including the YuleWalker equations and the Batch Least Squares
method. The spectral estimate is obtained by substituting these
estimates into equation 323. The Batch Least Squares method
identifies normal HRV data as a fifth order autoregressive process
(50), and spectral estimates using this method show remarkable
similarity to those obtained from interval spectra and the spectrum
of counts using stochastic methods.
Theoretical advantages of using time series techniques in HRV
analysis include the ability to plot a pole diagram of the process,
which sunmarises the relevant process dynamics succinctly and
allows direct comparison between processes. In addition, well
established techniques for forecasting and prediction may be used
as a method of arrhythmia detection and evaluation (49).
3.3 Spectral Analysis Using a Conmmication Model
37
A third method used to perform HRVSA views the heart beat
event series (HBES) as a modulated spike train, and attempts to
reconstruct the modulating signal ty effectively demodulating the
HBES. This method converts a discontinuous signal (i.e., the HBES)
into a continuous signal which represents the original modulating
signal, termed the autonomic activity signal (AAS), to which
standard techniques of spectral estimation may then be applied. In
this sense, the method resembles that used to estimate the spectrum
of counts by the Holden French algorithm, where lowpass filtering
is used to convert the discontinuous HBES into a continuous signal
which can then be analyzed by standard spectral analysis techniques
(i.e., regular sampling and application of the discrete finite
Fourier transform). However, in the Holden French algorithm,
demodulation is performed by lowpass filtering alone, while in a
communications model, demodulation may take a more complex form.
The design of a demodulation scheme requires a model of the method
of modulation. A model that has gained wide acceptance is the
integral pulse frequency modulator (IPFM), because of its
qualitative similarity to the known methods of spike generation in
the heart.
Chapter IV examines the physiologic justification for using
the IPFM model and explores some of its known properties, such as
frequency response, which are important in the design of a
demodulation scheme. Chapter V then presents original material
examining limitations of the model, such as nonlinear effects and
aliasing, which limit the ability to perform perfect demodulation.
38
These results are then incorporated into the design of a
demodulation scheme and experimental results are shown which
demonstrate the method's superiority over the stochastic methods of
spectral estimation discussed in the present chapter.
39
CHAPTER IV
THE INTEGRAL PULSE FREQUENCY MODULATOR
The integral pulse frequency modulator (IPFM), also known as
an integrate to threshold modulator or a voltage to frequency
converter, is a mathematical model which converts a continuous
input signal m(t) into a frequency modulated series of pulses p(t).
In coitmunications systems, this method of modulation is known as
pulse position modulation. The IPFM model has been widely used
(38,43,5360) in the investigation of HRV and in the related
problem of nerve fiber spike trains.
This chapter describes some previously derived properties of
the IPFM and discusses the physiologic justification for its use.
There has also been debate in the literature (6568) whether it is
more appropriate to use heart rate (i.e., instantaneous frequency)
or heart period (i.e, RR intervals), where the two are related as
demonstrated in equation 11, in analyzing HRV. Section 4.2
contains an original derivation which demonstrates that, based on
the IPFM model, the instantaneous heart rate (IHR) is linearly
related to the input signal m(t) over all modulation depths, while
40
the interbeat intervals are nonlinearly related to m(t). Thus, it
is shown that the IHR has advantages in recovering the original
modulating signal from the modulated spike train.
Section 4.3 analyzes the physiologic justification for the
use of the IPFM model. Section 4.4 derives the frequency response
of an ideal IPEM, which will be necessary in designing a
demodulation scheme. Section 4.5 presents the spectrum of a pulse
train modulated by a sinusoidal input signal, and demonstrates the
problems which occur in attempting to use stochastic point process
theory in deriving the spectrum of this pulse train. This will
then lead to Chapter V, where an original demodulation scheme,
based on first deriving an intermediate signal, is presented and
its properties evaluated.
4.1 Description of the Integral Pulse Frequency Modulator
In the most general case of the IPFM, an input signal m(t) is
integrated with respect to time until a threshold r(t) is reached;
when the integrated signal exceeds r(t), a pulse is generated and
the integrator is reset to zero, whereupon the cycle starts again.
The sequence of emitted pulses is the modulated pulse train p(t).
This system is shown schematically in Figure 41.
41
> INC
// //, '/ V/
TIME
tik
Figure 41
The integral pulse frequency modulator, demonstrating a block
diagram of the modulator and the mechanism of impulse modulation.
In the study of HRV, several restrictions are imposed on this
general model:
1) The threshold r(t) is assumed to be constant in time and
can be replaced by the constant r.
2) The integrator is assumed to be a perfect integrator;
specifically, the integrated value does not decay in
time (as it does in the socalled "leaky integrator").
3) The integrator is reset instantaneously to zero; i.e.,
there is no "dead time" or refractory period.
42
4) The spike train output p(t) consists of a train of delta
functions.
5) The input signal ro(t) can be decomposed into a constant
term mo and a time varying component mi(t); i.e., m(t) =
mo + mi(t).
4.2 Analysis of InputOutput Relationships of the IFFM
Using the model of Figure 41, any two pulse occurrence tiroes
ti and ti+i are related by:
In the special case where the input signal does not contain a
timevarying component (i.e., m(t) = mo), then equation 41 reduces
to
41)
42)
<
= 1)
and the unmodulated output frequency fo is given by
43)
in which a linear relationship is preserved between the pulse train
repetition frequency fo and the input signal mo.
43
In the more general case, however, where m(t) contains a
timevarying component, then equation 41 becomes
r
madt +
*! I
i
*!!
s J mx{t)dt
<
44)
Equation 44 can be manipulated to yield:
r
tt*i ~ t(
m0
45)
Defining f(i,i+l) as the instantaneous frequency (or, in the
case of heart rate variability, the instantaneous heart rate) over
the time interval (ti, ti+i), equation 45 becomes
/(i,i+ 1 )
1
t(i t
46)
Equation 46 demonstrates that the instantaneous frequency f(i,i+l)
is proportional to the DC value of m(t) plus a term proportional to
the mean value of mi(t) over the interval (ti, ti+i). Defining
mm(i,i+l) as the mean value of mi(t) in the time interval
(ti, ti+i) yields the result:
44
/(i,i+l)= /0 + 7'nm(Â£,i+ 1) ^(m.0 + mm(i,i+ 1)}
47)
If the instantaneous frequency f(i,i+l) is taken as the output
front a demodulator, then equation 47 demonstrates the desired
linear relationship between the input signal m(t) and the
demodulated output f(i,i+l). This is the discrete analog of the
instantaneous frequency in continuous FM systems, which also
maintains a linear relationship with the input signal. Figure 42
illustrates the relationship between m(t), p(t), and f(i,i+l),
demonstrated for a sinusoidal input.
M(I)
Figure 42
The relationship between the input m(t) and the output p(t) of an
IPFM, and the derived instantaneous frequency f(i,i+l). It is shown
in the text that f(i,i+l) is proportional to the mean value of ra(t)
over each time interval (ti, ti+i).
45
It should be noted that f(i,i+l) is associated not with ra(t) at
a specific point in time, but rather to its mean value over the
interval (ti, ti+i); this is a direct consequence of the averaging
effect of the integration. This suggests that f(i,i+l) should be
considered a property of the entire interval (ti, ti+l) and should
be associated with the entire interval, as demonstrated in Figure
42. To minimize errors in the recovery of m(t) from f(i,i+l), the
deviations of rai(t) from mm(i,i+l) must be small over each sampling
interval. A quantitative analysis of this requirement is given in
Chapter V.
If f(i,i+l) is taken to be the estimate of m(t) to be used in
estimating the spectrum of m(t), a fundamental problem arises in
that f(i,i+l) is not a widesense stationary signal so that,
strictly speaking, its spectum is undefined. In addition, f(i,i+l)
is not a continuous signal, and its spectrum is contaminated by
artifacts due to the discontinuities. Two methods of obtaining a
more wellbehaved function from f(i,i+l) for use in spectral
estimation are demonstrated in Chapter V.
Returning to equation 45, it can be seen that inverting both
sides of equation 45 yields:
1 48)
J o+ ?
demonstrating that the interbeat interval is nonlinearly related
to mm(i,i+l). Equation 48 can be expanded in a Taylor series,
yielding:
46
(i i + 1 )
m o
mnCi,i+ 1)
f f L g
49)
If mm(i,i+l) << mo (i.e., the timevarying component is much
smaller than the DC component), then second and higher order terms
can be ignored and
Therefore, to first order, both the instantaneous frequency and the
interbeat intervals are linearly related to rnm(i,i+l). If, however,
the deviations become large compared to the mean repetition rate,
then the instantaneous frequency retains its linear relationship
while the interbeat intervals take on a nonlinear relationship to
mm(i,i+l). This suggests that, if the IPFM is taken as a valid
model of impulse generation, then the instantaneous frequency is a
better choice of measured variable than is the interbeat interval.
Experimental evidence suggesting that this is, in fact, the case in
an actual physiologic system is presented in Section 4.4.
Chapter 2 described the process of temporal and spatial
integration that occurs at the cell membrane and their effects on
the spontaneous pacemaker rate. The IPFM model precisely reflects
this integrative process. The input signal m(t) is decomposed into
a steadystate component mo, which reflects the spontaneous
410)
4.3 Physiologic Justification of the IPFM Model
47
pacemaker potential, and a timevarying component mi(t), which is
integrated over time and influences the heart rate on a
beatbybeat basis.
Section 4.1 discussed some of the simplifying assumptions used
in this model; these will now be explored in more detail.
1) The integrator threshold r(t) is assumed to be a constant
value r. While the membrane threshold can be affected by changes in
ionic concentration gradients and temperature, these are either
longterm or pathological conditions. In the normal cell, the
membrane threshold appears to be a fixed quantity which is well
modeled by a constant value. A theoretical analysis has been
performed (52) which models r(t) as a random sequence: i.e., r(t)
is constant over each interval, but the interval values are the
realization of a random process. This analysis leads to a random
crossing problem, and the crossing problem between a random
function and a nonrandom function is not solved in general.
Therefore, for tractability and from known physiologic properties
of the membrane, r(t) is considered to be a constant value r.
2) The integrator is assumed to be a perfect integrator. This
implies that the time constant in equation 23 is long compared
to the mean interbeat interval. While is not known for the
complex interaction of autonomic nerve fibers at the sinus node,
measurements of from isolated animal muscle cell preparations
range from 35 to 75 msec (19), suggesting that this assumption may
48
not be valid. An alternative formulation to the perfect IPFM which
explicitly includes the effect of is known as the "leaky
integrator model" (53), which is described by
It can be shown, however, that the leaky integrator exhibits
certain characteristics, such as phaselocking, which are not
observed in HRV studies (53); in addition, the relevance of values
of observed in isolated nervemuscle cell preparations to the
much more complex problem of multiple nerve fibers ennervating the
sinus node is not clear. Thererfore, in the absence of a
theoretical model or reasonable experimental evidence to the
contrary, a perfect integrator model is assumed.
3) The integrator resets instantaneously. In reality,
electrically active cells demonstrate a "dead time" or refractory
period after impulse generation, during which time they are
unexcitable. This period lasts approximately 12 msec (16) and
consists of an absolute refractory period, during which time the
cell is absolutely unexcitable, and a relative refractory period,
during which time a stimulus of supernormal intensity may cause
impulse generation. While theoretical analyses have been performed
using a model with an exponentiallydecaying refractory period
(54), computer simulations have shown that the effect on impulse
generation is negligible if the refractory period is much less than
411)
49
the mean interbeat interval (55). Therefore, while this effect may
be significant in nerve fiber spike trains (with tens or hundreds
of spikes per second), the effect is negligible in heart beat
generation, where the mean interbeat interval is on the order of
8001000 msec.
4) The spike trains consists of delta functions. This
requirement is equivalent to stating that the time of a heart beat
can be measured to any desired accuracy. While the R wave duration
is typically 80 msec, an algorithm which detects the same point on
each R wave (e.g., the R wave peak) approximates this requirement.
While motion, such as breathing, can cause small shifts in the
relative time of the R wave peak with respect to the time of
impulse generation (56), these effects are small compared to the
sampling intervals typically used in data acquisition systems.
Therefore, errors in detecting the precise time of a heart beat
tend to be a function of the discrete sampling used rather than due
to any physiologic effect, and the output can be considered to be a
train of delta functions without loss of generality.
This analysis serves to indicate that while the perfect IPFM
may not be an exact model of heart beat impulse generation,
deviations from this model tend to be small and can probably be
ignored without serious consequence. The one potentially important
modification, the "leaky integrator," mast remain an open question
until experimental data concerning the time constant of the system
50
becomes available. Overall, the IPFM model is physiologically
plausible, mathematically tractable, and fits the known methods of
impulse generation in the heart.
4.4 System Properties of the IPEM
This section will derive the frequency response of the IPFM, as
knowledge of its frequency response will be important in the design
of an algorithm to recover the modulating signal m(t) from the
observed pulse train p(t). Because the IPFM is a nonlinear system,
the approach will follow that of Knight (54), which is to model the
IPFM as a linear system with a small perturbation. This will first
require the derivation of the frequency response of an intermediate
system: the sampling intergrator.
A sampling integrator converts a continuous input signal v(t)
into a set of discrete outputs a(tk) spaced at equal time intervals
T = tk tki. The output at time tk is proportional to the average
value of the input over the previous time interval; thus
The transduction from V(t) to a(tk) is strictly linear. To
determine the frequency rexponse, let V(t) = Vk e (jcot), yielding
412)
413)
51
It is convenient to express 413 in terms of the sampling frequency
fo = 1/T. The frequency response becomes
a(*) I '*
S.
0
B(cu)
414)
For >/,, B(to)" l and nulls are found at
The characteristics of the IPFM have been previously described
by equation 41, here modified to
* 415)
r = J V(t)dt
For V(t) = Vo (i.e., a constant), the output frequency fo has been
shown in equation 43 to be fo = Vo/r, and thus the system is
linear.
Now assume that the relevant values depart only slightly from
their steady state; i.e.,
tt ~ Â£*,o + f/t,i
Â£*1 = ffc1,0 +
vmV'+Vi (o 4~16
where the subscript o refers to the previously solved steady state
and the subscript 1 refers to a small departure. All three
departure terms will tend to alter the value of the integral 415.
However, since the integral must still total to the value r, these
departures must add to zero. Ignoring products of departures (which
will be small to second order), equation 415 yields
52
0 = Vra(t*lit*ili)+ j Vx(f)dt
tfcI.o
The alteration in the kth time span may then be written
417)
Tn,i t*,i ~ ifci.i
bringing equation 417 to the form
418)
*t..
0V0Tkil + J V^t)dt
,k.c~Tk.o
419)
This can be interpreted as showing that the time span must be
readjusted to just compensate a sampling integrator contribution
accumulated from Vl(t). Considering an oscillatory input perturba
tion
in equation 419 yields
T
1.1
To
v0
B(cv)V
420)
421)
or
T t.i
Vxitk)
To
Vo
S(co)
422)
53
which shows that the frequency response, using the interval
perturbation Tk,l as the output, has a gain proportional to the
mean interval and a 180 phase shift with respect to the input.
The frequency perturbation may be found from the interval
perturbation by using
/ = 1/7
1
7 + 7,
423)
in equation 421, yielding
_ J_ 424)
/*. i = i,i
* 0
From 421, this yields the desired frequency response
/ k.l 1 n, N
V it ) V T ^
vi[tk) Voi o 4_2S)
It must be noted that this analysis will hold to a desired
degree of accuracy only by choosing the input perturbation Vi(t)
sufficiently small. At zero input frequency, the inputoutput
relationship is strictly linear, and up to input frequencies a fair
fraction of it mean discharge rate, the preturbation may be made a
fair fraction of the mean input level.
54
Assuming, without loss of generality, that r = 1, 425 can be
manipulated to yield
fi(yo;) = eyB///sinc[/r///0] 426)
and
B(you) =sinc[/r///0] 427)
) = arg B(you)
where
I b(/co)  is the gain of the system
$(/>) is the phase shift
fo is the mean repetition rate.
Note that for this frequency response function, the gain is not
a linear function of the mean repetition rate, as it was for the
frequency response using the interbeat interval as the output
quantity (equation 422). In addition, there is no phase shift
between the input and output. This relationship was investigated
experimentally (57) in a frog nervemuscle preparation, where a
sinusoidally modulated pulse train was applied to the nerve and the
generated muscle tension was taken as the output variable. The
derived frequency response of this system was in agreement with
equations 427 but not equation 422. This important experimental
result confirms the concept that pulse frequency, rather than
interpulse interval, is the significant physiologic variable, which
lends credence to the derivation presented in Section 4.2. This
55
systems response analysis, therefore, lends significant support to
the use of the IPFM as a valid model and the choice of the
instantaneous frequency as the relevant variable.
4.5 The Spectrum of a SinusoidallyModulated Pulse Train
If a modulating signal
fm is the modulating frequency
9 is the phase of the modulating signed
is applied as the input of an ideal IPFM, then it can be shown (58)
that the output pulse train, which consists of delta functions, can
be expressed as
nij (t) = m.j cos(2/r/mt+ 0)
428)
where
m
429)
where
I is the impulse content
fo 1/T is the unmodulated pulse repetition frequency
mi/r is the modulation depth
is the phase of the modulating signal
56
is the time from the origin t=0 to the time of the first
impulse
Jn is a Bessel function of the first kind of order n.
The terms of equation 429 each represent a component of the
spectrum and are classified according to their frequencies of
occurrence in Figure 43.
Ifn
iWyrfni)
mil
w
MWMW(Wiffo)
1^1
Ko'V^Vrfni)
"fm 0
t
Figure 43
Spectrum of a sinusoidallymodulated pulse train produced by the
IPFM.
The first term, Ifo, is the zero frequency or DC component
and represents the average value of the pulse train. The second
57
term is the modulation or signal component, since it occurs at the
frequency fm of the modulation signal, and has magnitude mil/r,
which is proportional to the modulation depth. The remaining terms
are components and sidebands occurring at the carrier frequency and
its multiples fo, 2fo, 3fo, ...
The important features of this spectrum are:
1) The zero frequency component is accompanied by a single
component, without multiples or harmonics, at the
modulating frequency. Demodulation is therefore
possible by lowpass filtering without harmonic or
other distortion.
2) The spectrum contains components at all harmonics of the
carrier frequency, each accompanied by an infinite set
of adjacent components (or sidebands) which are removed
from these harmonics by multiples of the modulation
frequency.
These properties have led investigators to believe that the
modulating signal can be recovered from the modulated pulse train
by lowpass filtering, or equivalently, time averaging, with the
only requirement on the filter being that its averaging time be
short with respect to the modulation and long with respect to the
carrier period. It is clear that this is identical to the procedure
described in Section 3.1.5 in estimating the spectrum of counts by
the HoldenFrench algorithm, which uses an ideal digital lowpass
filter. A hardware filter using a cosine taper has also been
described and produces virtually identical results (59).
58
These results, however, are valid only for a single
modulating sinusoid. In the case of a modulating signal consisting
of the sum of two or more sinusoids, an expression for p(t,) cannot
be obtained in closed form. Simulation studies of a modulating
signal consisting of two sinusoids of different frequencies (60)
have consistently shown that the low frequency part of the spectrum
of counts, as well as the interval spectrum, contain components not
only at the two modulating frequencies but also at multiple sum and
difference frequencies between the carrier and modulating
frequencies. In fact, for the case of two sinusoids, the low
frequency (i.e., up to the carrier frequency) part of the spectrum
of counts and the interval spectrum contain eight distinct peaks,
and the situation worsens as more modulating frequencies are
introduced. In fact, it is difficult to know how to interpret the
spectrum of a pulse train generated by a modulating signal
consisting of a continuum of components, as presumably a
physiologic signal does.
It is not surprising that these spectra poorly reflect the
spectrum of the original modulating signal, since these spectra are
those of the frequency modulated pulse train p(t), and frequency
modulation is, in general, a nonlinear process (61). The quantity
of interest is the spectrum of the modulating signal m(t) and not
that of the pulse train p(t). Chapter V describes two methods of
obtaining an estimate of M(f), the desired spectrum of the
modulating signal m(t), by first deriving an intermediate signal,
the instantaneous frequency.
59
CHAPTER V
SPECTRAL ESTIMATION UTILIZING
THE INSTANTANEOUS FREQUENCY
The previous chapter demonstrated that the interval spectrum
and spectrum of counts of a modulated pulse train generated by an
IPFM differ from the spectrum of the original modulating signal by
the presence of spectral components at linear combinations of the
carrier and modulating frequencies. As the spectrum of interest is
that of the modulating signal, and the pulse train is of interest
only insofar as it carries information about the original
modulating signal, it is desirable to first derive a signal which
bears a linear relationship to the modulating signal m(t), the
spectrum of which can then be estimated. It was shown in Section
4.2 that the instantaneous frequency f(i,i+l) has this desirable
property of being linearly related to ra(t). The instantaneous
frequency, however, as illustrated in Figure 42, is a staircase
function and thus suffers from nonstationarity and discontinui
ties. Sections 5.1 and 5.2 demonstrate two original methods of
generating a continuous, stationary function from the instantaneous
frequency. Section 5.3 analyzes the sources of error in these
methods and describes a type of error that has not been previously
described in the literature. Section 5.4 gives a detailed analysis
60
of spectral estimation utilizing the techniques described in
Sections 5.1 and 5.2. Chapter VI then contains computer simulations
demonstrating the usefulness, as well as the limitations, of these
methods of estimating the spectrum.
5.1 LowPass Filtering
Section 4.3 demonstrated that the frequency response of an
IPFM, using the instantaneous frequency as the output quantity, is
given by
S(/) = e'"'"'Sinc[/r///0]
This suggests that if a filter with frequency response
H(/) = e
mtn.______I_____
sinc[rt///0]
51)
52)
is cascaded to an IPFM, then the Fourier transform of the output,
Y(f), is
YU) HU)BU)MU)
= e
*mtnt
1
mtn.
sinc[/r///0]
M{f) \ 0
sinc[/r///0] MU)
53)
which is precisely the desired spectrum of the modulating signal
m(t). The frequency range of the filter is limited to f < fo to
avoid a vanishing denominator (and subsequently infinite gain) in
H(f); in fact, there are additional, more stringent restriction on
61
the range of f. First, f can be restricted to f < fo/2 by the
Nyquist theorem. However, the derivation of B(f) in Section 4.4 was
based on modeling the perturbation in f as small compared to the
mean pulse repetition frequency; therefore, if f deviates
significantly from fo, B(f) no longer accurately describes the
frequency response. This will be examined in further detail in
Section 5.3.
While, in theory, a filter having frequency response H(f)
given in equation 52 could be designed, there are several
practical problems. First, because of the presence of fo in the
frequency response of H(f), the filter roust be adaptive; i.e., the
filter characteristics change for each set of data points.
Secondly, to perform filtering in the time domain, the Fourier
transform of sine1 [*f/fo] must be obtained, which is not
available in closed form. Nevertheless, suitable approximations
could be obtained and, in principle, such an adaptive filter could
be approximately realized. Since, however, the IPFM is a
mathematical abstraction which probably does not precisely describe
physiologic impulse generation, and, in addition, since B(f) is an
approximation which is only accurate for f << fo, it seems unwise
to expend inordinate effort on a complex filter designed to
precisely match its characteristics.
Recognizing the general principle, however, it appears that
lowpass filtering f(i,i+l) may yield a reasonable approximation of
M(f). Specifically, if f(i,i+l) is passed through a filter with
frequency response:
62
////.
e
0;
0
otherwise 54^
(i.e., a linear phase ideal lowpass filter), the Fourier transform
of the output Y(f) is
X(/) sinc[jr///0]tf(/) ; 0
0 ; otherwise 55)
which is easier to realize than that of equation 52, but still
requires an adaptive filter. Y(f) could then be normalized by
dividing the output spectrum by sine [nf/fo] and exact recovery of
M(f) for f << fo is still theoretically possible.
A final simplification is possible by noting that the resting
human heart rate is rarely slower than 60 beats per minute (or fo =
1 Hz). Therefore, passing f(i,i+l) through a filter
W(/) = e*'",/05 ; 0
0 ; otherwise 56)
then the output Y(f)
/(/)= sinc[/r//0.5] MU) : 0
0 ; otherwise 57)
satisfies the Nyquist criterion, and normalizing Y(f) by dividing
by sinc[2*f] allows errorfree recovery of M(f) in the range of
0 < f < 0.5 Hz and does not require adaptive filtering. The
63
disadvantage of this method is that if the mean heart rate is
greater than 1 beat/second, then spectral information is available
above 0.5 Hz which is not utilized.
Summarizing, it is seen that the spectrum M(f) of the
modulating signal m(t) can be recovered exactly for f << fo by
suitable lowpass filtering of f(i,i+l). An quantitative analysis
of the requirement that f be "much less" than fo is given in
Section 5.3.2. It is important to recognize that this method
differs from the lowpass filtering described in the estimation of
the spectrum of counts. In the case of the spectrum of counts, it
is the pulse train p(t) that is lowpass filtered, while in this
method, an intermediate signal f(i,i+l), which bears a linear
relationship to the original modulating signal, is first generated
and then lowpass filtered.
5.2 Polynomial Interpolation
While the previous section demonstrated that appropriate
lowpass filtering of f(i,i+l) provides an elegant method of
recovering the input signal m(t), sophisticated and timeconsuming
adaptive filtering was sometimes required, whereas, if adaptive
filtering is not used, then spectral information can potentially be
lost. Another method of estimating m(t) is to approximate f(i,i+l)
by a piecewise linear function and to use this function as an
estimate of m(t) for use in spectral estimation. In this context,
linear interpolation is a timevarying lowpass filter. While this
method is computationally simple, errors occur when m(t) deviates
64
significantly from linearity over a sampling interval (ti, ti+i);
however, as will be demonstrated in Section 5.3.2, deviations of
m(t) from linearity over a sampling interval occurs when it is no
longer true that f < < fo. Nevertheless, this method is computation
ally simple and fast and provides results similar to those obtained
by more sophisticated filtering techniques. Beutler (46) has shown
that the interpolation error in linear polygonal interpolation may
be smaller than that obtained from classical WienerKolmogorov
filtering, as the simple interpolators constitute timevarying
filters, whereas the WienerKolmogorov filter is constrained to be
timeinvariant.
Interpolation requires discrete sampled values at times
tk, k = l,2,...,n, while f(i,i+l) is a (discontinuous) function
defined at all time t. Some manipulation will be required to
convert f(i,i+l) into a set of discrete sampled values suitable for
interpolation. It was shown in Section 4.1 that f(i,i+l) is
linearly related to the mean of ra(t) over the time interval
(ti, ti+i), denoted rnm(i,i+l). While the mean value theorem from
basic calculus guarantees that there exists a time tm (where
ti < tm < ti+i) such that m(tm) = mm(i,i+l), there is, in general,
no way to determine tm. However, if m(t) is approximated by a
piecewise linear function, then it is, in fact, possible to
specify tm and, consequently, to specify m(tm), which is the
requirement for subsequent interpolation.
Let m(t) be approximated by a set of linear functions
{gi(t)}, defined as
65
g i(0 = att + bt
58)
where
ai is the slope and bi is the y intercept of the line
approximating m(t) over the interval (ti, ti+i).
Then the mean value of gi(t) over the interval (ti, ti+i) is
which, after algebraic manipulation, yields
gm.f(0
510)
Therefore, if tm,i is defined as (ti + ti+l)/2 (i.e., the midpoint
of the interval), then
9iltm.t)* rnm(i,i+1) 51.
The significance of this procedure is that now a set of discrete
times {tm,i} are defined, one in each sampling interval, for which
a set of sample values {mm(i,i+l)} are associated, as contrasted
with f(i,i+l), which associated nm(i,i+l) with the entire time
interval (ti, ti+i). This procedure is illustrated graphically in
Figure 51.
66
H(T)
Figure 51
Illustration of the derivation of the points mm(i,i+l): the set of
discrete data points used in linear interpolation. The justifica
tion for the association of mm(i,i+l) with the midpoints of the
instantaneous frequency signal f(i,i+l) is given in the text.
It is important to note that the choice of tm,i =
(ti+i + ti)/2 was based on the assumption that m(t) could be
approximated by a linear function over each interval (ti, ti + i). If
ro(t) deviates significantly from linearity over a given sampling
interval, then the value of tm, i will be incorrect and the sample
value will be erroneous, which will impact on the accuracy of the
interpolation to be described below.
67
Once the set of sample times {tm,i} and the associated set of
sample values {rnm(i,i+l)} are obtained, then linear interpolation
may be used to approximate m(t). This linearly interpolated
function will be denoted m~(t), as shown in Figure 52.
H(I)
Figure 52
Construction of m~(t): the linearly interpolated estimate of m(t).
M~(t) is derived by linearly interpolating between the points
mm(i,i+l), and the estimate can be compared to the original m(t) at
the top of the figure. tT (t) can be regularly sampled and its
spectrum M~(f) can be estimated as the discrete Fourier transform
of nT(t).
68
The spectrum of m(t) can be estimated by regularly sampling m~ (t)
and applying standard spectral estimation techniques to the sampled
values of rcT(t), as will be described in Section 5.4
5.3 Sources of Error: Aliasing and Nonlinear Effects
In order to analyze the errors which occur in the use of the
instantaneous frequency f(i,i+l) in the estimation of the input
signal m(t), it is instructive to more closely examine the
relationship between the two signals. The sources of error fall
into two general categories: aliasing and nonlinear effects.
5.3.1 Aliasing
Aliasing is a term which is generally used to describe the
effect of high frequency components being detected as a low
frequency component in a sampled data system due to an
insufficiently rapid sampling rate. Aliasing, however, may be
generalized to include any incorrectly detected signal due to
constraints in the detection system, and, in the IPFM, a unique
type of aliasing also occurs. Specifically, if rai(t) is an input
signal to an IPFM which generates fi(i,i+l) as the instantaneous
frequency output, then does there exist another signal m2(t) which
will also generate fi(i,i+l)? In general, the answer is yes, which
can be simply observed by noting that if fl(i,i+l) is the input to
an IPFM, the instantaneous frequency f(i,i+l) output will be the
identical fi(i,i+l). This occurs because the instantaneous
frequency output f(i,i+l) depends on the mean values of the input
69
signal over each sampling interval; therefore any set of input
signals which have the same set of mean values over each sampling
interval will generate the identical f(i,i+l) output. There are,
therefore, an infinite set of input signals which may have the same
set of mean values in each sampling interval and may therefore be
the input signal for any given output f(i,i+l).
This type of aliasing, which is unique to the IPFM, has not
been described in the literature, and is here denoted by the term
ftignftl aliasing.
The problem of signal aliasing does not have any simple
solution. While some possible input signals m(t) for a specific
output f(i,i+l) may be excluded by imposing certain constraints on
m(t), such as m(t) must be continuous, m(t) must be bandlimited to
a given frequency, etc., in the absence of direct experimentally
measured values of m(t) in vivo, these constraints are arbitrary
and are therefore suspect.
Sunmarizing, therefore, it is possible that an unexpected
input signal m2(t) may generate an observed f(i,i+l) which is then
interpreted as being generated by a different input mi(t), and
there is no way to determine which was the actual input signal.
While the potential for signal aliasing has been ignored in the
literature on HRV, it may turn out to be an important effect.
However, in order to determine its potential importance, it will be
necessary to measure m(t) (i.e., the autonomic activity signal)
directly in vivo by placing microelectrodes into the fibers of the
SA node of the heart and comparing the directly measured signal
70
with that determined by the instantaneous frequency method, and
this has never been done. For the time being, therefore, the
potential for signal aliasing must be kept in mind in interpreting
HRV data by the instantaneous frequency method.
In addition to signal aliasing, typical aliasing also occurs
in the IPFM. If fo is is the mean unmodulated pulse repetition
frequency of the IPFM, then fo can be taken as the sampling
frequency for the input signal. The Nyquist theorem then predicts
that if the input signal contains frequencies greater than fo/2,
then these higher frequencies will be aliased into lower
frequencies. An added complication, however, occurs due to the fact
that the sampling frequency is not constant, tut varies with the
input waveform (which is known as intrinsic sampling (45)). This
follows directly from equation 46, where it is shown that the
interspike intervals are a function of the mean value of the input
waveform. Therefore, as the heart rate increases, the interspike
intervals become shorter and the sampling frequency increases; the
converse occurs as the heart rate decreases. As will be seen in
Chapter VI, this causes unusual aliasing effects to occur.
5.3.2 Nonlinear Effects
The IPFM is, in general, a highly nonlinear device. While
its behavior was constrained to be linear in Section 4.4 by
constraining the deviations in the interspike intervals to be very
small compared with the mean interspike interval, it becomes
important to explore the situations in which its behavior becomes
71
nonlinear and what effects this nonlinearity will have in the
attempt to reconstruct m(t) from f(i,i+l). Analysis of the
nonlinear behavior of the IPFM has not previously appeared in the
literature on HRV, but it is of critical importance because, as
will be seen in Chapter VI, the IPFM operates in its nonlinear
range in most physiologic situations.
If the timevarying portion of the input signal m(t) is taken
to be the sum of two sinusoids; i.e.,
then any two contiguous spike occurrence times ti and ti+i are
related by
nri1(O = t,1sin(co1t + ^1)+62sin(a)2t + 02)
512)
513)
)+ j [b j sin (gu t + ,) + 62sin(cu 2t + 02)}dt
which, after integration, becomes
514)
72
It is possible to eliminate ti+i from equation 514, and thus have
a relationship that depends only on the starting time of the
interval, by recognizing that
Â£<1 = tt + tt.i 1(
= + (t{*i tj
t,+l/f(i,i+I) 5~15
Substitution of equation 515 into 514, and considerable algebraic
manipulation, then yields the relationship
.10
f(i
516)
Equation 516 is exact; i.e., there are no requirements of small
deviations or other aspects of linear behavior. Thus, f(i,i+l) is
seen to be a function of ti, ,, ua, bi, b2, and *s. However,
equation 516 is a nonlinear implicit transcendental equation
whose properties are not at all obvious from inspection.
If the frequencies , are small compared to the
instantaneous frequency f(i,i+l), then to first order
73
f(i.i+l)J f(i,i+l)
517)
f(i,i+l)J f(i.i+l)
gj 2 ^ o>2
and equation 516 reduces to
+ b1sin(a>1t(+01) + 62sin(oo2tj+02)
518)
f(i,i+l) =
r
which demonstrates a linear relationship between the instantaneous
frequency and the modulating signal, confirming the derivations in
Sections 4.1 and 4.2.
If, however, w, and a>a are not sufficiently small for the
linear approximations in equations 517 to be accurate, then to
second order
2f(i,i + l)z
cos
2f(i,i+l)2
519)
and equation 516 becomes
520)
74
which is a quadratic equation in f(i,i+l) with solution
itb, / , b 2 / \
f(M+l) ^r+^rsinlaj,i, + ^,) + sin(oj2ti + 02)
*\^w(a)1(aj2,^ ,02) + x(cu, +co2.0, + ^2) +y(2cu, ,2^, ,2a?2,2^2)
where
m: 2m
w(a>, .cu2.0, ,02) jÂ£[b,sin(cu1t, + ^ J + baSinfcOj^ + ^a)]
r r
52lb)
2
+ [b ,a>, cos(cu ,t, + ^,) + b2co2cos(a)2fl + ^2)]
x(w, +a>a,01 + *2)=^7icos((a>1 + cu2)t, + 0, + *J)cos((a),co2)t, + *1*a)]
52lc)
y(2
>?<
2rs
bl bi it \\ b? / / >\
521d)
which demonstrates that, to second order, f(i,i+l) is a function of
the additional parameters ,cua, *,+$, $,0a, 22coa, 2$,,
and 2*a.
It is instructive to examine the frequency ranges which
exhibit first, second, and higher order behavior. The first order
approximations in equations 517 are accurate to within 5% for
oo,/f(i.i+1) ranging from 0 to .314, and the second order
approximations in equations 519 are accurate to within 5% for
co,/f(ii+i) ranging from .314 to .84. Converting to cyclical
75
frequency, this implies that the system behaves approximately
linearly for 0 < f/f(i,i+l) <.0.1 and exhibits second order
behavior in the range 0.1 < f/f(i,i+l) <. 0.26. Above this range,
the system exhibits increasingly nonlinear behavior involving
higher order terms in a>. Therefore, the limit for useful signal
recovery is approximately half that of the Nyquist frequency, where
f/f(i,i+l) = 0.5.
Since the normal human heart beats in the range of 60 to 100
beats per minute, or an instantaneous frequency of 1.0 to 1.67 Hz,
this then suggests that for heart rate data the system operates
approximately linearly for 0 < f < 0.10 0.16 Hz, and demonstrates
second order behavior in the range of
0.10 0.16 < f < 0.26 0.43 Hz.
Referring to equations 520 and 521, if the input
frequencies are restricted to the ranges which demonstrate linear
and second order behavior, these equations suggest that f(i,i+l)
will have the following properties;
a) It will have a DC term on the order of mo/r.
b) It will have components on the order of bi/r, b2/r, ...,
bk/r at input frequencies ...co, for an input signal consisting
of k sinusoids.
c) It will have second order terms at linear combinations of
the input frequencies and phases.
d) Unlike the spectrum of counts in Section 4.5, it does not
contain components at multiples of the mean repetition frequency
fo mo/r.
76
As will be seen in Chapter VI, these predictions are verified
by experimental data showing computer simulations of the operation
of an IPFM and the estimation of the power spectrum of f(i,i+l).
Summarizing this section, it is seen that if the input signal
m(t) consists of frequencies which are small compared to the
instantaneous frequency (i.e., f/f(i,i+l) < 0.1), then there is a
linear relationship between m(t) and f(i,i+l). This is the same
requirement as that stated in Section 5.2 for accuracy in linear
interpolation, which stated that m(t) be piecewise linear over
each sampling interval. As the frequencies in the input signal
increase into the range of 0.1 < f/f(i,i+l) <0.26, the instanta
neous frequency was shown to be a function of additional and
uncontrollable variables, such as linear combinations of the input
frequencies and phases. This behavior, however, can be compared to
the spectrum of counts described in Section 4.5, which exhibits
nonlinear behavior over its entire frequency range. Finally, for
0.25 < f/f(i,i+l) < 0.5 (which is still below the Nyquist
frequency), the behavior becomes increasingly nonlinear and
difficult to characterize.
5.4 Spectral Analysis of the Instantaneous Frequency
This section describes the process of deriving the power
spectrum of the instantaneous frequency f(i,i+l). As discussed in
Sections 5.1 and 5.2, this signal is not analyzed directly, but is
smoothed using either a linear filtering algorithm or linear
interpolation, both of which serve to remove the discontinuities
77
from f(i,i+l). As the smoothed versions of f(i,i+l) were shown to
be estimates of the original input signal m(t) for f < fo/2, the
lowpass filtered version of f(i,i+l) will be denoted nff(t), and
the linearly interpolated version will be denoted m'i(t).
It is important to note that both nf Â£(t) and m*i(t) are
signals defined for all t, which were derived from the initial
point process p(t), and can therefore be analyzed using standard
spectral analysis techniques for continuous signals. The techniques
of spectral analysis for continuous signal are wellknown (63) and
will not be repeated here. The interested reader is referred to
standard references on the subject (63, 64). A summary of the
process is given below.
It must first be noted that the term estimate is used here in
two different contexts. First, nff(t) and m~l(t) are both estimates
of the original modulating signal m(t), with estimation accuracies
which depend on the factors discussed in Section 5.3. In addition,
the power spectrum estimate M"l(f) is an estimate of the power
spectrum of m'l(t) and the power spectrum estimate M*f(f) is an
estimate of the power spectrum of nf Â£(t). The way in which the term
estimate is used in the subsequent discussion should be clear from
the context.
In preparation for determining the power spectral estimates
M'f(f) and M'i(f), the estimates m~i(t) and nff(t) are first
regularly sampled at a sampling interval At, where At is chosen such
that fs = \tM is at least twice the highest frequency present in
m~i(t) and m'(f). However, as demonstrated in Section 5.3, due to
78
nonlinearities in the IPFM, the spectral estimates are reliable
only for f < fo/2, where fo is on the order of 1 to 2 Hz. It is
therefore unnecessary for fs to be greater than 2 Hz, or = 0.5
seconds. Therefore, ^t is chosen to be 0.5 seconds.
The power spectral estimates M~i(f) and M~f(f) are defined as
the discrete finite Fourier transforms of the discrete autocovari
ance function estimates Ci(k) and Cf(k) respectively, where
j Nk
ci(fc) = Z ('fti(u + fc)ml)(/fl,(u)/ril)
522a)
C,(k)
Nk
^^J^[rhr{u + k)m.f){rh,<;u)rn,)
522b)
and where
_ j " 522c)
mi=
** ILa 1
_ i " 522d)
rnr = T}Lrhriu)
IN Ul
are the sample means.
Using the symmetry of the autocovariance function, the power
spectrum estimates are
79
1&i(/)C,(0) + 2 Â£c,(k)l/(k)cos(27ifk)
it" 1
523a)
N
l5l/(/) = C/(0) + 2^C/(A:)l/(fc)cos(2^/A:)
ti
where W(k) is a smoothing lag window described below.
While these relationships are exact for deterministic
signals, it can be shown that for stochastic inputs, the power
spectral estimators yield unbiased but hot consistent estimates of
the spectrum (63, 64). (Refer to the discussion of the periodogram
in Section 3.1.4.) The variance of the power spectral estimates can
be reduced, at the cost of introducing bias into the estimates, by
smoothing the sample autocovariance functions with a lag window.
The choice of the width of the lag window is a critical one,
since the variance of the estimate decreases, while the bias of the
estimate increases, as the width of the lag window is widened. The
choice of the optimal width of the lag window is empirical, as the
window should be narrow enough to be able to distinguish distinct
spectral peaks, yet if it is overly narrow, instabilities appear in
the spectral estimate in the form of spurious oscillations (63).
While lag windows of various forms have been described (e.g.,
Bartlett, Parzen, Tukey, etc.), the width of the window is a more
critical parameter than the exact shape (63). For use here, a Tukey
window, which is a cosine tapered lag window, has been utilized.
Various lag window widths have been tested using clinical data
80
(Chapter VIII), and a window with bandwidth of approximately 0.04
Hz has been found to be optimal, distinguishing between the
spectral peaks present in clinical data while still minimizing
spurious instabilities in the spectrum. Therefore, all spectral
analysis performed herein smooths the estimated autocorrelation
function with a Tukey window with bandwidth 0.04 Hz prior to
computing the discrete Fourier transform.
An additional modification which has been found to be
important is the ability to remove the DC and very low frequency
components from the estimated signals mf(t) and mi(t). Due to the
finite length of the data record, the DC component of the spectrum
exhibits spurious leakage into the low frequency part of the
spectrum. In fact, due to the Tukey window bandwidth of 0.04 Hz,
the DC component will leak into frequencies up to 0.04 Hz, where
they mix with frequencies in the 0.04 to 0.08 Hz range. In
addition, slow trends in the mean of the signal cause the signal to
exhibit nonstationary behavior, which is undesirable. It has been
found that these slow trends in the mean of clinical data contain
frequencies which are limited to approximately 0.04 Hz. In order to
minimize both of these effects, the estimated signal mf(t) and
mi(t) are digitally highpass filtered prior to undergoing
autocorrelation. Two different types of filters have been used: a
recursive three pole 1 db Chebyshev HR filter, which frequently
exhibits undesirable ringing from the DC component and has a
nonlinear phase characteristic, and a nonrecursive linear phase
FIR filter, which exhibits excellent performance and has been used
81
exclusively in the signal processing which follows. For the sixty
term FIR filter used, the DC component is down 330 db from the
passband gain.
Finally, it is occasionally useful to work not with the power
spectra M"l(f) and Mf(f) directly, but rather with their
normalized power spectral density functions S~i(f) and Sf(f),
where
S,(/)W,(/)/C,(0) 524a)
S,(/) = 524b)
Summarizing, the spectrum of the input signal m(t) of an IPFM
is estimated by first deriving the instantaneous frequency f(i,i+l)
from the modulated pulse train p(t). The instantaneous frequency is
smoothed by either linear interpolation or lowpass filtering,
generating an estimate of m(t). This estimate is then regularly
sampled at 2 Hz, and highpass filtered using a sixtyterm FIR
filter with cutoff frequency 0.04 Hz to remove the DC and very low
frequency components. The autocorrelation function of the highpass
filtered function is then calculated. The autocorrelation function
is smoothed with a Tukey lag window with bandwidth 0.04 Hz, and the
discrete Fourier transform is then applied to the smoothed
autocorrelation function, generating an estimate of either the
power spectrum or the power spectral density function.
Examples of the application of this technique to simulated
IPFM data is demonstrated in Chapter VI and to clinical HRV data
Chapter VIII.
83
CHAPTER VI
SPECTRAL ANALYSIS OF SIMULATED IPFM DATA
In order to test the validity of the analysis presented in
Chapter IV on the system properties of the IPFM and the analysis in
Chapter V on methods of performing spectral analysis utilizing the
instantaneous frequency, an IPFM was prograrmed on a digital
computer to generate spike trains from known modulating signals.
This simulation allows an input signal ro(t) consisting of up to
three sinusoids of varying frequencies, amplitudes, and phases to
modulate an IPFM, which then produces a modulated pulse train p(t)
which can be analyzed according to the algorithm described in
Section 5.4. The unmodulated pulse rate of the IPFM was set to
fo = 1 Hz, corresponding to a mean heart rate of 60 beats per
minute.
Because the power spectrum of the input signed m(t) is
userselectable and therefore known, it is possible to compare the
input power spectrum with those derived utilizing the methods of
Chapter V. The goals of these simulations are to examine the
characteristics of the derived power spectrum of IPFM data in order
to evaluate the following questions:
a) Does the simulated data support the theoretical analysis
of the IPFM derived in Chapters IV and V?
84
b) Is linear interpolation or lowpass filtering the
instantaneous frequency f(i,i+l) more accurate in
deriving the power spectrum of IPFM data in the
physiologic range?
c) To what extent do the nonlinear effects described in
Section 5.3 distort the power spectrum, and how can they
best be minimized?
d) How does aliasing affect the derived power spectrum?
For each of the simulations presented in this Chapter, the
analysis is presented in a set of graphs which demonstrates the
exact input modulating signal m(t), the estimate of the AC portion
of the input signeil reconstructed by either linear interpolation or
lowpass filtering of the instantaneous frequency as described in
Sections 5.1 and 5.2, and the power spectral density estimate of
the reconstructed signal as described in Section 5.4. Each spectral
density function is normalized so that the largest spectrad
component has a magnitude of 1.0; the absolute AC power is noted on
each graph to adlow intergraph comparisons.
The simulations included in this chapter were chosen in order
to demonstrate important properties derived in Chapters IV and V,
and involve input signals which are similar to those seen in actual
clinical data. Sixty seconds of data is anadyzed in each
simulation; while more accurate spectral estimates could be
obtained by longer sampling periods, this time period corresponds
to that used in analyzing actual clinical data, where an updated
spectral estimate needs to be available every minute.
85
6.1 Sinulations of Linear and Nonlinear Frequency Response Effects
6.1.1 A Single Modulating Sinusoid
A single modulating sinusoidal input signal consisting of
m(t)60 61)
is first analyzed. Figure 61 (following page) demonstrates the
estimate obtained using linear interpolation and that obtained by
lowpass filtering. As the frequency of the modulating sinusoid is
0.1 Hz and therefore falls within the frequency range for linear
behavior of the IPFM derived in Section 5.3.2, it is seen that the
reconstructed signal is virtually identical for both methods of
estimating the input signal, and the spectral densities are also
identical. The width of each spectral component reflects the width
of the 0.04 Hz Tukey smoothing window.
86
Figure 61
The reconstructed signals and spectral density functions derived
from a single modulating sinusoid at 0.1 Hz. The spectral densities
are virtually identical, consisting of a single frequency component
at 0.1 Hz.
87
6.1.2 Two Modulating Sinusoids
Figure 62 (following page) demonstrates an input signal
composed of two sinusoids, where
m(t) = 60{ 1 + .lsin(2/r(.06)t)+.lsin(2/r(.I2)0> 62)
As the sinusoid at f = 0.12 Hz is slightly in the nonlinear
frequency range of the IPFM, second order effects are noticeable in
the spectral density function derived from lowpass filtering;
i.e., the peaks differ in magnitude by approximately 10%. The
linearly interpolated spectral density function does not show this
nonlinear effect for this particular signal.
Of great importance, note that there are no spectral
components noted at linear combinations of the input frequencies;
i.e., at 0.06 Hz (f2 fi) or at 0.18 Hz (f2 + fl), as there are in
the spectra derived by stochastic point process theory (Section
4.5). This confirms the derivations of Sections 5.1 and 5.2, which
predict that spectral components will appear at linear combinations
of the input frequencies only when nonlinear effects become
important.
