2.1.

Data Collection

The visual stimulation task was performed on a healthy, 36-year-old male, with no known history of neurological, psychiatric, or vision problems. NIRS measurements and simultaneous EKG and respiration recordings were collected during our study. Optical data was collected using NIRS1 (TechEn, Milford, Massachusetts), as used in previous studies.^{13, 37} We used two source lasers, at 682 and $830\mathrm{nm}$ , respectively, and two detectors (Hamamatsu C5460-01). The lasers were intensity modulated by an approximately 5-kHz square wave and an in-phase/quadrature-phase (IQ) detection scheme. The final analog output was low-pass-filtered at $3.4\mathrm{Hz}$ using a second-order RC low-pass filter. An EKG was measured using chest leads, and respiration rate was collected using a piezo-electric respiration sensor (Pneumotrace II, UFI, Morro Bay, California) tied at the upper abdomen, and the measurements were digitally low-pass filtered at $1\mathrm{Hz}$ (eighth-order IIR filter). Heart rate (HR) was calculated from the EKG recordings, done on a beat-to-beat basis from QRS waves. First the derivative of the EKG signal is calculated, and then a threshold is applied to the EKG derivative so that QRS waves can be detected. The threshold is manually selected—it can be either for the rising or falling edge of the QRS wave to get the best result. Errors from this straightforward QRS detection will be corrected manually, and the calculated beat-to-beat HR will then be digitally low-pass filtered at $1\mathrm{Hz}$ (sixth-order IIR filter). The optical signals, together with the EKG and respiration measurements, were simultaneously sampled at $200\mathrm{Hz}$ and transferred to a computer using a 16-bit A/D conversion (Measurement Computing, PCM-DAS16D/16); thus, all signal channels were strictly co-registered to enable accurate comparison of timings of event onsets.

The lights from two diode lasers were combined using a bifurcated glass fiber bundle with a nominal $2.7\text{-}\mathrm{mm}$ core diameter (Fiberoptics Technology, prototype). Detector fibers were glass fiber bundles, also with $2.7\text{-}\mathrm{mm}$ core diameter. The optical fiber bundles (optodes) were attached to a flexible piece of plastic that was attached to a black Velcro headband that fastened snugly around the head of the subject. Soft black Velcro was also used to absorb stray light beneath the probe. The combined light of two wavelengths were in contact with tissue in one source location [S in Fig. 1a ], and exiting light was collected from two detector locations, $1.5\mathrm{cm}$ and $4.5\mathrm{cm}$ from the source, respectively [D1 and D2 in Fig. 1a].

Fig. 1

(a) Probe configuration for the vision test. (b) Vision stimulation: alternating counter-phased radial checkerboard. (c) Rest period: uniform field with a central fixation cross.

During the visual task, the multichannel NIRS probe was positioned $2\mathrm{cm}$ lateral and parallel to the subject’s midline, over the primary/secondary visual cortex (confirmed by an existing MRI scan), with the source $1\mathrm{cm}$ inferior to the inion [Fig. 1a]. MRI images also show that the thickness of the scalp and skull at this location is approximately $2.1\mathrm{cm}$ , and with this probe configuration, photon propagation theory suggests that measurement from S-D1 will be sensitive to scalp and skull changes, while measurements from S-D2 will be sensitive to scalp, skull, CSF, vision cortex, and some white matter regions. The subject sat upright in a comfortable chair in a dark and quiet room facing a laptop screen. The stimulus consisted of a counterphased, 100% contrast radial checkerboard ( $5\text{-}\mathrm{Hz}$ luminance reversals) covering approximately $20\mathrm{deg}$ of the visual field [Fig. 1b]. This was alternated with a mean-luminance, uniform field with a central fixation cross [Fig. 1c]. One $8\text{-}\min$ run of 16 blocks was presented, where each block consisted of checkerboard oscillation for $15\mathrm{s}$ , followed by the fixation background for $15\mathrm{s}$ ; thus, the base visual stimulation frequency was $0.033\mathrm{Hz}$ .

2.2.

Data Analysis

The raw 200-Hz optical data were offset-corrected and further digitally low-pass filtered, in addition to the instrument filter, using a 1.25-Hz, sixth-order Butterworth filter. The optical measurements for each channel were converted to relative changes in the concentration of HHb and ${\mathrm{O}}_2\mathrm{Hb}$ using the modified Beer-Lambert law. ^{38, 39, 40, 41} Concentration measurements were then bandpass filtered using a sixth-order Butterworth low-pass filter at $1.25\mathrm{Hz}$ in series with a sixth order Butterworth high-pass filter at $0.004\mathrm{Hz}$ (to remove any slowly drifting signal components) and then down-sampled to $20\mathrm{Hz}$ . Based on our previous Monte Carlo simulations with common head structure and tissue optical properties,³⁶ we chose DPFs of 5.4 and 7.0 for S-D1 and S-D2 at $690\mathrm{nm}$ , and 5.1 and 6.6 for S-D1 and S-D2 at $830\mathrm{nm}$ . (We initially set all DPFs to 6.0 for the data analysis, and the results closely mirrored the results reported here, suggesting insensitivity to DPF errors.³⁶ The reported DPFs are based on more realistic Monte Carlo estimates of DPF differences across wavelengths and separations.) The time series of ${\mathrm{O}}_2\mathrm{Hb}$ concentration changes acquired from both source-detector pairs (S-D1 and S-D2) were fed into the adaptive filter. (HHb was filtered similarly but separately.) The adaptive filter used a finite impulse response (FIR) and transversal structure (tapped delay line), with 100 nodes (“weights”). The Widrow-Hoff least-mean-squared (LMS) adaptation algorithm was used to optimize the filter coefficients on a sample-by-sample basis. The details can be found in our previous report.³⁶ The convergence control parameter was $\mu =0.0001$ . In order to speed the convergence, we prenormalized the target and the reference data sets 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. If this initial guess were to remain constant, it would be similar to Saager and Berger’s approach, which assumes that the entire S-D1 time series can be least-squares fit to and subtracted from the S-D2 measurement.³⁵ Our approach differs in that (1) we assume a linear mapping between S-D1 and S-D2 and that this mapping is dependent on recent history (here, 100 weights instead of the entire time series, making the approach particularly suited to real-time applications), and (2) that this mapping is allowed to change over time, which in principle allows us to handle nonstationary time series (see Sec. 4.2).

As in previous NIRS studies, we performed block averaging using nonoverlapping portions of the data set triggered on stimulus onsets (averaging windows $2\mathrm{s}$ prior to checkerboard onset and $28\mathrm{s}$ after the onset of each block of visual stimulation).

After adaptive filtering, in order to retrieve the quantitative functional hemodynamic changes in the deep cortex layer in our previous Monte Carlo simulation study,³⁶ we calculated a sensitivity correction scaling factor. This is done by introducing a known perturbation into the blood content of the simulated gray matter layer only and then calculating the blood content changes using the surface measurements and the MBLL. The sensitivity correction factor is defined as the ratio of true blood concentration change to measured blood volume change. In the simulation with general physiological parameters, this sensitivity correction factor is calculated to be 27.5 for ${\mathrm{O}}_2\mathrm{Hb}$ and 26.4 for HHb, and we assumed the same correction factor in this human study.

Power spectral density (PSD) was used to identify the components and vision stimulation–associated changes in the hemodynamic responses or other physiological signals. PSD was estimated in MATLAB (based on Welch’s averaged, modified periodogram method, with a Hanning window for smoothing). In order to find peaks that are statistically significant, we chose a confidence interval of 0.95 and calculated the upper and lower bound of the power spectrum within this interval. Usually, we say a peak can be identified if the lower bound of the peak is above the averaged surrounding noise floor shown in the PSD.

With PSD, we can also perform a quantitative CNR analysis. For human subject data, we cannot totally separate signal and noise. Thus, we roughly calculate the signal power by integrating the power (in linear scale) from the PSD over the “signal bands,” which are the small frequency segments surrounding the stimulation frequency and its second harmonic. This way, the evoked brain activity peak will normally be included. We integrate the rest of the power spectrum (“noise bands”) to get the noise power. The unavoidable error in this calculation includes real signal power in the noise band and actual noise power in the signal bands. The CNR is then calculated as the square root of the ratio of signal power to noise power. It is worth mentioning that since the PSD is calculated using the whole data set (9674 points at $20\mathrm{Hz}$ sampling rate, no overlap in the PSD segmentation), the CNR calculated based on this PSD is the result of averaging the whole data set, not the CNR for a single stimulation.

3.1.

Origins of Global Interference: Comparing Optical, Respiration, and Cardiac Measurements

The respirometry and EKG were recorded simultaneously with optical measurements in order to help us determine the origins and understand the nature of the global interferences. The raw optical measurements (light intensity) at $830\mathrm{nm}$ collected from the “far” detector (S-D2 with $4.5\text{-}\mathrm{cm}$ source-detector separation) and from the “near” detector (S-D1 with $1.5\text{-}\mathrm{cm}$ source-detector separation) are shown in Figs. 2a and 2c , respectively. In order to more clearly view the details, such as the temporal correlation between different channels, we have shown here only $200\mathrm{s}$ of data. The associated power spectral densities are shown in Figs. 2b and 2d. The simultaneously recorded respiration signal (left $y$ axis with arbitrary unit, indicating the movement of the upper abdomen) and heart rate (right $y$ axis, calculated from the EKG) are presented in Fig. 2e, and the associated PSD is shown in Fig. 2f. The insets in Figs. 2b, 2d, 2f show the details of PSD at frequencies less than $0.2\mathrm{Hz}$ . Each PSD was estimated using all 96,736 data points so that more details can be seen at the low-frequency end. To avoid possible confusion, we need to point out again that the lower curve in Fig. 2e is the heart rate, not the EKG, and its spectrum in Fig. 2f is actually the spectrum for heart rate variation (HRV).

Fig. 2

(a) The raw optical measurement at $830\mathrm{nm}$ from S-D2 with $4.5\text{-}\mathrm{cm}$ source-detector separation (the inset of an enlarged data segment shows the heartbeats), and (b) its power spectral density magnitude (the inset shows the details at frequencies $<0.15\mathrm{Hz}$ ). (c) Raw optical measurement at $830\mathrm{nm}$ from S-D1 with 1.5-cm source-detector separation, and (d) its power spectral density magnitude. (e) The top curve is the normalized voltage output from the respiration sensor, and the bottom curve is the heart rate in beats per minute (bpm), calculated from EKG recordings, and (f) their power spectral density. The 1, 2, and 3 in (a), (c), and (e) indicate the sequence of an example wave in respiration, heart rate, and optical.

In Figs. 2a and 2c, we see that both optical measurements consist of a slow variation and a fast oscillation. The fast oscillation is obviously due to cardiac activity, it co-registers strictly with the EKG (not shown), and the inset in Fig. 2a shows an enlarged data segment to demonstrate the temporal details of these fast signal fluctuations. The PSD of the optical signals, shown in Figs. 2b and 2d, show that the cardiac activity peaks at $1.1\mathrm{Hz}$ , and its second harmonic at about $2.2\mathrm{Hz}$ is also very obvious.

From the respiration PSD [the inset in Fig. 2f], we see that the principal respiration frequency is at $0.16\mathrm{Hz}$ . However, for both optical measurements [seen in the insets in Figs. 2b and 2d], the major low-frequency peak is $\sim 0.1\mathrm{Hz}$ (with a $0.16\text{-}\mathrm{Hz}$ peak appearing about $10\text{-}\mathrm{dB}$ down, and among peaks that might be artifacts). The heart rate PSD [inset in Fig. 2f] helps in explaining the origin of this $0.1\text{-}\mathrm{Hz}$ peak in optical channels. The heart rate PSD shows $0.1\mathrm{Hz}$ and $0.16\mathrm{Hz}$ as the top two principal peaks, with the $0.1\text{-}\mathrm{Hz}$ peak as the highest and $2\mathrm{dB}$ above the $0.16\text{-}\mathrm{Hz}$ peak. (The latter is obviously due to the respiratory sinus arrhythmia.) This suggests that the $0.1\text{-}\mathrm{Hz}$ oscillation in the optical channels is physiological rather than artifact. Since $0.1\mathrm{Hz}$ is the prototypical LFO frequency, we suspect that this oscillation in the hemodynamics shown in the optical channels is vasomotor waves, and the $0.1\text{-}\mathrm{Hz}$ peak in the heart rate PSD is possibly an associated change in a vasomotor reflex. Full separation of component signals is difficult in this subject, because the subject’s respiration rate is relatively low, and we lack other more direct measurements such as invasive arterial pressure.

Optical measurements from far and near channels are very similar, actually a linear correlation analysis that shows that they are highly correlated ( $r=0.91$ , $p<{10}^{-20}$ ). The fact that the target (far) measurement and reference (near) measurement are highly correlated and that they all demonstrate major PSD peaks at heart beat and LFO frequencies indicates that the optical measurement from S-D2 is heavily contaminated by global interference (measured by S-D1), as was previously derived via Monte Carlo experiments.

We performed a T-test on the 16 averaged heart rates during vision stimulation periods $(67.0\pm 2.5\mathrm{bpm})$ and the 16 averaged heart rates during rest periods $(66.4\pm 2.5\mathrm{bpm})$ . The change was not significant [ $T\left(30\right)=-0.7$ , $p=0.24$ ], suggesting no obvious heart rate increase or associated hemodynamic interference due to the visual stimulation itself.

To summarize, the PSD in the optical and physiological measurements suggests that the systemic interference comes from mixed sources, with cardiac activity and LFO (Mayer waves or vasomotor waves, judging from the frequency) being dominate. The visual stimulation does not cause any heart rate increase.

3.2.

${\mathrm{O}}_2\mathrm{Hb}$ Changes during Visual Stimulation

Changes in ${\mathrm{O}}_2\mathrm{Hb}$ during visual stimulation are shown in Fig. 3 . Here, Figs. 3a and 3b represent the target measurement of ${\mathrm{O}}_2\mathrm{Hb}$ calculated from the raw “far” signal and its block average. Figures 3c and 3d show the reference ${\mathrm{O}}_2\mathrm{Hb}$ (“near” signal) and its block average. Although the target data set was expected to contain an increase in ${\mathrm{O}}_2\mathrm{Hb}$ following stimulus onset (and concomitant decrease following stimulus offset), neither the raw time series nor the block average show any obvious expected signal change. We hypothesized that the global interferences obscured the expected brain-related ${\mathrm{O}}_2\mathrm{Hb}$ change. The fact that the signal variations in the target ${\mathrm{O}}_2\mathrm{Hb}$ closely match those of the reference ${\mathrm{O}}_2\mathrm{Hb}$ (which should contain no visual response) suggested that indeed global interference may dominate the target data set.

Fig. 3

Adaptive filtering to remove global interference and 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.

The target data set after adaptive filtering is shown in Figs. 3g, 3h, 3i, 3j. Figures 3g and 3h show the raw adaptive filtered result with its block average, and 3(i) and 3(j) show the filtered result following sensitivity correction. Comparing Figs. 3g and 3h [or their sensitivity-corrected versions in 3i and 3j] with 3a and 3b, we see that after adaptive filtering, both LFO and cardiac interference are substantially reduced (actually, 80% of the signal variation was removed), while at the same time, the CNR of the hemodynamic responses to vision stimulation increased dramatically. Both the time series’ and the block averaged results show an expected increase following stimulation onset and return to baseline during rest periods, with temporal changes appropriately associated with the stimulation paradigm. For the detected visual response shown in Fig. 3j, we see a uniphase increase starting from approximately $2\mathrm{s}$ after stimulus onset, which plateaus at about $9\mathrm{s}$ . The response remains steady until after the visual stimulation ends (we suspect that the hump at about $17\mathrm{s}$ is artifact), and begins to decrease at about $20\mathrm{s}$ , recovering to baseline at about $25\mathrm{s}$ . Although unlike in our previous Monte Carlo simulation study, here we do not have a true ${\mathrm{O}}_2\mathrm{Hb}$ visual response value to compare with, the quantitative ${\mathrm{O}}_2\mathrm{Hb}$ changes match previous reports.^{9, 10}

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

The CNR improvement from adaptive filtering can be quantitatively estimated using power spectral density techniques. Figure 4a shows the PSD of the target ${\mathrm{O}}_2\mathrm{Hb}$ before adaptive filtering, in which we see two dominating peaks at about $1\mathrm{Hz}$ and $0.1\mathrm{Hz}$ , which correspond to cardiac activity and LFO. Figure 4b shows the details of the PSD of target ${\mathrm{O}}_2\mathrm{Hb}$ in 0 to $0.12\text{-}\mathrm{Hz}$ range (solid line), together with its upper and lower bound calculated at a confidence interval of 0.95. In Fig. 4b, using the average local noise level (horizontal dashed line; averaged PSD magnitude in the 0 to $0.12\mathrm{Hz}$ range) as the detection threshold (to compare with the lower bound of the PSD), we can see several peaks above this threshold (e.g., at $0.006\mathrm{Hz}$ , $0.012\mathrm{Hz}$ , $0.043\mathrm{Hz}$ , $0.066\mathrm{Hz}$ , $0.104\mathrm{Hz}$ , and $0.1\mathrm{Hz}$ ), the LFO peak being the most significant one, about $12\mathrm{dB}$ above. A peak of $0.033\mathrm{Hz}$ , the expected visual stimulation frequency, is actually among detected peaks; however, it is the second smallest compared with other detected peaks, only $7\mathrm{dB}$ above the threshold, and difficult to distinguish among the other detectable peaks as real or artifact. After adaptive filtering [Fig. 4c], the physiological interference due to LFO and cardiac activity peaks was significantly reduced. From the details presented in Fig. 4d, we can see that the noise floor (horizontal dashed line) was $10\mathrm{dB}$ less than the original, while the visual response peak remains [compared with Fig. 4b]. Now only two peaks can be detected ( $0.010\mathrm{Hz}$ and $0.033\mathrm{Hz}$ ), and the peak at $0.033\mathrm{Hz}$ , corresponding to vision stimulation, is the most significant peak and is $18\mathrm{dB}$ above the local noise floor.

Fig. 4

(a) and (b) present the power spectral density of ${\mathrm{O}}_2\mathrm{Hb}$ before adaptive filtering. (a) is the PSD of ${\mathrm{O}}_2\mathrm{Hb}$ acquired from S-D2 (4.5-cm source-detector separation) before adaptive filtering. (b) shows the details of the PSD in the 0 to 0.12-Hz range. The solid line is the PSD magnitude, and the dotted line is the upper and lower bound calculated at a confidence interval of 0.95. The horizontal dashed line indicates the average noise level (averaged PSD magnitude in the 0 to 0.12-Hz range). (c) and (d) present the filtered result in the same way as (a) and (b).

If we choose the signal bands of the vision response to be roughly 0.028 to $0.039\mathrm{Hz}$ (base band) and 0.061 to $0.072\mathrm{Hz}$ (second harmonic), and the rest of the spectrum in the 0 to $10\mathrm{Hz}$ range is the noise bands, our CNR analysis shown that adaptive filtering doubled the CNR (from 40.2% to 80.8%). The very low frequency of $\sim 0.01\mathrm{Hz}$ was also preserved. It remains to be determined whether signal variation at this frequency range is real or artifact.

3.3.

HHb Changes during Vision Stimulation

The adaptive filtering results for HHb are shown in Fig. 5 . As for Figs. 3i and 3j, Figs. 5i and 5j again utilized sensitivity correction. The PSD result is shown in Fig. 6 .

Fig. 5

Adaptive filtering to remove global interference and recover evoked brain activity. (a) Target HHb measurements from S-D2 with 4.5-cm source-detector separation, and (b) its block averaged result. (c) and (d) Reference 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.

Fig. 6

(a) and (b) present the power spectral density of HHb before adaptive filtering. (a) is the PSD of HHb acquired from S-D2 (4.5-cm source-detector separation). (b) shows the details of the PSD in the 0 to 0.12-Hz range. The solid line is the PSD magnitude, and the dotted line is the upper and lower bound calculated at a confidence interval of 0.95. The horizontal dashed line indicates the average noise level (averaged PSD magnitude in the 0 to 0.12-Hz range). (c) and (d) present the filtered result in the same way as (a) and (b).

One obvious difference in the HHb data as compared to ${\mathrm{O}}_2\mathrm{Hb}$ is that for HHb, the interference was relatively small and the hemodynamic response was visible before adaptive filtering [Figs. 5a and 5b]. In addition, the signal from the near detector (shallow layer hemodynamics) differs more substantially from the far detector signal (deep-layer hemodynamics)—that is, the far detector HHb signal was not dominated by synchronized global interference. We suspect that interferences here have local origins and thus are different at different detector locations. This explanation is supported by the PSD analysis. In the PSD of target HHb measurement [Figs. 6a and 6b], we actually see no obvious LFO, respiration, or cardiac activity peaks. In other words, the major global interference sources in HHb are not significant. In Sec. 5, we discuss why in this subject HHb did not demonstrate obvious global interference peaks, as previously observed in ${\mathrm{O}}_2\mathrm{Hb}$ .

Since the amount of interference in the target data set is relatively small and the visual response and the global interference are mostly independent of each other, the adaptive filtering did not substantially alter the signal. As seen in Figs. 5g and 5h, the overall amplitude of variations after adaptive filtering have only slightly changed (17% smaller), unlike the ${\mathrm{O}}_2\mathrm{Hb}$ adaptive filtering results, where 80% of the variations in the target data set were removed [Figs. 4g and 4h]. Comparing the results before and after adaptive filtering [Figs. 5a with 5g, respectively], we see that the visual response is clear in both, although adaptive filtering appears to slightly reduce the signal and increase noise. This can be more clearly seen in the PSD analysis in Figs. 5b and 5d, where we see a $5.9\text{-}\mathrm{dB}$ reduced peak at $0.033\mathrm{Hz}$ , and a $0.89\text{-}\mathrm{dB}$ elevated noise level. The previous CNR analysis shows a 35% reduction of CNR after adaptive filtering.

When comparing Figs. 5i and 5j with the low-pass filtering result shown in Figs. 5e and 5f (which use an eighth-order Butterworth filter with $0.125\text{-}\mathrm{Hz}$ bandwidth), we can see that the low-pass filtering result is smoother. The shape of the response in 5(j) (onset, rise, offset, and fall times) has less low frequency fluctuation and is slightly more temporally consistent with prototypical hemodynamic responses than that in Fig. 5f, but it is unclear if this is because of adaptive filtering.

For the filtered and sensitivity-corrected visual response shown in the block averaged and sensitivity corrected results in Fig. 5j, we see a uniphase decrease starting from approximately $1\mathrm{s}$ after the onset of the stimuli, which continues until about $8\mathrm{s}$ , when it reaches a valley. The signal remains at this value until after the visual stimulation ends, and it begins to increase at about $21\mathrm{s}$ , again recovering to baseline at about $25\mathrm{s}$ . These quantitative HHb changes due to vision stimulation also match those found in previous reports.^{9, 10}

The improvement from adaptive filtering on HHb in this case is not obvious. Although the shape of the filtered results [Fig. 5j] seems to agree more with expectations compared with Fig. 5b without adaptive filtering, it is unclear if this slight improvement is really from adaptive filtering, and the CNR estimated using PSD is actually reduced (from 66.5% to 43%). This suggests that in the cases where hemodynamic signals demonstrate obvious functional brain response and its PSD demonstrate no obvious global interference sources, adaptive filtering may be unnecessary.

4.1.

Stationary Comparison

In Fig. 7 , we compare our two raw concentration time series directly by mean-subtraction and dividing each time series by its own (temporal) standard deviation. Only a small segment is shown for clarity. Here, the normalized ${\mathrm{O}}_2\mathrm{Hb}$ concentration acquired from the near measurement (S-D1) is shown as the dashed line and the computed ${\mathrm{O}}_2\mathrm{Hb}$ global interference [Fig. 3a minus Fig. 3g] is the solid line. As shown, the two time series are highly overlapping. Actually, a linear regression of the two measurements [Fig. 8a ] shows that the correlation coefficient between these two time series is 0.96 $(p<{10}^{-20})$ . Although the estimated global interference may be different from real global interference, this correlation comparison demonstrates the similarity of the reference signal to the target signal after the visual response is removed and in turn explains why adaptive filtering removes the interference term very effectively.

Fig. 7

Comparison of normalized shallow-layer hemodynamics (dashed line) and estimated global interference (solid line) in ${\mathrm{O}}_2\mathrm{Hb}$ .

Fig. 8

(a) Linear regression between the global interference in ${\mathrm{O}}_2\mathrm{Hb}$ from S-D2 (vision response subtracted) and ${\mathrm{O}}_2\mathrm{Hb}$ from S-D1. (b) A continuous plot of one $20\text{-}\mathrm{s}$ segment. The arrow indicates the direction of change as a function of time.

4.2.

Dynamic Comparison

Our previous analysis suggests that there is reasonably a linear relationship between the shallow layer hemodynamics and global interference; however, this relationship changes with time. In Fig. 8b, one segment of global interference is plotted against the shallow layer hemodynamics (from $250\mathrm{s}$ to $270\mathrm{s}$ ; both are smoothed using a sixth-order $0.25\text{-}\mathrm{Hz}$ Butterworth filter), and the arrow shows how the data moves with time. The data did not strictly follow the line. Thus, we performed short-term linear regressions—that is, we performed linear fits using only 100 data points at a time $\left(5\mathrm{s}\right)$ and slid the data window through the entire data set, with the estimated global interference as the independent time series and the near ${\mathrm{O}}_2\mathrm{Hb}$ measurement as the dependent variable. The resulting linear fit slope and intercept values oscillate over time (Fig. 9 ). This dynamic feature can also be viewed in the adaptive filtering, where the node is adjusted at a point-by-point basis and is following the slow nonstationary variations of the time series to optimize the interference cancellation. The update of the first node as a function of time is presented in Fig. 10 .

Fig. 9

Short-term linear regression coefficients. The top curve is the linear regression slope, and the bottom curve is the intercept.

Fig. 10

Adaptive filtering node update (first node).

In this case study, the fact that the shallow layer hemodynamic changes correlated well with the global interference and that adaptive filtering in general follows the nonstationary changes in the time series explains why our method should work well in suppression of the global interference and improve the CNR in the evoked brain activity detection.