## 1.

## Introduction

Functional near-infrared spectroscopy (fNIRS) is a functional neuroimaging technique similar to functional magnetic resonance imaging (fMRI) and positron emission tomography (PET). It detects cerebral functional hemodynamic changes, but not the cerebral neural activation itself. fNIRS potentially has some advantages over other techniques: portability, inexpensive equipment, and availability for simultaneous use with other measurement devices including fMRI, PET, electroencephalography (EEG), and magnetoencephalography (MEG). However, undesirable artifacts are often observed in fNIRS measurements. Systemic physiological changes, physical movements, and even posture changes in a subject can generate such artifacts. These artifacts are observed because the detected light is influenced not only by the hemodynamic change in gray matter, but also by scattering and absorption changes in all tissues through which the light propagates. It is not known exactly how scattering and absorption changes occur in each tissue with a subject’s physiological change. Thus, various statistical methods, such as conventional filtering and multivariate data analysis, have been used to remove such artifacts from fNIRS measurements.

Artifacts induced by physiological activity, such as cardiac pulsation, respiration, and a blood pressure change, can be effectively removed by low-pass filtering or a moving average method. These techniques work well in a real-time environment because they are very simple to implement. However, they tend to overflatten the signal. Thus, for example, a fast-rising signal of brain activation may be deleted if they are applied. To solve this problem, Zhang, Brown, and Strangman^{1, 2} proposed an adaptive filtering technique based on low correlation between artifacts and cerebral functional activity. Saager and Berger^{3} applied the least square fit technique to this problem and conducted a detailed analysis of its performance. Kohno
^{4} used a powerful technique, independent component analysis, to separate a statistically independent signal from data contaminated with various artifacts.

These techniques and a conventional block averaging technique are useful for suppressing task-uncorrelated noises. However, if the data are contaminated with an artifact that is strongly correlated with a task sequence, these methods cannot work effectively. For example, Izzetoglu
^{5} pointed out that a subject’s head movement generates a blood flow change; thus, if the subject moves his head synchronously with a task, an artifact induced by this blood flow change will be observed. Since this kind of artifact is tightly correlated with the task sequence, the block averaging technique cannot remove it.

In general, to remove task-correlated artifacts, we adopt an experimental design including a reference task that is considered to generate the same signal as a target task except for a functional component. The artifacts can be canceled, but the functional component remains when we calculate the difference between the target and reference measurements. Such a design, however, usually requires a longer execution time, and it is very difficult to verify the equivalence of artifact components in target and reference tasks.

Another approach to this problem is multivariate data analysis. This approach is based on the empirical knowledge that the artifact signal is observed globally, and the functional signal is observed locally. Principal component analysis is used to extract the functional signal from multichannel measurement data. This method was reported to be effective.^{6, 7} However, the use of many optodes, most of which have to be positioned at areas unrelated to the target, is a big problem for practical measurement. In addition, if several activations at different areas occur simultaneously, separation of these signals will be difficult.

We do not know precisely how task-correlated artifacts are generated. Presumably they originate mainly from something other than activity in the cortical layer. Thus, it may be useful for functional detection if we can separate the signal of the cortical layer from that of other layers. Hence, several methods using multidistance probe arrangements for separating NIRS signals in lower layers from other layers have been proposed.

A NIRS algorithm using a multidistance probe arrangement, NIR spatially resolved spectroscopy (NIR-SRS), is known.^{8, 9} NIR-SRS was originally developed to measure the tissue oxygen index (TOI), and Al-Rawi, Smielewski, and Kirkpatrick reported that an extracarotid artery clamp had little influence on the TOI measured by NIR-SRS.^{10} This result may suggest that NIR-SRS can remove the interference of extracranial blood flow from brain functional activity signals. However, we cannot theoretically validate its ability to reduce artifacts of a superficial layer, because NIR-SRS is based on the theory of light propagation in a monolayered medium.^{11}

Improved approaches of multidistance probe arrangement based on two-layer models were reported. Choi
^{12} used an analytical solution of a two-layer model that contains an additional superficial layer on a basal semi-infinite medium. They obtained the hemoglobin oxygenation etc. of extra- and intracranial layers by a theoretical curve fitting to experimental data. Fabbri
^{13} conducted Monte Carlo simulations of a two-layered tissue model and showed that a multidistance probe arrangement is effective in detecting the absorption change in the lower layer (cortical layer), even when optical changes arise in both lower and upper layers. This is a viable technique for separating the signal at the cortical layer from artifacts, regardless of whether they are task correlated or uncorrelated.

Validity of these algorithms greatly depends on the optical models assumed. It should be examined by theoretical or simulation approaches whether such a simplified model of a few layers provides a sufficient condition to separate cortical signals from others. In addition, only absorption change in tissues has been considered in most algorithms based on the modified Lambert-Beer law. In a recent study, Tomita, Ohtomo, and Suzuki reported that hemodynamic changes could be accompanied by scattering changes caused by blood cell aggregation.^{14} Since scattering also causes optical attenuation changes, we think that analysis considering both absorption and scattering is important.

In this work, we present a theoretical analysis and Monte Carlo simulations of a five-layered tissue model in which both scattering and absorption changes occur, and show that a multidistance probe arrangement is effective in extracting absorption changes in the gray matter layer. The simulation result is used to determine the appropriate probe arrangement. The proposed method is compared with a conventional method through actual experiments on nonfunctional and functional tasks. Artifacts induced by body tilting, head nodding, and breath holding were clearly observed when a conventional method was used, and were appreciably reduced by the proposed method. Signals evoked by single-sided finger movements were observed globally at both hemispheres when we used a conventional method. On the other hand, localized signals at the primary motor area were observed by the proposed method. A statistically significant increase in oxygenated hemoglobin (HbO) and decrease in deoxygenated hemoglobin (HbR) were simultaneously observed only at each contralateral primary motor area.

## 2.

## Theory

## 2.1.

### Optical Multilayer Model

A typical probe arrangement of fNIRS is illustrated in Fig. 1 (source and detector 1). In this arrangement, light power is detectable only when scattering exists in a medium. Light from a source of power ${I}_{0,\lambda}$ propagates in a medium according to the spatial distribution of the scattering and absorption coefficients ( ${\mu}_{s,\lambda}^{\prime}$ and ${\mu}_{a,\lambda}$ ), and reaches a detector at a distance ${d}_{1}$ with power ${I}_{\lambda}$ . Here, $\lambda $ indicates the wavelength of the incident light. The optical attenuation is defined as ${A}_{d,\lambda}=-\mathrm{ln}({I}_{\lambda}\u2215{I}_{0,\lambda})$ . Since the input power ${I}_{0,\lambda}$ is difficult to measure, the optical attenuation change from a certain standard condition is used in an ordinary fNIRS measurement. When the change in scattering and absorption coefficients ( $\Delta {\mu}_{s,\lambda}^{\prime}$ and $\Delta {\mu}_{a,\lambda}$ ) from the standard condition is small, the optical attenuation change $\Delta {A}_{d,\lambda}$ can be represented as

## Eq. 1

$\Delta {A}_{d,\lambda}={\phantom{|}\frac{\partial {A}_{d,\lambda}}{\partial {\mu}_{a}}|}_{\mathrm{sv}}\Delta {\mu}_{a,\lambda}+{\phantom{|}\frac{\partial {A}_{d,\lambda}}{\partial {\mu}_{s}^{\prime}}|}_{\mathrm{sv}}\Delta {\mu}_{s,\lambda}^{\prime}.$If the second term of Eq. 1 is negligible, this equation corresponds to the modified Lambert-Beer law, which is the theoretical basis of a conventional fNIRS algorithm. In this case, the coefficient
$\text{\hspace{0.17em}}{\phantom{|}\partial {A}_{d,\lambda}\u2215\partial {\mu}_{a}|}_{\mathrm{sv}}$
is called the mean optical path length. In general, however, the second term cannot be neglected. For instance, Tomita, Ohtomo, and Suzuki experimentally demonstrated *in vitro* that the scattering property of a blood volume changes when its flow velocity changes.^{14} Thus, in this work we set the starting point of our theoretical discussion at Eq. 1.

We use a five-layered optical slab model of an adult head, which consists of layers of scalp (sc), skull (sk), cerebrospinal fluid (csf), gray matter (gm), and white matter (wm). Supposing the temporal change in scattering and absorption to be uniform in each layer, Eq. 1 is transformed to

## Eq. 2

$\Delta {A}_{d,\lambda}=\sum _{L}{l}_{d,\lambda}^{L}\Delta {\mu}_{a,\lambda}^{L}+\sum _{L}{m}_{d,\lambda}^{L}\Delta {\mu}_{s,\lambda}^{\prime L},$## Eq. 3

$\{\begin{array}{c}{l}_{d,\lambda}^{L}={\phantom{|}\frac{\partial {A}_{d,\lambda}}{\partial {\mu}_{a}^{L}}|}_{\mathrm{sv}}\\ {m}_{d,\lambda}^{L}={\phantom{|}\frac{\partial {A}_{d,\lambda}}{\partial {\mu}_{s}^{\prime L}}|}_{\mathrm{sv}}\end{array}\phantom{\}}.$## Eq. 4

$\{\begin{array}{l}{l}_{d,\lambda}^{L}=\{{A}_{d,\lambda}({\overline{\mu}}_{a}^{L}+{\mu}_{a}^{L},{\overline{\mathbf{\mu}}}_{o})-{A}_{d,\lambda}({\overline{\mu}}_{a}^{L},{\overline{\mathbf{\mu}}}_{o})\}\u2215{\mu}_{a}^{L}\\ {m}_{d,\lambda}^{L}=\{{A}_{d,\lambda}({\overline{\mu}}_{s}^{\prime L}+{\mu}_{s}^{\prime L},{\overline{\mathbf{\mu}}}_{o}^{\prime})-{A}_{d,\lambda}({\overline{\mu}}_{s}^{\prime L},{\overline{\mathbf{\mu}}}_{o}^{\prime})\}\u2215{\mu}_{s}^{\prime L}\end{array}\phantom{\}},$^{15}) and are shown in Table 1 . The program tMCimg, provided to the public by the Photon Migration Imaging Laboratory at Massachusetts General Hospital (see http://www.nmr.mgh.harvard.edu/PMI/index.htm), was used in the simulation.

## Table 1

Hemodynamic parameters, thickness, and optical properties for each layer of an adult head model.

Tissue type | Thickness[mm] | μa [mm−1] at 670∕690∕720∕750∕780∕830 [nm] |
---|---|---|

Scalp | 3 | $0.031\u22150.026\u22150.022\u22150.022\u22150.020\u22150.019$ |

Skull | 7 | $0.027\u22150.017\u22150.013\u22150.016\u22150.016\u22150.017$ |

CSF | 2 | $0.0044\u22150.0029\u22150.003\u22150.0045\u22150.0044\u22150.0056$ |

Gray matter | 4 | $0.048\u22150.039\u22150.036\u22150.038\u22150.036\u22150.041$ |

White matter | 34 | $0.024\u22150.018\u22150.014\u22150.016\u22150.016\u22150.018$ |

Tissue type | ${\mu}_{s}^{\prime}$ $\left[{\mathrm{mm}}^{-1}\right]$ at $670\u2215690\u2215720\u2215750\u2215780\u2215830$ [nm] | |

Scalp | $2.53\u22152.38\u22152.24\u22152.11\u22152.00\u22151.84$ | |

Skull | $2.33\u22152.13\u22151.92\u22151.79\u22151.66\u22151.47$ | |

CSF | $0.35\u22150.32\u22150.29\u22150.27\u22150.25\u22150.22$ | |

Gray matter | $2.69\u22152.57\u22152.46\u22152.39\u22152.31\u22152.10$ | |

White matter | $10.49\u221510.03\u22159.78\u22159.62\u22159.25\u22158.82$ |

Figure 2 shows a 2-D plot of optical attenuation changes when parameters ${\mu}_{a}^{L}$ and ${\mu}_{s}^{\prime L}$ increase by 10% from standard values in each layer $L$ . Five source-detector distances (10, 15, 20, 25, and $30\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$ ) were assumed, and an incident light of wavelength $780\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$ was used in the simulation. This figure shows the following.

• Scattering and absorption changes at the scalp and skull layers strongly affect the optical attenuation change. Interestingly, the attenuation change increases (skull layer) and decreases (scalp layer) proportionally if we increase the source-detector distance.

• The optical attenuation change induced by an absorption change in the csf layer is small, and the attenuation change due to a scattering change is comparable to that due to an absorption change in the gray matter layer. The observed attenuation change in the csf layer is thought to be negligible because that layer contains no blood cell, and thus scattering and absorption changes due to blood flow changes should be very small.

• The absorption change in the gray matter layer influences the attenuation change considerably. This change increases proportionally with the source-detector distance. Remarkably, the optical attenuation change caused by a scattering change in the gray matter layer is very small.

• The optical attenuation change induced by absorption and scattering changes in the white matter layer is very small.

## Eq. 5

$\Delta {A}_{d,\lambda}=({l}_{d,\lambda}^{\mathrm{sc}}\Delta {\mu}_{a,\lambda}^{\mathrm{sc}}+{m}_{d,\lambda}^{\mathrm{sc}}\Delta {\mu}_{s,\lambda}^{\prime \mathrm{sc}})+({l}_{d,\lambda}^{\mathrm{sk}}\Delta {\mu}_{a,\lambda}^{\mathrm{sk}}+{m}_{d,\lambda}^{\mathrm{sk}}\Delta {\mu}_{s,\lambda}^{\prime \mathrm{sk}})+{l}_{d,\lambda}^{\mathrm{gm}}\Delta {\mu}_{a,\lambda}^{\mathrm{gm}}.$## Eq. 6

$\Delta {A}_{d,\lambda}={l}_{d,\lambda}^{\mathrm{sp}}\Delta {\mu}_{a,\lambda}^{\mathrm{sp}}+{m}_{d,\lambda}^{\mathrm{sp}}\Delta {\mu}_{{s}^{\prime},\lambda}^{\mathrm{sp}}+{l}_{d,\lambda}^{\mathrm{gm}}\Delta {\mu}_{a,\lambda}^{\mathrm{gm}},$## Eq. 7

$\{\begin{array}{l}{l}_{d,\lambda}^{\mathrm{sp}}={l}_{d,\lambda}^{\mathrm{sc}}\\ {m}_{d,\lambda}^{\mathrm{sp}}={m}_{d,\lambda}^{\mathrm{sc}}\end{array}\phantom{\}},\phantom{\rule{1em}{0ex}}\text{if}\phantom{\rule{0.3em}{0ex}}\{\begin{array}{c}\Delta {\mu}_{a,\lambda}^{\mathrm{sk}}=0\\ \Delta {\mu}_{s,\lambda}^{\prime \mathrm{sk}}=0\end{array}\phantom{\}},$## Eq. 8

$\{\begin{array}{l}{l}_{d,\lambda}^{\mathrm{sp}}={l}_{d,\lambda}^{\mathrm{sc}}+{l}_{d,\lambda}^{\mathrm{sk}}\\ {m}_{d,\lambda}^{\mathrm{sp}}={m}_{d,\lambda}^{\mathrm{sc}}+{m}_{d,\lambda}^{\mathrm{sk}}\end{array}\phantom{\}},\phantom{\rule{1em}{0ex}}\text{if}\phantom{\rule{0.3em}{0ex}}\{\begin{array}{c}\Delta {\mu}_{a,\lambda}^{\mathrm{sk}}=\Delta {\mu}_{a,\lambda}^{\mathrm{sc}}\\ \Delta {\mu}_{s,\lambda}^{\prime \mathrm{sk}}=\Delta {\mu}_{s,\lambda}^{\prime \mathrm{sc}}\end{array}\phantom{\}}.$## 2.2.

### Artifact Cancellation by Multidistance Probe Arrangement

If we use two detectors at different source-detector distances ${d}_{1}$ and ${d}_{2}$ $({d}_{1}>{d}_{2})$ , the changes in optical attenuation are written as follows.

## Eq. 9

$\{\begin{array}{c}\Delta {A}_{{d}_{1},\lambda}={l}_{{d}_{1},\lambda}^{\mathrm{sp}}\Delta {\mu}_{a,\lambda}^{\mathrm{sp}}+{m}_{{d}_{1},\lambda}^{\mathrm{sp}}\Delta {\mu}_{s,\lambda}^{\prime \mathrm{sp}}+{l}_{{d}_{1},\lambda}^{\mathrm{gm}}\Delta {\mu}_{a,\lambda}^{\mathrm{gm}}\\ \Delta {A}_{{d}_{2},\lambda}={l}_{{d}_{2},\lambda}^{\mathrm{sp}}\Delta {\mu}_{a,\lambda}^{\mathrm{sp}}+{m}_{{d}_{2},\lambda}^{\mathrm{sp}}\Delta {\mu}_{s,\lambda}^{\prime \mathrm{sp}}+{l}_{{d}_{2},\lambda}^{\mathrm{gm}}\Delta {\mu}_{a,\lambda}^{\mathrm{gm}}\end{array}\phantom{\}}.$## Eq. 10

$\{\begin{array}{l}{k}_{a,\lambda}={l}_{{d}_{1},\lambda}^{\mathrm{sp}}\u2215{l}_{{d}_{2},\lambda}^{\mathrm{sp}}\\ {k}_{s,\lambda}={m}_{{d}_{1},\lambda}^{\mathrm{sp}}\u2215{m}_{{d}_{2},\lambda}^{\mathrm{sp}}\end{array}\phantom{\}}.$## Eq. 11

$\Delta {\mu}_{a,\lambda}^{\mathrm{gm}}=\frac{\Delta {A}_{{d}_{1},\lambda}-{k}_{\lambda}\Delta {A}_{{d}_{2},\lambda}}{{l}_{{d}_{1},\lambda}^{\mathrm{gm}}-{k}_{\lambda}{l}_{{d}_{2},\lambda}^{\mathrm{gm}}},$Figure 3 shows simulated values of ${k}_{a,\lambda}$ and ${k}_{s,\lambda}$ with various combinations of source-detector distances. In many cases of both the minimal and maximal conditions of the superficial layer, ${k}_{a,\lambda}$ and ${k}_{s,\lambda}$ had very similar values. However, when ${d}_{2}$ was $10\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$ , the values did not coincide sufficiently.

We simulated $\Delta {A}_{{d}_{1},\lambda}-{k}_{a,\lambda}\Delta {A}_{{d}_{2},\lambda}$ for various combinations of detector distances. The absorption coefficient at each activation volume of the gray matter layer along the line connecting a source and a detector is assumed to increase by 10% from its standard value. In fMRI studies on human adult brains, functional activations are often observed at a volume of several tens of ${\mathrm{mm}}^{3}$ . On the other hand, if we set a smaller activation volume (larger photon numbers); a longer computational time will be required to obtain sufficient quality of the simulation result. From this practical reason, we fixed a $5\times 5\times 4\phantom{\rule{0.3em}{0ex}}{\mathrm{mm}}^{3}$ volume as an activation volume in this simulation. The result is shown in Fig. 4, which shows the sensitivity distribution for detecting a local absorption change in the gray matter layer. The sensitivities obtained under the minimal and maximal skull layer conditions are very similar to each other. These tend to become high when ${d}_{1}$ is large and ${d}_{2}$ is small. Thus, ${d}_{1}$ should be large and ${d}_{2}$ should be small as much as possible to realize good sensitivity. However, if the source-detector distance is increased, the measurement noise would become large, since fewer photons can reach the detector. In addition, if vasculature and hemodynamic activity in the tissues are not sufficiently homogeneous, a larger distance between two detectors may violate our fundamental assumption that each tissue layer is a uniform slab.

Consequently, an effective probe arrangement should satisfy the following constraints.

• The source-detector distance combination should realize the condition ${k}_{a,\lambda}={k}_{s,\lambda}$ as much as possible.

• ${d}_{1}$ should be large and ${d}_{2}$ should be small to realize good sensitivity to local absorption changes in the gray matter layer.

• To keep the measurement noise small, ${d}_{1}$ should not be too large.

• The difference between ${d}_{1}$ and ${d}_{2}$ should be small enough to satisfy the assumption of uniform tissue layers.

## 2.3.

### Estimation of Partial Path Length Ratio $k$

To calculate the absorption change in the gray matter layer according to Eq. 11, we need to know
${k}_{\lambda}$
. We may be able to use a fixed value for
${k}_{\lambda}$
,^{16} however, it is better to estimate it before each experiment. We may be able to obtain it as
${k}_{\lambda}=\Delta {A}_{{d}_{1},\lambda}^{\text{rest}}\u2215\Delta {A}_{{d}_{2},\lambda}^{\text{rest}}$
,^{13} where
$\Delta {A}_{{d}_{1},\lambda}^{\text{rest}}$
and
$\Delta {A}_{{d}_{2},\lambda}^{\text{rest}}$
are optical attenuation changes at the resting state of a subject. Here, it is assumed that the absorption change in gray matter is small at the resting state. However, because the absolute value of
$\Delta {A}_{{d}_{2},\lambda}^{\text{rest}}$
could be very small, possibly zero, the ratio
$\Delta {A}_{{d}_{1},\lambda}^{\text{rest}}\u2215\Delta {A}_{{d}_{2},\lambda}^{\text{rest}}$
can be very unstable, even when the measurement noise is small. We proposed another estimation method that is stable and gives more accurate values in noisy conditions. The details of the method are described in Appendix { label needed for app[@id='xA'] } in Sec. 6.

We use a triple wavelength measurement, and ${k}_{\lambda}$ is estimated for each wavelength. Monte Carlo simulations under both the minimal and maximal conditions of the superficial layer indicate that the wavelength dependence of ${k}_{a,\lambda}$ and ${k}_{s,\lambda}$ is weak when ${d}_{1}=30\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$ and ${d}_{2}=20\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$ (Fig. 5 ). To statistically stabilize the estimation, we used the average $k$ of three estimated values when using Eq. 11.

## 2.4.

### Calculation of Changes in HbO and HbR

Based on measurement of the absorption changes $\Delta {\mu}_{a,\lambda}$ at several wavelengths, temporal hemoglobin changes $\Delta \mathrm{Hb}\mathrm{O}$ and $\Delta \mathrm{HbR}$ can be calculated as follows:

## Eq. 12

$\left(\begin{array}{c}\Delta \mathrm{Hb}\mathrm{O}\\ \Delta \mathrm{HbR}\end{array}\right)={\left(\begin{array}{cc}{\u03f5}_{\mathrm{Hb}\mathrm{O},{\lambda}_{1}}& {\u03f5}_{\mathrm{HbR},{\lambda}_{1}}\\ {\u03f5}_{\mathrm{Hb}\mathrm{O},{\lambda}_{2}}& {\u03f5}_{\mathrm{HbR},{\lambda}_{2}}\\ \vdots & \vdots \end{array}\right)}^{+}\left(\begin{array}{c}\Delta {\mu}_{a,{\lambda}_{1}}\\ \Delta {\mu}_{a,{\lambda}_{2}}\\ \vdots \end{array}\right),$^{17}

The proposed fNIRS method using a multidistance probe arrangement estimates the absorption change in the gray matter layer as follows [Eq. 11]:

## Eq. 13

$\Delta {\mu}_{a,\lambda}^{\mathrm{gm}}=\frac{\Delta {A}_{{d}_{1},\lambda}-k\Delta {A}_{{d}_{2},\lambda}}{{l}_{{d}_{1},\lambda}^{\mathrm{gm}}-k{l}_{{d}_{2},\lambda}^{\mathrm{gm}}}.$On the other hand, since a conventional fNIRS method assumes that tissues are a monolayer slab and no scattering change occurs, the absorption change measured by such a method is given as

## Eq. 14

$\Delta {\mu}_{a,\lambda}^{\text{mean}}=\frac{\Delta {A}_{{d}_{1},\lambda}}{{l}_{{d}_{1},\lambda}^{\text{mean}}},$Absolute values of mean optical path length
${l}_{{d}_{1},\lambda}^{\text{mean}}$
can be obtained by time-resolved spectroscopy. However, it is difficult to experimentally obtain partial path length
${l}_{{d}_{1},\lambda}^{\mathrm{gm}}$
and
${l}_{{d}_{2},\lambda}^{\mathrm{gm}}$
. In this work, we set
${l}_{{d}_{1},\lambda}^{\mathrm{gm}}-{k}_{\lambda}{l}_{{d}_{2},\lambda}^{\mathrm{gm}}=1$
and
${l}_{{d}_{1},\lambda}^{\text{mean}}=1$
for all wavelengths. As a result of this simplification, we cannot directly compare the magnitudes of
$\Delta \mathrm{Hb}\mathrm{O}$
and
$\Delta \mathrm{HbR}$
calculated by two different methods. Moreover, the simplification about wavelength dependence may lead to a distortion of hemodynamic changes (cross talk
^{18, 19, 20}). This distortion, however, is considered to be small, since a Monte Carlo simulation shows that the wavelength dependence of optical path length is not very strong (data not shown).

## 3.

## Methods

A healthy adult male volunteer participated in the experiment. This study was approved by the Institutional Review Board of National Institute of Advanced Industrial Science and Technology (Japan), and the participant gave his written informed consent.

For fNIRS measurement, we used a near-infrared oximeter with a multiprobe adapter system (NIRO-200 with C9866, Hamamatsu Photonics KK), which has two sources, each with wavelengths 776, 809, and $850\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$ , and eight detectors. The sampling frequency was $0.5\phantom{\rule{0.3em}{0ex}}\mathrm{Hz}$ . Each recorded signal was high-pass filtered off-line at a cut-off frequency of $0.01\phantom{\rule{0.3em}{0ex}}\mathrm{Hz}$ . Based on the simulation results, we set distances from the source to the detectors at ${d}_{1}=30\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$ and ${d}_{2}=20\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$ . Using an originally designed holder, the two sources and eight detectors were fixed on the scalp as shown in Figs. 6e and 6f . The detectors were numbered as shown in Fig. 6c. We implemented the conventional method using detectors 1, 3, 5, and 7 and the proposed method using paired detectors 1–2, 3–4, 5–6, and 7–8.

Paired detectors 1–2 and 7–8 were positioned on the activation area of the finger opposition task as follows. To locate the activation area during the task, T1-weighted images and echo planar images were measured using MR equipment (MR Signa 3.0T, GE Yokogawa Medical Systems KK). The subject was instructed to tap his thumb with his index finger at a frequency of $4\phantom{\rule{0.3em}{0ex}}\mathrm{Hz}$ . The block conditions were alternated using visual cues in this order: left-hand finger stimulation periods $\left(20\phantom{\rule{0.3em}{0ex}}\mathrm{s}\right)$ , rest periods $\left(20\phantom{\rule{0.3em}{0ex}}\mathrm{s}\right)$ , right-hand finger stimulation periods $\left(20\phantom{\rule{0.3em}{0ex}}\mathrm{s}\right)$ , and rest periods $\left(20\phantom{\rule{0.3em}{0ex}}\mathrm{s}\right)$ . A complete session consisted of an initial reference rest period $\left(20\phantom{\rule{0.3em}{0ex}}\mathrm{s}\right)$ and five trials of the block condition (rest, left finger task, rest, right finger task) without any interruption. T-contrast images of left- and right-hand finger opposition against rest periods were obtained using SPM5 (see http://www.fil.iion.ucl.ac.uk/spm/). The images of left-hand activity versus rest appear as colored pixels in Figs. 6a and 6b. The detector was located directly above the activation area, and the position of the activation area was identified by the overlying T1-weighted functional images. The Cz position of an EEG 10–20 system, which is defined as the midpoint of the medial sagittal surface curve connecting nasion and inion, was also identified using the T1-weighted image. The paired detectors 1–2 and 7–8 were located relative to the Cz position.

Optical attenuation changes at the resting state were measured for $10\phantom{\rule{0.3em}{0ex}}\mathrm{min}$ to estimate the ratio $k$ of each pair of detectors. Table 2 shows the correlation coefficients of the optical attenuation changes in each pair of detectors during the resting state. These values show that the attenuation changes at the paired detectors are highly correlated with each other. Based on the measured optical attenuation change, the $k$ value for each pair was calculated and is shown in the last column of Table 2.

## Table 2

Correlation coefficients of ΔAd,λ between adjacent detectors and k values.

Detector pair | Correlation coefficient | k | ||
---|---|---|---|---|

776nm | 809nm | 850nm | ||

1–2 | 0.905 | 0.912 | 0.919 | 1.37 |

3–4 | 0.885 | 0.909 | 0.928 | 1.36 |

5–6 | 0.939 | 0.954 | 0.963 | 1.23 |

7–8 | 0.979 | 0.983 | 0.982 | 0.89 |

To verify the effect of the proposed method on the cancellation of artifacts caused by posture change, body motion, and respiration, the subject, sitting in a chair, was instructed to perform three tasks: tilting his upper body forward by about $30\phantom{\rule{0.3em}{0ex}}\mathrm{deg}$ , nodding his head by about $30\phantom{\rule{0.3em}{0ex}}\mathrm{deg}$ with a frequency of $0.25\phantom{\rule{0.3em}{0ex}}\mathrm{Hz}$ , and holding his breath. In each task, the block conditions were alternated through task periods $\left(20\phantom{\rule{0.3em}{0ex}}\mathrm{s}\right)$ and rest periods $\left(20\phantom{\rule{0.3em}{0ex}}\mathrm{s}\right)$ using auditory cues. A complete session consisted of an initial reference rest period $\left(20\phantom{\rule{0.3em}{0ex}}\mathrm{s}\right)$ followed by a block of five task/rest sequences.

The finger opposition task was also performed to verify that the proposed method can detect a brain functional signal. To evaluate the performance of the proposed and the conventional (using a single detector and the modified Lambert-Beer law) methods for detecting brain functional activity in the finger opposition task, we analyzed statistical significance by the following procedure. The total task sequence comprises five trials of rest-task pairs. We averaged $\Delta \mathrm{Hb}\mathrm{O}$ and $\Delta \mathrm{HbR}$ at each task and rest periods (each has a duration of $20\phantom{\rule{0.3em}{0ex}}\mathrm{s}$ ). We adopted a temporal offset of $6\phantom{\rule{0.3em}{0ex}}\mathrm{s}$ for each period to take into account transient phases of hemodynamic change. Then we calculated the difference between the rest and task averages of five trials. We finally applied the paired $t$ -test to the resultant five differences.

## 4.

## Results

The conventional and proposed methods are compared for cases of body tilting, head nodding, and breath holding in Fig. 7 . $\Delta \mathrm{Hb}\mathrm{O}$ and $\Delta \mathrm{HbR}$ of five task/rest sequences are block-averaged. In all cases using the conventional method, $\Delta \mathrm{Hb}\mathrm{O}$ and $\Delta \mathrm{HbR}$ significantly changed during the task. The magnitude of these changes was comparable to or larger than that of the finger opposition task (shown in Fig. 8 ). The proposed method effectively reduced the magnitude of the changes.

The observed $\Delta \mathrm{Hb}\mathrm{O}$ and $\Delta \mathrm{HbR}$ for the finger opposition task using the two methods are shown in Fig. 8. In the conventional method, increases in $\Delta \mathrm{Hb}\mathrm{O}$ were observed at almost all detectors regardless of the side on which finger opposition was performed. The changes in $\Delta \mathrm{HbR}$ were smaller than those in $\Delta \mathrm{Hb}\mathrm{O}$ . On the other hand, significant simultaneous and opposing changes of $\Delta \mathrm{Hb}\mathrm{O}$ and $\Delta \mathrm{HbR}$ were observed only at the detector positions of 1-2 and 7-8 if we used the proposed method. These changes were observed at the contralateral sides in the finger opposition task. The changes in $\Delta \mathrm{Hb}\mathrm{O}$ were almost equal in magnitude to those in $\Delta \mathrm{HbR}$ .

The differences between the rest and task period averages of five trials are shown in the first (conventional method) and third (proposed method) rows of Fig. 9 . In this figure, averaged differences and their standard errors are given. The results of paired $t$ -tests are shown in the second (conventional method) and fourth (proposed method) rows of Fig. 9. Their numerical values are also given in Table 3 . The figure and table show the following. In the conventional method, the detector positions where $\Delta \mathrm{Hb}\mathrm{O}$ changed over a level of statistical significance of $p=0.05$ were not lateralized during finger opposition. $\Delta \mathrm{HbR}$ showed statistical significance for right finger opposition; however, it showed no significance for left finger opposition. These results suggest that the conventional method failed to detect a cortical functional activation. On the other hand, in the proposed method, the statistical significance of $\Delta \mathrm{Hb}\mathrm{O}$ and $\Delta \mathrm{HbR}$ were clearly lateralized. Statistical significances of the $\Delta \mathrm{Hb}\mathrm{O}$ increase and the $\Delta \mathrm{HbR}$ decrease were simultaneously observed for only two cases: 1 at the 1-2 detector position when left finger opposition was performed, and 2 at the 7-8 detector position when right finger opposition was performed. These positions correspond to locations directly above the activation area identified by fMRI. These results suggest that the proposed method works better than the conventional method in detecting cortical activation.

## Table 3

T -values of the paired t -test for the finger opposition task. Asterisks indicate statistical significance level of p=0.05 .

Detector | Conventional method | Detector | Proposed method | ||||||
---|---|---|---|---|---|---|---|---|---|

Left hand | Right hand | Left hand | Right hand | ||||||

ΔHbO | ΔHbR | ΔHbO | ΔHbR | ΔHbO | ΔHbR | ΔHbO | ΔHbR | ||

1 | ${6.528}^{*}$ | $-2.109$ | 1.527 | 0.873 | 1–2 | ${5.280}^{*}$ | $-{3.844}^{*}$ | $-1.174$ | $-0.277$ |

3 | 2.128 | 0.218 | ${3.691}^{*}$ | $-1.438$ | 3–4 | $-0.144$ | 0.130 | $-0.107$ | $-1.034$ |

5 | 2.563 | 0.833 | ${5.517}^{*}$ | $-2.027$ | 5–6 | $-0.337$ | $-0.026$ | ${3.714}^{*}$ | $-0.838$ |

7 | ${2.717}^{*}$ | 0.128 | ${6.287}^{*}$ | $-{5.630}^{*}$ | 7–8 | 0.279 | 0.130 | ${9.790}^{*}$ | $-{5.815}^{*}$ |

## 5.

## Discussion

Artifacts caused by physical and systemic physiological activities have been a serious problem in fNIRS measurement. We confirmed the existence of artifacts in conventional fNIRS measurement using three tasks: body tilting, head nodding, and breath holding. Their magnitude was comparable to or larger than that observed as a functional signal during finger opposition.

The instability of the optical contact between probes and scalp may have caused some of these artifacts. However, these kinds of signal changes are also observed even when we use another type of probe (A9782, Hamamatsu Photonics KK), which can be attached very stably on the forehead of a subject (data not shown here). Thus, we think that physical movements and systemic physiological activity changes strongly influence conventional fNIRS measurement and may generate serious artifacts, even if we can realize stable optical contacts. Experimental results showed that such artifacts can be eliminated effectively by the proposed method.

Some researchers have reported that NIRS signals were observed at both the contralateral and ipsilateral hemisphere when a subject performs a one-sided finger opposition.^{6, 21} We confirmed the existence of such a global change in the conventional fNIRS signal. Franceschini suggested that some systemic response accompanying the finger movement causes such global changes in
$\Delta \mathrm{Hb}\mathrm{O}$
and
$\Delta \mathrm{HbR}$
.^{6} This component cannot be easily removed by conventional data analysis, and it certainly interferes with detection of brain function. However, if these hemodynamic changes mainly occur in the scalp or skull layers, the interference can be reduced by the proposed method, and in fact, they were reduced well in our experiment (see Fig. 8).

There are some arguments about the correlation between the fNIRS signal and brain functional activity. Some studies claim that a strong correlation exists between
$\Delta \mathrm{HbR}$
decrease and the blood oxygen level dependent signal of fMRI.^{22} In the proposed method, the localization of
$\Delta \mathrm{HbR}$
decrease coincides well with the localization of fMRI signals. Therefore, this observation supports a consistent relationship between fNIRS and fMRI measurements.

Many researchers consider that an increase in
$\Delta \mathrm{Hb}\mathrm{O}$
corresponds to brain activity. In our result using the proposed method, the increase in
$\Delta \mathrm{Hb}\mathrm{O}$
was observed simultaneously with the decrease in
$\Delta \mathrm{HbR}$
; consequently, it is localized in the area where the fMRI signal was observed. On the other hand, in the conventional method, a certain global component of the
$\Delta \mathrm{Hb}\mathrm{O}$
increase at a superficial layer seems to be superposed on the
$\Delta \mathrm{Hb}\mathrm{O}$
increase in the gray matter. Apparently this component is not directly related to brain activity, because it is well known that the cerebral region associated with motor function is highly localized.^{23} If we consider such a contaminated signal of
$\Delta \mathrm{Hb}\mathrm{O}$
increase as a “functional” signal, an artificial correlation of
$\Delta \mathrm{Hb}\mathrm{O}$
to the task sequence may lead fNIRS into false-positive detection of brain activity.

In our multidistance probe arrangement, the distances between the source and each detector were set to 20 and
$30\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$
. However, some studies use another arrangement in which one detector was positioned close to the source to detect the signal from the scalp exclusively, while the other detector is positioned further away to detect the signal from gray matter. Leung, Elwell, and Delpy calculated the measurement accuracy and optimized distance conditions using various numbers of detectors.^{24} In the case using two detectors, they showed that a distance combination of 5 and
$50\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$
is optimal to detect attenuation change in brain. Saager and Berger also used this type of arrangement and showed by Monte Carlo simulations that it is possible to remove artifacts uncorrelated with a functional signal.^{3} Our method is different from their studies in the following three points. First, we intend to remove task-correlated or function-correlated artifacts. Thus, we cannot use the low correlation property between artifacts and functional signals. Second, we intend to remove artifacts caused by both the scattering and absorption changes in superficial layers. If a source-detector distance less than
$10\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$
is selected,
${k}_{a}={k}_{s}$
does not hold sufficiently (see Fig. 3), and it can counteract this type of artifact cancellation. Third, if we use their arrangements, the distance between two detectors becomes greater than
$20\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$
. This can violate the assumption of uniform slab in a two-layer model, because vasculature and hemodynamics at the actual superficial layer may not be sufficiently homogeneous.

We adopted the simplification that only functional hemodynamics occur in the gray matter layer. Under this simplification, detected signals from gray matter directly indicate a functional signal. However, systemic physiological activities can affect hemodynamics in the gray matter layer as well as in the superficial layer. Therefore, we discuss in Appendix { label needed for app[@id='xB'] } in Sec. 7 how systemic hemodynamics in gray matter, if it exists, affects the detection of brain function. In brief, the influence of systemic hemodynamics in gray matter can also be removed by the proposed method; see Appendix { label needed for app[@id='xB'] } in Sec. 7 for details.

In this work, data obtained from one single subject were shown. We implemented similar experiments on several subjects. Also in these cases, we confirmed effective artifact reduction and functional signal localization; however, it may require further validation with more subjects. Through further trials of the proposed method,
$\Delta \mathrm{Hb}\mathrm{O}$
and
$\Delta \mathrm{HbR}$
changes by cerebral activation of human adults, infants, and pathological cases may be revealed in the future. We expect that these can be adopted as cerebral hemodynamic response functions (HRF) in the statistical approaches based on the general linear model.^{25, 26} Such a statistical parameter mapping (SPM)-like analysis based on a precise HRF will increase the reliability of fNIRS measurement.

## Appendices

## { label needed for app[@id='xA'] }

### Appendix A: Procedure for Estimating the Path Length Ratio

If we have temporal absorption and scattering changes $\Delta {\mu}_{a,\lambda}\left(t\right)$ and $\Delta {\mu}_{s,\lambda}^{\prime}\left(t\right)$ , respectively, only at a superficial layer, the optical attenuation changes observed at detectors ${d}_{1}$ and ${d}_{2}$ are represented as

## Eq. 15

$\Delta {A}_{{d}_{1},\lambda}\left(t\right)={l}_{{d}_{1},\lambda}\Delta {\mu}_{a,\lambda}\left(t\right)+{m}_{{d}_{1},\lambda}\Delta {\mu}_{s,\lambda}^{\prime}\left(t\right)+{n}_{{d}_{1},\lambda}\left(t\right),$## Eq. 16

$\Delta {A}_{{d}_{2},\lambda}\left(t\right)={l}_{{d}_{2},\lambda}\Delta {\mu}_{a,\lambda}\left(t\right)+{m}_{{d}_{2},\lambda}\Delta {\mu}_{s,\lambda}^{\prime}\left(t\right)+{n}_{{d}_{2},\lambda}\left(t\right),$## Eq. 17

$\Delta {A}_{{d}_{1},\lambda}\left(t\right)={k}_{\lambda}\Delta {A}_{\lambda}+{n}_{{d}_{1},\lambda}\left(t\right),$## Eq. 18

$\Delta {A}_{{d}_{2},\lambda}\left(t\right)=\Delta {A}_{\lambda}+{n}_{{d}_{2},\lambda}\left(t\right),$To find ${k}_{\lambda}$ , we introduce the following criterion.

## Eq. 19

$J\left(k\right)=\frac{1}{\text{\hspace{0.17em}}}\sum _{t=1}^{N}{\Vert \Delta {A}_{{d}_{1},\lambda}\left(t\right)-k\Delta {A}_{{d}_{2},\lambda}\left(t\right)\Vert}^{2},$## Eq. 20

$=\frac{1}{N}\sum _{t=1}^{N}{\Vert ({k}_{\lambda}-k)\Delta {A}_{\lambda}+{n}_{{d}_{1},\lambda}\left(t\right)-k{n}_{{d}_{2},\lambda}\left(t\right)\Vert}^{2},$## Eq. 21

$={({k}_{\lambda}-k)}^{2}{\sigma}_{A,\lambda}^{2}+{\sigma}_{n,1}^{2}+{k}^{2}{\sigma}_{n,2}^{2},$We assume that the noise variance of detector ${d}_{1}$ ’s output is equal to that of detector ${d}_{2}$ ’s output, and we represent it as ${\sigma}^{2}$ . In general, noise variances ${\sigma}_{n,{d}_{1}}^{2}$ and ${\sigma}_{n,{d}_{2}}^{2}$ in the optical attenuation data differ because the log operation is used to calculate the optical attenuation. If we represent the average outputs of detectors ${d}_{1}$ and ${d}_{2}$ as ${I}_{{d}_{1}}$ and ${I}_{{d}_{2}}$ , respectively, we obtain the following approximate equations:

Then $J\left(k\right)$ becomes## Eq. 24

$J\left(k\right)={({k}_{\lambda}-k)}^{2}{\sigma}_{A,\lambda}^{2}+(\frac{1}{{I}_{{d}_{1}}^{2}}+\frac{{k}^{2}}{{I}_{{d}_{2}}^{2}}){\sigma}^{2}.$If $k={k}_{\lambda}$ , the first term of Eq. 24 becomes zero. Thus, the criterion $J\left(k\right)$ achieves the minimum value when $k={k}_{\lambda}$ if no noise exists $(\sigma =0)$ . However, since the second term of Eq. 24 works as a bias in the criterion, if noise exists, $k$ that minimizes $J\left(k\right)$ tends to be smaller than the correct value. To avoid this difficulty, we propose the following criterion ${J}^{\prime}\left(k\right)$ .

## Eq. 25

${J}^{\prime}\left(k\right)=\frac{J\left(k\right)}{1\u2215{I}_{{d}_{1}}^{2}+{k}^{2}\u2215{I}_{{d}_{2}}^{2}},$## Eq. 26

$=\frac{{({k}_{\lambda}-k)}^{2}}{1\u2215{I}_{{d}_{1}}^{2}+{k}^{2}\u2215{I}_{{d}_{2}}^{2}}{\sigma}_{A,\lambda}^{2}+{\sigma}^{2}.$## { label needed for app[@id='xB'] }

### Appendix B: Reduction of Systemic Hemodynamics in the Gray Matter Layer

If we assume that systemic physiological activities affect hemodynamics in the gray matter layer as well as in superficial layers, Eq. 9 changes as follows.

## Eq. 27

$\{\begin{array}{c}\Delta {A}_{{d}_{1},\lambda}={l}_{{d}_{1},\lambda}^{\mathrm{sp}}\Delta {\mu}_{a,\mathrm{sys},\lambda}^{\mathrm{sp}}+{m}_{{d}_{1},\lambda}^{\mathrm{sp}}\Delta {\mu}_{s,\mathrm{sys},\lambda}^{\prime \mathrm{sp}}+{l}_{{d}_{1},\lambda}^{\mathrm{gm}}(\Delta {\mu}_{a,\mathrm{sys},\lambda}^{\mathrm{gm}}+\Delta {\mu}_{a,\mathrm{func},\lambda}^{\mathrm{gm}})\\ \Delta {A}_{{d}_{2},\lambda}={l}_{{d}_{2},\lambda}^{\mathrm{sp}}\Delta {\mu}_{a,\mathrm{sys},\lambda}^{\mathrm{sp}}+{m}_{{d}_{2},\lambda}^{\mathrm{sp}}\Delta {\mu}_{s,\mathrm{sys},\lambda}^{\prime \mathrm{sp}}+{l}_{{d}_{2},\lambda}^{\mathrm{gm}}(\Delta {\mu}_{a,\mathrm{sys},\lambda}^{\mathrm{gm}}+\Delta {\mu}_{a,\mathrm{func},\lambda}^{\mathrm{gm}})\end{array}\phantom{\}}.$## Eq. 28

$\Delta {\mu}_{a,\mathrm{sys},\lambda}^{\mathrm{gm}}=\gamma \Delta {\mu}_{a,\mathrm{sys},\lambda}^{\mathrm{sp}}.$## Eq. 29

$\{\begin{array}{c}\Delta {A}_{{d}_{1},\lambda}=({l}_{{d}_{1},\lambda}^{\mathrm{sp}}+\gamma {l}_{{d}_{1},\lambda}^{\mathrm{gm}})\Delta {\mu}_{a,\mathrm{sys},\lambda}^{\mathrm{sp}}+{m}_{{d}_{1},\lambda}^{\mathrm{sp}}\Delta {\mu}_{s,\mathrm{sys},\lambda}^{\prime \mathrm{sp}}+{l}_{{d}_{1},\lambda}^{\mathrm{gm}}\Delta {\mu}_{a,\mathrm{func},\lambda}^{\mathrm{gm}}\\ \Delta {A}_{{d}_{2},\lambda}=({l}_{{d}_{2},\lambda}^{\mathrm{sp}}+\gamma {l}_{{d}_{2},\lambda}^{\mathrm{gm}})\Delta {\mu}_{a,\mathrm{sys},\lambda}^{\mathrm{sp}}+{m}_{{d}_{2},\lambda}^{\mathrm{sp}}\Delta {\mu}_{s,\mathrm{sys},\lambda}^{\prime \mathrm{sp}}+{l}_{{d}_{2},\lambda}^{\mathrm{gm}}\Delta {\mu}_{a,\mathrm{func},\lambda}^{\mathrm{gm}}\end{array}\phantom{\}}.$Here, we introduce ${k}_{a,\lambda}^{\prime}$ as

## Eq. 30

${k}_{a,\lambda}^{\prime}=\frac{{l}_{{d}_{1},\lambda}^{\mathrm{sp}}+\gamma {l}_{{d}_{1},\lambda}^{\mathrm{gm}}}{{l}_{{d}_{2},\lambda}^{\mathrm{sp}}+\gamma {l}_{{d}_{2},\lambda}^{\mathrm{gm}}}.$## Eq. 31

$\Delta {\mu}_{a,\mathrm{func},\lambda}^{\mathrm{gm}}=\frac{\Delta {A}_{{d}_{1},\lambda}-{k}_{\lambda}\Delta {A}_{{d}_{2},\lambda}}{{l}_{{d}_{1},\lambda}^{\mathrm{gm}}-{k}_{\lambda}{l}_{{d}_{2},\lambda}^{\mathrm{gm}}}.$## Acknowledgment

The authors thank the Photon Migration Imaging Laboratory at Massachusetts General Hospital for the Monte Carlo simulation code (available at http://www.nmr.mgh.harvard.edu/PMI/index.htm).

## References

**,” J. Biomed. Opt., 12 044014 (2007). https://doi.org/10.1117/1.2754714 1083-3668 Google Scholar**

*Adaptive filtering for global interference cancellation and real-time recovery of evoked brain activity: a Monte Carlo simulation study***,” J. Biomed. Opt., 12 064009 (2007). https://doi.org/10.1117/1.2804706 1083-3668 Google Scholar**

*Adaptive filtering to reduce global interference in evoked brain activity detection: a human subject case study***,” J. Opt. Soc. Am. A, 22 1874 –1882 (2005). https://doi.org/10.1364/JOSAA.22.001874 0740-3232 Google Scholar**

*Direct characterization and removal of interfering absorption trends in two-layer turbid media***,” J. Biomed. Opt., 12 062111 (2007). https://doi.org/10.1117/1.2814249 1083-3668 Google Scholar**

*Removal of the skin blood flow artifact in functional near-infrared spectroscopic imaging data through independent component analysis***,” IEEE Trans. Biomed. Eng., 52 934 –938 (2005). https://doi.org/10.1109/TBME.2005.845243 0018-9294 Google Scholar**

*Motion artifact cancellation in NIR spectroscopy using Wiener filtering***,” Psychophysiology, 40 548 –560 (2003). https://doi.org/10.1111/1469-8986.00057 0048-5772 Google Scholar**

*Hemodynamic evoked response of the sensorimotor cortex measured noninvasively with near-infrared optical imaging***,” J. Biomed. Opt., 11 054007 (2006). https://doi.org/10.1117/1.2363365 1083-3668 Google Scholar**

*Diffuse optical imaging of the whole head***,” Proc. SPIE, 2359 486 –495 (1995). https://doi.org/10.1117/12.209997 0277-786X Google Scholar**

*Absolute quantification methods in tissue near infrared spectroscopy***,” Proc. SPIE, 3597 582 –592 (1999). https://doi.org/10.1117/12.356862 0277-786X Google Scholar**

*A tissue oxygenation monitor using NIR spatially resolved spectroscopy***,” Stroke, 32 2492 –2499 (2001). https://doi.org/10.1161/hs1101.098356 0039-2499 Google Scholar**

*Evaluation of a near-infrared spectrometer (NIRO 300) for the detection of intracranial oxygenation changes in the adult head***,” Appl. Opt., 28 2331 –2336 (1989). https://doi.org/10.1364/AO.28.002331 0003-6935 Google Scholar**

*Time resolved reflectance and transmittance for the non-invasive measurement of tissue optical properties***,” J. Biomed. Opt., 9 221 –229 (2004). https://doi.org/10.1117/1.1628242 1083-3668 Google Scholar**

*Noninvasive determination of the optical properties of adult brain: near-infrared spectroscopy approach***,” Phys. Med. Biol., 49 1183 –1201 (2004). https://doi.org/10.1088/0031-9155/49/7/007 0031-9155 Google Scholar**

*Optical measurements of absorption changes in two-layered diffusive media***,” Neuroimage, 33 1 –10 (2006). https://doi.org/10.1016/j.neuroimage.2006.05.042 1053-8119 Google Scholar**

*Contribution of the flow effect caused by shear-dependent RBC aggregation to NIR spectroscopic signals***,” J. Biomed. Opt., 10 (1), 011015 (2005). https://doi.org/10.1117/1.1846076 1083-3668 Google Scholar**

*Wavelength dependence of crosstalk in dual-wavelength measurement of oxy- and deoxy-hemoglobin***,” J. Biomed. Opt., 14 (6), 064025 (2009). 1083-3668 Google Scholar**

*A Monte Carlo study of global interference cancellation by multidistance measurement of near-infrared spectroscopy***,” Anal. Biochem., 227 54 –68 (1995). https://doi.org/10.1006/abio.1995.1252 0003-2697 Google Scholar**

*Performance comparison of several published tissue near-infrared spectroscopy algorithms***,” Neuroimage, 22 583 –589 (2004). https://doi.org/10.1016/j.neuroimage.2004.02.023 1053-8119 Google Scholar**

*Separability and cross talk: optimizing dual wavelength combinations for near- infrared spectroscopy of the adult head***,” Neuroimage, 13 76 –90 (2001). https://doi.org/10.1006/nimg.2000.0674 1053-8119 Google Scholar**

*The accuracy of near infrared spectroscopy and imaging during focal changes in cerebral hemodynamics***,” J. Biomed. Opt., 14 (3), 034017 (2009). https://doi.org/10.1117/1.3147402 1083-3668 Google Scholar**

*A new crosstalk measure of near-infrared spectroscopy and its application to wavelength combination optimization***,” Med. Phys., 28 1108 –1114 (2001). https://doi.org/10.1118/1.1373401 0094-2405 Google Scholar**

*Wavelength dependence of the precision of noninvasive optical measurement of oxy-, deoxy-, and total-hemoglobin concentration***,” Neuroimage, 29 368 –382 (2006). https://doi.org/10.1016/j.neuroimage.2005.08.065 1053-8119 Google Scholar**

*A temporal comparison of BOLD, ASL, and NIRS hemodynamic responses to motor stimuli in adult humans***,” Brain, 120 141 –157 (1997). https://doi.org/10.1093/brain/120.1.141 0006-8950 Google Scholar**

*Localization of the motor hand area to a knob on the precentral gyrus—a new landmark***,” Phys. Med. Biol., 50 5783 –5798 (2005). https://doi.org/10.1088/0031-9155/50/24/002 0031-9155 Google Scholar**

*Estimation of cerebral oxy- and deoxy-haemoglobin concentration changes in a layered adult head model using near-infrared spectroscopy and multivariate statistical analysis***,” Neuroimage, 44 428 –447 (2009). https://doi.org/10.1016/j.neuroimage.2008.08.036 1053-8119 Google Scholar**

*NIRS-SPM: Statistical parametric mapping for near-infrared spectroscopy***,” J. Biomed. Opt., 12 064010 (2007). https://doi.org/10.1117/1.2804092 1083-3668 Google Scholar**

*Functional optical signal analysis: a software tool for near-infrared spectroscopy data processing incorporating statistical parametric mapping*