_{2}) and deoxygenated (HbR) hemoglobin are measured on the forehead with multichannel NIRS during hyper- and hypocapnia. PCA and ICA are used separately to identify and remove signal contribution from extracerebral tissue, and the resulting estimates of cerebral responses are compared to the expected cerebral responses. Both methods were able to reduce extracerebral contribution to the signals, but PCA typically performs equal to or better than ICA. The improvement in 3-cm signal quality achieved with both methods is comparable to increasing the source-detector separation from 3 to 5 cm. Especially PCA appears to be well suited for use in NIRS applications where the cerebral activation is diffuse, such as monitoring of global cerebral oxygenation and hemodynamics. Performance differences between PCA and ICA could be attributed primarily to different criteria for identifying the surface effect.

## 1.

## Introduction

Medical near-infrared spectroscopy (NIRS) is a noninvasive technique for measuring cerebral blood oxygenation and volume changes using light in the
$650\text{-}\phantom{\rule{0.3em}{0ex}}\text{to}\phantom{\rule{0.3em}{0ex}}950\text{-}\mathrm{nm}$
wavelength range.^{1, 2} A light source is used to deliver light into tissue, and light exiting the tissue is measured typically at distances of approximately
$1\phantom{\rule{0.3em}{0ex}}\text{to}\phantom{\rule{0.3em}{0ex}}5\phantom{\rule{0.3em}{0ex}}\mathrm{cm}$
from the source. Hemodynamic changes can then be estimated from their effects on light propagation in tissue,^{3} e.g., during anaesthesia and surgery,^{4, 5} natural sleep,^{6, 7} hyper- or hypocapnia,^{8} or various stimuli.^{9, 10}

Parameters measured with NIRS include changes in the concentrations of oxygenated
$\left(\left[\mathrm{Hb}{\mathrm{O}}_{2}\right]\right)$
and deoxygenated ([HbR]) hemoglobin and cytochrome c oxidase,
^{11, 12} blood flow and volume,^{13, 14} and tissue oxygenation.^{15, 16, 17} NIRS is currently the only noninvasive method available for long-term monitoring of cerebral oxygenation and hemodynamics, and has also the advantages of relatively low cost and high temporal resolution.^{18} Disadvantages include low spatial resolution when used for imaging purposes due to the strong light scattering properties of tissue,^{19} sensitivity to systemic hemodynamic fluctuations caused by heartbeat and respiration,^{18} and interference from hemodynamic changes in the skin and muscle tissue (surface tissue) surrounding the cranium, especially at measurement distances shorter than
$3\phantom{\rule{0.3em}{0ex}}\text{to}\phantom{\rule{0.3em}{0ex}}4\phantom{\rule{0.3em}{0ex}}\mathrm{cm}$
.^{20, 21} Since signal artifacts caused by changes in the source-tissue coupling affect all detector channels and are often indistinguishable from physiological interference from surface tissue, the terms surface effect and extracerebral interference are used in this work to signify the combined effect of physiological interference and coupling artifacts.

Cerebrovascular reactivity to changes in blood
${\mathrm{O}}_{2}$
or
$\mathrm{C}{\mathrm{O}}_{2}$
content has been previously studied with NIRS using relatively simple source-detector (SD) configurations and methodology that is vulnerable to extracerebral interference.^{22, 23} We have in this study compared principal (PCA) and independent component analysis (ICA) in removing the surface effect from
$\left[\mathrm{Hb}{\mathrm{O}}_{2}\right]$
and [HbR] signals measured during repetitive periods of hypercapnia (increased level of blood
$\mathrm{C}{\mathrm{O}}_{2}$
) and hypocapnia (reduced level of blood
$\mathrm{C}{\mathrm{O}}_{2}$
). Since PCA is often used as a preprocessing method for ICA, the latter is typically considered to be the more powerful method. However, the performance of ICA is often sensitive to noise, and means for identifying the surface effect differ between the two methods. One of the advantages of both PCA and ICA is that they do not rely on identifying a predetermined waveform of the surface effect or the cerebral response, so they can be applied in situations where the exact shape of the response or interference varies between individuals, or the location of the response in time is not known. Despite advances in multichannel optical tomography, many clinical applications of NIRS are based on monitoring the overall oxygenation of tissue at a single or few locations to detect unexpected oxygenation changes.

To our knowledge, neither PCA nor ICA has been previously applied to hyper- or hypocapnia NIRS data, and no comparison of the two methods has been presented in the context of NIRS measurements. PCA and ICA have been previously used in topographic activation studies for removing the surface effect from NIRS signals and extracting cerebral activation signals during motor tasks.^{24, 25, 26} Since these studies have applied PCA and ICA only to datasets collected using a large number of channels (SD pairs), we also investigated whether the methods could be successfully applied in monitoring applications using a relatively small number of channels.

## 2.

## Materials and Methods

## 2.1.

### Study Protocol

The data used in this study were acquired in a series of 11 hypercapnia and 11 hypocapnia measurements on 11 human volunteers, previously published in Ref. 27. Blood
$\mathrm{C}{\mathrm{O}}_{2}$
content is an important regulator of cerebral blood flow (CBF),^{28} so that an increase (decrease) in blood
$\mathrm{C}{\mathrm{O}}_{2}$
leads to an increase (decrease) in CBF. These CBF changes should lead to parallel changes in cerebral
$\left[\mathrm{Hb}{\mathrm{O}}_{2}\right]$
and opposite changes in [HbR] if cerebral blood volume and oxygen consumption stay unchanged.

Each measurement consisted of three $2\text{-}\mathrm{min}$ periods of hyper- or hypocapnia. Since the recovery of blood $\mathrm{C}{\mathrm{O}}_{2}$ to rest level takes some time, especially after hypocapnia, the hyper- or hypocapnia periods were separated by $4\text{-}\mathrm{min}$ intervals of rest. Hypercapnia was induced by having the subject breathe $\mathrm{C}{\mathrm{O}}_{2}$ -enriched air through a mask, hypocapnia was induced by hyperventilation, and during the rest periods the subject breathed normally. In the hypocapnia measurements, end-tidal $\mathrm{C}{\mathrm{O}}_{2}$ was 2.0 to 3.0% below the rest level (approximately 5%) at the end of each hyperventilation period. In the hypercapnia measurements, a similar increase was produced. The study protocol was approved by the ethical committee of Helsinki University Hospital, and informed consent was obtained from each subject before the study.

## 2.2.

### Instrumentation

The NIRS device used in the measurements was developed in the Laboratory of Biomedical Engineering (currently the Department of Biomedical Engineering and Computational Science) at Helsinki University of Technology.
^{29, 30} The device is based on frequency domain technology, but only attenuation changes of signals were considered in this study, so the results are comparable to those obtainable with continuous wave technology.

At the beginning of each measurement, two light source fibers and ten detector fiber bundles were arranged on the forehead of the subject at $1\text{-}\mathrm{cm}$ intervals in two parallel rows (Fig. 1 ). Signal amplification was adjusted independently for each SD pair to avoid phase-amplitude cross talk at short measurement distances while maximizing signal amplification. During the measurement, light pulses were time-multiplexed at two wavelengths (760 and $830\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$ ) and two source positions sequentially. Each pulse lasted for approximately $1\phantom{\rule{0.3em}{0ex}}\mathrm{sec}$ , resulting in a $4\text{-}\mathrm{sec}$ cycle $(2\phantom{\rule{0.3em}{0ex}}\text{sources}\times 2\phantom{\rule{0.3em}{0ex}}\text{wavelengths}\times 1\phantom{\rule{0.3em}{0ex}}\mathrm{sec})$ . The raw data were averaged over each pulse, and data from different wavelengths were temporally aligned by interpolating the data to double sampling frequency (approximately $0.47\phantom{\rule{0.3em}{0ex}}\mathrm{Hz}$ ) and time-shifting data collected with one wavelength by one sample.

Light attenuation in tissue is commonly modeled in NIRS using the modified Beer-Lambert law:^{3}

## Eq. 1

$A={\mathrm{log}}_{10}\phantom{\rule{0.2em}{0ex}}\frac{{I}_{0}}{I}=d\cdot \mathit{DPF}\cdot \sum _{i}{\alpha}_{i}\cdot {c}_{i}+G,$## Eq. 2

$\left(\begin{array}{c}\Delta \left[\mathrm{Hb}{\mathrm{O}}_{2}\right]\\ \Delta \left[\mathrm{HbR}\right]\end{array}\right)=\frac{{\left({\alpha}^{T}\alpha \right)}^{-1}{\alpha}^{T}}{d}\left(\begin{array}{c}\Delta {A}_{760\phantom{\rule{0.3em}{0ex}}\mathrm{nm}}\u2215{\mathit{DPF}}_{760\phantom{\rule{0.3em}{0ex}}\mathrm{nm}}\\ \Delta {A}_{830\phantom{\rule{0.3em}{0ex}}\mathrm{nm}}\u2215{\mathit{DPF}}_{830\phantom{\rule{0.3em}{0ex}}\mathrm{nm}}\end{array}\right),$^{11}Literature values are also available for $\mathit{DPF}$ , but since frequency domain technology allows estimation of $\mathit{DPF}$ directly from calibrated phase data, we used average $\mathit{DPF}$ values over all subjects for each source fiber and wavelength. The use of data-derived $\mathit{DPF}$ values instead of literature values in the estimation of $\Delta \left[\mathrm{Hb}{\mathrm{O}}_{2}\right]$ and $\Delta \left[\mathrm{HbR}\right]$ does not affect the applicability of our results to continuous wave measurements, since population averages of $\mathit{DPF}$ were used for individual measurements.

The modified Beer-Lambert law is based on the assumption of a homogeneous medium, whereas the brain is surrounded by layers of bone and scalp. The scalp and brain form two anatomically and physiologically separate compartments, resulting in a partial volume effect where a change in
$\left[\mathrm{Hb}{\mathrm{O}}_{2}\right]$
or [HbR] in either compartment registers as a parallel but smaller change in the
$\Delta \left[\mathrm{Hb}{\mathrm{O}}_{2}\right]$
or
$\Delta \left[\mathrm{HbR}\right]$
signal. The contribution of cerebral tissue to the NIRS signal increases with SD separation, but scalp contribution is present at all SD separations.
^{21, 31}

Scalp bloodflow is affected by changes in heartbeat, respiration, and blood pressure, as well as neuronal and hormonal regulation of vasodilation and vasoconstriction. Hemodynamic changes in the scalp are typically nonlocalized, leading to a surface effect that typically affects $\Delta \left[\mathrm{Hb}{\mathrm{O}}_{2}\right]$ and $\Delta \left[\mathrm{HbR}\right]$ estimates on all NIRS channels. Since changes in the source-tissue coupling will also result in uniform attenuation changes in all signals recorded using that particular source, we define the term surface effect to signify all signal changes that are present in most or all of the measurement channels and do not originate from cerebral hemodynamic changes. To accurately estimate cerebral hemodynamics, the surface effect in the $\Delta \left[\mathrm{Hb}{\mathrm{O}}_{2}\right]$ and $\Delta \left[\mathrm{HbR}\right]$ signals has to be minimised.

The use of two sources and ten detectors produced a total of 20 $\Delta \left[\mathrm{Hb}{\mathrm{O}}_{2}\right]$ and 20 $\Delta \left[\mathrm{HbR}\right]$ signals at SD separations of approximately $1\phantom{\rule{0.3em}{0ex}}\text{to}\phantom{\rule{0.3em}{0ex}}5\phantom{\rule{0.3em}{0ex}}\mathrm{cm}$ . We have analyzed the signals separately with both PCA and ICA to identify and remove the surface effect. As the modified Beer-Lambert law only gives concentration changes relative to an arbitrarily chosen baseline level, the signals were zero meaned before applying PCA or ICA.

## 2.3.

### Principal Component Analysis

PCA is a technique commonly used for reducing the dimensionality of a dataset.^{32} It is based on interpreting a set of
$n$
measured signals
${\mathbf{x}}_{k}$
as a linear combination of
$n$
uncorrelated components
${\mathbf{y}}_{k}$
:

## Eq. 3

$\mathbf{X}=\left(\begin{array}{c}{\mathbf{x}}_{1}\\ {\mathbf{x}}_{2}\\ \vdots \\ {\mathbf{x}}_{n}\end{array}\right)=\mathbf{W}\mathbf{Y}=({\mathbf{w}}_{1}{\mathbf{w}}_{2}\dots {\mathbf{w}}_{n})\left(\begin{array}{c}{\mathbf{y}}_{1}\\ {\mathbf{y}}_{2}\\ \vdots \\ {\mathbf{y}}_{n}\end{array}\right),$Since $\mathbf{Y}$ can be expressed as $\mathbf{Y}={\mathbf{W}}^{-1}\mathbf{X}$ , the $k$ ’th principal component can be removed from the data by removing ${\mathbf{y}}_{k}$ from $\mathbf{Y}$ :

## Eq. 4

${\mathbf{X}}_{r}=\mathbf{W}{\mathbf{Y}}_{r}=\mathbf{W}{({\mathbf{y}}_{1}\dots {\mathbf{y}}_{k-1}\phantom{\rule{0.2em}{0ex}}\mathbf{0}\phantom{\rule{0.2em}{0ex}}{\mathbf{y}}_{k+1}\dots {\mathbf{y}}_{n})}^{T}.$## 2.4.

### Independent Component Analysis

Like PCA, ICA is also based on representing
$n$
signals
${\mathbf{x}}_{k}$
as a combination of
$n$
components
${\mathbf{y}}_{k}$
[Eq. 3]. However, in ICA the mixing matrix
$\mathbf{W}$
is estimated so that the statistical independency of the components is maximized. ICA algorithms are based on various methods such as maximizing the nonGaussianity of components or maximum likelihood estimation for finding the optimal
$\mathbf{W}$
. They also commonly employ a preprocessing method such as PCA or singular value decomposition (SVD) for whitening the data.^{32}

In this study we have chosen to use the second-order blind identification (SOBI) ICA algorithm.^{33} Whereas many ICA algorithms treat the data as random variables with no time structure, the SOBI algorithm attempts to minimize the lagged covariances of the components, making it robust against noise in time-dependent signals such as NIRS data. The algorithm consists of three steps: whitening the data, e.g., by using the SVD of the data covariance matrix, computing time-delayed covariance matrices for the whitened data, and joint diagonalization of the covariance matrices. For the third step we used the joint diagonalization method presented in Ref. 34. A Matlab implementation of the method is available on the Internet.^{35} The number of time lags used in calculating the covariance matrices was empirically set to 100.

A common feature of most ICA algorithms is that they do not provide any information on the relevance of the individual components (compare
${\lambda}_{k}$
in PCA).^{32} However, the spatial distribution of the components can still be used to identify the surface component. In Ref. 25, a coefficient of spatial uniformity (CSU) was defined for component
$k$
as

## 2.5.

### Effect of Channel Configuration on Principal and Independent Component Analysis

Previous NIRS studies using PCA and ICA have used dozens of SD pairs covering a large portion of the head surface.^{24, 25, 26} Also, they have generally attempted to detect localized changes caused by activation of a specific area of the brain (e.g., motor cortex). In this case, the surface tissue signal is expected to affect all SD pairs, while the cerebral signal is present only in a small number of pairs. In contrast, the cerebral hemodynamic changes caused by hyper- or hypocapnia would be present to some extent in a majority of the SD pairs, which might make it difficult to automatically distinguish between the surface and cerebral components.

We tested whether PCA and ICA would work if a minimal number (two) of SD pairs was used, so that one pair represented primarily surface tissue and the other contained also a cortical component. This would correspond to the monitoring of cerebral oxygenation with simple equipment in a clinical setting. We used as input for PCA and ICA $\Delta \left[\mathrm{Hb}{\mathrm{O}}_{2}\right]$ and $\Delta \left[\mathrm{HbR}\right]$ values from source 1 and detectors 1 and 3. Using detector 3, which probes only a relatively small volume of cortical tissue, allows benchmarking of filtering results against unfiltered signals recorded at SD separations of $3\phantom{\rule{0.3em}{0ex}}\text{to}\phantom{\rule{0.3em}{0ex}}5\phantom{\rule{0.3em}{0ex}}\mathrm{cm}$ . For reference purposes, we also carried out the analyses using all 20 SD pairs.

## 3.

## Results

## 3.1.

### Original Measurements

Figure 2 shows means and STDs of unfiltered $\Delta \left[\mathrm{Hb}{\mathrm{O}}_{2}\right]$ and $\Delta \left[\mathrm{HbR}\right]$ from three SD pairs over all hypercapnia measurements. These three SD pairs together provide a representative view of the waveforms of signals measured at different distances from the sources, so signals from other detectors are omitted to simplify the graphical presentation. The response to hypercapnia in the $5\text{-}\mathrm{cm}$ signals corresponds to the expected cerebral response, whereas the $1\text{-}\mathrm{cm}$ signals do not exhibit any sign of such response. The $3\text{-}\mathrm{cm}$ signals show features from both the 1- and $5\text{-}\mathrm{cm}$ signals, but distinguishing the cerebrovascular $\Delta \left[\mathrm{Hb}{\mathrm{O}}_{2}\right]$ reaction to hypercapnia from other fluctuations in the $3\text{-}\mathrm{cm}$ signal would be difficult without knowledge of the hypercapnia periods.

Figure 3 shows corresponding averages over all hypocapnia measurements. The signal changes at $5\phantom{\rule{0.3em}{0ex}}\mathrm{cm}$ are consistent with the expected cerebrovascular reaction to a rapid hyperventilation-induced decrease and gradual recovery of blood $\mathrm{C}{\mathrm{O}}_{2}$ content. In contrast, although the $3\text{-}\mathrm{cm}$ signals also include a cerebral component, the hyperventilation periods could not be reliably identified from those signals alone.

In both hyper- and hypocapnia, the STD of the $1\text{-}\mathrm{cm}$ $\Delta \left[\mathrm{Hb}{\mathrm{O}}_{2}\right]$ signals is two to four times larger than the STD of the 3- and $5\text{-}\mathrm{cm}$ signals, and there are no consistent responses to hyper- or hypocapnia in the mean $1\text{-}\mathrm{cm}$ signals. Since the contribution of surface tissue to the signal decreases with increasing SD separation, these findings indicate that the surface effect is in practice independent of blood $\mathrm{C}{\mathrm{O}}_{2}$ content, and its time behavior varies greatly between individuals.

## 3.2.

### Filtered Signals

We used PCA and ICA for two different channel configurations to filter (i.e., remove the surface component from) the $3\text{-}\mathrm{cm}$ signals corresponding to source 1 and detector 3. Means and STDs of the original, unfiltered signals are shown in Figs. 2b and 3b. The filtering was first carried out with two SD pairs, with means and STDs of the filtered signals shown in Figs. 4a and 4c (hypercapnia) and Figs. 5a and 5c (hypocapnia). Both PCA and ICA enhance the cerebral $\Delta \left[\mathrm{Hb}{\mathrm{O}}_{2}\right]$ waveform in hypercapnia, but in the case of hypocapnia, PCA enhances the cerebral waveform of both $\Delta \left[\mathrm{Hb}{\mathrm{O}}_{2}\right]$ and $\Delta \left[\mathrm{Hb}\mathrm{R}\right]$ , while the improvement in signal quality from ICA is negligible.

When the filtering was repeated with 20 SD pairs, results from ICA were clearly inferior to PCA. We suspected that the performance of ICA might be related to the number of independent components estimated, so we reduced the number from 20 to 5 by reducing the dimensionality of the data in the whitening step of ICA. This is commonly done when the number of channels is much larger than the expected number of actual signals present in the data. The choice to use exactly five components was to some extent arbitrary, since there was in practice no difference between results obtained using three to ten components, and two components produced clearly inferior results with hypocapnia data. ICA performed slightly better after dimensionality reduction, producing in hypercapnia results equal to 20-channel PCA; but in hypocapnia, results from ICA remained inferior to PCA. The results of PCA filtering and five-component ICA filtering of the $\Delta \left[\mathrm{Hb}{\mathrm{O}}_{2}\right]$ and $\Delta \left[\mathrm{Hb}\mathrm{O}\mathrm{R}\right]$ signals corresponding to source 1 and detector 3 with 20 SD pairs are shown in Figs. 4c and 4d and Figs. 5c and 5d.

In most cases, signal changes and STDs are smaller in the filtered $3\text{-}\mathrm{cm}$ signals than in the unfiltered $5\text{-}\mathrm{cm}$ signals. This is an expected consequence of the partial volume effect in the modified Beer-Lambert law: since the $3\text{-}\mathrm{cm}$ signals sample a smaller volume of cerebral tissue than the $5\text{-}\mathrm{cm}$ signals, the magnitude of concentration changes in the $3\text{-}\mathrm{cm}$ signals after surface component removal should be smaller than in the original $5\text{-}\mathrm{cm}$ signals.

## 3.3.

### Quantitative Performance Evaluation

For quantitative performance evaluation, an estimate of the cerebral response during the measurement period is required. The $5\text{-}\mathrm{cm}$ $\Delta \left[\mathrm{Hb}{\mathrm{O}}_{2}\right]$ and $\Delta \left[\mathrm{HbR}\right]$ signals are affected by both the surface effect and cortical hemodynamics, but since there is no consistent response in surface tissue to hyper- or hypocapnia [Figs. 2a and 3a], averaging the $5\text{-}\mathrm{cm}$ signals over all measurements effectively cancels out the surface effect and provides the best available estimate of the cerebral response. We therefore used the averages of the unfiltered $5\text{-}\mathrm{cm}$ signals [Figs. 2c and 3c] to represent the expected cerebral response ${x}_{\mathrm{exp}}$ in hyper- and hypocapnia, respectively.

To provide a quantitative measure of error, we calculated the mean squared error (MSE) between ${x}_{\mathrm{exp}}$ and the PCA or ICA filtered $3\text{-}\mathrm{cm}$ signal $x$ for the two channel configurations:

## Eq. 6

$\mathrm{MSE}=\frac{1}{N}\sum _{i=1}^{N}{[x\left(i\right)-{x}_{\mathrm{exp}}\left(i\right)]}^{2},$Figure 6 shows the mean and STD of the MSE criterion for the filtered $3\text{-}\mathrm{cm}$ signals and the unfiltered 3- and $5\text{-}\mathrm{cm}$ signals in the two SD configurations and breathing conditions. To find out whether the MSEs of the filtered $3\text{-}\mathrm{cm}$ signals and the unfiltered $5\text{-}\mathrm{cm}$ signals were lower than the MSE of the unfiltered $3\text{-}\mathrm{cm}$ signals in a statistically significant sense, the corresponding MSE distributions were compared using a Student’s $t$ -test for unequal sample variances. For 20-channel ICA, MSE values are shown for the five-component case, but similar results were obtained with three to ten components (compare Sec. 3.2).

Overall, both PCA and ICA reduced extracerebral contribution to the $3\text{-}\mathrm{cm}$ signal, but PCA always produced lower MSEs than ICA. The lowest MSEs in both breathing conditions were given by two-channel PCA. Results from two-channel ICA were almost equal to two-channel PCA in hypercapnia, but were clearly inferior in hypocapnia. In the 20-channel case, results from PCA and ICA were almost similar after the performance of ICA was improved by dimensionality reduction. Two-channel ICA performed poorer than 20-channnel ICA in hypocapnia, but better in hypercapnia.

PCA produced a statistically significant MSE reduction in six out of eight cases at the $p<0.05$ level, and in eight out of eight at the $p<0.10$ level (Fig. 6). The corresponding figures for ICA were only one out of eight cases at $p<0.05$ and two out of eight at $p<0.10$ . The MSE difference between the unfiltered 3- and $5\text{-}\mathrm{cm}$ signals was statistically significant in only one out of four cases at the $p<0.10$ level. The average MSEs of the PCA filtered $3\text{-}\mathrm{cm}$ signals were in seven out of eight cases lower than those of the unfiltered $5\text{-}\mathrm{cm}$ signals, while ICA filtered $3\text{-}\mathrm{cm}$ signals had a lower MSE than the unfiltered $5\text{-}\mathrm{cm}$ signals in four out of eight cases. However, these MSE differences were typically not statistically significant.

Qualitative analysis supports the results of the quantitative MSE analysis. The shape of the cerebral response to hypercapnia is enhanced more in the two-channel case than in the 20-channel case, especially by PCA (Fig. 4). In hypocapnia, the best result is achieved with two-channel PCA, while ICA filtering completely fails to enhance the shape of the cerebral response, regardless of channel configuration (Fig. 5).

## 3.4.

### Comparison of the Coefficient of Spatial Uniformity and $\lambda $ Criteria

To determine how unambiguously the surface component could be identified with different methods, we examined the distributions of ${\lambda}_{k}$ and $\mathrm{CSU}\left(k\right)$ over the measurements. If the data originates from similar processes in all measurements, and if the mean value of ${\lambda}_{2}$ is much smaller than ${\lambda}_{1}$ , and STD of ${\lambda}_{2}$ is small compared to the mean, it is likely that the first PCA component originates from the same process in all measurements. Similar argumentation can be applied to the CSU criterion.

Means and STDs of ${\lambda}_{2}$ and CSU(2) in the two-channel configuration are shown in Table 1 . The mean values of both ${\lambda}_{2}$ and CSU(2) are less than half of ${\lambda}_{1}$ and CSU(1), respectively, and ${\lambda}_{2}$ is significantly smaller than CSU(2) ( $p<0.05$ using student’s $t$ -test for unequal variances).

## Table 1

Means and STDs of λ2 and CSU(2) in the two-channel configuration, when λ1 and CSU(1) have been normalized to unity.

CSU(2) | λ2 | |
---|---|---|

$\mathrm{Hb}{\mathrm{O}}_{2}$ , hypercapnia | $0.26\pm 0.28$ | $0.06\pm 0.04$ |

HbR, hypercapnia | $0.47\pm 0.23$ | $0.15\pm 0.13$ |

$\mathrm{Hb}{\mathrm{O}}_{2}$ , hypocapnia | $0.32\pm 0.26$ | $0.07\pm 0.05$ |

HbR, hypocapnia | $0.25\pm 0.22$ | $0.09\pm 0.09$ |

Mean and STD of the five largest $\mathrm{CSU}\left(k\right)$ and ${\lambda}_{k}$ in the case of 20 SD pairs are shown in Fig. 7 . The distributions of $\mathrm{CSU}\left(k\right)$ and ${\lambda}_{k}$ are fundamentally different, and CSU(2) and even CSU(3) are much closer to CSU(1) than ${\lambda}_{2}$ is to ${\lambda}_{1}$ . Also, the STDs are larger for $\mathrm{CSU}\left(k\right)$ than for ${\lambda}_{k}$ . $\mathrm{CSU}\left(k\right)$ are shown for ICA with the number of data dimensions reduced to five. Without dimensionality reduction, the mean value of CSU(2) was approximately 0.8 and that of CSU(3) approximately 0.6. This suggests that in the case of ICA, the CSU criterion for identifying the surface component is prone to error when the number of channels increases. On the other hand, in the case of PCA, the good filtering results and small mean value of ${\lambda}_{2}$ confirm the assumption that the first PCA component corresponds to the surface effect.

## 4.

## Discussion

Extracerebral contamination of signals is a fundamental problem of NIRS, and several methods have been developed to resolve the issue. The depth sensitivity of NIRS can be improved, e.g., by time domain spectroscopy, where rapid light pulses are delivered into tissue and detected with a time-gated system to screen for photons with longer flight times, as they are more likely to have probed deeper into the tissue.^{36} Time domain spectroscopy allows estimating hemodynamic changes in different tissue layers with a single-detector setup,
^{37} reducing the number of detectors in NIRS imaging applications,^{38} and improving the contrast-to-background ratio for activation measurements.^{39} However, most commercial NIRS devices are still based on continuous wave technology, most likely due to the higher costs and more complicated signal processing required by time or frequency domain technology.^{40}

Adaptive filtering has been proposed for removing noncerebral interference from NIRS signals. Such algorithms typically attempt to minimize the correlation of NIRS signals and reference signals such as blood pressure, respiration, or short-range NIRS measurements.
^{41, 42, 43, 44, 45} The advantages of adaptive filtering include computational simplicity, real-time application, and the ability to adapt to changes in the interfering signals. However, if the reference signals are also correlated with the cerebral signal, relevant information may be removed by the filter. For example, when using a short-range NIRS measurement for the reference signal, light can in some cases penetrate into the outer layers of the brain. This can occur, e.g., as a result of individual anatomical variations, especially in infants who have thin scalps and skulls. A second limitation that also affects PCA and ICA is the possible coupling of cerebral and noncerebral hemodynamics, e.g., in sensorimotor studies if the stimulus or task causes a change in blood pressure or respiration simultaneously with cerebral activation. Thus, the selection of the appropriate filtering method for each application requires careful testing.

Another approach to surface effect reduction is to model the absorption changes in different tissue layers separately. The simplest models feature only two slab-like layers, with the other one representing extracerebral and the other one cerebral tissue.
^{46, 47} However, these models require either *a priori* information on tissue anatomy, or frequency or time domain measurements to estimate the path length of light in tissue. Least-squares fitting may in some cases be used to circumvent such problems.
^{48} More sophisticated tissue models utilizing anatomical or physiological *a priori* information, e.g., from magnetic resonance imaging (MRI) or computed tomograhpy (CT), have been developed for tomographic NIRS imaging.
^{49, 50, 51, 52} However, obtaining MRI or CT images for routine monitoring is typically impractical or impossible, and a simple measurement geometry and ease of measurement are often more valuable criteria than high spatial resolution.

PCA has been used in NIRS to estimate and remove systemic interference from signals recorded during sensorimotor stimuli.^{24} The authors used relatively long SD separations, approximately
$3\phantom{\rule{0.3em}{0ex}}\text{to}\phantom{\rule{0.3em}{0ex}}4\phantom{\rule{0.3em}{0ex}}\mathrm{cm}$
, and assumed that the spatial pattern of the interference is stable in time. Baseline data were then collected prior to the stimuli to estimate the interference subspace. In another sensorimotor activation study, the long-distance
$\Delta \left[\mathrm{Hb}{\mathrm{O}}_{2}\right]$
signal, but not the
$\Delta \left[\mathrm{HbR}\right]$
signal, was shown to correlate strongly with extracerebral tissue hemodynamics.^{31} PCA has also been used with simulated data to determine the optimal SD configuration for NIRS measurements.^{53}

Applications of ICA to NIRS data can be divided into two categories. One approach is to identify a component corresponding to cerebral hemodynamics, typically the hemodynamic response to cortical activity. This requires either assumptions on the shape of the hemodynamic response,^{54} knowledge of its location in time,^{42} or manual screening of the components.^{26} The other approach, also taken in this study, is based on identifying the surface component, e.g., based on its spatial distribution, and removing it.
^{25}

ICA is a nonspecific term covering several different algorithms using different measures for statistical independence. To our knowledge, no comparative study of different ICA algorithms in NIRS data processing has been made. Our choice of the SOBI algorithm for this study is based on its general applicability to noisy time-dependent signals, and it has been found to perform better than PCA, e.g., in ocular artifact reduction from electroencephalographic (EEG) recordings.^{55} The SOBI algorithm has been used in a previous NIRS study for detection of motor cortex activity.^{42} We also tested the FastICA algorithm, which was used for NIRS data in Ref. 54, but preliminary results were inferior to SOBI.

The use of PCA or ICA for surface effect removal relies on the assumption that cerebral hemodynamics are uncorrelated with or independent of scalp tissue circulation, which may not hold in all cases. For example, in sensorimotor activation studies, the stimulus causing the activation may also cause a simultaneous change in heart or respiration rate, whereas in our study the hyperventilation task may have a similar effect. However, both PCA and ICA have been shown to perform well in activation studies,^{24, 26} and our results strongly suggest that they can be used to improve the detection of cerebrovascular responses to hyper- and hypocapnia. In particular, we did not observe any consistent extracerebral hemodynamic response to changes in blood
$\mathrm{C}{\mathrm{O}}_{2}$
level. Thus, the level of physiological coupling between surface and cerebral circulation appears to be in practice low enough to allow the use of PCA and ICA for surface effect removal in various experimental settings.

We compared PCA and ICA in removing surface tissue interference from NIRS signals measured during hypercapnia and hypocapnia. We carried out the analyses with two different channel configurations: a minimal two-channel configuration that might be used in a clinical monitoring setting, and a more extensive 20-channel configuration maximizing the amount of information gathered. We found that both PCA and ICA were able to reduce the MSE between the $3\text{-}\mathrm{cm}$ signal and the expected cerebrovascular response, modeled by the averaged $5\text{-}\mathrm{cm}$ signal, but PCA performed consistently equal to or better than ICA. Except for ICA in hypocapnia, the qualitative and quantitative improvement in $3\text{-}\mathrm{cm}$ signal quality achieved with both methods was typically comparable to increasing the SD separation from $3\phantom{\rule{0.3em}{0ex}}\text{to}\phantom{\rule{0.3em}{0ex}}5\phantom{\rule{0.3em}{0ex}}\mathrm{cm}$ , with best results obtained with two-channel PCA.

Since ICA typically contains PCA as a preprocessing stage and is considered statistically the more powerful method, the better performance of PCA is unexpected. We believe that this result is explained by the criteria used in PCA and ICA for identifying the surface component. The PCA method utilizes eigenvalues of the data covariance matrix, which contain information on the strength of the individual components in the data. Our results show that since the surface effect is present to some extent on all channels, the corresponding principal component is considerably stronger than any other components, even when the cerebral hemodynamic change is relatively strong and global, as in hyper- or hypocapnia. ICA does not inherently provide information on the significance of individual components, so the surface effect has to be identified based on its spatial distribution, which should be relatively uniform. The resulting CSU criterion, introduced in Ref. 25, appears to be less applicable for surface effect removal than the PCA eigenvalue criterion.

There are two possible reasons for the failure of the CSU criterion. First, the spatial pattern of cerebral hemodynamic changes caused by hyper- and hypocapnia is more diffuse than in the original study.^{25} From the 20 SD pairs recorded in this study, the cerebral signal is clearly present in perhaps 8 to 12, depending on the measurement. This could result in erroneously identifying the cerebral component as extracerebral interference in some measurements. Unfortunately, no statistical data on the distribution of CSU in different measurements was provided in Ref. 25. Therefore, it is difficult to determine to what extent the spatial distribution of cerebral activation might affect our results. Second, there are in practice only two signal sources present in our data: the cerebral response and the extracerebral interference. As the number of channels analyzed increases, the number of ICA components essentially estimating random noise also increases. Consequently, the probability that at least one of these random components has a spatial distribution more uniform than the extracerebral component increases, especially since the spatial distribution of the surface effect may be slightly skewed toward short SD separations due to the partial volume effect at long SD separations. In this case, removing the component with the highest CSU would not necessarily remove the surface effect. Reducing the dimensionality of our data both improved the performance of ICA and reduced the value of CSU(2) in relation to CSU(1). Therefore, the number of estimated components appears to be at least a contributing factor in the performance of ICA with the CSU criterion.

Since many ICA methods are parameter-dependent, finding the optimal set of parameters for each application could improve the results of ICA. For example, in our case, decreasing the dimensionality of data from 20 to 5 improved the performance of 20-channel ICA to the level of 20-channel PCA. Generally, our results suggest that increasing the number of channels decreases the performance of both PCA and ICA, and our tests with other channel configurations than the two reported here support this conclusion. For example, using source 1 and detectors 1 through 5 produced a larger MSE for both PCA and ICA than the two-channel configuration presented before, although including the $5\text{-}\mathrm{cm}$ signal could have been expected to facilitate distinguishing between the cerebral and surface components. One possible reason for this may be SD pair-specific variations in the surface effect resulting from changes in the detector-tissue coupling and differences in scalp vasculature under the detectors. As the number of channels and total area probed increase, such as in imaging applications, modeling the surface effect with only one PCA or ICA component becomes increasingly difficult compared to our application, where the monitoring can be confined to a comparatively small area. The performance of both PCA and ICA might be improved by manual selection of the cerebral or surface component(s), but this could also introduce selection bias and would increase the amount of manual work required as the number of channels increases.

## 5.

## Conclusions

Both PCA and ICA have been previously used for removing extracerebral interference from NIRS signals collected with several channels, but our results indicate that the methods are also suitable for applications where a small number of channels are available or desirable and the pattern of cerebral activation is diffuse, such as monitoring cerebral oxygenation or determining the cerebrovascular response to changes in blood $\mathrm{C}{\mathrm{O}}_{2}$ content. No benefit from choosing ICA over PCA is observed, and in all quantitative MSE comparisons, PCA performs better than ICA. Performance differences between the two methods could be attributed primarily to different criteria for identifying the surface effect.

## Acknowledgments

This work was financially supported by the Instrumentarium Foundation and the Research Foundation of Helsinki University of Technology. The original measurements were carried out at the BioMag Laboratory in Helsinki University Hospital.

## References

**,” Science, 198 1264 –1267 (1977). https://doi.org/10.1126/science.929199 0036-8075 Google Scholar**

*Noninvasive, infrared monitoring of cerebral and myocardial oxygen sufficiency and circulatory parameters***,” IEEE Signal Process. Mag., 18 57 –75 (2001). https://doi.org/10.1109/79.962278 1053-5888 Google Scholar**

*Imaging the body with diffuse optical tomography***,” Advances in Electromagnetic Fields in Living Systems, 4 Springer Science, New York (2005). Google Scholar**

*Diffuse optical imaging***,” Anesth. Analg. (Baltimore), 99 (5), 1365 –1375 (2004). https://doi.org/10.1213/01.ANE.0000134808.52676.4D 0003-2999 Google Scholar**

*Neurological monitoring for congenital heart surgery***,” J. Clin. Anesth., 17 426 –430 (2005). https://doi.org/10.1016/j.jclinane.2004.09.007 0952-8180 Google Scholar**

*A clinical evaluation of near-infrared cerebral oximetry in the awake patient to monitor cerebral perfusion during carotid endarterectomy***,” Chest, 109 916 –921 (1996). https://doi.org/10.1378/chest.109.4.916 0012-3692 Google Scholar**

*Changes in cerebral oxygenation and hemodynamics during obstructive sleep apneas***,” Brain Res., 866 313 –325 (2000). https://doi.org/10.1016/S0006-8993(00)02320-9 0006-8993 Google Scholar**

*Intracerebral hemodynamic probed by near infrared spectroscopy in the transition between wakefulness and sleep***,” Brain Res., 1107 206 –214 (2006). https://doi.org/10.1016/j.brainres.2006.06.002 0006-8993 Google Scholar**

*Cerebrovascular reactivity to hypercapnia in migraine patients measured with near-infrared spectroscopy***,” Neurosci. Lett., 316 75 –78 (2001). https://doi.org/10.1016/S0304-3940(01)02372-2 0304-3940 Google Scholar**

*Variation of temporal characteristics in human cerebral hemodynamic responses to electric median nerve stimulation: a near-infrared spectroscopic study***,” NeuroReport, 16 1373 –1377 (2005). https://doi.org/10.1097/01.wnr.0000175247.35837.15 0959-4965 Google Scholar**

*Bilateral hemodynamic responses to auditory stimulation in newborn infants***,” Univ. College of London, (1991). Google Scholar**

*The application of near infrared spectroscopy to non invasive monitoring of cerebral oxygenation in the newborn infant***,” J. Biomed. Opt., 12 (2), 024002 (2007). https://doi.org/10.1117/1.2718541 1083-3668 Google Scholar**

*Near-infrared spectroscopic quantification of changes in the concentration of oxidized cytochrome c oxidase in the healthy human brain during hypoxemia***,” Adv. Exp. Med. Biol., 317 235 –245 (1992). 0065-2598 Google Scholar**

*Measurement of cerebral blood flow in adult humans using near infrared spectroscopy–methodology and possible errors***,” Pediatr. Neurol., 23 323 –327 (2000). https://doi.org/10.1016/S0887-8994(00)00195-8 0887-8994 Google Scholar**

*Hemodynamic responses to photic stimulation in neonates***,” Ped. Anesthesia, 18 160 –166 (2008). Google Scholar**

*Comparison of different near-infrared spectroscopic cerebral oxygenation indices with central venous and jugular venous oxygenation saturation in children***,” (2003). Google Scholar**

*Method for spectrophotometric blood oxygenation monitoring***,” J. Biomed. Opt., 12 (6), 062106 (2007). https://doi.org/10.1117/1.2804911 1083-3668 Google Scholar**

*Functional near-infrared spectroscopy: current status and future prospects***,” Opt. Lett., 29 1506 –1508 (2004). https://doi.org/10.1364/OL.29.001506 0146-9592 Google Scholar**

*Improving the diffuse optical imaging spatial resolution of the cerebral hemodynamic response to brain activation in humans***,” Neuroimage, 8 69 –78 (1998). https://doi.org/10.1006/nimg.1998.0348 1053-8119 Google Scholar**

*A theoretical study of the signal contribution of regions of the adult head to near-infrared spectroscopy studies of visual evoked responses***,” J. Clin. Monit Comput., 14 353 –360 (1998). https://doi.org/10.1023/A:1009957032554 1387-1307 Google Scholar**

*Sensitivity of near infrared spectroscopy to cerebral and extra-cerebral oxygenation changes is determined by emitter-detector separation***,” Phys. Med. Biol., 50 R1 –R43 (2005). https://doi.org/10.1088/0031-9155/50/4/R01 0031-9155 Google Scholar**

*Recent advances in diffuse optical imaging***,” J. Biomed. Opt., 12 (6), 062113 (2007). https://doi.org/10.1117/1.2804705 1083-3668 Google Scholar**

*Age effects on brain oxygenation during hypercapnia***,” J. Biomed. Opt., 10 (1), 011014 (2005). https://doi.org/10.1117/1.1852552 1083-3668 Google Scholar**

*Eigenvector-based spatial filtering for reduction of physiological interference in diffuse optical imaging***,” 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***,” Neuroimage, 41 (Suppl. 1), c1 (2008). 1053-8119 Google Scholar**

*Blind ICA discrimination of evoked cortical responses in humans with DOT***,” Helsinki Univ. of Technology, (2004). Google Scholar**

*Instrumentation for diffuse optical imaging and near-infrared spectroscopy and their applications in human brain studies***,” Anesthesiology, 88 1365 –1386 (1998). https://doi.org/10.1097/00000542-199805000-00029 0003-3022 Google Scholar**

*Carbon dioxide and the cerebral circulation***,” Rev. Sci. Instrum., 73 3306 –3312 (2002). https://doi.org/10.1063/1.1497496 0034-6748 Google Scholar**

*Instrumentation for the accurate measurement of phase and amplitude in optical tomography***,” Rev. Sci. Instrum., 76 044302 (2005). https://doi.org/10.1063/1.1884193 0034-6748 Google Scholar**

*Instrumentation and calibration methods for the multichannel measurement of phase and amplitude in optical tomography***,” Neuroimage, 41 (Suppl. 1), e1 (2008). https://doi.org/10.1016/S1053-8119(08)70004-1 1053-8119 Google Scholar**

*Measurement of brain activation using near-infrared spectroscopy: comparison of principal components for signal changes between short and long source-detector spacings***,” IEEE Trans. Signal Process., 45 434 –444 (1997). https://doi.org/10.1109/78.554307 1053-587X Google Scholar**

*A blind source separation technique using second-order statistics***,” SIAM J. Matrix Anal. Appl., 17 161 –164 (1996). https://doi.org/10.1137/S0895479893259546 0895-4798 Google Scholar**

*Jacobi angles for simultaneous diagonalization***,” (2009) http://www.tsi.enst.fr/~cardoso/jointdiag.html Google Scholar**

*Joint diagonalization***,” Rev. Sci. Instrum., 70 193 –201 (1999). https://doi.org/10.1063/1.1149565 0034-6748 Google Scholar**

*Multichannel photon counting instrument for spatially resolved near infrared spectroscopy***,” Phys. Med. Biol., 46 879 –896 (2001). https://doi.org/10.1088/0031-9155/46/3/320 0031-9155 Google Scholar**

*Determining changes in NIR absorption using a layered model of the human head***,” J. Biomed. Opt., 10 064008 (2005). https://doi.org/10.1117/1.2136312 1083-3668 Google Scholar**

*Extraction of depth-dependent signals from time-resolved reflectance in layered turbid media***,” J. Biomed. Opt., 10 011013 (2005). https://doi.org/10.1117/1.1852553 1083-3668 Google Scholar**

*Improved sensitivity to cerebral hemodynamics during brain activation with a time-gated optical system: analytical model and experimental validation***,” J. Biomed. Opt., 12 062104 (2007). https://doi.org/10.1117/1.2804899 1083-3668 Google Scholar**

*Progress of near-infrared spectroscopy and topography for brain and muscle clinical applications***,” Phys. Med. Biol., 48 1491 –1504 (2003). https://doi.org/10.1088/0031-9155/48/11/301 0031-9155 Google Scholar**

*Time-series estimation of biological factors in optical diffusion tomography***,” Med. Biol. Eng. Comput., 42 92 –99 (2004). https://doi.org/10.1007/BF02351016 0140-0118 Google Scholar**

*Detection of fast neuronal signals in the motor cortex from functional near infrared spectroscopy measurements using independent component analysis***,” Neuroimage, 30 88 –101 (2006). https://doi.org/10.1016/j.neuroimage.2005.09.016 1053-8119 Google Scholar**

*Dynamic physiological modeling for functional diffuse optical tomography***,” 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 case study***,” 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***,” 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 spectroscopic approach***,” 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***,” Phys. Med. Biol., 44 2703 –2721 (1999). https://doi.org/10.1088/0031-9155/44/11/302 0031-9155 Google Scholar**

*Optical tomographic reconstruction in a complex head model using a priori region boundary information***,” Appl. Opt., 42 2881 –2887 (2003). https://doi.org/10.1364/AO.42.002881 0003-6935 Google Scholar**

*Monte Carlo prediction of near-infrared light propagation in realistic adult and neonatal head models***,” Appl. Opt., 44 2049 –2057 (2005). https://doi.org/10.1364/AO.44.002049 0003-6935 Google Scholar**

*Modeling anisotropic light propagation in a realistic model of the human head***,” Rev. Sci. Instrum., 77 114301 (2006). https://doi.org/10.1063/1.2364138 0034-6748 Google Scholar**

*Integrated measurement system for simultaneous functional magnetic resonance imaging and diffuse optical tomography in human brain mapping***,” 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***,” Med. Biol. Eng. Comput., 44 945 –958 (2006). https://doi.org/10.1007/s11517-006-0116-3 0140-0118 Google Scholar**

*Extraction of cognitive activity-related waveforms from functional near-infrared spectroscopy signals***,” Comput. Biol. Med., 38 348 –360 (2008). https://doi.org/10.1016/j.compbiomed.2007.12.001 0010-4825 Google Scholar**

*A comparative study of automatic techniques for ocular artifact reduction in spontaneous EEG signals based on clinical target variables: a simulation case*