Multi- and monofractal indices of
short-term heart rate variability
R. Fischer I
M. A k a y 2
P. Castiglioni 3
M. Di Rienzo 3
1Department of Biomedical Engineering, Rutgers University, Piscataway, USA
2Thayer School of Engineering, Dartmouth College, Hanover, USA
3LaRC, Centro di Bioingegneria, Fondazione Don C. Gnocchi, Milan, Italy
Abstract--Indices of heart rate variability (HRV) based on fractal signal models have
recently been shown to possess value as predictors of mortality in specific
patient populations. To develop more powerful clinical indices of HRV based on a
fractal signal model, the study investigated two HRV indices based on a monofractal
signal m o d e l called fractional Brownian motion and an index based on a multifractal
signal m o d e l called multifractional Brownian motion. The performance of the indices
was compared with an HRV index in c o m m o n clinical use. To compare the indices,
18 normal subjects were subjected to postural changes, and the indices were
compared on their ability to respond to the resulting autonomic events in HRV
recordings. The magnitude of the response to postural change (normalised by the
measurement variability) was assessed by analysis of variance and multiple comparison testing. Four HRV indices were investigated for this study: the standard
deviation of all normal R-R intervals; an HRV index c o m m o n l y used in the clinic;
detrended fluctuation analysis, an HRV index found to be the most powerful
predictor of mortality in a study of patients with depressed left ventricular function;
an HRV index developed using the m a x i m u m likelihood estimation (MLE) technique
for a monofractal signal model; and an HRV index developed for the analysis of
multifractional Brownian motion signals. The HRV index based on the MLE technique
was found to respond most strongly to the induced postural changes (95% CI). The
magnitude of its response (normalised by the measurement variability) was at least
25% greater than any of the other indices tested.
Keywords--Heart rate variability, Detrended fluctuation analysis, Fractional Brownian
motion, M a x i m u m likelihood estimation, Multifractal
\
Med. Biol. Eng. Comput., 2003, 41,543-549
1 Introduction
VARIABILITY IN the time interval between normal heartbeats
provides information about the modulation of autonomic
outflow to the heart. Mathematical metrics designed to characterise this variability have been shown to provide clinically
significant information regarding the status of autonomic cardiac
control and to be useful as a prognostic tool in a number of
clinical settings. The analysis of heart rate variability (HRV)
begins with a calculation of the time intervals between successive R-R intervals in a standard ECG record. The variability of
the intervals is then characterised, typically with some time- or
frequency-domain metric, in the time domain, for example, the
variability can be calculated with the standard deviation of
normal-to-normal (SDNN) interbeat time intervals, in the
frequency domain, the variability is often calculated as the
total spectral power within a predefined frequency band.
Correspondence should be addressed to Professor Metin Akay;
emaih metin.akay@dartmouth.edu
Paper received 17 October 2002 and in final form 8 April 2003
MBEC online number: 20033800
© IFMBE: 2003
Medical & Biological Engineering & Computing 2003, Vol. 41
J
in population studies, decreased HRV has been shown to
possess predictive value for mortality in both the general
population(TSUJI et al., 1994; 1996) and subjects with cardiovascular pathology. A relationship between decreased HRV and
mortality in post-myocardial infarction (post-Mi) patients has
been established in multiple studies (COPIE e t al., 1996;
KLEIGERet al., 1987; VAISHNAVet al., 1994). Moreover, in the
findings of the Multicenter Post-Infarction Project, HRV (using
the SDNN index) was found to have the strongest univariate
correlation with mortality of any Holter ECG variable
(KLEIGERe t al., 1987). Despite the strong predictive value of
decreased HRV, it may not have sufficient sensitivity and specificity to be used alone (TASK FORCEOF EUROPEAN SOCIETYOF
CARDIOLOGY AND NORTH AMERICAN SOCIETYOF PACING ~;
ELECTROPHYSIOLOGY,1996). in clinical work, an index of HRV
is often coupled with other risk factors to improve its clinical
utility. However, it may be possible to improve the sensitivity
and specificity of HRV assessment by developing more powerful
indices of HRV. This is a goal of the research presented here.
1.1 Characteristics of the HRV signal
The value of the HRV measurement has inspired research
efforts to associate aspects of the measurement with underlying
543
autonomic influences. Studies of the spectral content of HRV
identified two major harmonic components: one at frequencies
greater than 0.15 Hz, originally thought to be associated with the
modulation of parasympathetic influences on HRV, and a
component at or below 0.1 Hz, orignally thought to be associated
with modulation of sympathetic influences on HRV (MALLIANI
et al., 1991; SAUL et al., 1987). Broadband spectral power not
associated with these two components was typically disregarded
or believed to be measurement noise.
This interpretation of the spectral content of HRV was
reassessed however, in light of the discovery that HRV
possessed a strong 1/f component over a wide frequency
range KOBAYASHIand MUSHA, 1982). The 1/f component is
present over many orders of magnitude, from time periods of
10 s to 10 h, and accounts for approximately 85% of the spectral
power of resting human HRV (YAMAMOTO and HUGHSON,
1994). This finding has significant implications for researchers
seeking to interpret features observed in HRV frequency spectra.
First, the broad spectral aspects of the 1/fcomponent tend to bias
the magnitude of peaks in HRV spectra, especially in the lowfrequency region. Secondly, the 1/f component implies a nonstationary process, meaning that simple HRV metrics such as
SDNN will be influenced by the size of the data analysis window
(YAMAMOTOand HUGHSON, 1994).
Changes in the 1/f aspect of HRV apparently have physiological meaning as well. The LF coml~onent was found to be not
strictly l/f, but more generally l/f", where fl is the spectral
exponent, a parameter that determines the specific scaling
behaviour in the frequency domain. Estimators of the spectral
exponent were then derived, with the intention of providing
useful information about the dynamics of the autonomic control
of heart rate. Indeed, this is apparently the case; the spectral
exponent of HRV has been shown to change in various
physiological situations such as orthostatic challenge (BUTLER
et al., 1993) ageing (LIPSITZe t al., 1990) and physical exercise
(NAKAMURA et al., 1993). For example, experiments in which
subjects were tested at increasing levels of head-up tilt displayed
a progressive decrease in parasympathetic activity and a corresponding increase in spectral exponent (BUTLERe t al., 1993).
1.2 Fractal nature o f H R V
The non-stationary character of HRV and its scaling aspect
indicate that ffactal signal models may be appropriate for
describing and analysing HRV. Typically, in a fractal signal
analysis, we are interested in the calculation of certain scaling
exponents that fully characterise the scaling properties. There are
two general types of fractal signal model, monofractal and
multifractal.
1.2.1 Monofractal signals: A monofractal signal can be characterised with a single scaling exponent. One of the most
widely used monofractal signal models is the fractional
Brownian motion (FBM) model. The FBM is a statistical
signal model that includes many of the observed characteristics of HRV, including both the non-stationary aspect and
scaling behaviour in the time and frequency domains. To
characterise the variability of an FBM process, only a single
parameter is required, the global Hurst exponent H. if we
consider an HRV signal as the manifestation of an FBM
process, we can then use the parameter H as the HRV index.
indices of HRV based on the FBM model have been shown to
carry prognostic information about mortality not extractable
from traditional methods of HRV analysis in clinical studies of
at-risk populations (MiKIKALLIO et al., 1997; 1998; 1999). in
particular, an HRV index that uses a technique called
detrended fluctuation analysis (DFA) to calculate H was
544
shown to be the best univariate predictor of mortality in a
study of 159 patients with depressed left ventricular function
following acute myocardial infarction (MAKIKALLIDe t al.,
1999).
1.2.2 Multifractal signals: Multifractal signals have a highly
complex structure that requires multiple scaling exponents for
a complete description. A multifractal signal can be decomposed into many subsets, each characterised by different local
Hurst exponents h. The goal of a multifractal analysis then is
to calculate both the local exponents of the subsets and the
fractal dimension associated with the organisation of
the subsets. This information is described with a multifractal
spectrum D(h), where D(h) is the fractal dimension of the
subset of the signal having local Hurst dimension h. A
common way to perform the analysis is by the wavelet transform. In a multifractal formalism developed by ARNEODOet al.
(1991), signal singularities and the associated scaling exponents are recovered from the scaling of the wavelet transform
along the so-called modulus maxima line. The modulus
maxima line connects local maxima in the wavelet transform
across the multiple scales of the wavelet decomposition. The
scaling of the partition function Z(s, q) can be estimated from
the modulus maxima values. The partition function is defined
as the sum of the qth powers of the set of local modulus
maxima at scale s.
For small scales, the partition function has an approximate
relationship with the scaling function z(q) on the moments q
Z(s, q) ~_ s ~(q)
(1)
By linear regression on a log-log plot, the scaling function is
recovered and then used to derive the multifractal spectrum
through a Legendre transform
D(h) = qh - z(q)
(2)
For healthy subjects, it was shown (IVANOVet al., 1999) that a
multiffactal spectrum calculated in this manner on long (~6 h)
records of HRV includes a broad range of local scaling
exponents. However, for subjects with congestive heart
failure, the range of scaling exponents is narrow, which is
characteristic of a monofractal signal. These results imply that,
under healthy conditions, HRV has a multiffactal nature and that
a loss of multifractality may be associated with cardiac
pathology.
However, although this approach admits the presence of
variability in the local Hurst exponent, it is statistical in nature
and provides only a global estimate of scaling (STRUZlK,2000a).
in some applications, however, the scaling properties can be
non-stationary and can change in response to physiological
events. For these applications, we would like information
about local changes in scaling, and a single global estimate of
scaling would not be useful. Fortunately, recent research has
identified procedures that provide a stable local estimate of
scaling and are also consistent with the global multifractal
analysis approach (BENASSIe t al., 1998; STRUZlK, 2000a; b).
With this in mind, our interest lies in analysis techniques that
focus on local scaling information and permit investigation into
the temporal variation of the scaling, and yet are consistent with
a multifractal interpretation of HRV. Conventional multifractal
analysis could not be used to investigate local scaling owing to
the limited data available in local analysis data windows.
However, other techniques are available that provide a stable
estimate of the local scaling properties expected in a multifractal
signal. In one approach, a generalisation of fractional Brownian
motion called multifractional Brownian motion (MBM) has been
developed (AYACHEand VEHEL, 2000; BENASSI et al., 1998).
The major difference between FBM and MBM is that, for MBM,
Medical & Biological Engineering & Computing 2003, Vol. 41
the scaling exponent can vary with time. However, an important
attribute of the MBM model with regard to its analysis is that it
has been shown to behave locally around a given point as a
fractional Brownian motion (BENASSI et al., 1998; PELTIERand
VEHEL, 1995). Thus FBM estimators can be used for the
calculation of a local estimate of the Hurst exponent.
1.3 Research goals
In the research presented here, our goal was to identify HRV
indices that are consistent with a multifractal signal model for
HRV and that are appropriate for the study of short-term
autonomic influences on HRV. As stated previously, the use
of FBM estimators for a local estimate of scaling is not
inconsistent with a multiffactal interpretation of HRV. This is
fortunate: powerful estimators of FBM have been developed that
are appropriate for short time series (FISCHER and AKAY 1996;
PILGRAM and KAPLAN, 1998).
Fractal estimators can be classified into two categories,
non-parametric and parametric. Non-parametric estimators
typically depend on performing a linear regression fit of a
scale-dependent quantity against scale in a logarithmic representation. For real-world fractals, the scaling relationship does
not extend to arbitrarily large and small scales; consequently,
estimators using this approach must select a relevant range of
scales for analysis. Parametric estimators rely on a parametric
representation of the autocorrelation function or spectral density.
In this study, we utilised a non-parametric FBM estimator
called detrended fluctuation analysis (DFA), an estimator found
to be among the most powerful non-parametric methods in
testing on simulated data (TAQQUe t at., 1995). We also
investigated two parametric estimators. One estimator, developed using the MBM model, analyses the quadratic variations
of a signal at a given point, and the other is a maximum
likelihood estimator (MLE) for FBM. In the theory of statistical
estimation, an MLE is often regarded as the best possible
estimator, because it is asymptotically unbiased and efficient
(HAYKIN, 1991). An important goal of our research was to
determine whether the desirable properties of the MLE provided
the basis for a useful HRV index in a clinical setting.
The quantity produced by all these estimation techniques is h,
the local value of the Hurst exponent, it is this value that is used
as an index of heart rate variability. To test the effectiveness of
these estimators as indices of HRV, we used data having known
autonomic 'markers': clinical HRV recordings in which normal
subjects underwent posmral shift at a specified time. it has been
shown in previous studies that posmral tilt induces significant
autonomic change that can be detected by HRV analysis (FEI
et al., 1995). The autonomic change is thought to be the result of
vagal withdrawal and increased sympathetic outflow (VYBIRAL
et al., 1989). The HRV indices are then compared based on their
response to this marker, relative to the inherent measurement
variability.
2 Methods
2.1 H R V indices
2.1.1 Detrended fluctuation analysis estimator." The details of
the DFA algorithm employed in this study have been
described previously (PENG et al., 1995). The algorithm is a
method for characterising the self-similar (or fractal) nature of
a signal in the time domain. Briefly, in this algorithm, a time
series is divided into segments of a given length. Each
segment is 'detrended' by the subtraction of its linear trend.
The root-mean-square fluctuation of each segment is then
calculated and averaged over all segments. The average
Medical & Biological Engineering & Computing 2003, Vol. 41
fluctuation is recalculated at a succession of segment lengths
n and then plotted on a log-log scale against n. For a selfsimilar signal, this relationship will be linear, and the slope of
the plot can be used to characterise the ffactal component and
calculate a scaling exponent.
The implementation of the DFA algorithm used here is the
implementation developed by PENS et al. (1995), with one
exception. The data submitted to the algorithm were differenced
prior to analysis so that the scaling relationship was shifted by
one. The calculated scaling exponent is then the Hurst exponent.
2.1.2 Maximum likelihood estimator." The theory of maximum likelihood estimation is derived from statistical estimation theory. In this approach, we considered each interbeat
time interval to be the manifestation of a random process. The
random process was represented by a probabilistic data model;
this was the fractional Brownian motion model. As stated
previously, the FBM model is dependent on a single parameter
H. We used an MLE (described in Appendix 1) to estimate H
given the input HRV data.
2.1.3 Multifractional Brownian motion estimator." The multifractal estimator employed in this study was developed based
on a multifractional signal model called multifractional Brownian motion. The MBM model was developed independently
by BENASSI et at. (1997) and PELTIER and LEVY-VEHEL
(1995). The MBM model was developed to address a perceived shortcoming of the FBM model: the irregularity of an
FBM is the same along its entire path. For an MBM, we no
longer consider a global scaling exponent H for the entire
path; instead, we calculate a local h that can vary with time.
An estimator for MBM has been developed that calculates
quadratic variations of the signal in the neighbourhood of a
point t to produce an estimate of the scaling exponent, it is
described in detail in Appendix 2.
2.1.4 Standard deviation o f normal interbeat intervals: The
SDNN calculation was straightforwardly computed as the
sample standard deviation of the RR intervals for each data
window.
2.2 Data collection
In this study, we used posmral changes to elicit changes in
HRV that could be used as a basis for comparison of the HRV
indices. Posmral change is a commonly used method to alter the
balance between sympathetic and parasympathetic modulation
of heart rate. The test subjects were asked to lie in a supine
position, for approximately 10min in a quiet room with dim
lighting, while electrocardiographic (ECG) signals were
recorded. The subjects were then asked to sit upright as data
were recorded for another 10 min interval.
A posmral shift to an upright position is expected to elicit a
change in autonomic outflow, such that parasympathetic modulation of HR is suppressed and sympathetic modulation is
enhanced (VYBIRAL et al., 1989), increasing the value of h
measured by the HRV indices. To assess the ability of the indices
to respond to this posmral change, we used each estimator to
calculate a series of HRV estimates before and after the
autonomic event. The mean HRV estimates before and after
the autonomic event could then be calculated to characterise the
ability of the HRV index to respond to the underlying autonomic
event.
Although the magnitude of the index's response to the underlying autonomic event is important, it is also important to
consider the response in the context of the inherent variability
545
supine
upright =
the same data set, the chance of erroneously declaring a result
as statistically significant increases as we do more tests. The
Bonferroni error protection method was used here to insure that
the 'standard of proof' was made more stringent when we
performed multiple comparison testing using a single data set
(ZAR, 1996).
1.05E0.75~1.05
3 Results
rr
d- 0.75F
1.05~
0.75~samplenumber
Fig. 1 Typical HRV records for test subjects in upright and supine
conditions
of the measurement; in other words, to consider the signal-tonoise aspects of the HRV index. For a given HRV index, this was
accomplished by calculating the following dimensionless ratio
over the complete group of subjects:
mean(HR V ,pright)
SNR ----
-
mean(HR V s,pi,e )
~/var(HRV,pright) + var(HRVs,pi,e)
(3)
Data collection: Data were collected from 21 healthy, agematched subjects by ECG*, and finger arterial blood pressure
was measured by a non-invasive device t. Each recording was
digitised at 200 Hz and visually checked for artifacts. To
calculate the RR intervals, a parabolic interpolation algorithm
was used to refine the accuracy of the estimate of the time of
occurrence of each R-wave. Abnormal QRS complexes were
quite rare in the data recordings, occurring at a frequency of
approximately 0.0004%. Records with a large proportion of
artifacts were not considered; this resulted in three subjects being
removed from this study. Preceding the HRV analysis, each of
the 36 data records (18 subjects in two conditions: upright and
supine) was windowed to a 512-point record length. This was the
largest common window size that could be extracted from all
data records. Typical R-R interval data for subjects in the upright
and supine conditions are shown in Fig. 1.
Statistical analysis: For each subject, data records before and
after the postural shift were analysed with the four HRV indices.
To compute mean and variance for the estimates produced by
each index, a sequence of data analysis windows of 128-point
length were used to calculate the HRV estimates. The analysis
windows were overlapped by 50%, resulting in a total of 14
HRV estimates for each subject, seven before the postural shift
and seven after. The estimates were then used to calculate mean
and variance values for the estimates and to calculate SNR using
(3). One-way analysis of variance was then performed on the
signal-to-noise ratio (SNR) values for all subjects for the indices
to determine if there was a significant difference in the SNR
values produced by a particular HRV index. To determine which
indices had higher SNR, comparisons of the DFA, MLE and
MBM indices were performed using the SDNN index as control
at a confidence interval of 95%. in this situation, however, it is
important to adjust the error level to reflect the fact that, as
multiple comparison tests are performed on hypotheses from
*Cardioline Delta 1 Plus, Remco Italy
tFinapres, Ohmeda 2300, Ohmeda Monitoring
Englewood, CO, USA
546
The alteration in h induced by the postural change for all
subjects is shown graphically in Fig. 2. Each line segment in the
Figure represents the change in estimated HRV for a given
subject due to the postural change. Although the true value ofh is
unknown for these subjects, we can make some general observations. For the DFA and MLE indices, the value of h always
increased for the supine to upright transition, with the exception
of one subject. For the MBM index, however, the expected shift
in h was reversed for this subject as well as two others. From the
Figure, we can also compare the distribution of the h estimates
between the three indices. Relative to the MLE index, the
estimates produced by the DFA index are biased low. The
estimates produced by the MBM index are biased high,
beyond the expected range of 0.0-1.0.
For statistical testing, the SNR was calculated for the 18
subjects for the indices and then analysed with the one-way
analysis of variance. At ap-level of 0.014, the analysis found
that the SNR of the indices was not equal. Multiple contrast
testing (using the SDNN index as control) was then performed to
determine which indices had significantly greater SNR. Using
the Bonferroni error protection criterion at a confidence level of
95%, it was found that the MLE index responded more strongly
to the postural changes than the SDNN index. At this confidence
level, neither the DFA nor the MBM indices were found to differ
significantly from the SDNN index. The statistical analysis is
summarised in Table 1.
4 Discussion and conclusions
The discovery that the HRV signal has a strong fractal
component has inspired research into a variety of techniques
for characterising the ffactal aspect. HRV indices based on these
techniques have already achieved success in clinical applications. The success of these methods may lie in the fact that they
examine more than just the magnitude of HRV variability (they
characterise the correlation features of HRV at multiple time
scales) and perhaps are a better way to characterise the status of
the physiological control systems that are responsible for HRV.
Studies of relatively long records of HRV ( ~ 6 h ) have
demonstrated the multifractal nature of HRV. A complete
//~
0.8-
1.0-
0.60.40.20
Z//~
S
supine upright supineupright
supineupright
Fig. 2 Changes in Hurst exponent H induced by postural shift from
Systems,
supine to upright condition, as estimated by (0) DFA, (•)
MLE and (EB)MBM HRV indices
Medical & Biological Engineering & Computing 2003, Vol. 41
Table 1 One-way ANOVA analysis o f HRVestimator response to postural change. Analysis finds that
SNRs of various estimators" are not all equal (p = 0.014) and that only SNR for MLE estimator is
significantly larger than SDNN estimator
SNR by HRVmetric
MLE
DFA
SDNN
MBM
n
18
18
18
18
Source of variation
HRV estimator
Within cells
Total
SSq
256
155
181
Contrast
MLE against SDNN
DFA against SDNN
MBM against SDNN
Difference
1.60
1.07
0.513
Mean
2.52
1.99
0.917
1.43
SD
2.35
1.53
0.716
0.848
SE
0.555
0.361
0.169
0.200
DF
3
68
71
MSq
8.67
2.28
F
3.80
P
0.0140
Bonferroni 95% CI
0.37-2.84
(significant)
0.16-2.31
0.72-1.75
description of HRV at this time scale requires a spectrum of
scaling exponents. This spectrum reflects the fact that there is
variability among the local scaling exponents and it describes the
distribution of these exponents. In the research presented here,
we have addressed the scaling behaviour of HRV on a much
shorter time scale, in the order of a few minutes. Here, we are
concerned with tracking changes in the local scaling exponent
instead of describing the global distribution of scaling exponents. Our results show that an index based upon a monofractal
signal model is best at detecting short-term alterations to
autonomic balance. These results encourage further development of HRV indices based upon mono fractal signal models and
highlight the power of such indices for the analysis of physiological events in HRV.
Another outstanding question regards the accuracy o f a fractal
model for HRV. Although effective indices of HRV can be
derived from fractal models, it is by no means certain that a
fractal model gives a complete description of the HRV signal.
Certainly, HRV has other signal components, such as measurement noise, and harmonic components, such as respiratory sinus
arrhythmia. Signal models that include these components may
be the basis for more effective indices of HRV.
Mathematically, to find the ML estimate, we need to express
the probability density of the model as a function of the
parameters (the likelihood function) and then attempt to find
values of the parameters that maximise this function. By utilising
the increments of FBM, we create a process that is zero-mean,
normally distributed, strict-sense stationary and more amenable
to analysis. The probability density function for the increments
of FBM is apparently unimodal (LUNDAHL et al., 1986), and so,
by taking the logarithm (a monotonic function) of the likelihood,
we simplify our calculation without changing the location of the
maximum. From (2),
N
1
logp(x; H) = - ~ - l o g 2 n - ~loglRI
-
~xTR 1X
(4)
with the elements of R available from the autocorrelation
function as in (3).
The two parameters to be estimated are H and o-2. if we
decompose R as
R
= o-2R'
(5)
the likelihood function can be written as
N
Appendix 1
N
1
logp(x; H ) = - ~ - l o g 2 n - ~-log(o-2) - ~loglR '1
Maximum likelihood estimator f o r H
We sought to derive an algorithm (the estimator) that would
provide an estimate of H g i v e n the input HRV data. Further, we
required that the algorithm have certain desirable statistical
properties: estimates produced by an MLE can, under certain
circumstances, be unbiased and have minimum variance.
The derivation of the MLE starts with an expression for the
autocorrelation function for the increments of FBM (LUNDAHL
et al., 1986)
r[k]
v.
= -~- Llk + 112H -- 2lkl 2H -[- Ik -- 112H]
(1)
With the autocorrelation known, the probability density function
for the increments is straightforward (KAY, 1993)
1
{ ~
p(x; H) = (2n)N/21RI1/2 exp --
}
xrR
1x
(2)
1 xrR,_l x
(6)
2o-2
Equation (6) can be minimised over o-2 by taking the derivative
with respect to o-2 and letting the derivative go to zero
82
xTR ~ 1x
-- - -
(7)
N
inserting this in (6) leaves us with the final expression of the
likelihood that we will maximise over H
N
logp(x; H ) = - ~ - l o g 2 n - ~ -
1
2 loglR'l
N 1
[xTR~ Ix~
og~
~-
-)
N
(8)
where x is the dataset, R is the covariance matrix
e i j = r [ l i - jl]
(3)
and r[k] is the autocorrelation function given in (1).
Medical & Biological Engineering & Computing 2003, Vol. 41
To maximise this function over H, we need only search in the
interval [0, 1]. For the implementation of the MLE algorithm
used here, lower-upper (LU) decomposition of R was used to
547
find the inverse and determinant o f R, and Brent's algorithm was
used to maximise the likelihood (PRESS et al., 1992).
Appendix 2
Multifractional Brownian motion estimator
The estimator calculates quadratic variations of the signal in
the neighbourhood o f a point t to produce an estimate o f the
scaling exponent. The signal X is observed at sampling points
p / N , p = O, K, N. For any point o f time t during the observation
o f the signal, a neighbourhood of t is defined as
v~,N(t ) = { P - - t
_< ~}
(1)
w h e r e p = 0 , . . . ,N, tE (0, 1), e > 0 , N > 0 .
Local variations o f the signal are then calculated in the
neighbourhood as the sum o f discrete second derivatives o f
the signal
Ve,N(t)= ~
(x(P+l,~_2x(P)
+ x(p~_))2
(2)
This calculation is repeated again on a rescaled version o f the
signal for which we consider only every other point, V~,N/Z (t).
The local estimate o f the scaling exponent is then obtained as
h(t)
=
1 (log2 Ve'N/2(t) 1)
-~
~,N(t) ~-
(3)
Acknowledgem ents
This work was supported by NIH grant (HL 65732).
References
ARNEODO, A., BACRY,E., and MuzY, J. (1991): 'Wavelets and multifractal formalism for singular signals: Application to turbulence
data', Phys. Rev. Lett., 67, pp. 3515-3518
AYACHE, A., and VEHEL, J. (2000): 'The generalized multifractional
Brownian motion', Stat. Infe~ Stochast. Process., 3, pp. 7-18
BENASSI, A., COHEN, S., and ISTAS, J. (1998): 'Identifying the multifractional function of a Gaussian process', Statist. Prob. Lett., 39,
pp. 337-345
BENASSI, A., JAFFARD,S., and Roux, D. (1997): 'Gaussian processes
and pseudodifferential elliptic operators', Revista Mathernatica
iberoamericana, 13, pp. 19-89
BUTLER, G., YAMAMOTO, Y., XING, H., NORTHEY, D., and
HUGHSON, R. (1993): 'Heart rate variability and fractal dimension
during orthostatic challenges', J. Appl. Physiol., 75, pp. 2602-2612
CoPIg, X., HNATKOVA,K., STAUNTON,A., FEI, L., CAMM, A., and
MALIK, M. (1996): 'Predictive power of increased heart rate versus
depressed left ventriculax ejection fraction and heart rate variability
for risk stratification after myocardial infarction. Results of a twoyear follow-up study', J Am. Coll. Cardiol., 27, pp. 270-276
FISCHER, R., and AKAY, M. (1996): 'A comparison of analytical
methods for the study of fractional Brownian motion', Ann.
Biomed. Eng., 24, pp. 537-543
HAYKIN, S. (1991): 'Adaptive filter theory, 2nd edn' (Prentice Hall,
Englewood Cliffs, 1991)
IVANOV, E, AMARAL, L., GOLDBERGER, A., HAVLIN, S.,
ROSENBLUM, M., STRUZIK, Z., et al. (1999): 'Multifractality in
human heartbeat dynamics', Nature, 399, pp. 461-465
KAY, S. M. (1993): 'Fundamentals of statistical signal processing,
vol. h estimation theory' (Prentice-Hall, Upper Saddle River, 1993)
KLEIGER, R. E., MILLER, J. R, BIGGER, J. T., MOSS, A. J., and
MULTICENTER POST-INFARCTION RESEARCH GROUP (1987):
548
'Decreased heart rate variability and its association with increased
mortality after acute myocardial infarction', Am. J. Cardiol., 59,
pp. 256-262
KOBAYASHI,M., and MUSHA, T. (1982): ' l / f fluctuation of heaxtbeat
period', IEEE Trans. Biomed. Eng., 29, pp. 456-457
LIPSITZ, L. A., MIETUS, J., MOODY, G. B., and GOLDBERGER,A. L.
(1990): 'Spectral characteristics of heart rate variability before and
during postural tilt. Relations to aging and risk of syncope',
Circulation, 81, pp. 1803-1810
FEI, L., ANDERSON, M. H., STATTERS, B. M., MALIK, M., and
CAMM, A. J. (1995): 'Effects of passive tilt and submaximal
exercise on spectral heart rate variability in ventriculax fibrillation
patients without significant structural heart disease', Am. Heart
J, 129, pp. 285-289
LUNDAHL, T., OHLEY, W. J., KAY, S. M., and SIFFERT, R. (1986):
'Fractional Brownian motion: A maximum likelihood estimator and
its application to image texture', IEEE Trans. Med. Imag., 5, pp.
152-161
MALLIANI, A., PAGANI, E, LOMBARDI,E, and CERUTTI, S. (1991):
'Cardiovascular neural regulation explored in the frequency
domain', Circulation, 84, pp. 482-492
MA.KIKALLIO, T. H., SEPPA.NEN, T., AIRAKSINEN, K. E. J.,
KOISTINEN, J., TULPPO,M. P., PENG, C. K., et al. (1997): 'Dynamic
analysis of heart rate may predict subsequent ventricular tachycaxdia after myocaxdial infarction',Am. J. Cardiol., 80, pp. 779-783
M~.KIKALLIO,T. H., RISTIMA.E,T., AIRAKSINEN,K. E. J., PENG, C. K.,
GOLDBERGER, A. L., and HUIKURI, H. V (1998): 'Heart
rate dynamics in patients with stable angina pectoris and utility
of fractal and complexity measures', Am. J. Cardiol., 81,
pp. 27-31
MA.KIKALLIO,T. H., HOIBER, S., KOBER, L., TORP-PEDERSEN, C.,
PENG, C. K., GOLDBERGER,A. L., et al. (1999): 'Fractal analysis of
heart rate dynamics as a predictor of mortality in patients with
depressed left ventriculax function after acute myocardial infarction', Am. J. Cardiol., 83, pp. 836-839
NAKAMURA,Y., YAMAMOTO,Y., and MURAOKA,I. (1993): 'Autonomic
control of heart rate during physical exercise and fractal dimension
of heaxt rate variability', J. Appl. Physiol., 74, pp. 875-881
PELTIER, R., and VEHEL, J. L. (1995): 'Multifractional Brownian
motion: Definition and preliminary results'. Research Report
2645, INRIA National Institute for Research in Computer Science
& Control, France
PENG, C. K., HAVLIN,S., STANLEY,H. E., and GOLDBERGER,A. L.
(1995): 'Quantification of scaling exponents and crossover phenomena in nonstationary heartbeat time series', Chaos, 5, pp. 82-87
PILGRAM, B., and KAPLAN,D. T. (1998): 'A comparison of estimators
for 1/f noise', Physica D, 114, pp. 108-122
PRESS, W. H., TEUKOLSKY, S. A., VETTERLING, W. T., and
FLANNERY, B. E (1992): 'Numerical recipes in C, 2nd edn'
(Cambridge University Press, New York, 1992)
SAUL, J. R, ALBRECHT,E, BERGER, R. D., and COHEN, R. J. (1987):
'Analysis of long term heart rate variability: methods, 1/f scaling
and implications', Cornput. Cardiol., 14, pp. 419-422
STRUZlK, Z. R. (2000a): 'Determining local singularity strengths
and their spectra with the wavelet transform', Fractals, 8, pp.
163-179
STRUZlK, Z. R. (2000b): 'Revealing local variability properties of
human heartbeat intervals with the local effective Holder exponent'.
CWI Technical Report INS-R0015, National Research Institute for
Mathematics & Computer Science, The Netherlands
TAQQU, M. S., TEVEROVSKY, V, and WILLINGER, W. (1995):
'Estimators for long-range dependence: am emprical study',
Fractals, 3, pp. 785-798
TASK FORCE OF EUROPEAN SOCIETYOF CARDIOLOGYAND NORTH
AMERICAN SOCIETYOF PACING & ELECTROPHYSIOLOGY(1996):
'Heart rate variability. Standards of measurement, physiological
interpretation, and clinical use', Circulation, 93, pp. 1043-1065
TsuJI, H., LAWRENCE, M. G., VENDITTI, E J., MANDERS, E. S.,
EVANS, J. C., FELDMAN,C. L., et al. (1994): 'Reduced heart rate
variability and mortality risk in am elderly cohort. The Fraxningham
Heart Study', Circulation, 90, pp. 878-883
TsuJI, H., LAWRENCE, M. G., VENDITTI, E J., MANDERS, E. S.,
EVANS, J. C., FELDMAN, C. L., et al. (1996): 'Impact of reduced
heart rate variability on risk for cardiac events. The Fraxningham
Heart Study', Circulation, 94, pp. 2850-2855
Medical & Biological Engineering & Computing 2003, Vol. 41
VAISHNAV, S., STEVENSON, R., MARCHANT, B., LAGI, K.,
RANJADAYALAN,K., and TIMMIS, A. (1994): 'Relation between
heart rate variability early after acute myocardial infarction and
long-term mortality', Am. J. Cardiol., 73, pp. 653-657
VYBIRAL, T., BRYG, R. J., MADDENS,M. E., and BODEN, W. E. (1989):
'Effect of passive tilt on sympathetic and parasympathetic components of heart rate variability in normal subjects', Am. J. CardioL,
63, pp. 1117-1120
YAMAMOTO,Y., and HUGHSON, R. L. (1994): 'On the fractal nature of
heart rate variability in humans: effects of data length and
fl-adrenergic blockade', Am. d Physiol., 35, pp. R40-49
ZAR, J. H. (1996). 'Biostatistical analysis, 3rd edn' (Prentice Hall,
Upper Saddle River, 1996)
Medical & Biological Engineering & Computing 2003, Vol. 41
Author's biography
METIN AKAY is Associate Professor of Engineering, Psycology and
Brain Sciences, and Computer Science at Dartmouth College,
Hanover. He is author/coauthor of 14 books and is the founding
and current editor-chief of the Biomedical Engineering Book Series
published by the Wiley and IEEE Press. Dr. Akay's Neural Engineering and Informatics Lab is interested in investigating the respiratory
somatosensory (RSS) responses of patients with obstructive sleep
apnea syndrome (OSAS) and the effect of developmental abnormalities and maturation on the dynamics of respiration. He teaches
courses in the areas of Biomedical Informatics, Signal Processing,
Neural Networks and Computations.
549