Introduction
This page provides software and documentation related to our point process algorithms
for estimation of instantaneous measures of Heart Rate Variability
(HRV). Only a limited set of algorithms are presented here. The programs are fully functional, although they are primarily aimed at demonstrating the efficacy of the models, and guiding the user along the kind of results that our frameworks can achieve. Programs more tailored for research and statistical analysis use require a level of assessment of the procedures, as well as fine parameter tuning, which go beyond the capabilities of a simple automatic set of routines. You are welcome to
contact us for assistance on any questions you may have on the following programs and tutorials, as well as if interested in pursuing research projects with our point process models and algorithms.
Tutorial and Downloads
This tutorial goes over several steps to guide the user through Heart Rate Variability (HRV)
analysis by means of Point Process Models. It shows the use of the
software tools collectively called PPHRV (Point Process Heart Rate
Variability) provided in this website. Please always read LICENSE.TXT for terms of use.
Preliminary steps
Required software.
The software tools used
within this tutorial can be run either with
GNU
Octave or
Maltlab,
under GNU Linux or
Windows.
Specifically, they have been
tested with Octave-3.2.4-mingw and Matlab 2010a under Windows Vista
32bit and
octave-3.2.4 and Matlab 2008 under Ubuntu 10.10 32bit.
Installation of Octave.
For those who have
neither Octave nor Matlab, the former is a free software (free as in
beer and as in speech) and the Windows version can be downloaded from
here.
During the installation, when prompted for which components to
install, it is sufficient to accept the default options. After
installation, run
pkg rebuild -noauto oct2mat
at the octave prompt. Note that in some Windows versions the right-button cut-and-paste sequence commonly used to copy and run the instruction lines may not work. In such cases it is sufficient to right-click the title bar and choose Edit _ Paste.
Setup of the compiler for Matlab.
Those who have Matlab installed, should setup the Matlab C compiler by
following the following easy steps. First issue
mex -setup
at the Matlab prompt, then let mex locate the installed compilers
(press 'y') and accept the Lcc compiler (usually the first selection,
type '1'), finally confirm the proposed location of the compiler (press
'y').
Installation of the
PPHRV software tools.
Some
of the tools presented here are simple ".m" scripts. Other tools are
made of several ".m" functions and, possibly, some compiled shared
libraries (.dll or .so). For simplicity, we assume that all the files
and the datasets described below have been saved in the same working
folder, for example "C:\pphrv" under Windows or "~/pphrv" under Linux.
This way, in order to make them available to the interpreter, it
is sufficient to instruct
the interpreter (Octave or Matlab) to switch to that folder by issuing,
for example,
cd c:\pphrv
at the beginning of each session.
Datasets
Point
processes can be used as models for random events in time. A point
process can be described completely by the (random) time of occurrence
of the events. In the case of the analysis of heart rate and
heart
rate variability, the events of interest are the R waves of the
electrocardiogram. Therefore the input of our algorithms is the series
of the time of occurrence of the R waves. A point process analysis
cannot be conducted starting from the resampled tachogram because the
information about the time of occurrence of the beats is completely
lost. Instead, if the tachogram is available, the series of R events
can be reconstructed issuing
R = cumsum(RR);
where
RR
is the tachogram, i.e., the series of heart beat intervals. If only the
raw ECG is available, the series of R-events can be obtained through a
peak identification software (out ofthe scope of this tutorial).
Examples of R-event series can be found in
this
section of the
Physionet
archive.
In this tutorial we will make use of recordings from the
Meditation database.
RR series
1) Download the
Y2.pre
(right-click and select "save link as") file from Physionet and save it
in the working folder.
2) Load the
file and plot
the RR series:
data = load('Y2.pre');
R = data(:,1); % series of times of R-events [s]
figure
plot(R(2:end), 1000*diff(R)) % [ms]
3) Then, by visual inspection, select a segment that looks
approximately stationary, for example the heart beats with indices 75
to 300.
R = data(75:300,1); % subset of times of R-events [s]
figure
plot(R(2:end), 1000*diff(R)) % [ms]
xlabel('time [s]')
ylabel('RR [ms]')
to obtain a plot like the following.
History-dependent Inverse Gaussian regression
This section shows how to use the
regr_likel
tool to find
the parameters of a history-dependent inverse Gaussian (IG)
distribution by maximizing the likelihood. With the default options,
regr_likel
models
the mean of the IG as a regression of the last 8 observed RR intervals
plus a constant term. These parameters can be changed.
1) Download
regr_likel.zip
and extract it in your working folder.
2) Enter the following instructions to perform the regression:
[Thetap,Kappa,opt] = regr_likel(R);
3) If you want to see the set of parameters maximizing the likelihood,
enter
Thetap',Kappa,opt.Theta0
. For the segment
considered above you should see
Thetap=[0.770, ...,
-0.308]',
Kappa=1655.3, and
opt.Theta0=0.52257.
The variable
Thetap is a vector with the weight of the
previous RR intervals in the regression of the first moment (
µ)
of the IG while
Kappa is the shape parameter of the IG (it
is often called
λ but we use
Kappa to avoid
confusion with the lambda function that we will introduce later).
Spectral analysis
1) Download the
spectral.m
routine from
here
and save it in the working folder.
2) The power spectral density of the model found can plotted using the
following commands:
meanRR = mean(diff(R)); % average RR interval
Var = meanRR^3 / Kappa; % variance of an inverse Gaussian
Var = 1e6 * Var; % from [s^2] to [ms^2]
spectral(Thetap, Var, 1/meanRR);
% with Thetap: coeffs of IG regression; Var: variance;
1/meanRR: sampling freq
xlabel('f [Hz]')
ylabel('PSD [ms^2/Hz]')
In this simple form,
spectral.m
plots the power spectral
density of the RR series (black lines) and the components associated to
each real pole or pair of conjugate poles (coloured lines).

With more input and output arguments,
spectral
can give
information about the poles, their nominal frequency, their power, and
also plot the power spectral density of each individual pole in the
full range of frequency between -fsamp/2 and fsamp/2:
% [tot,comps,compsp,f,pole_freq,pole_pow,pole_res,poles] =
spectral(A, Var, fsamp, f, do_comps)
% with:
% A is a vector of length P with the regression terms of
the IG model (Thetap)
% Var is the variance of the IG
% fsamp is the sampling frequency
% f is the number of frequency points to use in [0 fsamp/2]
if f positive or
% [-fsamp/2
fsamp/2] if f negative
% do_comps can be
% 0: do not find
spectral components but only the total PSD
% 1: find
spectral components both for each individual pole (output
%
"comps") and for pair of conjugate poles (output "compsp")
% 2: as in the
case 1 and also plot "compsp" (default)
% 3: as in the
case 1 and also plot "comps"
% the output arguments are:
% tot is a vector with the total PSD (plotted in black)
% comps is a matrix with P rows, each row is the PSD of a
single pole
% (plotted in
colour with do_comps=3)
% compsp is a matrix with <=P rows, each row is the PSD
of a single real pole
% or of a pair of
conjugate poles (plotted in colour with do_comps=2)
% f is a vector with the frequency points where the PSD is
evaluated
% pole_freq is a vector with the nominal frequency of each
pole
% pole_pow is a vector with the power of each pole
% pole_res is a vector with the residual of each pole
% poles is a vector with the position of each pole in the
complex plane
3) Run
[tot,comps,compsp,f,pole_freq,pole_pow,pole_res,poles] =
spectral(Thetap, Var, 1/meanRR, -2048, 3);
to obtain the following plot:

and the following information:
pole_freq',pole_pow',pole_res,poles
% pole_freq' = 0.0689 -0.0689 0.118 -0.118 0.308 -0.308 0.466 -0.466
% pole_pow' = 239 239 303 303 174 174 20.5 20.5
% pole_res =
% 0.493531 + 0.372013i
% ....
% poles =
% 0.71300 + 0.30299i
% ....
4) To check that the total power in the range
[-fsamp/2,fsamp/2] is exactly the sum of the power of all the poles,
enter:
sum(tot)*(f(2)-f(1))
% ans = 1473.3 % [ms^2]
sum(pole_pow)
% ans = 1473.3 % [ms^2]
Point-process assessment of time-varying HRV indices through local
likelihood of IG distribution
In the example above, the RR series was assumed to be stationary.
It is possible to use a point-process approach to study a
non-stationary signal and track the time-varying parameters that, at
each moment in time, define the probability distribution of the time of
arrival of the next heart beat as a regression of the last observed RR
intervals.
In order to estimate the time-varying parameters of the inverse
Gaussian model download
pp_likel.zip,
then run the following:
R = data(:,1); % series of times of R-events [s]
[Thetap,Mu,Kappa,L,opt] = pplikel(R);
t = opt.t0 + (0:length(Mu)-1) * opt.delta;
Var = opt.meanRR.^3 ./ Kappa; % variance of an inverse Gaussian
Var = 1e6 * Var; % from [s^2] to [ms^2]
figure; hold on
plot(R(2:end), 1000*diff(R), 'r*')
plot(t, 1000*Mu)
legend('RR', 'First moment of IG distribution')
xlabel('time [s]')
ylabel('[ms]')
This plots the observed RR intervals. It also plots, at time
t,
the first moment (mean) of the IG distribution of the time of arrival
of the next R event, based on the history of RR intervals up to time
t.

Note that the first 90 seconds of the data do not have an estimate
because a window is used to initialize the estimation procedure.
The function
pplikel
can be invoked with additional
arguments
in order to change the settings used in the estimation.
For example
pplikel(R, 'delta', 0.001, 'P', 12, 'weight', 0.99)
uses a time resolution of 1ms, a regressive order 12, and a forgetting
factor 0.99
for the local likelihood.
The returned variable
Thetap is a matrix that, for each
sampling point (at time resolution
delta, default 5ms),
gives the
P weights of the regression of the first moment
of the IG from the last
P observed RR intervals. Similarly,
Mu and
Kappa are vectors with, respectively, the
first moment and the shape parameter of the IG at each sampling point.
All the information needed to compute the spectra along time are in the
matrix
Thetap, and the vectors
Mu
and
Kappa.
The spectrum can be easily compared at different time points; for
example to compare the spectra at time t1=900s and t2=1220s:
1) obtain the sample indices
corresponding to the two time points that you want
to compare
i1 = round((900 - opt.t0) / opt.delta);
i2 = round((1220 - opt.t0) / opt.delta);
2 ) use
spectral
again to plot the spectra:
spectral(Thetap(:,i1),
Var
(i1), 1/opt.meanRR(i1));
xlabel('f [Hz]'); ylabel('PSD [ms^2/Hz]')
spectral(Thetap(:,i2), Var(i2), 1/opt.meanRR(i2));
xlabel('f [Hz]'); ylabel('PSD [ms^2/Hz]')

From each spectrum, the spectral indices can be computed at each time
point.
Using the tool
hrv_indices
, the time-varying power of LF,
power of HF and LF/HF index (reflecting the sympathovagal balance) can
be estimated for the entire segment considered:
[powLF, powHF, bal] = hrv_indices(Thetap, Var
,
1./opt.meanRR);
and then plotted:
figure
plot(t, powLF)
xlabel('t [s]'); ylabel('powLF');
figure
plot(t, powHF)
xlabel('t [s]'); ylabel('powHF');
figure
plot(t, bal)
xlabel('t [s]'); ylabel('LF/HF'); ylim([-1, 20])
Goodness-of-fit assessment
The key point of a point process approach to the study of heart rate
variability is that goodness-of-fit methods based on the theory of
point processes can be used to
assess the agreement between human
heart beat series and model estimates of these series derived from the
algorithms. Using the time-rescaling theorem, the goodness-of-fit
between human
heart beat series and model estimates can be assessed by KS-plot
analysis.
To run the goodness-of-fit analysis you must have run the
pplikel
and used the computed hazard-rate function lambda (
L) as
input to the following steps.
1) Download the tool
goodness_of_fit.zip.
2) Plot the previously computed lambda function with the following
instructions:
figure
plot(t, L)
xlabel('t [s]'); ylabel('Lambda (hazard-rate) function');

3) Compute the KS-plot using
ks_plot.m
and running the
following line:
[KSdistance,Z] = ks_plot(R, L, opt.delta)

4) Compute the autocorrelation of the residuals using
check_corr.m
and running the following line:
[corr,threshold] = check_corr(Z);
Coming soon
Extended multivariate point process models.
Point process nonlinear analysis.
Ancillary tools.
Analysis of specific data sets.
Active message board.
Useful links related to cardiac
monitoring.