methods:hrv:hrv_stochastic_processes_and_resampling

*When measuring beat to beat intervals of the heart (the inverse of the heart rate), these samples are detected at intervals that are not equal in time. The series of samples is inherently asynchronous. We should consider whether this makes a big error or not. This page describes some background theory on how to process asynchronous signals like heart beat intervals.*

More details on the mathematical background of signal processing of stochastic signals can be found on the web at Wolfram MathWorld^{1)}.

The digital signal coming directly from a pulse detector (for example implemented in an Arduino or MX4 microcontroller) has the strange phenomenon that samples do not come in at equidistant periods: the signal is asynchronous. It is sketched in figure 1 as a series of samples $x_{i}$. Between every two subsequent data samples $x_{i}$ and $x_{i-1}$ there is a square surface because the height of the sample is equal to the time-spacing with the previous sample.

The signal is characterized by:

- Samples $x_{i}$ are discrete in amplitude: these are unsigned-long (4 bytes = 32 bits) numbers representing the time intervals in µs. Because the actual time resolution as specified by Arduino is $4 \mu s$, the two lower bits are not significant, meaning we have 30 bits effectively in the amplitudes of $x_{i}$.
- The lengths of the time intervals $\left [t_{i} - t_{i-1} \right ]$ are equal to the numbers represented in the samples, so we can say

\begin{equation} \left [ t_{i}-t_{i-1} \right ]=x_{i} \label{eq:TimeSeries} \end{equation}

However, the time intervals are transmitted over the USB/serial line with a very time-unreliable protocol. There is not enough significance in the timing information, and we should take the amplitude information as the source for signal processing.

Over a certain time period, we may assume the asynchronous signal to be “stationary”. This means that the characteristics (average, standard deviation, etc.) at $t_{1}$ are identical to the characteristics at $t_{N}$ for every $N$. The signal becomes non-stationary, for example, when we change the level of physical exercise or arousal during a single measurement. In that case, we have to work with a running average or an alternative representation technique.

**When we have a stationary situation, and the variation on the intervals $ \left [t_{i} - t_{i-1} \right ]$ is small with respect to the absolute interval, we can consider the data sequence $x$ as a quasi-synchronous signal. This means the mathematical average of the samples x is equal to the average heart-rate and equal to the nominal sample interval.**

Realizing this, we can think of the signal as a synchronous signal (meaning fixed periods between the samples), but with the real intervals in the amplitudes $x_{i}$. In mathematical writing:

\begin{equation} \left [ t_{i}'-t_{i-1}' \right ]=\mu_{1} \label{eq:IntervalIsAmplitude} \end{equation}

with $ \mu _{1}$ (first moment expectation value) the average of the intervals $x_{i}$ and $t_{i}'$ the equivalent of the sample times $t_{i}$, mapped onto an equidistant time grid. Two options are possible as shown in figure 2.

The Arduino heart pulse interval samples are the black sample peaks. They appear at asynchronous intervals. In the top graph, the samples are mapped onto an equidistant (= synchronous) sample peak series (red) with a spacing equal to the average heart rate $ \mu _{0}$. In the lower graph, we have converted to a synchronous sample peak series as well, but now with a subsampled higher data rate with arbitrary spacing. In the second case we assume the heart rate was equal to $x_{i}$ during the whole previous interval $t_{i-1}$ to $t_{i}.

The first option with $ \mu _{1}$ as the synchronous time base works well for local heart rate averages, standard deviations and even FFT’s (spectral analysis). The second solution becomes the only option when comparing two heart-rates (cross correlations). Satisfying results have been observed with re-sample intervals of $100ms$ up to $500ms$. This is a good compromise between numerical reliability and the increase of data points.

Characteristics in the time domain, can be described in terms of “moments” of the “expectation value”. The expectation value of the $k$-th moment is defined as

\begin{equation} \mu_{k}\equiv \left \langle x^{k} \right \rangle = \int_{-\infty }^{\infty} x^{k}f\left ( x \right )dx \end{equation}

for continuous functions. In this definition, $f(x)$ is the probability function of the “process” or “signal” $x$. For a finite discrete realisation (“ensemble”) $x_{1}..x_{N}$ this becomes

\begin{equation} \mu_{k}\equiv \left \langle x^{k} \right \rangle = \sum_{ i=1 }^{N} x_{i}^{k}f\left ( x \right ) \end{equation}

Note that a process where the probability function f(x) can be determined from each individual realisation $x_{1}..x_{N}$ is called “ergodic”. Ergodic processes make the mathematics much easier. For example, for our signal the probability function $f(x)$ is equal to $1/N$, because every sample has the same probability of occurrence and the process is ergodic.

The first order and second order moments have a well-known meaning. The first order moment ($k = 1$) is equal to

\begin{equation} \mu_{1}\equiv \left \langle x \right \rangle = \sum_{ i=1 }^{N} x_{i}f\left ( x \right ) \end{equation}

where we can include the probability function $f(x)=1/N$ and find

\begin{equation} \mu_{1} \equiv \left \langle x \right \rangle = \sum_{ i=1 }^{N} x_{i}f\left ( x \right ) = \sum_{ i=1 }^{N} \frac{x_{i}}{N} = \frac{1}{N} \sum_{ i=1 }^{N} x_{i} \end{equation}

which is known as the average of ensemble $x$. The second order moment is normally taken around the mean, and given by

\begin{equation} \mu_{2} \equiv \left \langle \left ( x-\mu_{1} \right )^{2} \right \rangle = \sum_{ i=1 }^{N} \left ( x_{i}-\mu_{1} \right ) ^{2} f\left ( x \right ) = \frac{1}{N} \sum_{ i=1 }^{N} \left ( x_{i}-\mu_{1} \right ) ^{2} = \sigma^{2} \end{equation}

The term $\sigma _{2}$ is the variance and is the square of the standard deviation $\sigma$.

Variance can be interpreted as the correlation of a signal with itself. Consider the expectation value of a signal $x$ on time t with the same signal $x$ on time $t+\tau$. We can write

\begin{equation} R_{xx}\left ( \tau \right ) \equiv \left \langle x\left ( t \right ) x\left ( t+\tau \right ) \right \rangle \label{eq:Rxx_Autocorrelation} \end{equation}

which is called the “autocorrelation function”. Note that

\begin{equation} R_{xx} \left ( 0 \right ) = \mu_{1} ^{2} + \sigma ^{2} \end{equation}

is the absolute maximum of $R_{xx}$ in $\tau = 0$. It is important to know that the autocorrelation function is the Fourier transform of the spectral density function $S_{xx}(f)$. So, with our discrete signals, we can simply use (inverse) Fast Fourier Transforms (FFT’s) to generate frequency spectra of our data. According to the Wiener-Khinchin Theorem, we can even generate the power spectrum immediately from the raw data without generating the cross-correlation function first.

In time critical applications, where we want to plot the variations in mean and/or standard deviation, there is normally no time to calculate the average of the whole dataset, or a series of subsets. In such cases, a moving average is taken. The algorithms to calculate a moving average have the advantage that we do not have to store and recall all the previous samples. There are many moving average options, but the most common is the “exponential moving average” EMA: the new one is simply added using a weighing factor. The equation is

\begin{equation} \mathtt {EMA}_{i} = \alpha x_{i} + \left ( 1 - \alpha \right ) \mathtt {EMA}_{i-1} \end{equation}

showing that we only have to store the previous EMA. The degree of weighing decrease is expressed as a constant smoothing factor $\alpha$, a number between $0$ and $1$. A higher α discounts older observations faster. Alternatively, α may be expressed in terms of $N$ time periods, where $\alpha = 2/(N+1)$. For example, $N = 19$ is equivalent to $\alpha = 0.1$. The half-life of the weights (the interval over which the weights decrease by a factor of two) is approximately $N/2.8854$.

Note that in our case of heartbeats, we have sufficient time (almost a second) between two subsequent samples. Therefore, moving averages are not in the standard representation methods of heart-rate variability.

The signal has a discrete time and a discrete amplitude. This means we can transform to the frequency domain using Fast Fourier Transforms (FFT). For a sampled sequence of $N$ points $x_{1}..x_{N}$, the FFT returns N frequency values $x_{1}..x_{N}$. This maps onto a horizontal frequency axis as

- $X_{1}$ is the power of the DC content, or in other words the power of the average heart rate. Because this is a huge peak in our case with respect to the more interesting information in the higher frequencies, it can be wise to remove the average DC signal from the signal $x_{1}..x_{N}$ first, and then do the FFT.
- $X_{N/2}$ is the power at the half the sampling frequency, so corresponding to the power at a frequency of ($2 \mu _{1})^{-1}Hz$. Note that due to the relation between amplitude and average asynchronous sample interval (equation \eqref{eq:IntervalIsAmplitude}), the spectrum has a relation between the frequency and power axis as well. This will not be the case when the synchronous subsampling of figure 2 is used. Half the sampling frequency is the highest frequency in the FFT result: the higher points are a symmetrical copy of the lower frequency points.
- The result is that the frequency resolution is equal to $(\mu _{1}N)^{-1}Hz$. So with $X_{1}$ at $0Hz$, $X_{2}$ is at $(\mu _{1}N)^{-1}Hz$.

A practical example is the following. When the Arduino board returns a nominal heart beat interval of $800ms$, this corresponds to a heart rate of $75 bpm$, or $0.125Hz$. An FFT of any number of samples would contain frequency information up to half the sampling frequency, so $0.625Hz$. With one minute of sampling, $60$ samples, the lowest frequency visible in the spectrum is $0.021Hz$. So, for observing variations in the heart rate lower than $0.021Hz$, we have to measure for at least one minute.

Besides the concept of autocorrelation, we can define cross-correlations between signals $x$ and $y$ (similar to the autocorrelation of equation \eqref{eq:Rxx_Autocorrelation}), using

\begin{equation} R_{xy}\left ( \tau \right ) \equiv \left \langle x\left ( t \right ) y\left ( t+\tau \right ) \right \rangle \end{equation}

The (inverse) Fourier transform of the cross-correlation function is called the cross-spectrum. For linear system we can find

\begin{equation} \left | S_{xy}\left ( f \right ) \right |^{2} = S_{xx} \left ( f \right ) S_{yy} \left ( f \right ) \end{equation}

a condition sometimes used to to check the linear relation between $x$ and $y$.

Wolfram Mathworld, http://mathworld.wolfram.com/

methods/hrv/hrv_stochastic_processes_and_resampling.txt · Last modified: 2015/11/29 00:35 (external edit)