Translator Disclaimer
1 August 2011 Estimating and validating the interbeat intervals of the heart using near-infrared spectroscopy on the human forehead
Author Affiliations +
In studies with near-infrared spectroscopy, the recorded signals contain information on the temporal interbeat intervals of the heart. If this cardiac information is needed exclusively and could directly be extracted, an additional electrocardiography device would be unnecessary. The aim was to estimate these intervals from signals measured with near-infrared spectroscopy with two novel approaches. In one approach, we model the heartbeat oscillations in these signals with a Fourier series where the coefficients and the fundamental frequency can continuously change over time. The time-dependent model parameters are estimated and used to calculate the interbeat intervals. The second approach uses empirical mode decomposition. The signal measured with near-infrared spectroscopy is empirically decomposed into a set of oscillatory components. The sum of a specific subset of them is an estimate of the pure heartbeat signal in which the diastolic peaks and consequential interbeat intervals are detected. We show in simultaneous electrocardiography and near-infrared spectroscopy measurements on 11 subjects (8 men and 3 woman with mean age 32.8 ± 8.1 yr), that the interbeat intervals (and the consequential pulse rate variability measures), estimated using the proposed approaches, are in high agreement with their correspondents from electrocardiography.



Near-infrared spectroscopy (NIRS) potentially measures changes in oxy- and deoxyhemoglobin1 caused by Mayer waves, breathing, cardiac activity, and brain activity. Figure 1a shows a NIRS signal featuring a heartbeat with a period of ≈1 s and a Mayer wave with a period of ≈10 s. The challenge is to extract the pure heartbeat and to estimate the heart rate variability (HRV) from NIRS signals.

HRV is (i) quantified by a set of statistical measures which result from time domain or frequency domain analysis of the time intervals between adjacent single heartbeats, called normal to normal (NN) intervals, in an electrocardiography (ECG) measurement, (ii) a quantitative and diagnostic marker of the autonomic nervous system's control on the heart rate, (iii) used in research and clinical studies,2 e.g., a relationship between the autonomic nervous system's activity and cardiovascular mortality3, 4, 5 has been shown.

The pulse rate variability (PRV) refers to the same set of statistical measures like HRV, but is extracted from photoplethysmography (PPG) signals. Since the operating principles of NIRS and PPG are similar, we use the term PRV when the measures are derived from NIRS signals; the term HRV is associated with ECG signals.

In NIRS signals, the heartbeat component is often present irrespective of the sensor's position. Thus, e.g., during functional studies, information on the NN intervals, and consequently on PRV, is already recorded. If this cardiac information is needed exclusively and would directly be extracted from NIRS, an ECG device would not be necessary.

In this paper, we focus on describing two approaches of how to estimate NN intervals from NIRS signals and showing proof of concept, which is a basis for future clinical studies. The approaches are validated by comparing their estimates with the corresponding ones from ECG and the resulting PRV and HRV measures, i.e., the standard deviation of the NN intervals (SDNN)

Eq. 1

[TeX:] \documentclass[12pt]{minimal}\begin{document} \begin{equation} S=\sqrt{\frac{1}{L} \sum _{l=1}^{L} (\chi _{l} - \overline{\chi })^2}, \end{equation} \end{document} S=1Ll=1L(χlχ¯)2,
with [TeX:] \documentclass[12pt]{minimal}\begin{document}$\overline{\chi }=\frac{1}{{L}} \sum _{l=1}^{L} \chi _{l}$\end{document} χ¯=1Ll=1Lχl and the square root of the mean squared differences of successive NN intervals (RMSSD)

Eq. 2

[TeX:] \documentclass[12pt]{minimal}\begin{document} \begin{equation} R=\sqrt{\frac{1}{L-1} \sum _{l=1}^{L-1} (\chi _{l+1} - \chi _{l})^2}, \end{equation} \end{document} R=1L1l=1L1(χl+1χl)2,
obtained from a set of L NN intervals χ1, …, χL. We omit frequency domain parameters, since they would not give additional insights concerning the agreement of NN intervals, SDNN, and RMSSD between NIRS and ECG.

One approach uses empirical mode decomposition (EMD); the other uses parameter estimation of a model for almost periodic signals (PEMAPS); both algorithms were designed to analyze nonstationary signals which do not necessarily contain strictly periodic components.



We tested the agreement between the NN intervals (and the resulting HRV measures) derived from ECG and the corresponding NN intervals (and the resulting PRV measures) estimated from NIRS. In 11 adult volunteers (8 men and 3 woman with mean age 32.8 ± 8.1 yr), ECG and NIRS were coregistered. During the experiment, each subject stood calmly, then sat calmly, and finally moved both arms and hands slowly while sitting; each of these three conditions took 5 min.



The ECG device was a MK3-ETA made by TOM Medical Entwicklungs GmbH; the NIRS device was a continuous wave MCPII.6

According to the user manual of the ECG device: electrode 1 was placed on the upper onset of the breastbone (sternal), electrode 2 was placed on the right lateral costal arch, and electrode 3 was placed submammary on the left.

The raw ECG signal is the sampled voltage difference between electrodes 1 and 3. Recommended ECG sampling rates are 250 to 500 Hz;2 we chose f ECG = 128 Hz to make ECG and NIRS signals, the latter sampled at a not alterable rate of 100 Hz, comparable.

NN intervals were extracted by detecting the peaks of the R waves in a detrended, but not further filtered, ECG signal. Detrending made peak detection more robust; we considered further denoising unnecessary, since the R waves in ECG signals are strong compared to the noise.

Electrode 2 was used for potential equalization.

The NIRS sensor is depicted in Fig. 2. Each source (light-emitting diode) sends light of constant intensity with wavelengths 750, 800, and 875 nm; each detector (photodiode) measures light intensity. The NIRS signal's sample values are proportional to the number of photons per time unit flowing through the photodiode, as well as the integral over the spectral sensitivity of the photodiode. MCPII was configured to drive 12 source/detector combinations, called “light paths,” each with 3 wavelengths, resulting in 36 data channels. The sampling rate was f NIRS = 100 Hz per data channel meaning that every 10 ms, 36 samples (1 sample per data channel) were acquired according to a time-multiplexed pattern.

Fig. 2

The NIRS sensor provides 4 light sources (circles) and 4 detectors (squares) and thus 16 light paths of which some are depicted as curved arrows.


NN intervals were estimated from one data channel which features a heartbeat component; the remaining 35 channels were ignored. The used light paths and wavelengths are stated in Table 1. The source/detector distances can be extracted from Fig. 2.

Fig. 8

The angulated line pairs represent R waves in a continuous-time ECG signal which is sampled (circles) at times 0, T, 2T, …. The positions of the R wave peaks are unknown. Each R peak (marked with an unfilled square) is enclosed by two samples; the higher one is detected.


Table 1

Cross correlation coefficients of NN intervals as defined in Eq. 14.

Subjectr ECG/PEMAPSr ECG/EMDSNRWave length [nm]Light- path

The NIRS sensor was placed on the right (from the subject's point of view) forehead, where usually in several data channels a clear heartbeat component is present.


Approach Based on PEMAPS

A periodic signal can be represented by a Fourier series. Specifically, one of the equidistant samples x 1, x 2, … of a real-valued periodic signal can be written as

Eq. 3

[TeX:] \documentclass[12pt]{minimal}\begin{document} \begin{equation} x_n = {\rm Re}\left(\sum _{k=0}^\infty A_k e^{jkn\Omega } \right) \end{equation} \end{document} xn=Rek=0AkejknΩ
with real coefficient A 0, complex coefficients A 1, A 2, …, and fundamental frequency [TeX:] \documentclass[12pt]{minimal}\begin{document}$\Omega \in \mathbb R$\end{document} ΩR .

We model the heartbeat component in NIRS signals as a trendless, almost periodic signal,7 i.e., A 1, A 2, … and Ω in Eq. 3 become time-dependent, and A 0 is discarded. Under this assumption, Eq. 3 changes to

Eq. 4

[TeX:] \documentclass[12pt]{minimal}\begin{document} \begin{equation} x_n = {\rm Re}\left(\sum _{k=1}^{K} A_{k,n} \cdot e^{j \Theta _n k}\right) \end{equation} \end{document} xn=Rek=1KAk,n·ejΘnk
with a single sample x n of the real-valued heartbeat component at discrete time n, time-dependent coefficients [TeX:] \documentclass[12pt]{minimal}\begin{document}$A_{1,n}, \ldots A_{K,n} \in {\mathbb C}$\end{document} A1,n,...AK,nC , time-dependent phase Θn ∈ [0, 2π], finite number of frequencies K, and

Eq. 5

[TeX:] \documentclass[12pt]{minimal}\begin{document} \begin{equation} A_{k,n+1} \approx A_{k,n}, \end{equation} \end{document} Ak,n+1Ak,n,

Eq. 6

[TeX:] \documentclass[12pt]{minimal}\begin{document} \begin{equation} \Theta _{n+1} = (\Theta _n + \Omega _n) \bmod 2\pi, \end{equation} \end{document} Θn+1=(Θn+Ωn)mod2π,

Eq. 7

[TeX:] \documentclass[12pt]{minimal}\begin{document} \begin{equation} \Omega _{n+1} \approx \Omega _n. \end{equation} \end{document} Ωn+1Ωn.
Equation 7 expresses the varying heart rate; Eq. 5 expresses the varying beat shape.

Let the NIRS signal be a noisy, trended version of the heartbeat component, i.e.,

Eq. 8

[TeX:] \documentclass[12pt]{minimal}\begin{document} \begin{equation} y_n = A_{0,n} + x_n + Z_n \end{equation} \end{document} yn=A0,n+xn+Zn
which means a single NIRS sample y n is a sum of the white Gaussian noise (sample Z n), the heartbeat (sample x n), and a slow trend (sample A 0, n). The latter models components slower than the heartbeat, i.e., Mayer waves, breathing, and brain activity.

Given a vector of N measured NIRS samples [TeX:] \documentclass[12pt]{minimal}\begin{document}$\bm y \break \stackrel{\scriptscriptstyle \bigtriangleup }{=}(y_1, \ldots, y_{N})$\end{document} y=(y1,...,yN) , the objective is to estimate the model parameter vector [TeX:] \documentclass[12pt]{minimal}\begin{document}$\bm \Theta \stackrel{\scriptscriptstyle \bigtriangleup }{=}(\Theta _1, \ldots, \Theta _{N})$\end{document} Θ=(Θ1,...,ΘN) and coefficient matrix

[TeX:] \documentclass[12pt]{minimal}\begin{document} \begin{equation*} \mathbf {A}=\left(\begin{array}{@{}c@{\quad}c@{\quad}c@{}}A_{0,1} & \ldots & A_{0,N} \\ \vdots & \ddots & \vdots \\ A_{K,1} & \ldots & A_{K,N} \end{array} \right), \end{equation*} \end{document} A=A0,1...A0,NAK,1...AK,N,
such that [TeX:] \documentclass[12pt]{minimal}\begin{document}$\textstyle \sum _{n=1}^{N} (y_n - x_n - A_{0,n})^2$\end{document} n=1N(ynxnA0,n)2 is minimal and reconstruct [TeX:] \documentclass[12pt]{minimal}\begin{document}$\bm x \stackrel{\scriptscriptstyle \bigtriangleup }{=}(x_1, \ldots, x_{N})$\end{document} x=(x1,...,xN) by applying the estimates in Eq. 4. We will use [TeX:] \documentclass[12pt]{minimal}\begin{document}$\bm A_{k,-}$\end{document} Ak, for the k'th row (harmonic index) of [TeX:] \documentclass[12pt]{minimal}\begin{document}$\mathbf {A}$\end{document} A .

Based on the measured NIRS samples [TeX:] \documentclass[12pt]{minimal}\begin{document}$\bm y$\end{document} y and a given estimate [TeX:] \documentclass[12pt]{minimal}\begin{document}$\hat{\mathbf {A}}$\end{document} Â (we mark estimates of parameters with a hat, e.g., [TeX:] \documentclass[12pt]{minimal}\begin{document}$\hat{\mathbf {A}}$\end{document} Â is an estimate of [TeX:] \documentclass[12pt]{minimal}\begin{document}$\mathbf {A}$\end{document} A ), [TeX:] \documentclass[12pt]{minimal}\begin{document}${{\bm \Theta}}$\end{document} Θ is estimated as

Eq. 9

[TeX:] \documentclass[12pt]{minimal}\begin{document} \begin{equation} \hat{\bm \Theta }={\mathop{\rm arg\, max}_{\bm \Theta \in [0, 2\pi ]^N}} f(\bm \Theta \,\, \vert \,\, \bm y, \hat{\mathbf {A}}) \end{equation} \end{document} Θ̂=argmaxΘ[0,2π]Nf(Θ|y,Â)
where the conditional probability density function f in Eq. 9 comprises the assumption in Eq. 8, the model 4, the relation 6, and the constraint 7. The latter is handled with adjustable strength by using prior knowledge of the upper and lower limits of Ωn. A heart rate has minimum H min and maximum H max values. Considering that [TeX:] \documentclass[12pt]{minimal}\begin{document}$\bm y$\end{document} y is sampled at f NIRS [Hz] and that during one single heartbeat the phases Θ1, Θ2, … in Eq. 4 traverse the interval [0, 2π], each heart rate H can be assigned to an angle growth Ω according to
[TeX:] \documentclass[12pt]{minimal}\begin{document} \begin{equation*} \Omega (H)=\frac{H \cdot 2\pi }{f_{\rm NIRS} \cdot 60}. \end{equation*} \end{document} Ω(H)=H·2πfNIRS·60.

The limits of Ωn are Ωmin = Ω(H min) and Ωmax = Ω(H max).

A k, n are estimated based on the measured NIRS samples y, estimates [TeX:] \documentclass[12pt]{minimal}\begin{document}${\hat{\bm A}}_{k-1,-}, \ldots, {\hat{\bm A}}_{0,-}$\end{document} Âk1,,...,Â0, , and [TeX:] \documentclass[12pt]{minimal}\begin{document}$\hat{\bm \Theta }$\end{document} Θ̂ from the previous iteration as

Eq. 10

[TeX:] \documentclass[12pt]{minimal}\begin{document} \begin{equation} \hat{A}_{k,n} = {\mathop{\rm arg\, max}_{A_{k,n} \in {\mathbb C}}}\ g(A_{k,n}\,\, \vert \,\, {\bm y}, {\hat{\bm{A}}}_{k-1,-}, \ldots, {\hat{\bm A}}_{0,-},{\hat{\bm \Theta }}) \end{equation} \end{document} Âk,n=argmaxAk,nCg(Ak,n|y,Âk1,,...,Â0,,Θ̂)
for increasing k. The function g in Eq. 10 comprises the assumption in Eqs. 8, 4, and the constraint 5. The latter is handled with adjustable strength by message damping as described in Ref. 7, Sec. 3.4 (“ [TeX:] \documentclass[12pt]{minimal}\begin{document}$\stackrel{\gamma }{=}$\end{document} =γ ”-node).

The probability density functions f in Eq. 9 and g in Eq. 10 are derived by using factor graphs and message passing algorithms described in Secs. 3.3 and 3.4 in Ref. 7.

The whole estimation algorithm is split into several building blocks whose interaction is depicted in Fig. 3.

Fig. 3

The building blocks of the PEMAPS algorithm.


Initially, the “ [TeX:] \documentclass[12pt]{minimal}\begin{document}${\bm A}_{0}$\end{document} A0 estimator” independently estimates the slow component [TeX:] \documentclass[12pt]{minimal}\begin{document}${\bm A}_{0,-}$\end{document} A0, by means of Eq. 10. Since for this [TeX:] \documentclass[12pt]{minimal}\begin{document}$\hat{\bm \Theta }$\end{document} Θ̂ is not needed, it is a one-time procedure based on [TeX:] \documentclass[12pt]{minimal}\begin{document}${\bm y}$\end{document} y only.

In the heartbeat component, most of the signal energy, apart from the noise, lies in the fundamental frequency coefficient [TeX:] \documentclass[12pt]{minimal}\begin{document}$\bm A_{1,-}$\end{document} A1, ; thus, a first rough estimate of the heartbeat component is a sinusoid. Its magnitude [TeX:] \documentclass[12pt]{minimal}\begin{document}$\tilde{A}_{1}$\end{document} Ã1 is calculated by the “initial A 1 estimator” block such that the sinusoid is of approximately the same energy as [TeX:] \documentclass[12pt]{minimal}\begin{document}${\bm y} - {\hat{\bm A}}_{0,-}$\end{document} yÂ0, .

Based on [TeX:] \documentclass[12pt]{minimal}\begin{document}${\hat{\bm A}}_{0,-}$\end{document} Â0, and [TeX:] \documentclass[12pt]{minimal}\begin{document}$\tilde{A}_{1}$\end{document} Ã1 , the “phase estimator” calculates [TeX:] \documentclass[12pt]{minimal}\begin{document}${\hat{\bm \Theta}}$\end{document} Θ̂ which is used by the “coefficient estimator” to calculate the full set of coefficient estimates [TeX:] \documentclass[12pt]{minimal}\begin{document}$\hat{{\bm A}}$\end{document} Â . For the specific algorithms of the phase estimator and the coefficient estimator refer to Ref. 7, Secs. 3.3 and 3.4.

At this point, it is possible to enter entries of [TeX:] \documentclass[12pt]{minimal}\begin{document}${\hat{\bm \Theta}}$\end{document} Θ̂ and [TeX:] \documentclass[12pt]{minimal}\begin{document}$\hat{\bm {A}}$\end{document} Â directly in Eq. 4 (done in the “reconstruction of [TeX:] \documentclass[12pt]{minimal}\begin{document}$\hat{x}$\end{document} x̂ ”-block) and obtain a first estimate [TeX:] \documentclass[12pt]{minimal}\begin{document}$\hat{\bm x}$\end{document} x̂ of the heartbeat component. The result can be improved by iterating in the coefficient estimator–phase estimator loop (Fig. 3) before calculating [TeX:] \documentclass[12pt]{minimal}\begin{document}$\hat{\bm x}$\end{document} x̂ .

Figure 1 illustrates how NN intervals are estimated from NIRS signals using PEMAPS. To convert the intervals Q in plot (e) of Fig. 1 to units of time, i.e., [s],

Eq. 11

[TeX:] \documentclass[12pt]{minimal}\begin{document} \begin{equation} t_Q = \frac{Q}{f_s} \end{equation} \end{document} tQ=Qfs
with sampling rate f s [Hz] can be used.

Fig. 1

(a) shows a raw NIRS signal. (b) shows the heartbeat component resulting from applying [TeX:] \documentclass[12pt]{minimal}\begin{document}${\hat{\bm A}}_{1,-}, \ldots, {\hat{\bm A}}_{K,-}$\end{document} Â1,,...,ÂK, and [TeX:] \documentclass[12pt]{minimal}\begin{document}${\hat{\bm \Theta}} = (\hat{\Theta }_1, \ldots, \hat{\Theta }_{1000})$\end{document} Θ̂=(Θ̂1,...,Θ̂1000) in Eq. 4. (c) shows [TeX:] \documentclass[12pt]{minimal}\begin{document}${\hat{\bm \Theta}}$\end{document} Θ̂ . (d) depicts [TeX:] \documentclass[12pt]{minimal}\begin{document}$\hat{\Theta }^{\prime }_1, \ldots, \hat{\Theta }^{\prime }_{999}$\end{document} Θ̂1,...,Θ̂999 with [TeX:] \documentclass[12pt]{minimal}\begin{document}$\hat{\Theta }_{n}^{\prime }=\hat{\Theta }_{n+1} - \hat{\Theta }_{n}$\end{document} Θ̂n=Θ̂n+1Θ̂n . Since with increasing n, the values [TeX:] \documentclass[12pt]{minimal}\begin{document}$\hat{\Theta }_n$\end{document} Θ̂n increase monotonically with respect to the modulo 2π function [dictated by, Eq. 6], [TeX:] \documentclass[12pt]{minimal}\begin{document}$\hat{\Theta }_{n}^{\prime }$\end{document} Θ̂n must be positive for all n except for the transition indices, where the modulo operator takes effect [peaks in (d)]. Finally, the samples between adjacent peaks are counted (e).



Approach Based on EMD

EMD decomposes a signal into a finite number of oscillatory modes, called intrinsic mode functions (IMFs), by their characteristic time scales. The IMFs are derived empirically from the measured signal without any prior knowledge or model.

In Ref. 8, Sec. 5 an IMF is formally defined, and the decomposition procedure, called “the sifting process,” is described in detail. The latter can be summarized as follows:

  1. Set i = 1.

  2. Set [TeX:] \documentclass[12pt]{minimal}\begin{document}$\bm u \stackrel{\scriptscriptstyle \bigtriangleup }{=}\bm y$\end{document} u=y , with measured samples [TeX:] \documentclass[12pt]{minimal}\begin{document}$\bm y=y_1, y_2, \ldots$\end{document} y=y1,y2,... .

  3. Set the initial residual [TeX:] \documentclass[12pt]{minimal}\begin{document}$\bm r_{1} \stackrel{\scriptscriptstyle \bigtriangleup }{=}\bm y$\end{document} r1=y (relevant in step 8).

  4. Find the local extrema of [TeX:] \documentclass[12pt]{minimal}\begin{document}$\bm u$\end{document} u .

  5. Fit a cubic spline through all maxima which is an upper envelope [TeX:] \documentclass[12pt]{minimal}\begin{document}$\bm e_{u}$\end{document} eu of [TeX:] \documentclass[12pt]{minimal}\begin{document}$\bm u$\end{document} u , analogously through all minima, which is a lower envelope [TeX:] \documentclass[12pt]{minimal}\begin{document}$\bm e_{l}$\end{document} el of [TeX:] \documentclass[12pt]{minimal}\begin{document}$\bm u$\end{document} u .

  6. The mean of the two envelopes [TeX:] \documentclass[12pt]{minimal}\begin{document}$\bm m_{i}=(\bm e_{l}+\bm e_{u})/2$\end{document} mi=(el+eu)/2 is subtracted from [TeX:] \documentclass[12pt]{minimal}\begin{document}$\bm u$\end{document} u : [TeX:] \documentclass[12pt]{minimal}\begin{document}$\bm h_{i}=\bm u-\bm m_{i}$\end{document} hi=umi .

  7. Treat [TeX:] \documentclass[12pt]{minimal}\begin{document}$\bm h_{i}$\end{document} hi as the new input signal [TeX:] \documentclass[12pt]{minimal}\begin{document}$\bm u$\end{document} u and go to step 4. Steps 4 – 7 are repeated until [TeX:] \documentclass[12pt]{minimal}\begin{document}$\bm h_{i}$\end{document} hi becomes a curve with as many zero crossings as extrema and the upper and lower envelopes of [TeX:] \documentclass[12pt]{minimal}\begin{document}$\bm h_{i}$\end{document} hi become symmetric with regard to the zero line. Then [TeX:] \documentclass[12pt]{minimal}\begin{document}$\bm h_{i}$\end{document} hi is the i'th IMF denoted as [TeX:] \documentclass[12pt]{minimal}\begin{document}$\bm c_{i}$\end{document} ci .

  8. Let [TeX:] \documentclass[12pt]{minimal}\begin{document}$\bm u \stackrel{\scriptscriptstyle \bigtriangleup }{=}\bm r_{i+1}=\bm r_{i} - \bm c_{i}$\end{document} u=ri+1=rici .

  9. Increase i and go to step 4.

The sifting process stops when [TeX:] \documentclass[12pt]{minimal}\begin{document}$\bm r_{i+1}$\end{document} ri+1 in step 8 is a constant, a monotonic slope, or a function with only one extremum. The original signal is given as

[TeX:] \documentclass[12pt]{minimal}\begin{document} \begin{equation*} \bm y=\sum _{i=1}^{N} \bm c_i + \bm r_{N+1}. \end{equation*} \end{document} y=i=1Nci+rN+1.

The IMFs with lower indices i represent fast and those with higher indices slow oscillations.

By examining the IMFs by eye, one can recognize (based on their mean period length) which of them are related to the heartbeat component. The sum of this subset of IMFs represents an estimate [TeX:] \documentclass[12pt]{minimal}\begin{document}$\hat{\bm x}$\end{document} x̂ of the heartbeat component, of which an example is illustrated in Fig. 4.

Fig. 4

An estimate of the heartbeat component computed using EMD. All diastolic maxima are marked with circles, some of the nondiastolic maxima are marked with squares.


In NIRS, near-infrared light penetrates tissue. The more blood the tissue contains, the more light is absorbed, i.e., the less light reaches the detector. Consequently, the light intensity at the detector is highest during diastole. We use the diastolic maxima (circles in Fig. 4) to estimate the NN intervals.

[TeX:] \documentclass[12pt]{minimal}\begin{document}$\hat{x}_n$\end{document} x̂n is assumed to be the m'th diastolic maximum P m with entry index I(P m) = n in [TeX:] \documentclass[12pt]{minimal}\begin{document}${\hat{\bm x}}$\end{document} x̂ , if it is the m'th entry of [TeX:] \documentclass[12pt]{minimal}\begin{document}${\hat{\bm x}}$\end{document} x̂ for which the conditions

Eq. 12

[TeX:] \documentclass[12pt]{minimal}\begin{document} \begin{eqnarray} \hat{x}_n - \hat{x}_{n-1} > 0, \nonumber \\ \hat{x}_n - \hat{x}_{n+1} > 0, \nonumber \\ \hat{x}_n > \xi _1, \end{eqnarray} \end{document} x̂nx̂n1>0,x̂nx̂n+1>0,x̂n>ξ1,

Eq. 13

[TeX:] \documentclass[12pt]{minimal}\begin{document} \begin{eqnarray} I(P_m) - I(P_{m-1}) > \xi _2 \end{eqnarray} \end{document} I(Pm)I(Pm1)>ξ2
hold. With appropriate ξ1 (e.g., [TeX:] \documentclass[12pt]{minimal}\begin{document}$\xi _1 \stackrel{\scriptscriptstyle \bigtriangleup }{=}45$\end{document} ξ1=45 in Fig. 4), the condition 12 discards the nondiastolic maxima (marked by squares in Fig. 4). Sometimes a diastolic maximum is lower than a nondiastolic one; then there is no ξ1, such that the diastolic maximum is detected and the nondiastolic one is omitted. Choosing ξ1 low enough such that all diastolic (and as well some nondiastolic) maxima are detected, using Eq. 13 to find neighboring maxima which are too close to each other, and discarding the smaller maximum solves the problem.

The NN intervals can now be estimated by counting the samples between the detected adjacent diastolic maxima. Finally, Eq. 11 can be used to convert the intervals to units of time.


Validation: Data Analysis

All recorded NIRS and ECG signals were evaluated offline. Segments with movement artifacts (sudden, typically between 0.5 and 2 s long changes, during which the standard deviation of the signal increases more than 100%) were excluded from the evaluation. In all our signals, such changes are distinct, and we recognized them easily by eye. Alternatively, the algorithm in Ref. 9 can be used to detect excessive values in the standard deviation of a NIRS or ECG signal, which are larger than an intuitively chosen threshold (see Ref. 9, Sec. 2.2.2). The corresponding NIRS and ECG samples are then assumed to be distorted by movement artifacts. Since in our case the latter are distinct without edge cases, this algorithm would exclude the same segments as we recognized by eye.

One subject was excluded from the evaluation due to technical problems.

From NIRS, NN intervals were estimated from a data channel with a clear heartbeat component; the remaining data channels were ignored.

From a detrended ECG signal, NN intervals were estimated by detecting the peaks of the R waves with the same procedure as detecting diastolic peaks in the heartbeat component estimated from NIRS signals using EMD (procedure described at the end of Sec. 5). For every subject/ECG signal, ξ1 in Eq. 12 was chosen individually. For all subjects/ECG signals, ξ2 = 61 samples in Eq. 13.

For all subjects/NIRS signals, PEMAPS was set up with K = 3 in Eq. 4, two iterations in the phase estimator/coefficient estimator loop (Fig. 3), message damping factors for harmonic indices k = 0: γ = 0.936, k = 1: γ = 0.98, k = 2: γ = 0.985 and k = 3: γ = 0.989 (see Ref. 7, Sec. 3.4, “ [TeX:] \documentclass[12pt]{minimal}\begin{document}$\stackrel{\gamma }{=}$\end{document} =γ ”-node) and Ωmin = 0.026 rad and Ωmax = 0.131 rad (see Sec. 4). [TeX:] \documentclass[12pt]{minimal}\begin{document}${\hat{\bm A}}_{0,-}$\end{document} Â0, was smoothed by feeding it back to the [TeX:] \documentclass[12pt]{minimal}\begin{document}$\bm A_0$\end{document} A0 estimator (Fig. 3) as the new input signal and re-estimating [TeX:] \documentclass[12pt]{minimal}\begin{document}$\bm A_{0,-}$\end{document} A0, . For every subject, this procedure was performed twice successively.

For the EMD-based approach, ξ1 was chosen for every subject/NIRS signals individually; ξ2 = 48 samples, for all subjects/NIRS signals.


Validation: Results

Assuming that, for the j'th subject we evaluated a set of W j NN intervals [TeX:] \documentclass[12pt]{minimal}\begin{document}$\bm \chi ^{(j)}_{\rm ECG} \stackrel{\scriptscriptstyle \bigtriangleup }{=}(\chi ^{(j)}_{E,1}, \ldots, \chi ^{(j)}_{E,W_j})$\end{document} χECG(j)=(χE,1(j),...,χE,Wj(j)) derived from ECG and [TeX:] \documentclass[12pt]{minimal}\begin{document}$\bm \chi ^{(j)}_{\rm NIRS} \stackrel{\scriptscriptstyle \bigtriangleup }{=}(\chi ^{(j)}_{N,1}, \ldots, \chi ^{(j)}_{N,W_j})$\end{document} χNIRS(j)=(χN,1(j),...,χN,Wj(j)) estimated from NIRS using PEMAPS and EMD, Table 1 shows, for each subject, i.e., for j = 1, …, 11, the cross-correlation coefficient

Eq. 14

[TeX:] \documentclass[12pt]{minimal}\begin{document} \begin{equation} r \stackrel{\scriptscriptstyle \bigtriangleup }{=}\frac{1}{W_j-1}\frac{\sum\limits _{w=1}^{W_j}(\chi ^{(j)}_{{\rm E},w} - \overline{\chi }^{(j)}_{\rm ECG}) \cdot (\chi ^{(j)}_{{\rm N},w} - \overline{\chi }^{(j)}_{\rm NIRS})}{\hat{\sigma }^{(j)}_{\chi _{\rm E}} \cdot \hat{\sigma }^{(j)}_{\chi _{\rm N}}} \end{equation} \end{document} r=1Wj1w=1Wj(χE,w(j)χ¯ECG(j))·(χN,w(j)χ¯NIRS(j))σ̂χE(j)·σ̂χN(j)
with empirical means [TeX:] \documentclass[12pt]{minimal}\begin{document}$\overline{\chi }^{(j)}_{\rm ECG} = {1}/{W_j} \sum _{w=1}^{W_j} \chi ^{(j)}_{{\rm E},w}$\end{document} χ¯ECG(j)=1/Wjw=1WjχE,w(j) and [TeX:] \documentclass[12pt]{minimal}\begin{document}$\overline{\chi }^{(j)}_{\rm NIRS} = {1}/{W_j} \sum _{w=1}^{W_j} \chi ^{(j)}_{{\rm N},w}$\end{document} χ¯NIRS(j)=1/Wjw=1WjχN,w(j) , and empirical standard deviations [TeX:] \documentclass[12pt]{minimal}\begin{document}$\hat{\sigma }^{(j)}_{\chi _{\rm E}} = \smash{\sqrt{{1}/{(W_j-1)} \sum _{w=1}^{W_j} (\chi ^{(j)}_{{\rm E},w} - \overline{\chi }^{(j)}_{\rm ECG})^2}}$\end{document} σ̂χE(j)=1/(Wj1)w=1Wj(χE,w(j)χ¯ECG(j))2 and [TeX:] \documentclass[12pt]{minimal}\begin{document}$\hat{\sigma }^{(j)}_{\chi _{\rm N}} = \smash{\sqrt{{1}/{(W_j-1)} \sum _{w=1}^{W_j} (\chi ^{(j)}_{{\rm N},w} - \overline{\chi }^{(j)}_{\rm NIRS})^2}}$\end{document} σ̂χN(j)=1/(Wj1)w=1Wj(χN,w(j)χ¯NIRS(j))2 . In addition, for each subject, the used wavelength, light path, and an estimate of the signal-to-noise ratio (SNR)

Eq. 15

[TeX:] \documentclass[12pt]{minimal}\begin{document} \begin{equation} {\rm SNR}=\frac{\sum _{n=1}^{N} \hat{x}_n^2}{\sum _{n=1}^{N} (y_n - \hat{x}_n - \hat{A}_{0,n})^2} \end{equation} \end{document} SNR=n=1Nx̂n2n=1N(ynx̂nÂ0,n)2
with reconstructed heartbeat sample [TeX:] \documentclass[12pt]{minimal}\begin{document}$\hat{x}_n$\end{document} x̂n in Eq. 4, NIRS sample y n, and estimated sample [TeX:] \documentclass[12pt]{minimal}\begin{document}$\hat{A}_{0,n}$\end{document} Â0,n of the slow trend in Eq. 8 are given in Table 1.

Table 2 shows, per subject, how SDNN changes, from standing to sitting and, during sitting from resting to moving hands. Table 3 shows the same for RMSSD. In both tables, the physiological adaptation phases, i.e., the first 40 s during sitting and the first 10 s during the exercise, were excluded from the analysis, and SDNN/RMSSD values of all 11 subjects are compared between ECG and NIRS using cross-correlation coefficients r and p-values (t-test); the usability of the t-test has been approved by a Lilliefors test.

Table 2

Change of SDNN (in percent) related to condition changing; cross-correlation coefficients r and p-values (t-test) give a comparison of all 11 SDNN values between ECG and NIRS.

From standing to sittingFrom resting to exercise
p4.61 × 10−141.22 × 10−113.46 × 10−135.1 × 10−12

Table 3

Change of RMSSD (in percent) related to condition changing; cross-correlation coefficients r and p-values (t-test) give a comparison of all 11 RMSSD values between ECG and NIRS.

From standing to sittingFrom resting to exercise
p6.37 × 10−114.28 × 10−068.14 × 10−093.59 × 10−05

Let [TeX:] \documentclass[12pt]{minimal}\begin{document}$\bm \chi _{\rm ECG}$\end{document} χECG be a vector of NN intervals derived from ECG of all 11 subjects; let [TeX:] \documentclass[12pt]{minimal}\begin{document}$\bm \chi _{\rm PEMAPS}$\end{document} χPEMAPS and [TeX:] \documentclass[12pt]{minimal}\begin{document}$\bm \chi _{\rm EMD}$\end{document} χEMD be the same, but estimated from NIRS using PEMAPS and EMD. The values of [TeX:] \documentclass[12pt]{minimal}\begin{document}$\bm \chi _{\rm ECG}$\end{document} χECG , [TeX:] \documentclass[12pt]{minimal}\begin{document}$\bm \chi _{\rm PEMAPS}$\end{document} χPEMAPS , and [TeX:] \documentclass[12pt]{minimal}\begin{document}$\bm \chi _{\rm EMD}$\end{document} χEMD at a given index are three estimates of the same NN interval. Figure 5 shows the agreement between [TeX:] \documentclass[12pt]{minimal}\begin{document}$\bm \chi _{\rm ECG}$\end{document} χECG and [TeX:] \documentclass[12pt]{minimal}\begin{document}$\bm \chi _{\rm PEMAPS}$\end{document} χPEMAPS on the one hand, and [TeX:] \documentclass[12pt]{minimal}\begin{document}$\bm \chi _{\rm ECG}$\end{document} χECG and [TeX:] \documentclass[12pt]{minimal}\begin{document}$\bm \chi _{\rm EMD}$\end{document} χEMD on the other hand. In all plots in this section, the lack of agreement between two measures is summarized by the mean [TeX:] \documentclass[12pt]{minimal}\begin{document}$\hat{\mu }$\end{document} μ̂ and the standard deviation [TeX:] \documentclass[12pt]{minimal}\begin{document}$\hat{\sigma }$\end{document} σ̂ of their difference points.

Fig. 5

Agreement between NN intervals in all subjects derived from ECG and those estimated from NIRS using PEMAPS (upper plot), and using EMD (lower plot); [TeX:] \documentclass[12pt]{minimal}\begin{document}$\hat{\mu }_1=0.00008$\end{document} μ̂1=0.00008 s, [TeX:] \documentclass[12pt]{minimal}\begin{document}$\hat{\sigma }_1=0.00755$\end{document} σ̂1=0.00755 s, [TeX:] \documentclass[12pt]{minimal}\begin{document}$\hat{\mu }_2\break=0.00008$\end{document} μ̂2=0.00008 s and [TeX:] \documentclass[12pt]{minimal}\begin{document}$\hat{\sigma }_2=0.01133$\end{document} σ̂2=0.01133 s. These types of plots, in conjunction with testing agreement, was proposed in Ref. 10. The size of a single point in the plot is proportional to the number of occurrences of the corresponding entry pairs.


Due to the sampling, the entries of [TeX:] \documentclass[12pt]{minimal}\begin{document}$\bm \chi _{\rm ECG}$\end{document} χECG , [TeX:] \documentclass[12pt]{minimal}\begin{document}$\bm \chi _{\rm PEMAPS}$\end{document} χPEMAPS , [TeX:] \documentclass[12pt]{minimal}\begin{document}$\bm \chi _{\rm EMD}$\end{document} χEMD , and thus the differences [TeX:] \documentclass[12pt]{minimal}\begin{document}$\bm \chi _{\rm ECG}$\end{document} χECG [TeX:] \documentclass[12pt]{minimal}\begin{document}$\bm \chi _{\rm PEMAPS}$\end{document} χPEMAPS and [TeX:] \documentclass[12pt]{minimal}\begin{document}$\bm \chi _{\rm ECG}$\end{document} χECG [TeX:] \documentclass[12pt]{minimal}\begin{document}$\bm \chi _{\rm EMD}$\end{document} χEMD , are discrete. Hence, a regular raster can be recognized in Fig. 5, and often the pair of values at a given index in [TeX:] \documentclass[12pt]{minimal}\begin{document}$\bm \chi _{\rm ECG}$\end{document} χECG and [TeX:] \documentclass[12pt]{minimal}\begin{document}$\bm \chi _{\rm PEMAPS}$\end{document} χPEMAPS or [TeX:] \documentclass[12pt]{minimal}\begin{document}$\bm \chi _{\rm ECG}$\end{document} χECG and [TeX:] \documentclass[12pt]{minimal}\begin{document}$\bm \chi _{\rm EMD}$\end{document} χEMD is not unique. Thus, the size of a single point in Fig. 5 is proportional to the number of occurrences of the corresponding entry pairs.

The mean of NN intervals in all subjects (derived from ECG signals), i.e., the mean of x-coordinates of all points in Fig. 5, is 0.84 s. An NN interval estimated with PEMAPS or EMD differs from this mean on average by [TeX:] \documentclass[12pt]{minimal}\begin{document}${\hat{\sigma }_1}/{0.84 {\rm s}}=0.9$\end{document} σ̂1/0.84s=0.9 %, or [TeX:] \documentclass[12pt]{minimal}\begin{document}${\hat{\sigma }_2}/{0.84 {\rm s}}\break =1.35$\end{document} σ̂2/0.84s=1.35 %, respectively.

Figure 6 shows SDNN, calculated from ECG using Eq. 1, S ECG, versus its difference to S PEMAPS and S EMD. Figure 7 shows the same for RMSSD. The mean of SDNN values in all subjects (derived from ECG), i.e., the mean of x-coordinates of all points in Fig. 6, is 0.08 s. A SDNN value, derived with PEMAPS or EMD, differs from this mean on average by [TeX:] \documentclass[12pt]{minimal}\begin{document}${\hat{\sigma }_3}/{0.08{\rm s}}=0.4$\end{document} σ̂3/0.08s=0.4 % or [TeX:] \documentclass[12pt]{minimal}\begin{document}${\hat{\sigma }_4}/{0.08{\rm s}}=0.88$\end{document} σ̂4/0.08s=0.88 %, respectively. The average differences for RMSSD are [TeX:] \documentclass[12pt]{minimal}\begin{document}${\hat{\sigma }_5}/{0.038{\rm s}}=2.55$\end{document} σ̂5/0.038s=2.55 % and [TeX:] \documentclass[12pt]{minimal}\begin{document}${\hat{\sigma }_6}/{0.038{\rm s}}=5.16$\end{document} σ̂6/0.038s=5.16 %.

Fig. 6

Agreement between SDNN derived from ECG and SDNN estimated from NIRS using PEMAPS (upper plot), and EMD (lower plot); [TeX:] \documentclass[12pt]{minimal}\begin{document}$\hat{\mu }_3=-0.00084$\end{document} μ̂3=0.00084 s, [TeX:] \documentclass[12pt]{minimal}\begin{document}$\hat{\sigma }_3=0.00032$\end{document} σ̂3=0.00032 s, [TeX:] \documentclass[12pt]{minimal}\begin{document}$\hat{\mu }_4=-0.0014$\end{document} μ̂4=0.0014 s, and [TeX:] \documentclass[12pt]{minimal}\begin{document}$\hat{\sigma }_4\break=0.0007\,{\rm s}$\end{document} σ̂4=0.0007s . Every point was derived from the entire experiment with one subject.


Fig. 7

Agreement between RMSSD derived from ECG and RMSSD estimated from NIRS using PEMAPS (upper plot), and EMD (lower plot); [TeX:] \documentclass[12pt]{minimal}\begin{document}$\hat{\mu }_5=-0.00247$\end{document} μ̂5=0.00247 s, [TeX:] \documentclass[12pt]{minimal}\begin{document}$\hat{\sigma }_5=0.00097$\end{document} σ̂5=0.00097 s, [TeX:] \documentclass[12pt]{minimal}\begin{document}$\hat{\mu }_6=-0.00452$\end{document} μ̂6=0.00452 s, and [TeX:] \documentclass[12pt]{minimal}\begin{document}$\hat{\sigma }_6=0.00196$\end{document} σ̂6=0.00196 s. Every point was derived from the entire experiment with one subject.


In Figs. 6, 7, the differences between the two measures are biased, i.e., [TeX:] \documentclass[12pt]{minimal}\begin{document}$\hat{\mu }_3, \ldots, \hat{\mu }_6$\end{document} μ̂3,...,μ̂6 are different from zero. This is caused by an error when detecting R peaks (ECG) and diastoles (NIRS) in discrete-time signals. The issue is explained in the appendix of this paper.


Discussion and Conclusion

As stated at the end of Sec. 7, the average discrepancies between SDNN from ECG and SDNN from NIRS are 0.4% (PEMAPS) and 0.88% (EMD). In Ref. 3, the risk of mortality was compared between two groups of subjects; in one group SDNN < 0.05 s, in the second group SDNN [TeX:] \documentclass[12pt]{minimal}\begin{document}$>$\end{document} > 0.1 s, which is [TeX:] \documentclass[12pt]{minimal}\begin{document}$>$\end{document} > 100% higher than in the first group. Conclusively, SDNN derived from NIRS with the proposed approaches could be sufficiently accurate to derive the risk of mortality.

As shown in Ref. 8, EMD can be applied to signals from various sources in nature and does not depend on amplitude or frequency of the oscillations to be reconstructed. The version of PEMAPS used in this work is limited to signals with a strong (initial estimate [TeX:] \documentclass[12pt]{minimal}\begin{document}$\tilde{A}_{1}$\end{document} Ã1 , Sec. 4) and rather high (≈0.4 Hz) fundamental frequency.

The NIRS instrumentation should capture at least the fundamental frequency of the heartbeat oscillation, i.e., the sampling rate of NIRS should at least be 2H max = 6 Hz if H max = 180 bpm is the maximal heart rate. We tested this successfully by interpolating a 6 Hz NIRS signal to 100 Hz from which the NN intervals were derived with the proposed approaches. The correlation coefficients were as high as deriving the NN intervals from a real 100 Hz NIRS signal. Interpolating was necessary, since low sampling rates induce an error when deriving NN intervals. In Ref. 11, for example, the influence of this error on the power spectrum of the NN intervals is quantified. We address the influence of this error on SDNN and RMSSD in the appendix of this paper.

Both methods can be used with arbitrary signal lengths. The only limiting factor is RAM. If the amount of RAM does not suffice to analyze a 24 h recording as one block, the signals can be split into several blocks, each being analyzed separately. The implementations we used in this work require for a 15 min long NIRS signal ≈1.3 GB (PEMAPS) and ≈60 MB (EMD) of RAM.

Between subjects, the used wavelengths and SNR Eq. 15 vary considerably, which has no high impact on the correlation coefficients in Table 1, e.g., compare subjects 1 and 8 concerning SNR and subjects 5 and 11 concerning wavelength.

The higher SDNN, the larger the variability between all NN intervals in the evaluated NIRS or ECG signal. The higher RMSSD, the larger the variability between successive NN intervals in the evaluated NIRS or ECG signal. According to Ref. 12, RMSSD relates mainly to the parasympathethic activity, whereas SDNN relates to sympathetic and parasympathetic activity. In all subjects (except for subject 8) in Table 3, RMSSD from ECG considerably increases from standing to sitting and decreases during sitting from resting to performing the exercise, which is also detectable from NIRS signals using PEMAPS and EMD. On the contrary, SDNN shows no such uniform changes over all subjects. In each subject in Table 2, the changes of SDNN from NIRS agree with their ECG correspondent more than the changes of RMSSD. In addition, the more the RMSSD and SDNN values differ from 0, the better they agree between ECG and NIRS.

In Refs. 13, 14, 15, 16, NN intervals and HRV measures derived from ECG were compared to their correspondents from PPG. In Ref. 13 correlations [TeX:] \documentclass[12pt]{minimal}\begin{document}$1>r>0.97$\end{document} 1>r>0.97 were found between HRV and PRV in 44 subjects. Frequency domain and time domain measures were calculated. All signals were sampled at 400 Hz. In Ref. 14 where 10 subjects participated, these correlations are in the range [TeX:] \documentclass[12pt]{minimal}\begin{document}$0.99985>r>0.962$\end{document} 0.99985>r>0.962 . All signals were sampled at 400 Hz. In Ref. 15 the median correlation coefficient of the NN intervals derived from ECG and their PPG correspondents in 10 subjects is r = 0.97 (compared to the medians 0.995355 and 0.988709 in Table 1). All signals were sampled at 1 kHz. In Ref. 16, the same coefficient was derived from 42 subjects as r = 0.91 (including an outlier). The ECG signals were sampled at 200 Hz; PPG signals were sampled at 100 Hz. We conclude that our results show a slightly higher agreement, although we used considerably lower sampling rates.

Many factors can affect the agreement between NN intervals derived from ECG and the corresponding ones estimated from NIRS using PEMAPS and EMD, e.g., external forces on the arterial vessels, pathologies, methodical problems, small and unnoticed movement artefacts, and variability in time which the pulse pressure waveform takes to propagate through the arterial tree. Our results show that with healthy subjects who are not exposed to mechanical forces and behave calmly during the experiment, the impact of these factors is negligible.

PEMAPS estimates correspond closer to ECG than the EMD estimates, i.e., [TeX:] \documentclass[12pt]{minimal}\begin{document}$\hat{\sigma }_1<\hat{\sigma }_2$\end{document} σ̂1<σ̂2 in Fig. 5, [TeX:] \documentclass[12pt]{minimal}\begin{document}$|\hat{\mu }_3|<|\hat{\mu }_4|$\end{document} |μ̂3|<|μ̂4| and [TeX:] \documentclass[12pt]{minimal}\begin{document}$\hat{\sigma }_3<\hat{\sigma }_4$\end{document} σ̂3<σ̂4 in Fig. 6, and [TeX:] \documentclass[12pt]{minimal}\begin{document}$|\hat{\mu }_5|<|\hat{\mu }_6|$\end{document} |μ̂5|<|μ̂6| and [TeX:] \documentclass[12pt]{minimal}\begin{document}$\hat{\sigma }_5<\hat{\sigma }_6$\end{document} σ̂5<σ̂6 in Fig. 7.

Compared to EMD, the version of PEMAPS used in this work is (i) slower and requires more RAM, (ii) less error-prone, and (iii) set up in a more general way, i.e., the input arguments of the PEMAPS implementation are the same for all subjects. When using EMD, the user must examine the intrinsic mode functions (see Sec. 5) by eye and set ξ1 in Eq. 12 manually for every subject.

According to Ref. 17 the clinical outlook of near-infrared techniques is noninvasive (i) brain imaging by providing functional and metabolic maps of the activated brain cortex, (ii) measurement of changes and absolute values in oxy- and deoxyhemoglobin, (iii) measurement of blood pressure changes, and (iv) measurement of respiratory rate. In all these applications, NIRS coregisters information on NN intervals; this could give new physiological insights. Furthermore, in multimodal measurement setups with strong electromagnetic fields, e.g., as caused by magnetic resonance imaging or computed tomography, ECG may be disturbed, while NIRS would function properly.


Appendix A: SDNN and Sampling

In this section, based on the error induced by the sampling, [TeX:] \documentclass[12pt]{minimal}\begin{document}$\hat{\mu }_3 < 0$\end{document} μ̂3<0 and [TeX:] \documentclass[12pt]{minimal}\begin{document}$\hat{\mu }_4 < 0$\end{document} μ̂4<0 in Fig. 6 are motivated, and a lower bound for [TeX:] \documentclass[12pt]{minimal}\begin{document}$\hat{\sigma }_1$\end{document} σ̂1 and [TeX:] \documentclass[12pt]{minimal}\begin{document}$\hat{\sigma }_2$\end{document} σ̂2 in Fig. 5 is derived.

Let I l be the unknown (continuous-time) position of the l'th maximum, and let [TeX:] \documentclass[12pt]{minimal}\begin{document}$\hat{I}_l$\end{document} Îl be the (discrete-time) position of the detected maximum (see Fig. 8). The continuous-valued discrepancy [TeX:] \documentclass[12pt]{minimal}\begin{document}$D_l=I_l - \hat{I}_l$\end{document} Dl=IlÎl is assumed to be uniformly distributed over the interval [TeX:] \documentclass[12pt]{minimal}\begin{document}$[-\frac{T}{2}, \frac{T}{2}]$\end{document} [T2,T2] with variance [TeX:] \documentclass[12pt]{minimal}\begin{document}$\sigma _{D}^{2}={T^2}/{12}$\end{document} σD2=T2/12 .

The discrepancy between the l'th real and the l'th detected NN interval is given as

Eq. 16

[TeX:] \documentclass[12pt]{minimal}\begin{document} \begin{eqnarray} \Lambda _{l} &=& I_{l+1} - I_{l} - (\hat{I}_{l+1} - \hat{I}_{l})\nonumber \\ &=& D_{l+1} - D_{l}. \end{eqnarray} \end{document} Λl=Il+1Il(Îl+1Îl)=Dl+1Dl.
It follows that the difference of two uniformly distributed random variables Λl must be triangularly distributed over the interval [ − T, T] with variance [TeX:] \documentclass[12pt]{minimal}\begin{document}$\sigma _{\Lambda }^{2}=2\sigma _{D}^{2}={T^2}/{6}$\end{document} σΛ2=2σD2=T2/6 . In our case, this variance is [TeX:] \documentclass[12pt]{minimal}\begin{document}$\sigma _{\Lambda _{\rm E}}^{2}={T_{\rm ECG}^2}/{6}={1}/({6f_{\rm ECG}^2}) \approx 0.00001017{\rm s}^2$\end{document} σΛE2=TECG2/6=1/(6fECG2)0.00001017s2 (ECG) and [TeX:] \documentclass[12pt]{minimal}\begin{document}$\sigma _{\Lambda _{\rm N}}^{2}={1}/({6f_{\rm NIRS}^2}) \approx 0.00001667{\rm s}^2$\end{document} σΛN2=1/(6fNIRS2)0.00001667s2 (NIRS).

By rearranging Eq. 16, it follows that

Eq. 17

[TeX:] \documentclass[12pt]{minimal}\begin{document} \begin{equation} \hat{I}_{l+1} - \hat{I}_{l}=I_{l+1} - I_{l} - \Lambda _{l}. \end{equation} \end{document} Îl+1Îl=Il+1IlΛl.
Let the l'th entry of [TeX:] \documentclass[12pt]{minimal}\begin{document}$\bm \chi _{\rm ECG}$\end{document} χECG in Fig. 5 be modeled, according to Eq. 17, as

Eq. 18

[TeX:] \documentclass[12pt]{minimal}\begin{document} \begin{equation} \chi _{{\rm E},l} \stackrel{\scriptscriptstyle \bigtriangleup }{=}I_{l+1} - I_{l} - \Lambda _{{\rm E}, l} \end{equation} \end{document} χE,l=Il+1IlΛE,l
with discrepancy ΛE, l between the l'th real and the l'th detected NN interval; analogously

Eq. 19

[TeX:] \documentclass[12pt]{minimal}\begin{document} \begin{equation} \chi _{{\rm N},l} \stackrel{\scriptscriptstyle \bigtriangleup }{=}I_{l+1} - I_{l} - \Lambda _{{\rm N}, l} \end{equation} \end{document} χN,l=Il+1IlΛN,l
in the case of NIRS. Let [TeX:] \documentclass[12pt]{minimal}\begin{document}$\sigma ^2_{\chi _{E}}$\end{document} σχE2 be the variance of χE, l and let [TeX:] \documentclass[12pt]{minimal}\begin{document}$\sigma ^2_{\chi }$\end{document} σχ2 be the variance of I l + 1I l. With respect to Eq. 18,

Eq. 20

[TeX:] \documentclass[12pt]{minimal}\begin{document} \begin{equation} \sigma ^2_{\chi _{E}}=\sigma ^2_{\chi } + \sigma _{\Lambda _{E}}^{2}. \end{equation} \end{document} σχE2=σχ2+σΛE2.
Likewise, let [TeX:] \documentclass[12pt]{minimal}\begin{document}$\sigma ^2_{\chi _{N}}$\end{document} σχN2 be the variance of χN, l. With respect to Eq. 19,

Eq. 21

[TeX:] \documentclass[12pt]{minimal}\begin{document} \begin{equation} \sigma ^2_{\chi _{N}}=\sigma ^2_{\chi } + \sigma _{\Lambda _{N}}^{2}. \end{equation} \end{document} σχN2=σχ2+σΛN2.
From [TeX:] \documentclass[12pt]{minimal}\begin{document}$\sigma _{\Lambda _{E}} < \sigma _{\Lambda _{N}}$\end{document} σΛE<σΛN , Eqs. 20, 21, it follows that [TeX:] \documentclass[12pt]{minimal}\begin{document}$\sigma _{\chi _{E}} < \sigma _{\chi _{N}}$\end{document} σχE<σχN ; this motivates [TeX:] \documentclass[12pt]{minimal}\begin{document}$\hat{\mu }_3 < 0$\end{document} μ̂3<0 and [TeX:] \documentclass[12pt]{minimal}\begin{document}$\hat{\mu }_4 < 0$\end{document} μ̂4<0 in Fig. 6, since S ECG is an estimator of [TeX:] \documentclass[12pt]{minimal}\begin{document}$\sigma _{\chi _{E}}$\end{document} σχE and S EMD and S PEMAPS are estimators of [TeX:] \documentclass[12pt]{minimal}\begin{document}$\sigma _{\chi _{N}}$\end{document} σχN

In Fig. 5, according to Eqs. 18, 19, the quantity on the y-axes can be modeled as χE, l − χN, l = ΛN, l − ΛE, l. The distribution of χE, l − χN, l has standard deviation [TeX:] \documentclass[12pt]{minimal}\begin{document}$\smash{\sqrt{\sigma _{\Lambda _{E}}^{2} + \sigma _{\Lambda _{N}}^{2}}} \break \approx 0.00518$\end{document} σΛE2+σΛN20.00518 s. The latter may be seen as a lower bound for [TeX:] \documentclass[12pt]{minimal}\begin{document}$\hat{\sigma }_1$\end{document} σ̂1 and [TeX:] \documentclass[12pt]{minimal}\begin{document}$\hat{\sigma }_2$\end{document} σ̂2 , since the error in estimating NN intervals is not only caused by the sampling.

Appendix B: RMSSD and Sampling

In this section, based on the error induced by the sampling, [TeX:] \documentclass[12pt]{minimal}\begin{document}$\hat{\mu }_5 < 0$\end{document} μ̂5<0 and [TeX:] \documentclass[12pt]{minimal}\begin{document}$\hat{\mu }_6 < 0$\end{document} μ̂6<0 in Fig. 7 are motivated.

By expanding the squared term in Eq. 2, it follows that the expectation value

Eq. 22

[TeX:] \documentclass[12pt]{minimal}\begin{document} \begin{eqnarray} E[R^2]=2\sigma ^2 + 2\mu ^2 - \frac{2}{L-1}\sum _{l=1}^{L-1}E\big [\chi _{l+1} \cdot \chi _{l}\big ] \end{eqnarray} \end{document} E[R2]=2σ2+2μ22L1l=1L1Eχl+1·χl
with variance σ2 and mean μ of χl.

Let [TeX:] \documentclass[12pt]{minimal}\begin{document}$\mu _{\chi _{E}}$\end{document} μχE be the mean of χE, l in Eq. 18; let [TeX:] \documentclass[12pt]{minimal}\begin{document}$\mu _{\chi _{N}}$\end{document} μχN be the mean of χN, l in Eq. 19; let μχ be the mean of I n + 1I n. Since the means of ΛE, l and ΛN, l are 0,

Eq. 23

[TeX:] \documentclass[12pt]{minimal}\begin{document} \begin{eqnarray} \mu _{\chi _{E}}=\mu _{\chi _{N}}=\mu _{\chi }. \end{eqnarray} \end{document} μχE=μχN=μχ.

After replacing χl by χE, l, and thus σ2 and μ by [TeX:] \documentclass[12pt]{minimal}\begin{document}$\sigma ^2_{\chi _{E}}$\end{document} σχE2 and [TeX:] \documentclass[12pt]{minimal}\begin{document}$\mu _{\chi _{E}}$\end{document} μχE , in Eq. 22 and using Eqs. 20, 23, the expectation value of the square of R ECG in Fig. 7 is given as

Eq. 24

[TeX:] \documentclass[12pt]{minimal}\begin{document} \begin{eqnarray} E\big [R_{\rm ECG}^2\big ] &=& 2\big(\sigma ^2_{\chi } + \sigma _{\Lambda _{E}}^{2}\big) + 2\mu ^2_{\chi _{E}} - \ldots \nonumber \\ &&\frac{2}{L-1}\sum _{l=1}^{L-1}E\big [\chi _{{E},l+1} \cdot \chi _{{E},l}\big ]. \end{eqnarray} \end{document} ERECG2=2σχ2+σΛE2+2μχE2...2L1l=1L1EχE,l+1·χE,l.
By (i) substituting χE, l and χE, l + 1 in Eq. 24 according to Eq. 18, (ii) then expanding their product, (iii) assuming that ΛE, l is independent of I l and I l + 1 and (iv) using EE, l] = 0, it follows that

Eq. 25

[TeX:] \documentclass[12pt]{minimal}\begin{document} \begin{eqnarray} E\big [R_{\rm ECG}^2\big ] &=& 2\big(\sigma ^2_{\chi } + \sigma _{\Lambda _{E}}^{2}\big) + 2\mu ^2_{\chi _{E}} - \ldots \nonumber \\ && \frac{2}{L-1}\sum _{l=1}^{L-1}E\big [(I_{l+2} - I_{l+1})(I_{l+1} - I_{l})\big ].\nonumber\\ \end{eqnarray} \end{document} ERECG2=2σχ2+σΛE2+2μχE2...2L1l=1L1E(Il+2Il+1)(Il+1Il).
Analogously, it follows from Eqs. 21, 23 that

Eq. 26

[TeX:] \documentclass[12pt]{minimal}\begin{document} \begin{eqnarray} E\big [R_{\rm NIRS}^2\big ] &=& 2\big(\sigma ^2_{\chi } + \sigma _{\Lambda _{N}}^{2}\big) + 2\mu ^2_{\chi _{N}} - \ldots \nonumber \\ && \times\,\frac{2}{L-1}\sum _{l=1}^{L-1}E\big [(I_{l+2} - I_{l+1})(I_{l+1} - I_{l})\big ].\nonumber\\ \end{eqnarray} \end{document} ERNIRS2=2σχ2+σΛN2+2μχN2...×2L1l=1L1E(Il+2Il+1)(Il+1Il).
With respect to Eq. 23, the difference between Eqs. 25, 26 is
[TeX:] \documentclass[12pt]{minimal}\begin{document} \begin{eqnarray*} E\big [R_{\rm ECG}^2\big ] - E\big [R_{\rm NIRS}^2\big ] &= 2\sigma _{\Lambda _{E}}^{2} - 2\sigma _{\Lambda _{N}}^{2}. \nonumber \end{eqnarray*} \end{document} ERECG2ERNIRS2=2σΛE22σΛN2.
From [TeX:] \documentclass[12pt]{minimal}\begin{document}$\sigma _{\Lambda _{E}}^{2} < \sigma _{\Lambda _{N}}^{2}$\end{document} σΛE2<σΛN2 [paragraph after Eq. 16], it follows that
[TeX:] \documentclass[12pt]{minimal}\begin{document} \begin{eqnarray*} E\big [R_{\rm ECG}^2\big ] - E\big [R_{\rm NIRS}^2\big ] < 0\nonumber \end{eqnarray*} \end{document} ERECG2ERNIRS2<0
which motivates [TeX:] \documentclass[12pt]{minimal}\begin{document}$\hat{\mu }_5 < 0$\end{document} μ̂5<0 and [TeX:] \documentclass[12pt]{minimal}\begin{document}$\hat{\mu }_6 < 0$\end{document} μ̂6<0 in Fig. 7.


We thank all subjects for their cooperation and time, Patrick Flandrin and Gabriell Rilling for the fast EMD implementation ( we used in this work, Swiss National Science Foundation (SNF) for funding this project (App. Nos. K-23K0-120727 and 405740-113370), and Christoph Reller from the Signal and Information Processing Laboratory (ISI) at ETH Zurich for constructive discussions.



I. Tachtsidis, C. E. Elwell, T. S. Leung, C. -W. Lee, M. Smith, and D. T. Delpy, “Investigation of cerebral haemodynamics by near infrared spectroscopy in young healthy volunteers reveals posture dependent spontaneous oscillations,” Physiol. Meas., 25 437 –445 (2004). Google Scholar


Task Force of The European Society of Cardiology and The North American Society of Pacing and Electrophysiology, “Heart rate variability, standards of measurement, physiological interpretation, and clinical use,” Eur. Heart J., 17 354 –381 (1996). Google Scholar


R. E. Kleiger, J. P. Miller, J. T. Bigger, and A. J. Moss, “Decreased heart rate variability and its association with increased mortality after acute myocardial infarction,” Am. J. Cardiol., 59 256 –262 (1987). Google Scholar


M. Malik, T. Farrell, T. Cripps, and A. J. Camm, “Heart rate variability in relation to prognosis after myocardial infarction: selection of optimal processing techniques,” Eur. Heart J., 10 1060 –1074 (1989). Google Scholar


J. T. Bigger, J. L. Fleiss, R. C. Steinman, L. M. Rolnitzky, R. E. Kleiger, and J. N. Rottman, “Frequency domain measures of heart period variability and mortality after myocardial infarction,” Circulation, 85 164 –171 (1992). Google Scholar


D. Haensse, P. Szabo, D. Brown, J.-C. Fauchère, P. Niederer, H.-U. Bucher, and M. Wolf, “A new multichannel near infrared spectrophotometry system for functional studies of the brain in adults and neonates,” Opt. Exp., 13 4525 –4538 (2005). Google Scholar


I. Trajkovic, C. Reller, M. Wolf, and H. -A. Loeliger, “Modelling and filtering almost periodic signals by time-varying Fourier series with application to near infrared spectroscopy,” 632 –636 (2009). Google Scholar


N. E. Huang, Z. Shen, S. R. Long, M. C. Wu, H. H. Shih, Q. Zheng, N. Yen, C. C. Tung, and H. H. Liu, “The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis,” Proc. R. Soc. Lond. A, 454 903 –995 (1996). Google Scholar


F. Scholkmann, S. Spichtig, T. Muehlemann, and M. Wolf, “How to detect and reduce movement artifacts in near-infrared imaging using moving standard deviation and spline interpolation,” Physiol. Meas., 31 (5), 649 –662 (2010). Google Scholar


J. M. Bland and D. G. Altman, “Statistical methods for assessing agreement between two methods of clinical measurement,” Lancet, 327 (8476), 30+7–310 (1986). Google Scholar


M. Merri, D. C. Farden, J. G. Mottley, and E. L. Titlebaum, “Sampling frequency of the electrocardiogram for spectral analysis of the heart rate variability,” IEEE Trans. Biomed. Eng., 37 (1), 99 –106 (1990). Google Scholar


P. K. Stein, M. S. Bosner, R. E. Kleiger, and B. M. Conger, “Heart rate variability: a measure of cardiac autonomic tone,” Ame. Heart J., 127 (5), 1376 –1381 (1994). Google Scholar


R. Rauh, R. Limley, R. Bauer, M. Radespiel-Troger, and M. Mueck-Weymann, “Comparison of heart rate variability and pulse rate variability detected with photoplethysmography,” Proc. SPIE, 5474 115 –126 (2004). Google Scholar


S. Lu, H. Zhao, K. Ju, K. Shin, M. Lee, K. Shelley, and K. H. Chon, “Can photoplethysmography variability serve as an alternative approach to obtain heart rate variability information?,” J. Clin. Monito. Comput., 22 (1), 23 –29 (2008). Google Scholar


N. Selvaraj, A. Jaryal, J. Santhosh, K. K. Deepak, and S. Anand, “Assessment of heart rate variability derived from finger-tip photoplethysmography as compared to electrocardiography,” J. Med. Eng. Technol., 32 (6), 479 –484 (2008). Google Scholar


G. Lu, F. Yang, J. A. Taylor, and J. F. Stein, “A comparison of photoplethysmography and ECG recording to analyse heart rate variability in healthy subjects,” J. Med. Eng. Techn., 33 (8), 634 –641 (2009). Google Scholar


M. Wolf, M. Ferrari, and V. Quaresima, “Progress of near-infrared spectroscopy and topography for brain and muscle clinical applications,” J. Biomed. Opt., 12 062104 (2007). Google Scholar
©(2011) Society of Photo-Optical Instrumentation Engineers (SPIE)
Ivo Trajkovic, Felix Scholkmann, and Martin Wolf "Estimating and validating the interbeat intervals of the heart using near-infrared spectroscopy on the human forehead," Journal of Biomedical Optics 16(8), 087002 (1 August 2011).
Published: 1 August 2011


Back to Top