1.

Introduction

Over the past 15 years, near-infrared spectroscopy (NIRS) and diffuse optical imaging (DOI) have been used to detect hemodynamic or neuronal changes associated with functional brain activity in a variety of experimental paradigms. ^{1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16} Compared with existing functional methods (e.g., fMRI, PET, EEG, and MEG), the advantages of NIRS and DOI for studying brain function include good temporal resolution, measurement of both oxygenated hemoglobin $\left({\mathrm{O}}_2\mathrm{Hb}\right)$ and deoxygenated hemoglobin $\left(\mathrm{HHb}\right)$ , nonionizing radiation, portability, and low cost.^{4, 6} Disadvantages include modest spatial resolution, limited penetration depth, potential sensitivity to hair absorption and motion artifacts, and global interference (also called systemic physiological interference).

Global interference can arise from at least two spatial sources: 1. in the superficial layers (such as scalp and skull), and 2. inside brain, due to factors such as heart activity, respiration, and spontaneous low frequency oscillations [i.e., low frequency oscillations (LFOs) and very low frequency oscillations (VLFOs)]. ^{17, 18, 19, 20, 21} In empirical studies of brain function using NIRS and DOI, the amount of global interference varies from subject to subject and from time to time. In some cases, the amount of interference is small and evoked brain activity can be seen in the raw measurement; other times the amount of interference is too large for the evoked brain activity to be detected without signal processing.¹⁸ Several methods have been explored for the removal of global interference and improvement of evoked brain activity measurements. Low pass filtering is the most common and straightforward, as it is highly effective in removing cardiac oscillations.^{22, 23} However, for physiological variations such as respiration, LFOs, and VLFOs, there is a significant overlap between their frequency spectra and that of the hemodynamic response to brain activity. Frequency-based removal of these interferences can therefore result in large distortion and inaccurate timing for the recovered brain activity signal. Other methods for improving the contrast-to-noise ratio (CNR, mean signal during task performance minus mean signal during baseline, divided by the noise) for NIRS-based brain function measurements include adaptive average waveform subtraction,²⁴ direct subtraction of a “nonactivated” NIRS waveform,²² state space estimation,^{25, 26, 27} and principal components analysis.²⁸

Recently, Morren adopted the technique of adaptive filtering to remove cardiac oscillations, using signals acquired from pulse oximetry as a reference.²⁹ Adaptive filtering has been widely used in interference cancellations,^{30, 31} and has great potential in the removal of global interferences and recovery of evoked brain activity detection, as demonstrated in EEG and MRI studies.³² The advantage of adaptive filtering includes its capability of following the signal’s nonstationary changes and its simple implementation with low computational overhead. Since biomedical signals are generally nonstationary and real-time features are desired for most NIRS applications, adaptive filtering has the potential to be a good fit for NIRS applications. Morren showed that this method effectively removes the cardiac-related signal variation in the optical measurements. However, in the application of evoked brain activity detection, global interference includes not only cardiac oscillations but also other physiological variations such as vasomotor waves and respiration, which are not represented by pulse oximetry signals. Moreover, pulse oximetry is often measured from fingers or toes, far from the head, and acquires measurements at different wavelengths. Thus, this reference signal is less representative of the NIRS signal observed during head measurements and hence is not optimal for reducing head/brain-based interferences.

In evoked brain activity detection, a good reference measurement used in adaptive interference cancellation should be highly correlated to the global interference. Ideally, it measures directly the interference while avoiding any sampling of the evoked response. Previously, researchers have used optical measurements with short interoptode distances for monitoring superficial layer hemodynamics.^{17, 20} For example, McCormick have attempted to measure interference from superficial layers using short source-detector separation to correct cerebral oxygen delivery monitoring.²⁰ No detailed algorithm for interference correction was presented in their publication, and the superficial layer interference measurement was used simply to visually compare measurements from far source-detector separations rather than for interference correction. The methodology we developed combines a multiseparation probe for data collection and adaptive filtering for signal processing, and it can be used in conjunction with the existing methods such as low pass filtering. This method has low computational requirements, and hence can be implemented in real time, an important feature needed for potential applications for real-time brain function localization procedures, biofeedback, and potential neuroprosthetic devices.

2.

Multidistance Optode Configuration and Adaptive Filtering Algorithm to Remove Global Interference

According to a photon transport theory,³³ photons propagating through a highly scattering tissue travel along a zig-zag path before they are detected. The collective photon propagation follows a roughly banana-shaped pattern (formulated by three-point Green’s function^{34, 35}) when reflection geometry is used, as in most applications of NIRS in the measurement of neuronal activity (see Fig. 1 ). With appropriate source and detector placement, we can make one channel primarily sensitive to the shallow layer hemodynamic changes (S-D1, with close separation of source and detector; Fig. 1) and another channel sensitive to hemodynamic changes, both in the shallow layer (unavoidably) and on the cerebral cortex (S-D2, with far separation of source and detector). In adaptive interference cancellation, measurements from S-D2 can be used as a target signal channel, and with measurements from S-D1 used as reference channel. This process is equivalent to assuming a linear mapping between the shallow layer hemodynamics, acquired from S-D1, and the global interference in the target measurement from S-D2. Optimized cancellation is then achieved by a point-to-point optimization of this linear mapping. This cancellation, and hence the improvement in CNR, is maintained even when hemodynamic changes in the superficial layer are nonstationary, so long as the changes are relatively slow compared to the adaptive filter convergence rate.

Fig. 1

Multisource-detector separation approach. Superficial layer hemodynamics acquired from S-D1 will be used to estimate the global interference presented in the target measurement from S-D2, which is then canceled using adaptive filtering.

Our study used the continuous wave (CW) NIRS method, where two wavelengths (690 and $830\mathrm{nm}$ ) of light are shone into the head, detected at the scalp’s surface, and then are converted to relative changes in the concentration of deoxyhemoglobin $\left(\mathrm{HHb}\right)$ and oxyhemoglobin $\left({\mathrm{O}}_2\mathrm{Hb}\right)$ using the modified Beer-Lambert law. ^{5, 34, 36, 37} First we calculate the changes in the absorption coefficient by:

An adaptive filter with a finite impulse response (FIR) and transversal structure (tapped delay line) is used in our global interference cancellation.³⁸ The filter output signal $e_i$ is given by:

Here the $\Delta \left[\mathrm{HHb}\right]$ (or $\Delta \left[{\mathrm{O}}_2\mathrm{Hb}\right]$ ) acquired from far source-detector separations (S-D2 in Fig. 1), which contains evoked brain hemodynamic changes, is used as the target measurement (the signal channel), denoted by $y_i$ . The subscript $i$ is the index of the time point. The $\Delta \left[\mathrm{HHb}\right]$ (or $\Delta \left[{\mathrm{O}}_2\mathrm{Hb}\right]$ ) acquired from short source-detector separations (S-D1 in Fig. 1) is used as reference measurement (the reference channel), denoted by $x_i$ . $M$ represents the order of the filter and the ${\omega }_k$ are the filter coefficients, where $k$ is the coefficient index. Since the coefficients are adjusted by the filter output $e_i$ on a sample by sample basis, we use ${\omega }_{k,i}$ to denote the $k$ ’th coefficient at time $i$ . Coefficients were updated via the Widrow-Hoff least mean square (LMS) algorithm.³⁹ This algorithm is simple and fast, features needed for real-time applications, especially for applications such as diffuse optical imaging where hundreds of channels would have to be processed together in real time. The LMS algorithm for optimization is:

The processing steps proceed as follows. First, we calculate the $\Delta \left[\mathrm{HHb}\right]$ (or $\Delta \left[{\mathrm{O}}_2\mathrm{Hb}\right]$ ) time series for both close and far source-detector separations, $x_i$ and $y_i$ , respectively. These time series become the inputs to the adaptive filter. The adaptive filter converts $x_i,$ the hemodynamic and oxygenation variation associated with the superficial layers, to an estimate of the global interferences embedded in $y_i$ . Finally, this estimate is then subtracted from the original time series. The transfer function of the adaptive filter is optimized dynamically, via LMS [Eq. 4], to ensure the best quality of cancellation and account for variations of the living tissue. To expedite the convergence of the LMS algorithm, we can normalize the two time series $x_i$ and $y_i$ , so that both have standard deviations close to one. In real application on human subjects, such normalization could be achieved by collecting (e.g.) relatively short (30 to $60\sec$ ) pretest recordings prior to running the main experiment. Pretests will also allow pretraining of the adaptive filter to acquire good initial FIR filter coefficients.

The performance of the filter is controlled by the order of the FIR filter $M$ and the step size used in updating the FIR nodes $\mu$ . Note that if $M=1$ , the adaptive filter becomes a straight subtraction:

3.

Monte Carlo Simulation Study of Evoked Brain Activity Recovery

3.2.

Simulating Global Interference and Evoked Functional Brain Hemodynamics

Generally, multiple Monte Carlo simulations have to be performed to acquire the simulated time series. However, using the tMCimg Monte Carlo program (implemented in C),^{40, 41} pathlengths through the simulated tissue are recorded for each photon detected. So for different time points where only absorption changes, the result can be calculated without rerunning the whole simulation. We launched 100,000,000 photons from the source for the simulated data collection at each time point. As shown in Fig. 2 , the size of the simulated tissue is $150\times 100\times 50{\mathrm{mm}}^3$ . The thickness and scattering properties selected for each layer can be found in the first column of Table 1 . Because the source and detectors are in the middle of the simulated tissue surface, the boundary effect was ignored. The index of refraction mismatch was also ignored in this study (both tissue and air were given a refraction index of 1).

Fig. 2

Geometry of the Monte Carlo simulation for removal of global interferences and recovery of evoked brain activity.

Table 1

Hemodynamic parameters used in the Monte Carlo simulation.

Head layers	Blood content	Baseline concentration (μM)	A (μM) , respiration	B (μM) , heartbeat	C (μM) , evoked response
Scalp, $7\mathrm{mm}$ $({\mu }_s^{\prime }=10{\mathrm{cm}}^{-1})$	${\mathrm{O}}_2\mathrm{Hb}$	39	2	1.1	0
	$\mathrm{HHb}$	16	0.13	0.074	0
Skull, $7\mathrm{mm}$ , $({\mu }_s^{\prime }=12{\mathrm{cm}}^{-1})$	${\mathrm{O}}_2\mathrm{Hb}$	39	2	1.2	0
	$\mathrm{HHb}$	16	0.13	0.078	0
CSF, $1\mathrm{mm}$ $({\mathrm{\mu }}_s^{\prime }=0.1{\mathrm{cm}}^{-1})$	${\mathrm{O}}_2\mathrm{Hb}$	11.7	0.2	0.12	0
	$\mathrm{HHb}$	4.8	0.01	0.006	0
Gray matter, $3\mathrm{mm}$ $({\mu }_s^{\prime }=5{\mathrm{cm}}^{-1})$	${\mathrm{O}}_2\mathrm{Hb}$	56	2	1.2	15
	$\mathrm{HHb}$	20	0.13	0.076	$-4$
White matter, $33\mathrm{mm}$ $({\mu }_s^{\prime }=7{\mathrm{cm}}^{-1})$	${\mathrm{O}}_2\mathrm{Hb}$	56	2	1.1	0
	$\mathrm{HHb}$	20	0.13	0.072	0

The hemodynamic changes in the scalp, skull, CSF, and gray and white layers were simulated as a combination of cardiac fluctuation $c\left(t\right)$ , respiratory fluctuations $r\left(t\right)$ , and functional hemodynamic responses $v\left(t\right)$ . In a real experiment, there would be a certain amount of uncorrelated changes in the superficial layers compared with deep layers. For example, the skin may sweat, which will only produce variations in the scalp layer. To simulate this phenomenon, we also introduced a slow varying random time series $g\left(t\right)$ in the scalp layer response. In summary:

Eq. 6

$$f_{\mathrm{HHb}}^1\left(t\right)={\mathrm{HHb}}_0^1+g\left(t\right)[A_{\mathrm{HHb}}^{1}c\left(t\right)+B_{\mathrm{HHb}}^{1}r\left(t\right)+C_{\mathrm{HHb}}^{1}v\left(t\right)],$$

Eq. 7

$$f_{{\mathrm{O}}_2\mathrm{Hb}}^1\left(t\right)={\mathrm{O}}_2{\mathrm{Hb}}_0^1+g\left(t\right)[A_{{\mathrm{O}}_2\mathrm{Hb}}^{1}c\left(t\right)+B_{{\mathrm{O}}_2\mathrm{Hb}}^{1}r\left(t\right)+C_{{\mathrm{O}}_2\mathrm{Hb}}^{1}v\left(t\right)],$$

Eq. 8

$$f_{\mathrm{HHb}}^{2,3,4,5}\left(t\right)={\mathrm{HHb}}_0^{2,3,4,5}+A_{\mathrm{HHb}}^{2,3,4,5}c\left(t\right)+B_{\mathrm{HHb}}^{2,3,4,5}r\left(t\right)+C_{\mathrm{HHb}}^{2,3,4,5}v\left(t\right),$$

Eq. 9

$$f_{{\mathrm{O}}_2\mathrm{Hb}}^{2,3,4,5}\left(t\right)={\mathrm{O}}_2{\mathrm{Hb}}_0^{2,3,4,5}+A_{{\mathrm{O}}_2\mathrm{Hb}}^{2,3,4,5}c\left(t\right)+B_{{\mathrm{O}}_2\mathrm{Hb}}^{2,3,4,5}r\left(t\right)+C_{{\mathrm{O}}_2\mathrm{Hb}}^{2,3,4,5}v\left(t\right),$$

where $f_{\mathrm{HHb}}^1\left(t\right)$ , $f_{{\mathrm{O}}_2\mathrm{Hb}}^1\left(t\right)$ , $f_{\mathrm{HHb}}^{2,3,4,5}\left(t\right)$ , and $f_{{\mathrm{O}}_2\mathrm{Hb}}^{2,3,4,5}\left(t\right)$ represent the concentration of deoxygenated and oxygenated hemoglobin in each layer as a function of time, with the superscripts 1 to 5 indicating the layer index for scalp, skull, CSF, and gray and white matters, respectively. ${\mathrm{HHb}}_0$ and ${\mathrm{O}}_2{\mathrm{Hb}}_0$ represent the average or baseline concentrations. The coefficients $A$ , $B$ , and $C$ with layer index as superscript and $\mathrm{HHb}$ or ${\mathrm{O}}_2\mathrm{Hb}$ as subscript are the hemodynamic variation amplitude control parameters. They are used to adjust the magnitude of variations of deoxygenated and oxygenated hemoglobin concentrations in each layer due to cardiac pulsation, respiration, and evoked brain response. The values of the previously described parameters, including the baseline concentration and the variation amplitude control parameters $A$ , $B$ , $C$ for deoxygenated and oxygenated hemoglobin in each layer, can be found in Table 1. Since this study focuses on signal and hemodynamic variations, the constant absorption from water and other background chromophores is not considered. The parameters used are based on Choi ⁴² and others, ^{34, 43, 44, 45} together with our own human subject data.

The cardiac and respiratory oscillations $c\left(t\right)$ and $r\left(t\right)$ were both simulated as amplitude- and frequency-varying sinusoidal oscillations:

Eq. 10

$$c\left(t\right)=M_{\mathrm{heart}}\sin \left(2\pi f_{\mathrm{heart}}t\right),$$

Eq. 11

$$r\left(t\right)=M_{\mathrm{resp}}\sin \left(2\pi f_{\mathrm{resp}}t\right),$$

where $M_{\mathrm{heart}}$ , $f_{\mathrm{heart}}$ , $M_{\mathrm{resp}}$ , and $f_{\mathrm{resp}}$ are all random variables, and were generated by low pass filtered Gaussian white noise (with offset added). The average value and standard deviation of the previous variables and the bandwidth of the low pass filter (fourth order, Butterworth) used to filter the Gaussian white noise are listed in Table 2 . These values were chosen based on our human subject data.

Table 2

Simulation parameters for amplitude and frequency of interference oscillations.

Parameters	Baseline (average) value	Standard deviation	Bandwidth of the low pass filter (Hz)
$M_{heart}$ (a.u.)	1	0.32	0.8
$f_{heart}$ (Hz)	1.1	0.003	0.2
$M_{resp}$ (a.u.)	1	0.18	0.2
$f_{resp}$ (Hz)	0.18	0.002	0.1
$g\left(t\right)$ (a.u.)	1	0.1	0.01

The evoked hemodynamic response $v\left(t\right)$ was defined as the convolution of the stimulation paradigm $s\left(t\right)$ , [ $s\left(t\right)=0$ for rest and 1 for stimulation] and a prototypical hemodynamic impulse response $h\left(t\right)$ ⁴⁶:

Eq. 12

$$s\left(t\right)=\{\begin{array}{c}0,t∊\mathrm{rest}\hfill \\ 1,t∊\mathrm{stimulation}\hfill \end{array},$$

Eq. 13

$$h\left(t\right)={\left(\frac{t}{bc}\right)}^b\exp (b-\frac{t}{c});b=8.6,c=0.547,$$

Eq. 14

$$v\left(t\right)=h\left(t\right)\otimes s\left(t\right),0\le t\le 200.$$

Independent scalp variation $g\left(t\right)$ is also generated by biased and low pass filtered Gaussian white noise; its related information can be found in Table 2.

The hemoglobin concentration changes at each layer were converted to reduced absorption coefficients using a linear transform with specific extinction coefficients at 690 and $830\mathrm{nm}$ . Monte Carlo simulations were performed using the prior parameters (about $15\mathrm{h}$ of compute time on a Gateway 450 laptop), and the simulated fluence rate for a 2000 point time series ( $200\sec$ of data at a sampling rate of $10\mathrm{Hz}$ ) was calculated without rerunning the Monte Carlo for each time point (the calculation for one entire time series took about $32\mathrm{h}$ in Matlab). Noise was then added to the acquired simulated optical measurements. Simulated electronic noise was generated by low pass filtered white noise to $3\mathrm{Hz}$ , with a standard deviation of 1/100 standard deviation caused by respiration and heart beat.

3.4.

Simulated Optical Measurements

Figure 3 shows the raw simulated optical measurement at $690\mathrm{nm}$ and its power density spectrum, acquired from the 1.5- and 4.5-cm source-detector separations. The simulated measurements qualitatively match human subject data.⁴⁷ The gray bars indicate blocks of evoked stimulation. As from the raw time series and from the power density spectrum, neither the simulated measurements from S-D1 nor from S-D2 show any obvious evoked response, because the signal variation is dominated by global interference, i.e., respiration and heartbeat.

Fig. 3

Simulated optical measurements at $690\mathrm{nm}$ from S-D1 and S-D2 (a) and (c), together with their power density spectrum (b) and (d).

3.5.

Adaptive-Filtering-Based Removal of Global Interference and Recovery of Simulated Hemodynamic Response

After the simulated optical measurements are acquired, the data analysis described in Sec. 2.1 is applied. For the filtering parameters, considering both convergence speed and filter stability, we chose $\mu =0.0001$ and we use $M=100$ , about twice as large as the average period of the respiration fluctuation. To speed up the convergence, we prenormalized the target and the reference datasets by dividing by their estimated standard deviation, so that both time series have a standard deviation of approximately 1. After the adaptive filtering, the quantitative values of the filtered results were recovered by multiplying back the standard deviation. The initial guess of the adaptive filter weights was $[10\dots 0]$ , in other words we assume that the shape of the global interference and the S-D1 measurements are identical at the beginning of the adaptive filtering. $\left[{\mathrm{O}}_2\mathrm{Hb}\right]$ and $\left[\mathrm{HHb}\right]$ were filtered separately and are presented in Figs. 4 and 5 , respectively. In these figures, we show both raw time series of calculated blood concentration (the first column) and their block averaged results (the second column). For the block average results, we average from $5\sec$ prior through $25\sec$ after the onset of each simulated stimulation. The gray bars indicate the period of active stimulation.

Fig. 4

Adaptive filtering to remove global interference and to recover evoked brain activity. (a) target ${\mathrm{O}}_2\mathrm{Hb}$ measurements from S-D2 with 4.5-cm source-detector separation and (b) its block averaged result. (c) and (d) Reference measurements from S-D1 with 1.5-cm source-detector separation. (e) and (f) Low pass filtering result for the target measurements. (g) and (h) Adaptive filtering result for the target measurement. (i) and (j) Adaptive filtering result with sensitivity correction (solid line), together with the true evoked brain activity (dashed line) used for quantitative comparison and evaluation of the performance of the adaptive filtering method.

Fig. 5

Adaptive filtering to remove global interference and to recover evoked brain activity. (a) Target $\mathrm{HHb}$ measurements from S-D2 with 4.5-cm source-detector separation and (b) its block averaged result. (c) and (d) Reference $\mathrm{HHb}$ measurements from S-D1 with 1.5-cm source-detector separation. (e) and (f) Low pass filtering result for the target measurements. (g) and (h) Adaptive filtering result for the target measurement. (i) and (j) Adaptive filtering result with sensitivity correction (solid line), together with the true evoked brain activity (dashed line) used for quantitative comparison and evaluation of the performance of the adaptive filtering method.

In Fig. 4, where the ${\mathrm{O}}_2\mathrm{Hb}$ result is presented, the target dataset for adaptive filtering is the calculated $\left[{\mathrm{O}}_2\mathrm{Hb}\right]$ from S-D2 with 4.5-cm source-detector separation [Figs. 4a and 4b], and the reference dataset is the calculated $\left[{\mathrm{O}}_2\mathrm{Hb}\right]$ from S-D1 with 1.5-cm source-detector separation [Figs. 4c and 4d]. In theory, the target dataset should contain the evoked hemodynamic response; however, neither its raw time series [Fig. 4a] nor its block average [Fig. 4b] show visible signal change correlated with the stimulation paradigm. In fact, its similarity (both in the raw time series and the block averages) to the reference dataset [Figs. 4c and 4d] indicates that it is dominated by global interferences.

After adaptive filtering, as seen from the filtered result [Figs. 4g and 4h], approximately 80% of the signal variation has been removed. This is another indication that the signal variations in the target dataset are predominantly global interferences that correlate well with the superficial layer responses. After we multiplied the sensitivity correction factor back into the filtered result, we compared it with the real evoked brain hemodynamic changes in the gray matter layer. The comparison can be seen in Figs. 4i and 4j, with the solid line as the filtered result recovered after sensitivity correction, and the dashed line as the true evoked brain response used in the simulation. Although obvious interferences still remain, the contrast-to-noise ratio is dramatically improved for evoked hemodynamic response detection. We also see, quantitatively, good agreement between the filtered result (after sensitivity correction) and the true value. Cardiac and respiration variations were significantly reduced, even without the application of any bandpass filtering. The CNR of the evoked response can be further improved by block averaging [Fig. 4j], since the global variations are not temporally correlated (due to the random factors in the frequency and amplitude of heartbeat and respiration), and can be canceled and thus suppressed significantly by block averaging.

Comparing the recovered evoked hemodynamics [solid line in the block average results, Fig. 4j] with the real underlying hemodynamics [dashed line in the block average result, Fig. 4j], we noticed several types of errors. 1. Residual interference: the recovered evoked response is not as smooth as the real hemodynamics, due to the leftover global interference after filtering. 2. Quantitative error: the amplitude of the recovered evoked response is about 85% of the real value. 3. Timing error or phase lag: the rise and fall of the evoked response has a small amount of time delay compared with the real hemodynamic response. Residual interference is expected, as no technique will completely remove truly random interference. In our adaptive filtering, the error may be because the random factor $g\left(t\right)$ in superficial layers reduced the correlation between the global interference in the superficial layers and the global interference in the brain layers. In addition, the convergence of the LMS algorithm is relatively slow, and the results may improve when another adaptive filtering algorithm with faster convergence is used.

We also compared the adaptive filtering result with a traditional low pass filtering result. $\left[{\mathrm{O}}_2\mathrm{Hb}\right]$ acquired from S-D2 [Fig. 4a] was low pass filtered with an eighth-order Butterworth filter with 0.5-Hz bandwidth, and the result is shown in Figs. 4e and 4f. From Fig. 4e, we can see that this low pass filter effectively removes the cardiac oscillations; however, systemic variations due to respiration remain. After block average, seen in Fig. 4f, the systemic interference was further suppressed; however, it is still too large for the evoked brain activity to be detected. Comparing Fig. 4f with the final adaptive filtering block averaged result [Figs. 4h and 4j], we can see that adaptive filtering removes the global interference and recovers the evoked brain activity much more effectively.

In Fig. 5, we show the equivalent of Fig. 4 but for $\left[\mathrm{HHb}\right]$ . The evoked response is visible, although noisy, after block averaging without adaptive filtering [Fig. 5b]. This is because the amount of global interference in $\left[\mathrm{HHb}\right]$ is relatively small. After adaptive filtering, the global interference is significantly removed, CNR increased, and the evoked response can be clearly seen [Figs. 5g and 5h]. The performance of adaptive filtering for interference removal and evoked brain activity detection can be further demonstrated with sensitivity correction [Figs. 5i and 5j]. As for ${\mathrm{O}}_2\mathrm{Hb}$ , we can see there are still interferences left; however, cardiac and respiration variations were significantly reduced, again without bandpass filtering. The quantitative value of the functional blood concentration change was recovered to about 80% of the expected value. When comparing adaptive filtering result with the low pass filtering result shown in Figs. 5e and 5f using an eighth-order Butterworth filter with 0.5-Hz bandwidth, we can see that this low pass filter again effectively removes the cardiac oscillations; however, systemic variations due to respiration remain. After block averaging [Fig. 5f], the evoked brain activity can be clearly seen, and systemic interference was further suppressed; however, the residuals are still quite obvious. In comparison, the final adaptive filtering block averaged result [Figs. 5h and 5j] has much less residual interferences and improved CNR.

Both $\left[{\mathrm{O}}_2\mathrm{Hb}\right]$ and $\left[\mathrm{HHb}\right]$ demonstrated approximately 20% quantitative error after adaptive filtering [Figs. 4j and 5j]. We hypothesized that there was a small sensitivity to brain tissue at a 1.5-cm separation, and the evoked brain activity detected by the reference channel was removed from the signal channel after adaptive filtering, resulting in reduced filtered results. This hypothesis was supported by the fact that when we re-ran the simulations using a reference separation of $1.0\mathrm{cm}$ in place of $1.5\mathrm{cm}$ , we were able to recover 95% of the expected value of functional blood concentration change (result not shown). Generally, using a reference that is partly correlated with the expected hemodynamic response will result in a loss in quantitative accuracy when filtering. Optimization of the selection of reference channels for the best filtering performance will be studied and reported in the future.

Together, our Monte Carlo simulations show that, in principle, our methodology significantly removes global interference and improves the CNR in evoked brain activity detection.

4.

Resistance to Positional Errors and Differential Pathlength Factor Errors

Errors in the location of optical sources and detectors and in the assumed DPFs are two major error sources in NIRS using CW measurements and MBLL. Both affect the quantitative values of the result significantly, but less on the relative variations or the shape of the $\left[{\mathrm{O}}_2\mathrm{Hb}\right]$ and $\left[\mathrm{HHb}\right]$ time series. Since this adaptive-filtering-based method uses mostly the shape of the time series rather than the quantitative values of the time series, it is resistant to the errors in source and detector position measurements or in DPFs. One example demonstrating this method’s resistance to DPF error is shown in Fig. 6 . The optical measurements are the same as described in Sec. 3.2; however, when calculating the $\left[{\mathrm{O}}_2\mathrm{Hb}\right]$ and $\left[\mathrm{HHb}\right]$ time series, the DPFs used for S-D1 and S-D2 are 4.3 and 5.6 at $690\mathrm{nm}$ , and 3.6 and 4.6 at $830\mathrm{nm}$ (i.e., introduction of 25% error). From the result we see that the quantitative values of $\left[{\mathrm{O}}_2\mathrm{Hb}\right]$ and $\left[\mathrm{HHb}\right]$ changed; however, the global interferences are still effectively removed and evoked responses successfully detected, although the quantitative values of the recovered evoked response changed due to DPF error. Generally speaking, when DPFs are underestimated, the reconstructed quantitative values will be overestimated and vice versa. This effect, however, will be systematic throughout an experiment, such that all measurements will be over- or underestimated in a similar way, preserving comparability in terms of relative change.

Fig. 6

Adaptive filtering results using inaccurate DPFs (25% error).

An example demonstrating our method’s resistance to positional error is shown in Fig. 7 . In this test, we introduced random positional error (generated in Matlab using the function “randn,” with standard deviation of $5\mathrm{mm}$ ). Again, global interferences are effectively removed and evoked responses are still successfully detected, although the quantitative values of the recovered evoked response changed due to source and detector positional error. Generally, when source-detector separation estimates are smaller than the real separations, the calculated blood concentration will be larger than the results using the correct separations, and vice versa.

Fig. 7

Adaptive filtering results using inaccurate source-detector separation ( $\pm 5\mathrm{mm}$ error).

5.

Comparison of Shallow Layer Hemodynamics and Global Interference

The performance of the adaptive-filtering-based removal of the global interference depends on how well the global interference (in the target measurement from S-D2) and the reference measurement (from S-D1) are correlated. In principle, global interference in the target dataset comes from all layers, since respiration, heart beat, and related hemodynamics could exist in all tissue types. Our reference measurement, however, is predominantly sensitive to shallow layer hemodynamics, and hence changes there may differ from those in deeper layers. Two major reasons conspire to make the reference measurements and the global interference in the target measurements correlate well: 1. sensitivity dominance of the superficial layer, and 2. physiological connection between superficial and deeper layers.

5.1.

Sensitivity Preeminence of Superficial Layers

First, superficial layer variation is the predominant component in the total global interference in the noninvasive surface measurement. According to photon transport theory, even for measurements with large source-detector separation like S-D2, the surface measurements are more sensitive to superficial layer changes.^{34, 35} Although the total global interferences consist of superficial layer hemodynamic variations and deep-layer hemodynamic variations, the superficial layer hemodynamic variation comprises a much larger portion of the total global interferences. Hence, a measurement of the superficial layer hemodynamics is a very good choice for a reference measurement in a noninvasive reflectance probe geometry.

Our Monte Carlo simulation allows us to look at the total global interference and its components (i.e., from superficial layers and from cerebral layers) separately. This is done as follows. To acquire total systemic interference, we performed simulation with the same parameters as those we described in Sec. 3.2, except that there is no evoked hemodynamic response added to the gray layer. Then, to acquire the component of global interference from superficial layers, we performed simulation with the respiration and cardiac hemodynamics restricted to scalp and skull layers only. Finally, to acquire the component of global interference from deep layers, we restricted the respiration and cardiac to gray and white matters layers only. Again, those simulation time series are calculated based on recorded photon pathlength information without rerunning individual Monte Carlos at each time point. Simulated optical measurements were acquired from S-D2 with 4.5-cm source-detector separation, and then blood content was calculated using MBLL.

The $\left[{\mathrm{O}}_2\mathrm{Hb}\right]$ result is shown in Fig. 8 , where the solid line represents the total global interference in ${\mathrm{O}}_2\mathrm{Hb}$ , acquired when respiration and cardiac hemodynamics exists in all layers; the dotted line represents the global variations from scalp and skull; and the dashed line represents the variations from deep brain layers. As Fig. 8 shows, the total global interference and its components look very similar. Quantitatively, the interference from scalp and skull is the major component, and it is about 80% of the total global interference, which is approximately the same proportion of the signal removed during adaptive filtering. This decomposition of total global interference clearly demonstrates the sensitivity preeminence of superficial layers for optical measurements with source-detector separation of $4.5\mathrm{cm}$ in our simulation.

Fig. 8

Total global interference and its components. Solid line (labeled 1) is the total global interference in the $\left[{\mathrm{O}}_2\mathrm{Hb}\right]$ from S-D2 in the simulation; dotted line (labeled 2) is the interference from superficial scalp and skull layers; and dashed line (labeled 3) is the interference from deep layers of gray and white matters.

Table 3

Calculated DPFs using CW and RF methods.

	690nm		830nm
Methods	S-D1 (1.5cm)	S-D2 (4.5cm)	S-D1 (1.5cm)	S-D2 (4.5cm)
CW method	5.4	7.0	5.1	6.6
RF method	5.5	7.1	5.2	6.7

6.

Real-Time Implementation of the Methodology

The LMS adaptive filtering algorithm we discussed is particularly suited to real-time implementation. We have implemented the algorithm in $\mathrm{C}++$ and tested it (Intel Pentium M processor, $2000\mathrm{MHz}$ , 400-MHz external bus), and it is capable of processing $833\mathrm{K}$ data points (64 bit data type) per second with adaptive filtering $(M=100)$ . For a typical DOI imaging system with 50 sources, 50 detectors (2500 channels), and 100-Hz sampling rate, the time required to process all data collected in $1\sec$ with adaptive filtering is $0.3\sec$ , or 30% of the total time (70% can be used for data acquisition and other processing). Therefore, a real-time removal of global interference and recovery of evoked brain activity is practical on even modestly current computer hardware, an important feature required for applications such as brain activity localization, biofeedback, and potential neuroprosthetic devices.

7.

Discussion

Using Monte Carlo simulation, we demonstrated that one should be able to combine a multiseparation probe configuration with adaptive filtering to dramatically reduce global interference and improve the CNR in evoked brain activity detection. We further demonstrated that superficial layer interference is the major component of the total global interference, making superficial layer hemodynamics a particularly good estimate for the total global interference. An added advantage of our method is its suitability for real-time use.

In previous studies,^{23, 34, 48} we utilized simple, traditional signal conditioning methods (high and low pass filtering). The implicit assumption was that the overlying layers, even if they do interfere with the brain measurement, do not correlate with the stimulation protocol and hence can be reduced or eliminated by block averaging or statistical methods. Across studies, however, we have consistently observed notable intersubject variability in 1. the robustness of evoked responses and 2. the magnitude of nonevoked variability (especially cardiac and respiratory oscillations). Based on our findings here—particularly the observation that approximately 80% of our raw signal is removed during filtering—we suggest that the robustness of noninvasively recorded evoked brain responses may, at least in some cases, be much more severely affected by interference from overlying tissue layers than previously assumed. Preliminary application of this technique to human data supports this suggestion.

While we tried to make the simulation as realistic as we could, it is not possible to include every potential physiological variable. In the Monte Carlo simulation work, we chose to include respiration and cardiac oscillations, with respiration being representative of lower frequency oscillations (and indeed we selected a relatively low frequency for respiration of $0.18\mathrm{Hz}$ ). In the human subject case report, the major interferences were found to be Mayer waves and cardiac, and it was demonstrated that they were also effectively removed by adaptive filtering. In principle, with appropriate reference signal and adaptive filter design, adaptive filtering helps to remove interferences with different types.

In this study we performed adaptive filtering after converting raw signal to blood content. Roughly speaking, the raw optical signal is an exponential function of the blood concentration or absorption coefficient; the total measurement is a product of contributions from global interference and evoked brain hemodynamics. Linear adaptive filtering does not work well in this case, since here the output signal is total signal minus (not divided by) estimated interference. Put another way, there is a nonlinear operation (logarithm) between raw intensity and $\Delta \left[{\mathrm{O}}_2\mathrm{Hb}\right]$ and $\Delta \left[\mathrm{HHb}\right]$ , and nonlinear operations cannot be transposed to get the same answer. However, after the raw optical measurements are converted to absorption coefficients (i.e., after taking the logarithm), then there should be minimal difference switching the linear procedures of blood conversion (using the absorption coefficient) and adaptive filtering. Since the result will be presented as hemodynamics (blood concentration variations) at the end, the DPF error cannot be avoided. Conceptually it is relatively easier to explain the physiology and discuss the interference from ${\mathrm{O}}_2\mathrm{Hb}$ and $\mathrm{HHb}$ separately. In addition, in many of our human subject studies we found that the interference in ${\mathrm{O}}_2\mathrm{Hb}$ and $\mathrm{HHb}$ are different. While information about ${\mathrm{O}}_2\mathrm{Hb}$ and $\mathrm{HHb}$ are merged in the raw signals at each wavelength, converting the raw signal to blood components first gives us the flexibility of choosing which blood species we will perform the adaptive filtering on to get the greatest benefit.

Three limitations of the method need to be mentioned. First, in this study the stimulation paradigm and the global interference were uncorrelated. In some experiments, however, the evoked brain activity and the global interference may be correlated to some degree. For example, tasks requiring exertion may increase respiration and heart rates in concert with task performance, thereby making the global interference and brain activity coupled. How well our adaptive filtering will work for this scenario is still under investigation. Generally, error will occur if the reference time series and the signal to detect are correlated, and this error depends on the magnitude of correlated components in the two time series. Second, good performance of our method depends on measurements from S-D1 with short source-detector separation correlating well with the global interference in measurements from S-D2. In our simulation the head is assumed to be layered homogenous, while those layers in an actual human head can be very heterogeneous, e.g., due to the heterogeneous distribution of large blood vessels. We are susceptible to independent optical changes around the location of D1, or movement of D1 independent of D2. In principle, adaptive filtering effectively removes only interferences that are contained both in the target measurement and in the reference channel; any independent changes or noise in the reference measurement can be added (with inverted sign) into the target measurement by the subtraction step in the filtering process. Thus, the reference measurement needs to be as free from artifact as possible. One possible solution might be to collect more than one short source-detector separation reference measurement and to combine or select these references as needed. Third, as compared with the true evoked hemodynamic response, the response recovered via adaptive filtering still has leftover interferences, as well as quantitative error and phase distortion. The focus of this study is to preliminarily demonstrate the methodology, and solutions to the previously mentioned limitations are under investigation.

Overall, the combination of a multidistance probe and adaptive filtering appears quite promising for removing a major source of error in NIRS for evoked brain activity detection. Such filtering may be particularly useful in studies of the optical fast signal,^{2, 16, 29} where signals are small and lowpass filtering removes the signal of interest along with the noise. In depth discussion about the optimization of the adaptive filter parameters and comparison of different adaptive filter structures will be reported in the future. In our tests, we have observed that the best parameters ( $M$ and $\mu$ ) in the adaptive filter are relatively robust and do not vary substantially from case to case. By optimizing the probe configuration and algorithm, for example choosing the best source-detector separations for the reference measurement and target measurement, or using adaptive filtering algorithms with faster convergence (e.g., normalized LMS, recursive least squares, adaptive RLS), we may be able to further improve the performance of this methodology.