Since it was proposed by Jobsis,1 functional near-infrared spectroscopy (fNIRS) has emerged as a low-cost and noninvasive technique for quantifying brain functional activities.2 This neuroimaging technique emits near-infrared light on a subject’s head using a source optode and then detects the scattered and diffused light through the head using a detector optode near the source optode. The attenuation of the detected light is linearly related to the concentration of chromophores, i.e., oxy- and deoxyhemoglobin in the subject’s brain. This relationship can be expressed by a modified Beer–Lambert law.3,4
fNIRS data reflect neuronal activity indirectly because this activity is accompanied by the metabolic demand for the oxygen carried by hemoglobin whose concentrations are calculated from fNIRS data.5,6 Because it allows for the noninvasive investigation of neuronal activities, fNIRS studies have now been performed in clinical situations, psychiatric evaluations, studies of cognition, and investigations of brain development.22.214.171.124.12.13.–14 These applications require that the fNIRS data be processed effectively and provide adequate information about the hemodynamics of a subject’s brain. However, obtaining the desired hemodynamic information from fNIRS data is challenging because fNIRS mixes a lot of physiological noise, including cardiac pulsation, respiration, and mean arterial blood pressure variations,1516.17.–18 with the desired data. Moreover, these physiological noises are not organized in an orderly fashion but are mixed in a complicated way. They are periodic and easily correlated with each other or with the desired hemodynamic response.19,20 In order to investigate activated areas in a subject’s brain during a task, it is necessary to separate the hemodynamic response from the physiological noises. Due to the difficulty of accurate separation, activation detection is one of the challenges in fNIRS, although it is a necessary step before analyzing brain function.
Two kinds of methods have been attempted for activation detection in fNIRS. The first, which was a model-driven method in which a variety of regressors were utilized to represent distinct physiological signals, adopted a general linear model (GLM) or adaptive GLM.21 For example, a polynomial basis function was used to describe instrumental and physiological drifts, and Fourier basis functions have been used to represent periodic signals, such as the respiration signal.2223.–24 Short-distance channels have also been utilized as one of the regressors to reduce physiological noises in fNIRS data, but this required extra measurements by pairs of source and detector in a short distance.2526.–27 The second type of method estimated the hemodynamic response function (HRF) at each channel by first using a deconvolution algorithm. Then the peak amplitude was used to identify the activated channels.16,28 The GLM method was comparatively easy to implement, but in practice, the estimation used a linear single scale model and the remaining error was difficult to simulate as Gaussian noise. The HRF-based method was more difficult to use because it was necessary to solve an ill-posed deconvolution problem to estimate HRF in each channel. The ill-posed estimation of HRF required an iterative way to obtain a solution, which made it difficult to get an accurate solution.
Because wavelet transform (WT) could decompose signals into a time-frequency domain, the transform has been utilized in fNIRS to process data, such as to detrend data29 and remove motion artifacts.30 The wavelet coherence was proposed to detect the correlation between two signals in the dual domain.31 Typical applications by the coherence analysis in fNIRS were to analyze interpersonal coherence in NIRS-based hyperscanning,32 to quantify the impact of physiological processes on fNIRS signals20 and to evaluate prefrontal functional connectivity between two hemispheres of patients’ brains. These applications took advantage of the coherence, which was investigated in the time-frequency domain. In this paper, we describe an activation detection method based on wavelet coherence. This method calculated the coherence between the predicted response and the concentrations of oxy- or deoxyhemoglobin, as obtained by fNIRS, during a task. The coherence was distributed in the time-frequency domain and decomposed across multiple scales after WT. This differed from the GLM method, which analyzes signals in a single domain and single scale. The predicted response was estimated by convolving the stimuli with the canonical HRF. Because coherence is a local measure of the correlation between two signals, the time-frequency coherence reflects the local correlation between the predicted response and the concentrations of hemoglobin at different times and frequencies. Therefore, our method was able to be used to detect activation by an elaborate analysis of the correlation. Its performance was also demonstrated by simulated and experimental data, as reported in the third section of the paper.
Optimally localizing the signal in the time-frequency domain plays a key role in signal analysis. The level of the localization can be used to evaluate the performance of a function in analyzing a signal in detail. If an appropriate balance between the duration in the frequency domain and the duration in the time domain does not exist, a function will fail to achieve the minimum of the uncertainty product, , and the function will not be able to analyze the details of the signal in either domain.
A wavelet is a function with a mean of zero. It is characterized by limited duration (localization) in both the frequency and time domains. Unlike the wavelet, the Fourier basis is only localized in the frequency domain, while it expands to infinity in the time domain.33 According to the Heisenberg uncertainty principle, a tradeoff always exists between localizations in the two domains. Wavelets have been studied to optimize the uncertainty product. Therefore, it can analyze signals optimally in both the time and frequency domains.34 Here, we employed Morlet WT to analyze signals measured by fNIRS. This wavelet includes a plane wave modulated by a Gaussian function as described by
The WT can decompose a signal in the time-frequency domain and indicate the power spectrum of a signal in terms of time and frequency. Because of its dual localizations, a wavelet can divide a given function into different frequency and time components and then investigate them individually. In the transform, a wavelet is stretched in time by adjusting its scale(s), , and is normalized to have unit energy. Hence, this transform has the advantage of low time resolution at low frequencies and high-time resolution at high frequencies. The WT of a time series (, ) with uniform time steps is defined as the convolution of with a scaled and normalized wavelet,35
Accordingly, the wavelet power is defined as .
The WT decomposes a signal into the time-frequency domain and its power spectrum describes the energy distribution of a signal in the dual time-frequency domain. But in some cases, we need to deal with both signals simultaneously and are interested in the correlation between them in the time-frequency domain. Hence, wavelet coherence was proposed. It reveals areas in the time-frequency domain and these areas are quantified with the coherence between them.35 The time period and frequency band covered by these areas indicate that the two signals are coherent or correlated in both the time period and the frequency band.
Given two time series, and , with WTs, and , the wavelet coherence is defined as , where * denotes complex conjugation. In the same manner, the cross wavelet power will be . Torrence and Grinsted defined the wavelet coherence of two time series as31,3637 The complex argument of describes the phase angle between and in the time-frequency domain. It is defined as
Wavelet coherence analysis provides the coherence between two signals, but the quantity is normalized between 0 and 1. This requires an additional metric that can assess the coherence confidence compared to that which occurs between two random signals, i.e., white Gaussian noise in our case.
In order to determine the confidence level of the wavelet spectrum, the first step is to select an appropriate background spectrum and then to compare the wavelet spectrum of a signal against the background one.35 Matteau-Pelletier studied noise in fNIRS and proved that it was a pink noise (a.k.a., noise), i.e., the power increased proportionally with a decrease in frequency.16 Kaulakys suggested that the noise can be analyzed by an autoregressive model.38 Referring to the theory in the Fourier analysis, the normalized Fourier power spectrum is given by6) and then the 95% confidence contour lines can be drawn by comparing the local wavelet power spectrum against . The contour can be used to identify the region where a significant power of the fNIRS channel data is located.
Activation Detection by Wavelet Coherence
fNIRS provides the oxy and deoxyhemoglobin concentrations, which provide functional information about certain underlying brain sensorimotor or cognitive tasks. Hence, their concentrations should reflect correlations between brain neurophysiology and behavioral paradigms. The correlations between these have been confirmed by both functional magnetic resonance imaging and fNIRS.39 In our case, the coherence, i.e., the time-locked correlation in the frequency domain, is used to evaluate the level of each channel correlated with the task paradigm.
Suppose that an experimental paradigm is expressed as , then a local hemodynamic response can be thought of as the output of a linear time-invariant system.40 In such a system, a hemodynamic response has a limited duration because of the brevity of the sensory or cognitive stimulus. Moreover, the response is independent of time, that is to say, the same response is replayed as long as the same stimulus is exerted, regardless of when it is implemented in the system. In addition, the responses to many successive inputs can be assumed to approach linearity. Hence, we can predict the hemodynamic response by the convolution of with an impulse response, ,41
As stated before, fNIRS measures oxy- and deoxyhemoglobin concentrations. The two chromophores reflect local consumption of oxygen in a region of the cortex since oxygen is necessary to sustain neuronal firing in response to sensory or cognitive stimuli. Oxygen consumption, in turn, is physiologically coupled to local hemodynamics. Accordingly, the fNIRS data are related to the hemodynamic response. However, the physiological noise that is mixed into the fNIRS data is likely to be periodic but not always synchronized with the predicted response. Moreover, the hemoglobin concentrations obtained by fNIRS and the predicted response are correlated in both the time and frequency domains. Therefore, wavelet coherence works well to identify the correlation in both domains.
Given channels that are used to collect the hemodynamic response of different regions of a subject’s brain, , , the WTs of the predicted response and the detected signals are denoted by and , respectively. Equation (3) can then be used to calculate the coherence between the fNIRS channel data and the predicted response. Since there are channels collecting hemodynamic data for a subject’s brain, we can obtain coherence distributions of channel data with the predicted response, , .
These distributions all covered the same time period and the specified frequency band [0.01, 5] Hz. The upper band limit was determined by the sampling rate of our fNIRS system. It sampled fNIRS data at 10 Hz, which was able to capture signals at a frequency lower than 5 Hz according to the sampling theorem. The lower limit was set following normal usages in fNIRS data postprocessing. Signals under 0.01 Hz in fNIRS did not contain meaningful physiological information but instead were composed of systematic drifts or noises. It has been agreed with in fNIRS community. So people would like to filter out the part by a high-pass filter with a cutoff frequency at 0.01 Hz. The hemodynamic information recorded by fNIRS reflected changes in concentrations of hemoglobin in the blood vessels of the brain, plus superficial tissue, which vary in their deoxy- and oxyhemoglobin concentrations in response to neuronal activity and various physiological processes. The frequency of the information we collected covered the band [0.01, 5] Hz. A mean power was calculated for each channel in the region of interest in the two domains. Then every distribution was represented by a mean power value, which indicates the level of the coherence of the corresponding channel with the predicted response. Finally, these values were sorted and activation in every channel was quantified by the mean power. These adjacent activated channels were linked together to represent the regions that were activated during the task.
Two experiments were conducted to demonstrate the performance of the activation detection method. In the first experiment, we adopted simulated data to verify that the activated channels were accurately detected. In the second experiment, data obtained while using visual stimuli were used to analyze the performance of the method in detecting the activation area.
The fNIRS measurements were collected using a continuous wave system (ETG-4000, Hitachi Medical Co., Japan). The system generated two wavelengths of near-infrared light (690 and 830 nm) and collected the hemoglobin concentration at a 10 Hz sampling rate. optode probe sets (consisting of seven detectors and eight light emitters with 3.0 cm of source-detector separation) were used to produce 22 measurement channels to allow for the measurement of the brain.
Experiment on Simulated Data
The experiment using simulated data with 22 channels evaluated the accuracy of the detection method. The data were measured during a resting state for 300 s. We recruited one person to participate in the experiment. The subject sat quietly on a chair while we recorded the fNIRS data at 22 channels. These data were used as the background signals. We then added a hemodynamic response into several channels and determined whether the method was able to detect these activated channels.
We introduced the hemodynamic response data based on having a task that started at 70 s with a trial that stayed “on” for 10 s and then was “off” for 10 s. There were three trials in the task at 70, 110, and 190 s, respectively. The task protocol is shown in Fig. 1(b). In order to simulate an activated channel, we predicted a response by convolving the canonical HRF with the set of stimuli using Eq. (7), as shown in Fig. 1(a). We contaminated the response with Gaussian white noise to reduce the signal-to-noise ratio to 20 dB. In addition, cardiac noise at 0.8 Hz, Mayer wave noise at 0.1 Hz,42 and low frequency oscillations at 0.01 Hz were added into the predicted response.43 Finally, the predicted response with physiological noise was added into the data measured in the resting state at channels 6, 17, and 19. These data were used as simulated data to test our detection method.
A scalogram is dedicated for wavelet analysis since it decomposes a signal in terms of time and scale. The scalogram is a visual method to display the decomposition with respect to the time (, axis) and scale (, axis). Its counterpart is the spectrogram, which is frequently used in Fourier transform to display the power spectrum along with the time. The scalogram of the channel 17 is shown in Fig. 2(a). In this figure, the upper left line plot presents the power spectral density of the channel; the lower right line plot is the data in the time domain; and the upper right color map is the scalogram of the WT. The color bar in Fig. 2(a) represents the confidence level. From the wavelet scalogram, we were able to determine the energy distribution of channel 17 with a simulated predicted response in the time-frequency domain. Three frequency bands clearly indicated the three noises involved: cardiac noise at 0.8 to 1.5 Hz, Mayer wave noise at 0.1 to 0.4 Hz, and low frequency oscillations at 0.01 Hz. The desired hemodynamic response was hidden in these noises. The scalogram presents regions where the confidence level was over 95%. This means that the power in these regions was significantly greater than Gaussian white noise. We used the region in the scalogram with a confidence level over 95% as masks because this confidence level guaranteed a mask that was the least affected by other physiological signals. Although the resulting mask may sometimes have been larger than necessary because of using this method, the result showed that it was still effective to locate critical coherence.
The expected response was added into the measured resting-state data in channels 6, 17, and 19 with amplitudes of , , and , respectively. The data from four channels—6, 7, 17, and 18—are shown in Fig. 3. Here, the data from channels 7 and 18 are shown for comparison with the data containing the desired responses in channels 6 and 17. As these plots reflect, we were not able to detect any differences between the data with the desired response and the data without the response.
Figure 4 shows the wavelet coherence results for the two channels, 17 and 18, with the predicted response, respectively. The response is the result of the convolution of the stimuli responses with the HRF. The wavelet coherence is employed to investigate the correlation between channel data and the predicted response. The color bar is scaled with respect to the result by Eq. (3). Arrows in the coherence show the phase lag between the measured signal and the predicted response. Arrows pointing to the right mean that there is no lag between these signals and those pointing to the left mean that there is a phase lag between them.
Comparing these two coherence results, we cannot see any significant coherence at channel 18 in the period, [70, 240] s. Neither the coherence nor the phase angle shows any significant correlation between the predicted response and this channel data. The coherence at channel 17 presents significant red areas surrounded by black contours. One of them is located at the frequency band, [0.04, 0.08] Hz and has arrows pointing to the right. By referring to the scale of the color bar, red areas indicate where significant coherence takes place. Point-right arrows illustrate that no phase lag exists between the two signals. Significant red coherence areas plus point-right arrows confirm that data in channel 17 are closely related with the predicted response. The remaining one red area is located in a majority at around [0.01, 0.03] Hz and with arrows pointing down.
We used the mask in Fig. 2(b) to delineate the region at which the desired response can be expected to appear. Because these noises were periodic and might be related to activation signals, the results of wavelet coherence of channel data with a predicted response might present a false coherence. We used the mask to indicate areas with the coherence of interest. Therefore, one of the purposes of using the mask was to unveil the desired coherence. Meanwhile, the mask was used to filter out coherence caused by physiological noises. These regions were set as the regions of interest to indicate where we should calculate the coherence. Then the coherence in every channel was selected and summed together. The normalized integral coherence at every channel is shown in Fig. 5(a). The normalized coherence of channels 6, 17 and 19 is significantly larger than those of the other channels. We have drawn an activation map using these normalized coherences in Fig. 5(b). These activated areas also emerged at channels 6, 17 and 19. The experiment was also performed by the GLM method. The result is shown in Fig. 5(c). In a comparison of the two results in Figs. 5(b) and 5(c), GLM was able to detect an apparent activated channel but failed to detect channels with a weak response. The method by wavelet coherence was able to detect these active channels.
Experiment on Real Data
An experiment on real data was performed to validate the detection method. All experimental procedures were approved by the Beijing Normal University Institutional Review Board. Research was carried out according to the principles of the Declaration of Helsinki, and the experiments were conducted with the permission and written consent of each participant.
In the experiment, we recruited 21 subjects and they conducted a visual stimulation task. The task was intended to detect activations in the visual cortex under stimuli of oriented gratings. It has been well explored that ocular dominance columns and orientation preference maps coexist and overlap across cortical tissue in the primary visual cortex (V1).44 The experiment was to verify whether the active regions detected by fNIRS data were consistent with the biophysiological mechanism. Figure 6(a) shows the experimental paradigm in a trial and Fig. 6(b) presents the probe positioned on the participant’s scalp to collect fNIRS data. In Fig. 6(b), red dots denote sources and blue dots denote detectors. Numbers indicate the channels in the configuration.
In the first 120 s, the subjects were told to stay quiet and the data were measured in the resting state. After this initial 120 s, the subjects were periodically exposed to a visual stimulus containing oriented cardinal or oblique gratings in 0, 45, and 90 deg in a random order as shown in Fig. 6(a). A block included visual stimuli for 20 s and then a resting state for 20 s. There were a total of five blocks. The predicted response to the stimuli paradigm was estimated by convolving the experimental paradigm with the canonical HRF and is shown in the lower subfigure in Fig. 7(a). The WT spectrum for the response in the time-frequency domain is shown in the upper right of Fig. 7(a). In order to illustrate the frequency bands of significant power, we show the spectrum in terms of the frequency in the upper left subfigure in Fig. 7(a). From the subfigure, we are able to recognize the main band at 0.025 Hz. The black contour delineates the significant region from the noise. The region was set as a mask for the following analysis, as shown in Fig. 7(b).
We applied the wavelet coherence to analyze the coherence between the predicted response and each channel data point. The coherences for two channel data are drawn in Fig. 8 and its subfigures (a) and (b), which correspond to channels 10 and 11, respectively. There is no any meaningful coherence at the resting state from the beginning to 120 s in the two coherence distributions. After starting to perform the task at 120 s, coherence primarily appeared at 0.08 and 0.025 Hz. At the frequency band of interest, around 0.025 Hz, the coherence in channel 10 presents a rather weak significant correlation, whereas the coherence at channel 11 showed significant correlation. In order to reduce the effect of noise at the other frequency bands, we used the mask shown in Fig. 7(b) to delineate regions at which the desired response should be expected to be detected. All the coherences in the region were summed and the integral coherence was used to evaluate the level of the activation.
Because there were 21 subjects involved in the experiment, the normalized coherence across subjects is shown in Fig. 9. Channels 11 and 12 present significantly larger coherence than the other channels. We mapped the normalized coherence of the channels into the probe as we measured the fNIRS data. The map is shown in Fig. 9(b). The probe was positioned on the occipital lobe as shown in Fig. 6(b). Referring to the subfigure, channels 11 and 12 were projected to the V1 visual cortex of subjects. Therefore, the activation area detected by our method is consistent with the biophysiological mechanism.
The present detection method makes use of wavelet coherence to explore an activated area which is based on the WT decomposing channel data into the time-frequency domain across multiple scales. The analysis tool effectively separates physiological noise from the signal of interest in channel data, given the experimental paradigm of recording the data. Nevertheless, in some special cases, physiological noises may be synchronized completely with functional activation under certain tasks so that activation by the method will not be detected. We will recommend recording physiological noises with dedicated tools and then regressing them from the measured signals.
In general cases, it is essentially reasonable to assume that signals originating from functional activity are not synchronized with these physiological noises, because these noises’ sources—lung, blood flow, heart, etc.—are independent from the cortex. These sources operate in different mechanisms from that of the neurohemodynamics. Although their signals are mixed in a complicated manner in fNIRS data, these signals still have to follow their own mechanisms. In addition, the cross wavelet correlates the data with the experimental time-locked paradigm and works well to locally evaluate the correlation and indicate the parts that are consistent with the paradigm. The most highly correlated portion reveals the level of activation in a channel. Two experiments on simulated and real data demonstrated that the coherence effectively performed the detection.
We appreciate the assistance that Rhoda and Ed Perozzi provided in proofreading this work. This work was partially supported by the National Key Basic Research and Development Program (973) (Grant No. 2011CB707800), the National Key Scientific Instrument and Equipment Development Project of China (2012YQ120046), the National High Technology Research and Development Program of China (863) (Grant No. 2012AA011600), and the Natural Science Foundation of China (Grant No. 81101082).
Xin Zhang is an associate professor at the Institute of Automation, Chinese Academy of Sciences. He received his BS and MS degrees in biomedical engineering from the Capital Medical University in 2002 and Tsinghua University in 2006, respectively, and his PhD degree in electric and electronic engineering from the University of Hong Kong in 2010. His current research interests include neurophotonics and neurohemodynamics. He is a member of SPIE.
Jian Yu is a graduate student at the University of Electronic Science and Technology of China. He is a visiting student at the Institute of Automation, Chinese Academy of Science since August 2013. His work focuses on the research and development of functional near-infrared spectroscopy instruments.
Ruirui Zhao is a graduate student at the University of Electronic Science and Technology of China and her major is biomedical engineering. She is a visiting student at the Institute of Automation, Chinese Academy of Sciences. She received her BS degrees in communication engineering from Henan Normal University in 2008.
Wenting Xu is an assistant engineer at the Institute of Automation, Chinese Academy of Sciences. He received his BS and MS degrees in biomedical engineering from the Taishan Medical University in 2011 and University of Electronic Science and Technology of China in 2014, respectively. His current work includes neuroimaging and neurohemodynamics.
Haijing Niu is a Project 985 researcher at the Beijing Normal University. She received her BS and MS degrees in physics optics from Yanbian University in 2002 and Tianjin University in 2006, respectively, and her PhD degree in information optics from Beijing Normal University in 2009. Her current research interests include neurophotonics and fNIRS-based brain connectome methodology. She is a member of SPIE, OSA, and IEEE.
Yujin Zhang is an assistant professor at the Institute of Automation, Chinese Academy of Sciences. She received her BS degree in communication engineering from Tianjin University in 2007, and her PhD degree in cognitive neuroscience from Beijing Normal University in 2013. Her current research interests include functional neuroimaging, near-infrared spectroscopy imaging, and neurohemodynamics.
Nianming Zuo is an associate professor at the Institute of Automation, Chinese Academy of Sciences (CASIA). He received his BS degree in mathematics from Shandong University in 2002 and PhD degree in medical imaging from Chinese Academy of Sciences in 2007, respectively. Before joining the Brainnetome Center, CASIA, he was a research scientist in Kodak Health Group (Shanghai, China). His current research interests include computational neuroscience and medical equipment development. He is a member of OHBM.
Tianzi Jiang is a professor of Brain Imaging and Cognitive Disorders, Institute Automation of Chinese Academy of Sciences (CASIA), and a professor of Queensland Brain Institute, University of Queensland. He is the director of Brainnetome Center of CASIA. His research interests include neuroimaging, brainnetome, imaging genetics. He is the author or coauthor of over 200 reviewed journal papers in these fields. He is an associate editor of IEEE Transactions on Medical Imaging and several other journals.