Power spectrum analysis of human heart rate variability

Material Information

Power spectrum analysis of human heart rate variability a quantitative measure of cardiovascular function and control
Riff, Kenneth Mark
Publication Date:
Physical Description:
167 leaves : illustrations ; 28 cm


Subjects / Keywords:
Heart beat -- Testing ( lcsh )
Cardiology -- Apparatus and instruments -- Testing ( lcsh )
bibliography ( marcgt )
theses ( marcgt )
non-fiction ( marcgt )


Includes bibliographical references (leaves 161-167).
General Note:
Submitted in partial fulfillment of the requirements for the degree, Master of Science, Department of Electrical Engineering, Department of Computer Science and Engineering.
Statement of Responsibility:
by Kenneth Mark Riff.

Record Information

Source Institution:
University of Colorado Denver
Holding Location:
Auraria Library
Rights Management:
All applicable rights reserved by the source institution and holding location.
Resource Identifier:
19719722 ( OCLC )
LD1190.E54 1988m .R53 ( lcc )

Full Text
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

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

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 non-equivalent spectra, the spectrum of counts and the
interval spectrum, and methods of estimating these spectra are
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, input-output relation-
ships, and analysis of its non-linear 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 non-linear behavior of

the IPFM. Two methods of converting the instantaneous frequency
signal into a stationary continuous signal suitable for spectral
analysis, linear interpolation and low-pass filtering, are derived
and analyzed. The sources of error in the use of the instantaneous
frequency, non-linearities 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 non-invasive 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.

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
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
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

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
4.1 Description of the Integral Pulse
Frequency Modulator.................................40
4.2 Analysis of Input-Output Relationships................42
4.3 Physiologic Justification.............................46
4.4 System Properties.....................................50
4.5 The Spectrum of a Sinusoidally-Modulated Pulse Train..55
THE INSTANTANEOUS FREQUENCY.............................59
5.1 Low-Pass Filtering....................................60
5.2 Polynomial Interpolation..............................63
5.3 Sources of Error: Aliasing and Non-linear Effects.....68
5.3.1 Aliasing........................................68
5.3.2 Non-linear Effects..............................70
5.4 Spectral Analysis of the Instantaneous Frequency......76
6.1 Simulations of Linear and Non-linear
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

6.1.5 Two Sinusoids in the Second Order
Non-linear Range..............................93
6.1.6 Two Sinusoids in Higher Order Non-linear 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 Non-classical Aliasing Effects.................100
7.1 Data Acquisition System..............................102
7.1.1 Patient Interface..............................102
7.1.2 Analog Signal Pre-Processing...................103
7.1.3 Analog-to-Digital 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

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
9.1 Discussion...........................................151
9.2 Summary.............................................155
9.3 Directions for Future Research..................... 159

2-1 Electrical model of the passive
cell membrane properties...............................12
2-2 The cell pacemaker potential...............................16
2-3 Block diagram of the autonomic nervous system's
influences on heart rate variability...................18
2-4 Control model of the baroregulatory system.................21
4-1 The integral pulse frequency modulator.....................41
4-2 The input-output relationship of the integral pulse
frequency modulator and the derivation of the
instantaneous frequency........................................44
4- 3 Spectrum of a sinusoidally-modulated 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
6-2 IPFM simulation: two modulating sinusoids...............88
6-3 IPFM simulation: three modulating sinusoids of
equal phases...........................................90
6-4 IPFM simulation: three modulating sinusoids of
different phases.......................................92
6-5 IPFM simulation: two modulating sinusoids, one in the
second order non-linear range of the IPFM.............94
6-6 IPFM simulation: two modulating sinusoids, one in the
higher order non-linear range of the IPFM.............96
6-7 IPFM simulation: classical aliasing effect..............99

6-8 IPEM simulation: non-classical aliasing effect.............101
8-1 Instantaneous heart rate of a healthy subject
resting in the supine position.........................117
8-2 Filtered instantaneous heart rate of
supine resting heart rate data.........................117
8-3 Compressed spectral array of supine resting
heart rate data........................................118
8-4 Instantaneous heart rate of a healthy subject
attaining the upright posture..........................121
8-5 Filtered instantaneous heart rate of healthy
subject attaining the upright posture.................121
8-6 Compressed spectral array of healthy subject
attaining the upright posture........................122
8-7 Total AC power response of healthy subject
attaining the upright posture........................124
8-8 Mean power frequency response of healthy subject
attaining the upright posture........................124
8-9 Standard error of the mean power frequency response of
healthy subject attaining the upright posture.........125
8-10 Midband power percentage response of healthy subject
attaining the upright posture........................125
8-11 Maximum power frequency response of healthy subject
attaining the upright posture........................125
8-12 Instantaneous heart rate of subject with autonomic
dysfunction attaining the upright posture.............128
8-13 Filtered instantaneous.heart rate of subject with
autonomic dysfunction attaining the upright posture....128
8-14 Compressed spectral array of subject with autonomic
dysfunction attaining the upright posture..............129
8-15 Total AC power response of subject with autonomic
dysfunction attaining the upright posture..............131
8-16 Mean power frequency response of subject with autonomic
dysfunction attaining the upright posture..............131

8-17 Standard error of the mean power frequency response of
subject with autonomic dysfunction attaining
the upright posture...........................................132
8-18 Midband power percentage response of subject with
autonomic dysfunction attaining the upright posture...132
8-19 Maximum power frequency response of subject with
autonomic dysfunction attaining the upright posture...133
8-20 Instantaneous heart rate response to exercise...........135
8-21 Filtered instantaneous heart rate
response to exercise................................ 135
8-22 Compressed spectral array to exercise...................136
8-23 Total AC power response to exercise.....................138
8-24 Mean power frequency response to exercise...............138
8-25 Standard error of the mean power frequency
response to exercise..................................139
8-26 Midband power percentage response to exercise...........139
8-27 Instantaneous heart rate response to standing in
untreated sympathicotonic orthostatic hypotension.....142
8-28 Filtered instantaneous heart rate response
to standing in untreated sympathicotonic
orthostatic hypotension.......................................142
8-29 Compressed spectral array on standing in
untreated sympathicotonic orthostatic hypotension.....143
8-30 Total AC power response to standing in
untreated sympathicotonic orthostatic hypotension.....144
8-31 Mean power frequency response to standing in
untreated sympathicotonic orthostatic hypotension.....145
8-32 Midband power percentage response to standing in
untreated sympathicotonic orthostatic hypotension.....145
8-33 Instantaneous heart rate response to standing in
treated sympathicotonic orthostatic hypotension.......147
8-34 Filtered instantaneous heart rate response
to standing in treated sympathicotonic
orthostatic hypotension

8-35 Compressed spectral array on standing in treated
sympathicotonic orthostatic hypotension..............148
8-36 Total AC power response on standing in treated
sympathicotonic orthostatic hypotension..............149
8-37 Mean power frequency response to standing in treated
sympathicotonic orthostatic hypotension..............150
8-38 Midband power percentage response to standing in
treated sympathicotonic orthostatic hypotension.......150

When the normal human heart rate is accurately measured, it
is found to be not quite regular, but rather undergoes spontaneous
beat-to-beat 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 "R-R 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 12-15 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 (2-4). Corresponding oscillations in continuously
recorded blood pressure measurements have been termed "Traube
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 continuously-recorded 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
intra-arterial 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 R-R
intervals, or the instantaneous heart rate (IHR) derived from
these, where the IHR is defined as

IHR (beats/minute) = 60/R-R interval (seconds)
Factors affecting the decision to use measurements of period (i.e.,
R-R 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 R-R intervals,
where the SD is defined in the usual fashion:
While the SD yields a measure of the spread of the R-R
intervals around the mean R-R 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 R-R
intervals, where
N- 1

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
(Max-Min HR) during the sampling period, where
Max-Min HR IHRmax IHRmin 1-4)
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 R-R interval (Max/Min
R-R) during the sampling period, where
Max/Min R-R RRmax/RRmin 1-5)
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.
1.1.2 Power Spectrum Analysis of Heart Rate Variability
Analysis of the power spectrum of a time-varying 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
(12-14), 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 non-invasive cardiovascular monitoring in a variety of

clinical situations. Unfortunately, careful examination of these
papers reveals that most authors have not considered some basic
problems in performing heart rate variability spectral analysis
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
As an overview of this thesis, Chapters II, III, and IV are a
wide-ranging 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,

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 non-linear 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.

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
Chapter IX discusses the results demonstrated in Chapter VIII
and suggests clinical implications and some potential uses for the

2.1 Basic Klectrophysiology
Living cells are surrounded by a cell membrane: a complex
layer of lipids and proteins which act as a semi-permeable 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] 2-1)
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

The relationship between the total membrane potential Em and
the concentrations of Na+ K+, and Cl- is given by the
Goldman-Hodgkin-Katz equation (17)
ftT P,[jr]. + P,.[A;q-). + fc.[Ct-], *-
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 2-2 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.

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
resistor-capacitor network similar to that used in cable theory and
depicted in Figure 2-1, and the passive electrical properties of
the cell are known as its cable properties.

Figure 2-1
Electrical model of a cable applicable to the resting cell
membrane. See text for physiologic equivalents of electrical
The behavior of a cable such as that in Figure 2-1 is
described by the following differential equation (17):

+ E-Er-0
^ m
f ifi^ m
is known as the membrane length constant, and
is known as the membrane time constant.

If current is injected into the cable at point X=0 and time
t=0 in Figure 2-1, then the voltage Em at point X=0 at time t is
given by
E n = E r + AE 2~4>
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 steady-state distribution of potential has been
reached, the potential at any point x along the cable is given by
Em = Er + AEexp ( -x/k) 2-5)
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
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

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 all-or-none
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

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 2-2. 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).

Figure 2-2
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 excitation-contraction coupling.
Normally, each heart beat is initiated by the depolarization of a
specialized group of pacemaker cells known as the sino-atrial (SA)
node. The wave of depolarization then spreads throughout the heart

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

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 non-linear, negative-feedback control system. A
simplified block diagram of the physiological interactions which
affect heart rate variability is demonstrated in Figure 2-3.
Figure 2-3
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))

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
well-understood example of the negative-feedback control maintained
by the ANS. Blood pressure is generally modelled try a simple Ohm's
Law relationship:
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 non-linear 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.

As an example of the complexity of the ANS, other controlled
variables can also influence heart rate and blood pressure, as seen
in Figure 2-3. 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
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.15-0.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 2-4. 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 non-linear switch

Linear filter
Srltdi eleiBit
Itfamce 7 'N ^ . M %
- + -M
cure time delay
Output jf(t)
diaractert sties
of uiuulli nuscle
embroiling psipheral
vascular resistance
or afferent
Figure 2-4
Control model of the baroreceptor system which displays non-linear
characteristics such as spontaneous oscillations and frequency
entrainment. The pire time delay in the linear filter gives rise to
spontaneous oscillations, while the non-linear 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>)

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
22.0(1 + 12 j(v)e~,a>T
(1 +92;o;)(l + 2y'u>)
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 2-9 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 non-linear switch element is added to
the control loop. Therefore, the model of Figure 2-4 gives
excellent agreement with experimental results.
Thus the ANS is seen to be a complex non-linear 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

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 2-3. 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
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.

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 wide-ranging 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 short-term 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
(26-28), 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.

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
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

relevant but non-equivalent spectra of random point
processes: the interval spectrum and the spectrum of
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
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

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) 3-1)
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

process are equivalent only through their complete distributions
(equation 3-1). If a statistical analysis is based only on first
and second order properties, the implication of equation 3-1 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, 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 3-3 then becomes
var (X)
k = ..., -1, 0, 1, ... 3-2)
k = ..., -1, 0, 1, ... 3-3)
k = ..., -1, 0, 1, ... 3-4)

and it is clear that the p*'s are the Fourier coefficients of f(<).
The inverse relationship is

2 it
1 + 2 ^ ptcos(fcco)
-n < uo < n
Since p = p- for real sequences {Xi}, and, from equation 3-5, f(to) =
f(-to), a spectral density for positive to can be defined:
/.() =2/(co)
= ^/ 1 + 2 ^ ptcos(fcu;)
0 < to £ n
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 (35-37).

m,(t) = A/^r)
= lim
prob{event in t + x ,t+ x + At | event at t}
and defining m as the mean rate of the process, then
m = 1 / E {X} = lim
AT-* o
Using the definitions in equations 3-7 and 3-8, it can be
shown (15) that the covariance density of the differential process
y + (T) = lim
cov{AN ttAN t.t)
= m[m- m + 6(z)}
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 J-m
+ f [m.{r)-m)e~/a>x dx
2n 2nJy n ' 1
Again, it is convenient to define a spectral density for
non-negative to by:

g + (co) = 2 g(o>)
e',wtdr; cu £ 0
Since the stationary processes are real, mf(*) = mf(-r) and mf(t)
tends to m as t-*, An alternative form of equation 3-11, therefore,
where mj(s) is the one-sided 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.
Also define:


Anoj(V') + jB^joo)
i a.

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 non-zero covariance into the estimate.
3.1.5 Estimation of the Spectrum of Counts
0 < ou < it
0 < ao < n
co i co 2 3-18)

The spectrum of counts can be estimated by a straight-forward
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 tk is the time of the kth occurrence of an event. For equal
intervals, (i.e., I = tk tk-i), 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 low-pass filter with
cut-off 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 low-pass filtered event series (38). An efficient
algorithm for performing this convolution was described by Holden
and French for neuronal spike train data (39-42). The so-called
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

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
R-R 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 (44-47) that if the mean sampling frequency exceeds
the Nyquist rate, then the signal can, in theory, be completely

reconstructed. Interpolation schemes have been described, using
simple polynomial interpolators, that provide better performance
than a Wiener-Kolmogorov 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, R-Ri 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
where p is the order of the autoregressive process and e(k) is a
zero mean white noise process with variance i.e.,
E{e(fc)> = 0
E = o26,k
The system function H(z) is the Z transform of 3-20:

H (z)
1 +
The power density spectrum S(f) is then determined as
S(/) = G7-At | 1 +

Various methods may be used to estimate the set of
autoregressive coefficients {ai} and the white noise variance a*,
including the Yule-Walker equations and the Batch Least Squares
method. The spectral estimate is obtained by substituting these
estimates into equation 3-23. 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

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 low-pass 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 low-pass 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 non-linear effects and
aliasing, which limit the ability to perform perfect demodulation.

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.

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,53-60) 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 (65-68) whether it is
more appropriate to use heart rate (i.e., instantaneous frequency)
or heart period (i.e, R-R intervals), where the two are related as
demonstrated in equation 1-1, 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

the interbeat intervals are non-linearly 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 4-1.

// //, '/ V/
Figure 4-1
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 so-called "leaky integrator").
3) The integrator is reset instantaneously to zero; i.e.,
there is no "dead time" or refractory period.

4) The spike train output p(t) consists of a train of delta
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 Input-Output Relationships of the IFFM
Using the model of Figure 4-1, any two pulse occurrence tiroes
ti and ti+i are related by:
In the special case where the input signal does not contain a
time-varying component (i.e., m(t) = mo), then equation 4-1 reduces
= -1|)
and the unmodulated output frequency fo is given by
in which a linear relationship is preserved between the pulse train
repetition frequency fo and the input signal mo.

In the more general case, however, where m(t) contains a
time-varying component, then equation 4-1 becomes
madt +
*! I
s J mx{t)dt
Equation 4-4 can be manipulated to yield:
tt*i ~ t(
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 4-5 becomes
/(i,i+ 1 )
t(i t|
Equation 4-6 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:

/(i,i+l)= /0 + 7'nm(£,i+ 1) ^(m.0 + mm(i,i+ 1)}
If the instantaneous frequency f(i,i+l) is taken as the output
front a demodulator, then equation 4-7 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 4-2
illustrates the relationship between m(t), p(t), and f(i,i+l),
demonstrated for a sinusoidal input.
Figure 4-2
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).

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
4-2. 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 wide-sense 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 well-behaved function from f(i,i+l) for use in spectral
estimation are demonstrated in Chapter V.
Returning to equation 4-5, it can be seen that inverting both
sides of equation 4-5 yields:
1 4-8)
J o+ ?
demonstrating that the interbeat interval is non-linearly related
to mm(i,i+l). Equation 4-8 can be expanded in a Taylor series,

(i i + 1 )
m o
mnCi,i+ 1)
f f L g
If mm(i,i+l) << mo (i.e., the time-varying 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 non-linear 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 steady-state component mo, which reflects the spontaneous
4.3 Physiologic Justification of the IPFM Model

pacemaker potential, and a time-varying component mi(t), which is
integrated over time and influences the heart rate on a
beat-by-beat 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
long-term 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 2-3 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

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 phase-locking, which are not
observed in HRV studies (53); in addition, the relevance of values
of observed in isolated nerve-muscle 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 1-2 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 exponentially-decaying refractory period
(54), computer simulations have shown that the effect on impulse
generation is negligible if the refractory period is much less than

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
800-1000 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

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 non-linear 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 tk-i. 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

It is convenient to express 4-13 in terms of the sampling frequency
fo = 1/T. The frequency response becomes
a(*) I '*
For >/,, B(to)" l and nulls are found at The characteristics of the IPFM have been previously described
by equation 4-1, here modified to
* 4-15)
r = J V(t)dt
For V(t) = Vo (i.e., a constant), the output frequency fo has been
shown in equation 4-3 to be fo = Vo/r, and thus the system is
Now assume that the relevant values depart only slightly from
their steady state; i.e.,
tt ~ £*,o + f/t,i
£*-1 = ffc-1,0 +
vm-V'+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 4-15.
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 4-15 yields

0 = Vra(t*lit*-ili)+ j Vx(f)dt
The alteration in the kth time span may then be written
Tn,i t*,i ~ ifc-i.i
bringing equation 4-17 to the form
0-V0Tkil + J V^t)dt
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-
in equation 4-19 yields
T t.i

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
7 + 7,
in equation 4-21, yielding
_ J_ 4-24)
/*. i = i,i
* 0
From 4-21, 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 input-output
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.

Assuming, without loss of generality, that r = 1, 4-25 can be
manipulated to yield
fi(yo;) = eyB///sinc[/r///0] 4-26)
|B(you) |=sinc[/r///0] 4-27)
) = arg B(you)
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 4-22). In addition, there is no phase shift
between the input and output. This relationship was investigated
experimentally (57) in a frog nerve-muscle 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 4-27 but not equation 4-22. 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

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 Sinusoidally-Modulated 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)
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

is the time from the origin t=0 to the time of the first
Jn is a Bessel function of the first kind of order n.
The terms of equation 4-29 each represent a component of the
spectrum and are classified according to their frequencies of
occurrence in Figure 4-3.


"fm 0


Figure 4-3
Spectrum of a sinusoidally-modulated pulse train produced by the
The first term, Ifo, is the zero frequency or DC component
and represents the average value of the pulse train. The second

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 low-pass 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
These properties have led investigators to believe that the
modulating signal can be recovered from the modulated pulse train
by low-pass 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 Holden-French algorithm, which uses an ideal digital low-pass
filter. A hardware filter using a cosine taper has also been
described and produces virtually identical results (59).

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 non-linear 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.

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 4-2, is a staircase
function and thus suffers from non-stationarity 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

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 Low-Pass 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
is cascaded to an IPFM, then the Fourier transform of the output,
Y(f), is
= e
M{f) \ 0 sinc[/r///0] MU)
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

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 5-2 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 sine-1 [*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
low-pass 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:


0 otherwise 5-4^
(i.e., a linear phase ideal low-pass filter), the Fourier transform
of the output Y(f) is
X(/)- sinc[jr///0]tf(/) ; 0 0 ; otherwise 5-5)
which is easier to realize than that of equation 5-2, 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*'",/0-5 ; 0 0 ; otherwise 5-6)
then the output Y(f)
/(/)= sinc[/r//0.5] MU) : 0 0 ; otherwise 5-7)
satisfies the Nyquist criterion, and normalizing Y(f) by dividing
by sinc[2*f] allows error-free recovery of M(f) in the range of
0 < f < 0.5 Hz and does not require adaptive filtering. The

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 low-pass 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 low-pass 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 low-pass 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 low-pass filtered.
5.2 Polynomial Interpolation
While the previous section demonstrated that appropriate
low-pass filtering of f(i,i+l) provides an elegant method of
recovering the input signal m(t), sophisticated and time-consuming
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 piece-wise 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 time-varying low-pass filter. While this
method is computationally simple, errors occur when m(t) deviates

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 Wiener-Kolmogorov
filtering, as the simple interpolators constitute time-varying
filters, whereas the Wiener-Kolmogorov filter is constrained to be
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
piece-wise 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

g i(0 = att + bt
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
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) 5-1.
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 5-1.

Figure 5-1
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 mid-points 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.

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 5-2.
Figure 5-2
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).

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 Non-linear 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 non-linear 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

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

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 4-6, 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 Non-linear Effects
The IPFM is, in general, a highly non-linear 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

non-linear and what effects this non-linearity will have in the
attempt to reconstruct m(t) from f(i,i+l). Analysis of the
non-linear 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 non-linear
range in most physiologic situations.
If the time-varying 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)
)+ j [b j sin (gu t + ,) + 62sin(cu 2t + 02)}dt

which, after integration, becomes


It is possible to eliminate ti+i from equation 5-14, 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 5-15 into 5-14, and considerable algebraic
manipulation, then yields the relationship


Equation 5-16 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 5-16 is a non-linear 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

f(i.i+l)J f(i,i+l)

f(i,i+l)J f(i.i+l)
gj 2 ^ o>2
and equation 5-16 reduces to
+ b1sin(a>1t(+01) + 62sin(oo2tj+02)
f(i,i+l) =
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 5-17 to be accurate, then to
second order
2f(i,i + l)z
and equation 5-16 becomes

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)
m: 2m
w(a>, .cu2.0, ,02)- j£[b,sin(cu1t, + ^ J + baSinfcOj^ + ^a)]
r r
+ [b ,a>, cos(cu ,t, + ^,) + b2co2cos(a)2fl + ^2)]
x(w, +a>a,01 + *2)=^7i|cos((a>1 + cu2)t, + 0, + *J)-cos((a),-co2)t, + *1-*a)]
y(2 >?<
bl bi it \\ b? / / >\
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 5-17 are accurate to within 5% for
oo,/f(i.i+1) ranging from 0 to .314, and the second order
approximations in equations 5-19 are accurate to within 5% for
co,/f(ii+i) ranging from .314 to .84. Converting to cyclical

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 non-linear 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 5-20 and 5-21, 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, 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.

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 piece-wise 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
non-linear 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 non-linear 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

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
low-pass 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 well-known (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

non-linearities 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 N-k
ci(fc) = Z ('fti(u + fc)-ml)(/fl,(u)-/ril)
^-^J^[rhr{u + k)-m.f){rh,<;u)-rn,)
and where
_ j " 5-22c)
** ILa 1
_ i " 5-22d)
rnr = T}Lrhriu-)
IN U-l
are the sample means.
Using the symmetry of the autocovariance function, the power
spectrum estimates are

1&i(/)-C,(0) + 2 £c,(k)l/(k)cos(27ifk)
it" 1
l5l/(/) = C/(0) + 2^C/(A:)l/(fc)cos(2^/A:)

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

(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 non-stationary 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 high-pass 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
non-linear phase characteristic, and a non-recursive linear phase
FIR filter, which exhibits excellent performance and has been used

exclusively in the signal processing which follows. For the sixty
term FIR filter used, the DC component is down -330 db from the
pass-band 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),
S,(/)-W,(/)/C,(0) 5-24a)
S,(/) = 5-24b)
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 low-pass filtering,
generating an estimate of m(t). This estimate is then regularly
sampled at 2 Hz, and high-pass filtered using a sixty-term FIR
filter with cut-off frequency 0.04 Hz to remove the DC and very low
frequency components. The autocorrelation function of the high-pass
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.

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
Because the power spectrum of the input signed m(t) is
user-selectable 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?

b) Is linear interpolation or low-pass 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 non-linear 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
low-pass 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 inter-graph 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.

6.1 Sinulations of Linear and Non-linear Frequency Response Effects
6.1.1 A Single Modulating Sinusoid
A single modulating sinusoidal input signal consisting of
m(t)-60 6-1)
is first analyzed. Figure 6-1 (following page) demonstrates the
estimate obtained using linear interpolation and that obtained by
low-pass 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.

Figure 6-1
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.

6.1.2 Two Modulating Sinusoids
Figure 6-2 (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> 6-2)
As the sinusoid at f = 0.12 Hz is slightly in the non-linear
frequency range of the IPFM, second order effects are noticeable in
the spectral density function derived from low-pass filtering;
i.e., the peaks differ in magnitude by approximately 10%. The
linearly interpolated spectral density function does not show this
non-linear 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 non-linear effects become