1.

Introduction

Near-infrared spectroscopy (NIRS) potentially measures changes in oxy- and deoxyhemoglobin¹ 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,² e.g., a relationship between the autonomic nervous system's activity and cardiovascular mortality^{3, 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)

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.

2.

Measurement

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.

3.

Instrumentation

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

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;² 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.

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.

4.

Approach Based on PEMAPS

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

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

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

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} $\mathit{y}\stackrel{△}{=}(y_1,...,y_N)$ , 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} $\mathit{\Theta }\stackrel{△}{=}({\Theta }_1,...,{\Theta }_N)$ and coefficient matrix

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

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} ${\widehat{\mathit{A}}}_{k-1,-},...,{\widehat{\mathit{A}}}_{0,-}$ , and [TeX:] \documentclass[12pt]{minimal}\begin{document}$\hat{\bm \Theta }$\end{document} $\widehat{\mathit{\Theta }}$ from the previous iteration as

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} ${\mathit{A}}_0$ estimator” independently estimates the slow component [TeX:] \documentclass[12pt]{minimal}\begin{document}${\bm A}_{0,-}$\end{document} ${\mathit{A}}_{0,-}$ by means of Eq. 10. Since for this [TeX:] \documentclass[12pt]{minimal}\begin{document}$\hat{\bm \Theta }$\end{document} $\widehat{\mathit{\Theta }}$ is not needed, it is a one-time procedure based on [TeX:] \documentclass[12pt]{minimal}\begin{document}${\bm y}$\end{document} $\mathit{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} ${\mathit{A}}_{1,-}$ ; 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} ${\tilde{A}}_1$ is calculated by the “initial A ₁ 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} $\mathit{y}-{\widehat{\mathit{A}}}_{0,-}$ .

Based on [TeX:] \documentclass[12pt]{minimal}\begin{document}${\hat{\bm A}}_{0,-}$\end{document} ${\widehat{\mathit{A}}}_{0,-}$ and [TeX:] \documentclass[12pt]{minimal}\begin{document}$\tilde{A}_{1}$\end{document} ${\tilde{A}}_1$ , the “phase estimator” calculates [TeX:] \documentclass[12pt]{minimal}\begin{document}${\hat{\bm \Theta}}$\end{document} $\widehat{\mathit{\Theta }}$ 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} $\widehat{\mathit{A}}$ . 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} $\widehat{\mathit{\Theta }}$ and [TeX:] \documentclass[12pt]{minimal}\begin{document}$\hat{\bm {A}}$\end{document} $\widehat{\mathit{A}}$ directly in Eq. 4 (done in the “reconstruction of [TeX:] \documentclass[12pt]{minimal}\begin{document}$\hat{x}$\end{document} $\widehat{x}$ ”-block) and obtain a first estimate [TeX:] \documentclass[12pt]{minimal}\begin{document}$\hat{\bm x}$\end{document} $\widehat{\mathit{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} $\widehat{\mathit{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],

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} ${\widehat{\mathit{A}}}_{1,-},...,{\widehat{\mathit{A}}}_{K,-}$ and [TeX:] \documentclass[12pt]{minimal}\begin{document}${\hat{\bm \Theta}} = (\hat{\Theta }_1, \ldots, \hat{\Theta }_{1000})$\end{document} $\widehat{\mathit{\Theta }}=({\widehat{\Theta }}_1,...,{\widehat{\Theta }}_{1000})$ in Eq. 4. (c) shows [TeX:] \documentclass[12pt]{minimal}\begin{document}${\hat{\bm \Theta}}$\end{document} $\widehat{\mathit{\Theta }}$ . (d) depicts [TeX:] \documentclass[12pt]{minimal}\begin{document}$\hat{\Theta }^{\prime }_1, \ldots, \hat{\Theta }^{\prime }_{999}$\end{document} ${\widehat{\Theta }}_1^{\prime },...,{\widehat{\Theta }}_{999}^{\prime }$ with [TeX:] \documentclass[12pt]{minimal}\begin{document}$\hat{\Theta }_{n}^{\prime }=\hat{\Theta }_{n+1} - \hat{\Theta }_{n}$\end{document} ${\widehat{\Theta }}_n^{\prime }={\widehat{\Theta }}_{n+1}-{\widehat{\Theta }}_n$ . Since with increasing n, the values [TeX:] \documentclass[12pt]{minimal}\begin{document}$\hat{\Theta }_n$\end{document} ${\widehat{\Theta }}_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} ${\widehat{\Theta }}_n^{\prime }$ 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).

5.

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:

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

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} $\widehat{\mathit{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} ${\widehat{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} $\widehat{\mathit{x}}$ , if it is the m'th entry of [TeX:] \documentclass[12pt]{minimal}\begin{document}${\hat{\bm x}}$\end{document} $\widehat{\mathit{x}}$ for which the conditions

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.

6.

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, ξ₁ in Eq. 12 was chosen individually. For all subjects/ECG signals, ξ₂ = 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} $\stackrel{\gamma }{=}$ ”-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} ${\widehat{\mathit{A}}}_{0,-}$ was smoothed by feeding it back to the [TeX:] \documentclass[12pt]{minimal}\begin{document}$\bm A_0$\end{document} ${\mathit{A}}_0$ estimator (Fig. 3) as the new input signal and re-estimating [TeX:] \documentclass[12pt]{minimal}\begin{document}$\bm A_{0,-}$\end{document} ${\mathit{A}}_{0,-}$ . For every subject, this procedure was performed twice successively.

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

7.

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} ${\mathit{\chi }}_{\mathrm{ECG}}^{\left(j\right)}\stackrel{△}{=}({\chi }_{E,1}^{\left(j\right)},...,{\chi }_{E,W_j}^{\left(j\right)})$ 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} ${\mathit{\chi }}_{\mathrm{NIRS}}^{\left(j\right)}\stackrel{△}{=}({\chi }_{N,1}^{\left(j\right)},...,{\chi }_{N,W_j}^{\left(j\right)})$ estimated from NIRS using PEMAPS and EMD, Table 1 shows, for each subject, i.e., for j = 1, …, 11, the cross-correlation coefficient

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.

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.

Let [TeX:] \documentclass[12pt]{minimal}\begin{document}$\bm \chi _{\rm ECG}$\end{document} ${\mathit{\chi }}_{\mathrm{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} ${\mathit{\chi }}_{\mathrm{PEMAPS}}$ and [TeX:] \documentclass[12pt]{minimal}\begin{document}$\bm \chi _{\rm EMD}$\end{document} ${\mathit{\chi }}_{\mathrm{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} ${\mathit{\chi }}_{\mathrm{ECG}}$ , [TeX:] \documentclass[12pt]{minimal}\begin{document}$\bm \chi _{\rm PEMAPS}$\end{document} ${\mathit{\chi }}_{\mathrm{PEMAPS}}$ , and [TeX:] \documentclass[12pt]{minimal}\begin{document}$\bm \chi _{\rm EMD}$\end{document} ${\mathit{\chi }}_{\mathrm{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} ${\mathit{\chi }}_{\mathrm{ECG}}$ and [TeX:] \documentclass[12pt]{minimal}\begin{document}$\bm \chi _{\rm PEMAPS}$\end{document} ${\mathit{\chi }}_{\mathrm{PEMAPS}}$ on the one hand, and [TeX:] \documentclass[12pt]{minimal}\begin{document}$\bm \chi _{\rm ECG}$\end{document} ${\mathit{\chi }}_{\mathrm{ECG}}$ and [TeX:] \documentclass[12pt]{minimal}\begin{document}$\bm \chi _{\rm EMD}$\end{document} ${\mathit{\chi }}_{\mathrm{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} $\widehat{\mu }$ and the standard deviation [TeX:] \documentclass[12pt]{minimal}\begin{document}$\hat{\sigma }$\end{document} $\widehat{\sigma }$ 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} ${\widehat{\mu }}_1=0.00008$ s, [TeX:] \documentclass[12pt]{minimal}\begin{document}$\hat{\sigma }_1=0.00755$\end{document} ${\widehat{\sigma }}_1=0.00755$ s, [TeX:] \documentclass[12pt]{minimal}\begin{document}$\hat{\mu }_2\break=0.00008$\end{document} ${\widehat{\mu }}_2=0.00008$ s and [TeX:] \documentclass[12pt]{minimal}\begin{document}$\hat{\sigma }_2=0.01133$\end{document} ${\widehat{\sigma }}_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} ${\mathit{\chi }}_{\mathrm{ECG}}$ , [TeX:] \documentclass[12pt]{minimal}\begin{document}$\bm \chi _{\rm PEMAPS}$\end{document} ${\mathit{\chi }}_{\mathrm{PEMAPS}}$ , [TeX:] \documentclass[12pt]{minimal}\begin{document}$\bm \chi _{\rm EMD}$\end{document} ${\mathit{\chi }}_{\mathrm{EMD}}$ , and thus the differences [TeX:] \documentclass[12pt]{minimal}\begin{document}$\bm \chi _{\rm ECG}$\end{document} ${\mathit{\chi }}_{\mathrm{ECG}}$ – [TeX:] \documentclass[12pt]{minimal}\begin{document}$\bm \chi _{\rm PEMAPS}$\end{document} ${\mathit{\chi }}_{\mathrm{PEMAPS}}$ and [TeX:] \documentclass[12pt]{minimal}\begin{document}$\bm \chi _{\rm ECG}$\end{document} ${\mathit{\chi }}_{\mathrm{ECG}}$ − [TeX:] \documentclass[12pt]{minimal}\begin{document}$\bm \chi _{\rm EMD}$\end{document} ${\mathit{\chi }}_{\mathrm{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} ${\mathit{\chi }}_{\mathrm{ECG}}$ and [TeX:] \documentclass[12pt]{minimal}\begin{document}$\bm \chi _{\rm PEMAPS}$\end{document} ${\mathit{\chi }}_{\mathrm{PEMAPS}}$ or [TeX:] \documentclass[12pt]{minimal}\begin{document}$\bm \chi _{\rm ECG}$\end{document} ${\mathit{\chi }}_{\mathrm{ECG}}$ and [TeX:] \documentclass[12pt]{minimal}\begin{document}$\bm \chi _{\rm EMD}$\end{document} ${\mathit{\chi }}_{\mathrm{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} ${\widehat{\sigma }}_1/0.84\mathrm{s}=0.9$ %, or [TeX:] \documentclass[12pt]{minimal}\begin{document}${\hat{\sigma }_2}/{0.84 {\rm s}}\break =1.35$\end{document} ${\widehat{\sigma }}_2/0.84\mathrm{s}=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} ${\widehat{\sigma }}_3/0.08\mathrm{s}=0.4$ % or [TeX:] \documentclass[12pt]{minimal}\begin{document}${\hat{\sigma }_4}/{0.08{\rm s}}=0.88$\end{document} ${\widehat{\sigma }}_4/0.08\mathrm{s}=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} ${\widehat{\sigma }}_5/0.038\mathrm{s}=2.55$ % and [TeX:] \documentclass[12pt]{minimal}\begin{document}${\hat{\sigma }_6}/{0.038{\rm s}}=5.16$\end{document} ${\widehat{\sigma }}_6/0.038\mathrm{s}=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} ${\widehat{\mu }}_3=-0.00084$ s, [TeX:] \documentclass[12pt]{minimal}\begin{document}$\hat{\sigma }_3=0.00032$\end{document} ${\widehat{\sigma }}_3=0.00032$ s, [TeX:] \documentclass[12pt]{minimal}\begin{document}$\hat{\mu }_4=-0.0014$\end{document} ${\widehat{\mu }}_4=-0.0014$ s, and [TeX:] \documentclass[12pt]{minimal}\begin{document}$\hat{\sigma }_4\break=0.0007\,{\rm s}$\end{document} ${\widehat{\sigma }}_4=0.0007\mathrm{s}$ . 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} ${\widehat{\mu }}_5=-0.00247$ s, [TeX:] \documentclass[12pt]{minimal}\begin{document}$\hat{\sigma }_5=0.00097$\end{document} ${\widehat{\sigma }}_5=0.00097$ s, [TeX:] \documentclass[12pt]{minimal}\begin{document}$\hat{\mu }_6=-0.00452$\end{document} ${\widehat{\mu }}_6=-0.00452$ s, and [TeX:] \documentclass[12pt]{minimal}\begin{document}$\hat{\sigma }_6=0.00196$\end{document} ${\widehat{\sigma }}_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} ${\widehat{\mu }}_3,...,{\widehat{\mu }}_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.

8.

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} ${\tilde{A}}_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} ${\widehat{\sigma }}_1<{\widehat{\sigma }}_2$ in Fig. 5, [TeX:] \documentclass[12pt]{minimal}\begin{document}$|\hat{\mu }_3|<|\hat{\mu }_4|$\end{document} $|{\widehat{\mu }}_3|<|{\widehat{\mu }}_4|$ and [TeX:] \documentclass[12pt]{minimal}\begin{document}$\hat{\sigma }_3<\hat{\sigma }_4$\end{document} ${\widehat{\sigma }}_3<{\widehat{\sigma }}_4$ in Fig. 6, and [TeX:] \documentclass[12pt]{minimal}\begin{document}$|\hat{\mu }_5|<|\hat{\mu }_6|$\end{document} $|{\widehat{\mu }}_5|<|{\widehat{\mu }}_6|$ and [TeX:] \documentclass[12pt]{minimal}\begin{document}$\hat{\sigma }_5<\hat{\sigma }_6$\end{document} ${\widehat{\sigma }}_5<{\widehat{\sigma }}_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 ξ₁ 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.