# Filtering beat-to-beat heart rate data

From an experiment, I have a dataset of beat-to-beat heart rate data: a list of the time between each heart beat in [ms]. The data is measured using an infrared optic sensor at the finger tip. The sensor frequently misinterprets a slight movement of the finger as an heart beat. Data therefore often looks somewhat like this:

{1000, 1000, 1000, 1000, 500, 500, 1000, 1000, 1000, 600, 400, 1000}


In this example, one can easily see that the 5th and 6th element should be one; same for the 10th and 11th. However, in real life the data looks more like this:

data = {981, 870, 1099, 1105, 650, 397, 920, 917, 1015, 1085, 210, 344, 457,
950}


where the 5-6 (650, 397) and 11-12-13 (210, 344, 457) should be taken together. It is easy to just delete the incorrect data by using something like:

DeleteCases[data,
x_ /; x < Mean[data] - StandardDeviation[data] ||
x > Mean[data] + StandardDeviation[data]]


...but I want to make a function that recognizes when multiple elements should be added to one.

One could just add every two, three or four (=length) elements and select the Cases where the result lies (for example) in the range Mean[data]±StandardDeviation[data]:

length = 2;
Position[Total[data[[# ;; # + length]]] & /@
Range[Length[data] - length],
x_ /; x > Mean[data] - StandardDeviation[data] &&
x < Mean[data] + StandardDeviation[data]]


Result:

{{5}, {11}, {12}}


This gives me an idea of where the incorrect data is. Unfortunately, after having this result, I don't know what to do with it... For example, I get confused by the fact that elements 11-12-13 return 2 cases of incorrect data when I use length=1. And maybe there are more (simple) ways to filter this data.

Question: can anyone give me a kick-start?

Edit: You can download an example of actual data here. Just Flatten[Import[filename,"Table"]]

-
Perhaps fortuitously Split works on your small example: Split[data, Abs[#2 - #1] < 300 &]. However, more data is required to test how robust a simple, single threshold to Split is... –  KennyColnago Nov 14 '12 at 21:58
Any chance this experiment attempts to measure heart rate variability? Just looking into something along these lines myself. –  Jagra Nov 14 '12 at 22:09
@KennyColnago: i don't quite understand what your code does, but it works for the small example. When the mean of the data shifts to (say) 1500, i have to change the 300 to something else. I think it's difficult to implement it in a function that handles lots of different data. And i added a link to some actual data. –  A. Goossens Nov 14 '12 at 22:14
@Jagra: indeed, doing a little research on heart rate variability and training load! I learned that almost everyone in this field just deletes the incorrect data, but to me that seems not te be the way to go... –  A. Goossens Nov 14 '12 at 22:17
What if one of the subjects has arrhythmia? (Speaking from experience, the beat can occasionally get weirdly irregular.) –  Szabolcs Nov 14 '12 at 22:54

You are looking for ideas, so I will venture a partial solution in the hope it might inspire something useful. The idea presented here is to exploit a statistical model of the heartbeat intervals as a way to test the goodness of any attempted clustering of the data.

The approach is general, because it is based on maximum likelihood methods, but I will illustrate it for a particularly simple, tractable, and (possibly) applicable case. (I hold out no hope that the following model is particularly correct; I only maintain it may be a good enough approximation to afford some insight and perhaps improve solutions to the problem.)

Let us suppose that the heartbeat intervals are independent and identically distributed according to some member $f_\theta$ of a parametric family, such as a Normal distribution (where $\theta=(\mu, \sigma)$ gives the mean and standard deviation). Let us also suppose that, independent of that process, there is another parameterized discrete random variable with values in the nonnegative integers giving the number of places where which each interval is broken. Let its probability function be $p_\lambda$: presumably, its expectation is close to $0$. Finally, let us suppose that conditional on these two processes, a third distribution governs the proportions in which the intervals are broken. Let its probability function be $u_\phi$.

When $f$ is Normal, $p$ is Poisson, and $u$ is uniform, the problem becomes analytically tractable. (Such tractability is not necessary, but it speeds the execution of the program and can provide insight into what's going on.) To write the log likelihood--which we intend to maximize--suppose that the data have already been grouped into sequences corresponding to a single beat. The contribution to the log likelihood from such a group $(x_i)$, $i=1,2,\ldots,k$, is

$$\Lambda = \Bigl[-\frac{(\sum x_i - \mu)^2}{2 \sigma^2} - \frac{1}{2} \log(2 \pi \sigma^2)\Bigr] + [-\lambda + k \log(\lambda) - \log(k!)] + [\log(k!)].$$

From left to right, the terms are grouped according to the contributions from the Normal distribution, the Poisson, and the uniform. Contingent on this grouping, closed formulas for values of $\mu, \sigma$, and $\lambda$ which maximize the log likelihood can be derived (and are simple): $\hat{\mu}$ is the mean of the group sums, $\hat{\lambda}$ is $1$ greater than the ratio of number of observations to number of groups, and $\hat{\sigma}$ is the standard deviation of the group sums about $\hat{\mu}$.

To enable automatic identification of heartbeats, it may be best to represent the groups of data by means of an indicator (0-1) vector $y$: the cumulative sums of $y$ designate the groups. Thus, $y$ is $1$ at the start of each group of readings and is $0$ for intermediate readings. Equivalently, $y$ is the indicator of the actual beats. Let $\Lambda^*(y)$ designate the maximized log likelihood for the grouping given by $y$.

What we need is a way to find such a $y$ for which $\Lambda^*(y)$ is as large as possible. This is a nonlinear binary optimization program. It seems a little tricky, but there's a natural underlying structure to the problem making it feasible. Generally, many heartbeats can be clearly identified using heuristics such as those suggested in other answers to the question. All we have to do is vary $y$ slightly around the endpoints of some of the more problematic groupings.

Let's see how this might play out in practice. Begin with the sequence of $215$ observations:

data = Flatten[Import["https://dl.dropbox.com/u/1989758/beattobeatdata.txt", "Table"]];
ListPlot[data, PlotRange -> {Full, Full}, AxesOrigin -> {0, 0},
PlotStyle -> PointSize[0.0125], AxesLabel -> {"Index", "Interval"}]


A break in its histogram suggests a starting threshold for identifying overly short intervals:

Histogram[data]


Somewhere around $550$ will do. We take this as an initial indicator and, preparatory to the next step, save it:

y = ConstantArray[1, Length[data]];
y[[ Flatten[Position[data, a_ /; a < 550]]]] = 0;
x = y;


Before going on, we need to be able to maximize $\Lambda(y)$. This function returns an array given the best value of the log likelihood and the associated parameter values, $\{\Lambda^*(y), \{\{\hat{\mu}, \hat{\sigma}^2\}, \hat{\lambda}\}$.

maxLogL[data_, x_] := Module[{logL, k = Total[x], n = Length[x], \[Lambda], \[Mu], \[Sigma]2},
\[Lambda] = n/k - 1;
\[Mu] = Total[data]/k;
\[Sigma]2 = Total[(Total /@ (First /@ # & /@
GatherBy[{data, Accumulate[x]}\[Transpose], Last]) - \[Mu])^2] / k;
{f[data, {x, {\[Mu], \[Sigma]2}, \[Lambda]}], {{\[Mu], \[Sigma]2}, \[Lambda]}}
];


It calls f to compute $\Lambda^*(y)$:

f[data_, {x_, {\[Mu]_, \[Sigma]2_}, \[Lambda]_}] := Module[{logL},
logL = Function[{z}, -(1/2.) (Total[z] - \[Mu])^2 / \[Sigma]2 - (1/2.) Log[2. \[Pi] \[Sigma]2]
- \[Lambda] + (Length[z] - 1) Log[\[Lambda]]];
Total[logL /@ (First /@ # & /@  GatherBy[{data, Accumulate[x]}\[Transpose], Last])]
];


The maximum log likelihood associated with the initial value of $y$ equals $-1293.79$:

maxLogL[data, y] // N


{-1293.79, {{898.461, 12047.1}, 0.0539216}}

Let's try to improve this by changing just one entry in $y$ at a time:

For[y = x; lMax = maxLogL[data, y] // First; i = 1, i <= Length[x], i++,
y[[i]] = 1 - x[[i]];
lMax0 = First[maxLogL[data, y]];
{lMax, y[[i]]} = If[ lMax0 < lMax, {lMax, x[[i]]}, {lMax0, y[[i]]}]
];
lMax


-1172.63

That's a tremendous improvement in $\Lambda^*$ from $-1293.79$ to $-1172.63$. (Increases of $2$ or so are often considered "significant." Remember, these numbers are natural logarithms.) By inspection of the data (using plots and tables), we can do a little better still:

y[[147]] = 1; y[[148]] = 0; maxLogL[data, y] // N


{-1142.2, {{889.738, 2671.79}, 0.0436893}}

The large increase of $30.4$ in $\Lambda^*$ demonstrates this change in the clustering was worthwhile.

Let's take a look at the results by plotting the current estimate of the heartbeats represented by $y$:

ListPlot[Total /@ (First /@ # & /@
GatherBy[{data, Accumulate[y]}\[Transpose], Last]),
PlotRange -> {Full, Full}, AxesOrigin -> {0, 0},
PlotStyle -> PointSize[0.0125], AxesLabel -> {"Beat", "Interval"}]


It is wise to explore the log likelihood function. The mean of $889.738$ looks solid, so let's draw contours of the log likelihood by varying the other two parameters, $\sigma$ and $\lambda$:

ContourPlot[f[data, {y, {889.738, \[Sigma]^2}, \[Lambda]}],
{\[Lambda], 0.02, 0.08}, {\[Sigma], 46, 58}],
Epilog -> {Red, PointSize[0.015], Point[{0.0436893, Sqrt[2671.79]}]}


The contour interval is $1$. We can therefore expect that this plot shows, approximately, a $95$% contour ellipse for $(\lambda, \sigma)$. The conclusion is that:

• The mean heartbeat interval is $889.7$. It is likely between $882$ and $898$ (as indicated by similar contour plots of $\mu$ versus $\sigma$ and $\mu$ versus $\lambda$, not shown here).

• The standard deviation of the intervals is $51.7 = \sqrt{2671.79}$ and is likely between $46$ and $58$ or so.

• The rate at which beats are interrupted is $\hat{\lambda} = 0.044$ and is likely between $0.02$ and $0.08$ or so.

There are some outlying data around indexes $50$ to $70$ in the raw data, as the previous list plot shows: this is a clear violation of the parametric assumptions of these calculations. Nevertheless, it appears they have done a good job of characterizing the data overall and of guiding us to a good solution of the original problem of clustering the data and classifying them into heartbeats and non-heartbeats.

There are obvious improvements to be made. The first is that this probability model is obviously deficient: any anomalies in the data are likely to be positively correlated; the variation of heartbeat intervals is not Normal; the mean interval might change over time. Introducing parameters to handle these complications is a matter of computational complexity but creates no new conceptual or computational problems.

The second is that I have not proposed an automated way to find $y$. As far as I can tell, Mathematica has no built-in procedure to handle this kind of optimization. It should succumb easily to a genetic algorithm, simulated annealing, or perhaps even a dynamic program, but I haven't the time to do that coding.

(My code is awfully quick and dirty too, but at this stage that is of little import.)

I have offered this account because the idea of using a likelihood maximizer to help evaluate possible solutions ought to be applicable no matter what solution method is ultimately adopted. By using a probability model it will be less ad hoc than most attempts and can even provide some insight into the heartbeat process itself.

-
Wow! That's a bit more than i expected, and more than my mind can grasp... I don't understand most of your statistics, but i think i get the basic idea. Thanks! I will probably reread this a couple of times in the next days... –  A. Goossens Nov 15 '12 at 10:21

Using your values for data, this seems to work:

error = StandardDeviation[data];
data //. {a___, b_, c__, d___} /; Abs[b + c - Median[data]] < error :> {a, b + c, d}


I used StandardDeviation[data] because that's what you used, but you can put in whatever error bound you think is best there. Also note that I replaced Mean[data] with Median[data], upon reviewing Rahul N.'s comment because I remembered that Median is more representative of the center of a skewed distribution than the Mean.

Rahul N. also suggested an "error" value which is the following:

error[b_, c_] := Min[Abs[b - Median[data]], Abs[c - Median[data]]]
data //. {a___, b_, c__, d___} /; Abs[b + c - Median[data]] < error[b,c] :> {a, b + c, d}


This minimizes deviation from the median by combining terms.

I'd like to note that using a constant error value speeds up the whole process twofold (when using your linked data set). However, if speed was a concern (which you didn't mention it was), we'd have to get less elegant and avoid pattern matching.

-
Replace the condition with Abs[b + c - Median[data]] < Min[Abs[b - Median[data]], Abs[c - Median[data]]] and you don't even need to choose an error bound. –  Rahul Nov 15 '12 at 1:48
Thanks a lot! Simple method and very good results. Could you explain why for this data Median is better than Mean? I understand that the data is skewed, but that's due to the errors, am i right? 'Perfect' data (without) errors will have a normal distribution. –  A. Goossens Nov 15 '12 at 10:24
My guess is Mean will average in lowliers whereas, assuming they are not terribly frequent, Median is likely to skip over them. –  Daniel Lichtblau Nov 15 '12 at 16:07
@A.Goossens I just remembered that from stats class :). Median is approx. equal to mean for normal distributions anyway, so it's safe to go with the median. –  VF1 Nov 15 '12 at 18:03

The simplest thing would probably be to put all heart beats in a list and perform a Fourier transform:

beats = ConstantArray[0., Total[data]];
beats[[Accumulate[data]]] = 1;
beats = GaussianFilter[beats, 100, Padding -> "Cyclic"];


Which gives you an array like this (I've picked a section of your data with a few extra pulses):

ft = Fourier[GaussianFilter[beats, 100, Padding -> "Cyclic"]];
ft[[1]] = 0;
maxFrequency = First[Flatten[Position[Abs[ft], Max[Abs[ft]]]]]-1;
ListPlot[Abs[ft[[;; 300]]], PlotRange -> All, Filling -> 0]


The fourier transform looks like this, with a clear peak at 38, so there are probably 38 "true" heart beats:

You can visualize the FT result:

phase = ft[[maxFrequency + 1]]/Norm[ft[[maxFrequency + 1]]];
Show[
Plot[Im[phase*Exp[I*(k - 1)*(maxFrequency)/Length[ft]*2 \[Pi]]], {k, 1, Length[beats]}],
Graphics[{Red, Point[{#, 1} & /@ Accumulate[data]]}],
AspectRatio -> 1/10, ImageSize -> 800]


The phase is not perfect, probably because the heart rate is not perfectly constant over the whole time.

The advantage of this method is that it is very simple and very reliable if there are extra beats or missing (undetected) beats. Maybe you can use it as a first step to get a reliable estimate of the heart rate, then use that to filter extra/missing pulses.

-
Fourier transformations are a good way to grasp the main properties of an heart beat signal (a good part of the heart rate variability analysis in science is done using fourier transforms), but i think using it to filter' or estimate the data removes to many of the characteristics of the signal. Can you endorse this? –  A. Goossens Nov 15 '12 at 11:21
I was trying to do something similiar but you beat me. However, you can save a line and using a SparseArray to make the beats array: beats = SparseArray[Thread[Accumulate@data -> 1]] –  s0rce Nov 16 '12 at 14:21

A bit crude, but it works. Idea: Find lowliers, then merge those that are contiguous.

lowposns = Flatten[Position[data, aa_ /; aa < 4/5*Max[data]]];
groups = Split[lowposns, #2 - #1 == 1 &];
regroup =
SortBy[Join[
Transpose[{Complement[Range[Length[data]], Flatten[groups]]}],
groups], #[[1]] &]

(* {{1}, {2}, {3}, {4}, {5, 6}, {7}, {8}, {9}, {10}, {11, 12, 13}, {14}} *)

newdata = Map[Total, Map[data[[#]] &, regroup]]

(* {981, 870, 1099, 1105, 1047, 920, 917, 1015, 1085, 1011, 950} *)


I'm sure there are better ways to code this.

-
Good approach, now i understand how to implement my original idea. Any idea why it doesn't give good results in my actual data? Could it be the definition of the 'lowliers'? –  A. Goossens Nov 15 '12 at 10:31