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.