1.

Introduction

In functional near-infrared spectroscopy (fNIRS) studies, cerebral blood volume and blood oxygenation changes are studied as indicators of brain activity. Hemodynamics are quantified based on localized changes in the concentrations of oxyhemoglobin (HbO) and deoxyhemoglobin (HbR). Three different classes of NIRS instrumentation are used to measure cerebral hemodynamics, and each have certain advantages and disadvantages. Time-domain NIRS (TD) and frequency domain NIRS (FD) are capable of providing absolute quantification of hemodynamics.^1–3 TD provides superior depth resolution by tracking the time of flight of photons.⁴ However, the complexity of TD and FD systems make them more costly for scale up to high-density measurement probes⁵ and whole head neuroimaging.⁶ The instrument of choice for high-density neuroimaging applications with NIRS is continuous-wave (CW).

With CW-NIRS systems, the attenuation of light is measured at multiple wavelengths between source and detector optodes that are separated by a few centimeters on the scalp. The change in light absorption or optical density $\mathrm{\Delta }{\mathrm{OD}}_{\lambda }$ for a given wavelength, $\lambda$, is the negative logarithm of the ratio of the detected light at a given time $t_2$ relative to the detected light at a reference time $t_1$

Changes in HbO and HbR can be calculated by solving a set of linear equations based on the modified Beer-Lambert Law (MBLL):⁷

Unlike TD or FD, CW-NIRS cannot be used to determine ${\mathrm{DPF}}_{\lambda }$. Typically, the assumed ${\mathrm{DPF}}_{\lambda }$ values used in CW-NIRS analysis are obtained from prior studies of healthy subjects and range from three to six.^2,8 However, ${\mathrm{DPF}}_{\lambda }$ can vary between subjects due to differences in the anatomical structure and tissue composition of the brain. Morphological changes in the brain induced by aging or neurological disorders can affect the optical properties of brain tissue, altering ${\mathrm{DPF}}_{\lambda }$ values.⁹ Furthermore, ${\mathrm{DPF}}_{\lambda }$ used in the MBLL is assumed constant for a specific wavelength and cannot account for partial volume effects that arise from assuming global changes in the concentrations of HbO and HbR in a brain volume consisting of multiple layers of tissue.¹⁰ Studies have shown that skull thickness and the cerebrospinal fluid (CSF) in the head can affect the mean optical path length and hence ${\mathrm{DPF}}_{\lambda }$ values.^11,12 The reason this is a problem is that systematic errors in ${\mathrm{DPF}}_{\lambda }$ lead to systematic errors in the quantification of cerebral oxygenation levels. In this work, we propose a method for calculating changes in the concentration of HbO and HbR while continuously correcting for systematic errors in ${\mathrm{DPF}}_{\lambda }$.

2.

Theoretical Formalism

CW-NIRS is an increasingly important tool for clinical monitoring and for diagnosis of neurological disorders. However, ${\mathrm{DPF}}_{\lambda }$ values cannot be known in advance in such patients. Efforts have been made to estimate mean optical path lengths in infants¹³ as well as in patients undergoing cardiopulmonary bypass.¹⁴ Having knowledge of ${\mathrm{DPF}}_{\lambda }$ values significantly improves the accuracy of calculated cerebral hemodynamics, which are important in the clinical evaluation of patients with stroke or brain injury.^15,16

When inaccurate ${\mathrm{DPF}}_{\lambda }$ values are used in the MBLL, the calculated concentration changes in $\mathrm{\Delta }\mathrm{HbO}$ and $\mathrm{\Delta }\mathrm{HbR}$ become inaccurate in two ways. The first type of inaccuracy is a loss of the absolute scale on the results; the second type of inaccuracy is called cross talk, which is when $\mathrm{\Delta }\mathrm{HbO}$ become contaminated by $\mathrm{\Delta }\mathrm{HbR}$ and vice versa. From a mathematical perspective, the errors $C$ in $\mathrm{\Delta }\mathrm{HbO}$ and $\mathrm{\Delta }\mathrm{HbR}$ are given by:

Pure scale errors result when $S\ne 1$ and pure cross talk errors occur when $S=1$ and $\mathrm{\Delta }{\mathrm{DPF}}_{\lambda }\ne 0$ for all but one of the wavelengths. In NIRS spectroscopic calculations involving multiple wavelengths, the errors $C$ in $\mathrm{\Delta }\mathrm{HbO}$ and $\mathrm{\Delta }\mathrm{HbR}$ can be computed as:

Several methods based on numerical and analytical approaches are already available for calculating ${\mathrm{DPF}}_{\lambda }$ in simulation. These methods typically rely on modeling light transport in the head using the radiative transport equation (RTE). The diffusion equation, which is derived from the $P_1$ approximation to the RTE^17,18 can accurately describe photon transport in a highly scattering medium. However, in a low scattering media like CSF, the diffusion approximation does not hold. Analytical solutions to the diffusion equation, which can be obtained for simple geometries¹⁹ are difficult to solve in complex geometries like the head. Instead, numerical approaches like the finite element method (FEM) are preferred. Another numerical approach is the FEM based Monte Carlo (MC) technique.^20–22 Although computationally expensive,²³ FEM based MC has the advantage that it can handle complex geometries like the head and can be used to solve the full RTE. The major disadvantages of simulation approaches to calculating ${\mathrm{DPF}}_{\lambda }$ values are inaccuracies in FEM representations of head anatomy and a lack of knowledge of the true optical properties in the mesh elements. Therefore, ${\mathrm{DPF}}_{\lambda }$ values obtained from simulation methods are likely to be different from the actual values, which can lead to cross-talk error in the calculated $\mathrm{\Delta }\mathrm{HbO}$ and $\mathrm{\Delta }\mathrm{HbR}$.^10,22,24 This motivates the need for experimental methods of calculating relative ${\mathrm{DPF}}_{\lambda }$ values from CW-NIRS data.

In the present work, we focus on calculating relative offsets in differential path length factors (${\mathrm{DPF}}_{\lambda }$) at wavelengths of 690, 785, and 830 nm with respect to a ${\mathrm{DPF}}_{\lambda }$ value at 808 nm that is held fixed ($\mathrm{\Delta }{\mathrm{DPF}}_{808}=0$). We apply the extended Kalman filter (EKF) to calculate the $\mathrm{\Delta }{\mathrm{DPF}}_{\lambda }$ values and the changes in concentration of HbO and HbR. Since $\mathrm{\Delta }{\mathrm{DPF}}_{\lambda }$ is a relative correction, the calculated $\mathrm{\Delta }\mathrm{HbO}$ and $\mathrm{\Delta }\mathrm{HbR}$ are corrected relative to one another but not corrected to absolute scales. The absolute accuracies of the calculated $\mathrm{\Delta }\mathrm{HbO}$ and $\mathrm{\Delta }\mathrm{HbR}$ values are determined by the a priori value of ${\mathrm{DPF}}_{\lambda }$ at 808 nm.

3.

Dynamic State Estimation Using the EKF Algorithm

The EKF algorithm operates on a state-space dynamical system that consists of an observation model and a process model. In our case, the observation model is based on the MBLL and is derived by relating $\mathrm{\Delta }{\mathrm{OD}}_{\lambda }$ to the state variables $\mathrm{\Delta }{\mathrm{DPF}}_{\lambda }$, $\mathrm{\Delta }\mathrm{HbO}$, and $\mathrm{\Delta }\mathrm{HbR}$. The subscripts $j$ and $k$ are introduced in the equations below to specify the source-detector pair and the discrete time step, respectively. The MBLL-based observation model in matrix form is:

In Eq. (10), ${\mathbf{X}}_{j,k}$ is a vector containing state variables that we want to calculate using the EKF algorithm:

The measurement transition matrix ${\mathbf{H}}_{j,k}$ determines the dynamic characteristics of the model and is formed from the coefficients of $\mathrm{\Delta }\mathrm{HbO}$ and $\mathrm{\Delta }\mathrm{HbR}$ in the MBLL equation [Eq. (2)] with the addition of $\mathrm{\Delta }{\mathrm{DPF}}_{{\lambda }_i,j,k}$ terms that correspond to the relative changes in the differential path length factor at wavelength ${\lambda }_i$ for source-detector pair $j$ and time step $k$. $\mathrm{\Delta }{\mathrm{DPF}}_{{\lambda }_i,j,k}$ is updated at each time step for ${\lambda }_i$ ($i=690$, 785, 830) and serves as a correction term for the assumed ${\mathrm{DPF}}_{{\lambda }_i,j}$. Three additional rows are also included in ${\mathbf{H}}_{j,k}$ so that we can specify statistical priors on $\mathrm{\Delta }{\mathrm{DPF}}_{{\lambda }_i,j,k}$:

We assume that the expected values of $\mathrm{\Delta }{\mathrm{DPF}}_{{\lambda }_i,j,k}$ are zero. These priors are incorporated in the measurement vector ${\mathbf{Y}}_{j,k}$ along with $\mathrm{\Delta }{\mathrm{OD}}_{{\lambda }_i,j,k}$ obtained for $i=690$, 785, 808, 830:

The observation model [Eq. (10)] contains an observation noise term ${\mathbf{v}}_{j,k}$, which is assumed to be zero mean white Gaussian noise (WGN) with covariance ${\mathbf{R}}_j$:

The process model describes the dynamics of the state variables and is assumed to be stochastic:

Equation (15) contains an innovation term, also known as process noise, defined in vector ${\mathbf{w}}_{j,k}$, that is assumed to be drawn from a zero mean multivariate normal distribution with covariance ${\mathbf{Q}}_j$:

The state transition matrix ${\mathbf{F}}_j$ in Eq. (15) is used to describe the relationship between the state variables $\mathbf{X}$ at time step $k$ and $k+1$. In our process model, ${\mathbf{F}}_j$ is a 5 by 5 identity matrix.

The input and output state transition matrices ${\mathbf{F}}_j$ and ${\mathbf{H}}_j$, along with the output vector ${\mathbf{Y}}_{j,k}$ are used in the EKF algorithm to estimate the state variables in ${\mathbf{X}}_{j,k}$. The EKF performs state estimation through an iterative prediction-correction scheme. In the prediction step, the most recent estimated state variables ${\widehat{\mathbf{X}}}_{j,k-1}$ and error covariance matrix ${\mathbf{P}}_{j,k-1}$ are projected forward in time to compute a predicted or a prior estimate ${\widehat{\mathbf{X}}}_{j,k}^{-}$ at the current time step $k$:

The predicted error covariance matrix is given by:

In the correction stage, a posterior estimate of the states ${\widehat{\mathbf{X}}}_{j,k}$ is calculated based on a linear combination of the prior estimate ${\widehat{\mathbf{X}}}_{j,k}^{-}$ and a weighted difference between the predicted measurement ${\mathbf{H}}_j{\widehat{\mathbf{X}}}_{j,k}^{-}$ and the actual measurement ${\mathbf{Y}}_{j,k}$:

The optimal Kalman gain weighting matrix, ${\mathbf{K}}_{j,k}$, is given by:

The computed ${\widehat{\mathbf{X}}}_{j,k}$ and ${\mathbf{P}}_{j,k}$ are then used in the next predictor step.

4.

Methods

5.

Results

The EKF and WLSQ solutions for the simulated NIRS data are shown in Fig. 1. The solutions are displayed after low-pass filtering with a 2.5 Hz cutoff frequency using a Butterworth filter of order 4. The filter was applied both in the forward and reverse direction to compensate for phase lags introduced by the filter. Figure 1(a) shows the EKF and WLSQ solutions for $\mathrm{\Delta }\mathrm{HbO}$, which agree closely with the known $\mathrm{\Delta }\mathrm{HbO}$ sequence as evidenced by $R$² statistics of 0.99 for EKF and 0.96 for WLSQ. The residual errors between the known $\mathrm{\Delta }\mathrm{HbO}$ sequence and the solution obtained using EKF and WLSQ are shown in Fig. 1(b). Similarly, Fig. 1(c) and 1(d) shows $\mathrm{\Delta }\mathrm{HbR}$ computed from EKF and WLSQ and their corresponding residual errors in $\mathrm{\Delta }\mathrm{HbR}$. It is evident in Fig. 1(c) that $\mathrm{\Delta }\mathrm{HbR}$ calculated using EKF agrees almost perfectly with the known sequence ($R^2=0.99$). In case of WLSQ, there are large discrepancies between the calculated $\mathrm{\Delta }\mathrm{HbR}$ and the known sequence ($R^2=0.66$). We can observe from the residual error plots for EKF that the magnitude of the error decreases rapidly during the first 2 s. This is approximately the amount of time it takes for the $\mathrm{\Delta }{\mathrm{DPF}}_{\lambda }$ values to converge to the known values in the simulation [Fig. 1(e)]. In case of WLSQ, the residual errors do not decrease over time [Fig. 1(b) and 1(d)]; they tend to show some degree of serial correlation and contain oscillatory components.

Fig. 1

EKF and WLSQ solutions for simulated NIRS data: (a) $\mathrm{\Delta }\mathrm{HbO}$ calculated using EKF and WLSQ and the true known sequence; (b) residual errors in $\mathrm{\Delta }\mathrm{HbO}$ for EKF and WLSQ; (c) $\mathrm{\Delta }\mathrm{HbR}$ calculated using EKF and WLSQ; (d) residual errors for $\mathrm{\Delta }\mathrm{HbR}$; (e) $\mathrm{\Delta }{\mathrm{DPF}}_{\lambda }$ terms calculated using EKF. The true $\mathrm{\Delta }{\mathrm{DPF}}_{\lambda }$ errors are shown in dotted lines.

In the simulated data, we had explicit prior knowledge of the measurement noise variance ${\sigma }_{\mathrm{\Delta }\mathrm{OD},{\lambda }_i}^2$ based on the simulated WGN that was introduced in $\mathrm{\Delta }{\mathrm{OD}}_{\lambda }$. The process noise variances ${\sigma }_{\mathrm{\Delta }\mathrm{HbO}}^2$ and ${\sigma }_{\mathrm{\Delta }\mathrm{HbR}}^2$ were found to be approximately $1\times {10}^{-12}{(\mu \mathrm{M})}^2$. Results from tuning the process noise variance in $\mathrm{\Delta }{\mathrm{DPF}}_{\lambda }$ are shown in Fig. 2. A minimum value in the percent mean RSS error occurs at ${\sigma }_{\mathrm{\Delta }\mathrm{DPF},Q}^2=1\times {10}^{-6}$. This value of ${\sigma }_{\mathrm{\Delta }\mathrm{DPF},Q}^2$ was used in the EKF applied to both the simulated and experimental NIRS data.

Fig. 2

Percent mean RSS error between known $\mathrm{\Delta }\mathrm{HbO}$ and $\mathrm{\Delta }\mathrm{HbR}$ and the solutions computed using EKF for nine evaluated levels of ${\sigma }_{\mathrm{\Delta }\mathrm{DPF},Q}^2$. The minimum RSS error corresponds to a process noise variance of $1\times {10}^{-6}$.

Figure 3(a) shows the percent cross-talk error in the solutions obtained from EKF and WLSQ. The percent cross-talk error in $\mathrm{\Delta }\mathrm{HbO}$ is significantly smaller in EKF relative to that found for WLSQ across all frequency bands (paired $t$-test, $p<0.05$). Although the cross-talk error in $\mathrm{\Delta }\mathrm{HbR}$ is slightly higher for EKF in the HF and CARDIAC band, the corresponding values in WLSQ are much higher at 45.9% and 35.1%, respectively.

Fig. 3

Simulated NIRS data residual analysis: (a) percent cross-talk error in $\mathrm{\Delta }\mathrm{HbO}$ and $\mathrm{\Delta }\mathrm{HbR}$ solution obtained using EKF ad WLSQ. Error bars show standard error; (b) mean sum squared cross-correlation coefficients for LF, HF, and CARDIAC filtered residuals for EKF and WLSQ.

Figure 3(b) shows the mean of the sum squared cross correlation coefficients obtained for the LF, HF, and CARDIAC frequency bands for EKF and WLSQ applied to simulated NIRS data. We can observe that the EKF has reduced mean sum squared cross-correlation coefficient values in the LF band relative to WLSQ. In case of the HF and CARDIAC bands, the mean squared cross-correlation coefficient values for EKF and WLSQ are comparable.

The results from EKF and WLSQ applied to experimental NIRS data segments are shown in Fig. 4. $\mathrm{\Delta }\mathrm{HbO}$ and $\mathrm{\Delta }\mathrm{HbR}$ computed using EKF agree closely with the WLSQ results [Fig. 4(a) and 4(b)]. Figure 4(c) shows the $\mathrm{\Delta }{\mathrm{DPF}}_{\lambda }$ values calculated for wavelengths of 690, 785, and 830 nm. Large deviations in the $\mathrm{\Delta }{\mathrm{DPF}}_{\lambda }$ values occur between 125 and 140 s. During this time interval, the difference between the EKF and WLSQ solution for $\mathrm{\Delta }\mathrm{HbO}$ and $\mathrm{\Delta }\mathrm{HbR}$ is also large [Fig. 4(d)]. The $R^2$ values for the plotted time interval are 0.97 for $\mathrm{\Delta }\mathrm{HbO}$ and 0.98 for $\mathrm{\Delta }\mathrm{HbR}$.

Fig. 4

EKF and WLSQ results obtained for experimental NIRS data collected from subject 3 (a) $\mathrm{\Delta }\mathrm{HbO}$ estimated using WLSQ and EKF (b) corresponding $\mathrm{\Delta }\mathrm{HbR}$ estimated using WLSQ and EKF; (c) EKF estimates of $\mathrm{\Delta }{\mathrm{DPF}}_{\lambda }$ at 690, 785, and 830 nm; (d) difference between WLSQ and EKF hemodynamic signals.

Figure 5 shows $\mathrm{\Delta }\mathrm{HbO}$ and $\mathrm{\Delta }\mathrm{HbR}$ results in the LF, HF and CARDIAC bands for EKF and WLSQ. The differences between the EKF and WLSQ solutions are also plotted.

Fig. 5

Band pass filtered $\mathrm{\Delta }\mathrm{HbO}$ (upper row) and $\mathrm{\Delta }\mathrm{HbR}$ (lower row) sequences from EKF and WLSQ for subject 3: (a) and (d) LF (0.04 to 0.15 Hz); (b) and (e) HF (0.15 to 0.4 Hz); (c) and (f) CARDIAC (0.4 to 2.0 Hz).

Figure 6 shows the residuals obtained for WLSQ and EKF and their autocorrelation sequences for subject 3. The WLSQ and EKF residuals shown in Fig. 6(a) and 6(b) appears to contain physiological signal components such as cardiac oscillations in the EKF residuals. There also appear to be low frequency components in both the WLSQ and EKF residuals. The autocorrelation structures shown in Fig. 6(c) and 6(d) shows reduced serial correlation at larger lags in EKF compared to WLSQ. The autocorrelation values for the 690 nm residual appear to be substantially smaller in EKF compared to WLSQ.

Fig. 6

(a) WLSQ residuals for wavelengths 690, 758, 808, and 830 nm; (b) EKF residuals for the four wavelengths; (c) autocorrelation coefficients for residuals from WLSQ with the 95% confidence interval shown in solid gray lines (d) autocorrelation coefficient of residuals from EKF with the 95% confidence interval shown in solid gray lines.

Table 1 summarizes the mean percentages of points in the autocorrelation functions for WLSQ and EKF that fall outside of the 95% confidence intervals. The percentages of points outside the confidence bounds for EKF are consistently smaller than for WLSQ.

Table 1

Mean percentages of points that fall outside of the 95% confidence intervals for the EKF and WLSQ autocorrelation functions.

Figure 7 compares the mean squared cross-correlation coefficient values obtained in the analysis of EKF and WLSQ residuals. All three subjects have smaller mean squared correlation coefficient values for EKF relative to those found for WLSQ (see Table 2 for ANOVA). The lowest mean squared cross-correlation coefficients were found for the HF band in EKF.

Fig. 7

Mean squared cross-correlation coefficients of residuals from the analysis of experimental data from three subjects. Results are shown for LF, HF and CARDIAC frequency bands with standard error bars.

Table 2

NOVA of mean squared cross-correlation values shown in Fig. 7.

ANOVA results are summarized in the Table 2. Algorithm (EKF versus WLSQ) contributes significantly ($p<0.001$) in all three frequency bands. Therefore, we can reject the null hypothesis that the mean squared cross-correlation coefficient values found for EKF and WLSQ are the same at those frequency bands. It is interesting to note that the wavelength factor also contributes significantly to the variability in the mean squared cross-correlation coefficient values indicating that the EKF corrections vary significantly with wavelength.

Figure 8 shows the mean percent deviations between EKF and WLSQ solutions in the LF, HF and CARDIAC frequency bands. For each frequency band, the mean percent deviations in $\mathrm{\Delta }\mathrm{HbO}$ and $\mathrm{\Delta }\mathrm{HbR}$ are significantly different from zero (two-tailed $t$-test: $p<0.001$). The mean deviation over all frequency bands was found to be 4.04% for $\mathrm{\Delta }\mathrm{HbO}$ and 4.92% for $\mathrm{\Delta }\mathrm{HbR}$. In the LF band, the mean percent deviations are comparable for both $\mathrm{\Delta }\mathrm{HbO}$ and $\mathrm{\Delta }\mathrm{HbR}$. The mean percent signal deviation is smaller in $\mathrm{\Delta }\mathrm{HbO}$ compared to $\mathrm{\Delta }\mathrm{HbR}$ in the HF and CARDIAC bands ($p<0.05$).

Fig. 8

Mean percent deviations between EKF and WLSQ solutions filtered in the LF, HF, and CARDIAC frequency bands shown with 95% confidence intervals from a paired $t$-test.

6.

Discussion and Conclusion

The proposed EKF algorithm is able to calculate $\mathrm{\Delta }\mathrm{HbO}$ and $\mathrm{\Delta }\mathrm{HbR}$ from simulated NIRS data by continuously updating $\mathrm{\Delta }{\mathrm{DPF}}_{\lambda }$ values. We observed in simulation that the residual errors in $\mathrm{\Delta }\mathrm{HbO}$ and $\mathrm{\Delta }\mathrm{HbR}$ decrease rapidly within the first 2 s as the EKF tracked the errors introduced in ${\mathrm{DPF}}_{\lambda }$. In the analysis of experimental data, we observed reduced autocorrelation of residuals for EKF compared to WLSQ. We believe that the superior performance of EKF is achieved due to weighting of the state variable innovations by the optimal Kalman gain at each time step. The application of the Kalman gain in this way helps to minimize the posterior error covariance over time. In WLSQ, the measurement data is weighted by constant factors over all time steps; so even though the solution is optimal in the least squares sense, WLSQ lacks the flexibility of EKF to account for time variation in ${\mathrm{DPF}}_{\lambda }$. As a result, EKF outperforms WLSQ in the estimation of $\mathrm{\Delta }\mathrm{HbO}$ and $\mathrm{\Delta }\mathrm{HbR}$.

We observed reduced cross-correlations in the residuals in wavelength pairs for EKF in the LF, HF and CARDIAC bands compared to WLSQ, which is consistent with the reduction of cross-talk error. Convergence of the EKF state estimates to the true states is dependent on prior knowledge of the process noise and measurement noise variance statistics. In the proposed method, we demonstrate how to obtain prior estimates of those statistics based on minimizing the residual error in EKF.

We also quantified the differences in the hemodynamic signals calculated using EKF and WLSQ in physiologically relevant frequency bands. We found significant differences in the $\mathrm{\Delta }\mathrm{HbR}$ sequences in the LF and HF bands. The average of differences between EKF and WLSQ solutions across all frequency bands was 3.75% for $\mathrm{\Delta }\mathrm{HbO}$ and 4.57% for $\mathrm{\Delta }\mathrm{HbR}$. We believe that this difference is due to a reduction of cross-talk error in $\mathrm{\Delta }\mathrm{HbO}$ and $\mathrm{\Delta }\mathrm{HbR}$ obtained using EKF. This is supported by the observation that the EKF solution obtained from simulated NIRS data had reduced cross-talk error relative to the WLSQ solutions. In the simulation, the errors in $\mathrm{\Delta }\mathrm{HbO}$ and $\mathrm{\Delta }\mathrm{HbR}$ from the known sequences were significantly lower in the EKF results compared to the WLSQ results.

The proposed EKF method offers several advantages over other methods that may be used to calculate ${\mathrm{DPF}}_{\lambda }$. Unlike other studies where ${\mathrm{DPF}}_{\lambda }$ is estimated based on the assumption that it is constant over time,^2,25,31 the proposed EKF algorithm continuously updates ${\mathrm{DPF}}_{\lambda }$ and it can account for variation in ${\mathrm{DPF}}_{\lambda }$ due to differences in anatomical structure and tissue composition of the brain. In cases where relative calibration of $\mathrm{\Delta }\mathrm{HbO}$ and $\mathrm{\Delta }\mathrm{HbR}$ is desired, the proposed method eliminates the need to measure ${\mathrm{DPF}}_{\lambda }$ using TD or FD instruments or applying simulation techniques discussed in Sec. 2.

The proposed EKF algorithm has potential applications in the clinical investigation of brain injury. A number of recent studies have used NIRS to monitor patients during stroke rehabilitation and recovery.^32–34 A limitation in these studies is the inability to accurately quantify cerebral hemodynamics in a region of interest. In the above referenced stroke studies, ${\mathrm{DPF}}_{\lambda }$ is assumed to be constant over time and the studies did not account for possible variation due to the stroke. However, ${\mathrm{DPF}}_{\lambda }$ could deviate substantially in stroke because of the high absorption coefficient of hemoglobin. As a result, the calculated $\mathrm{\Delta }\mathrm{HbO}$ and $\mathrm{\Delta }\mathrm{HbR}$ may not have accurately reflected the hemodynamics in those patients. The proposed algorithm could also be applied to other clinical diagnostic applications of NIRS such as breast cancer detection.³⁵ The optical properties of cancerous tissue are different from those of healthy tissue in the breast. Accurate assessment of hemodynamic changes in breast tissue also relies on using accurate estimates of DPF in the NIRS spectroscopic calculations. Further investigation is required to determine if the proposed EKF algorithm is able to correct path length errors in NIRS applications in stroke, cancer, and other clinical conditions.