2.1.

Wavelet Transform

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 $\mathrm{\Delta }\omega$ and the duration in the time domain $\mathrm{\Delta }t$ does not exist, a function will fail to achieve the minimum of the uncertainty product, $\mathrm{\Delta }\omega ·\mathrm{\Delta }t$, 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.³³ 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.³⁴ 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), $\mu =st$, 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 ($x_n$, $n=1,\dots ,N$) with uniform time steps is defined as the convolution of $x_n$ with a scaled and normalized wavelet,³⁵

Accordingly, the wavelet power is defined as ${|{\mathcal{W}}_{x_n}(n,s)|}^2$.

2.2.

Wavelet Coherence

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.³⁵ 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, $x_n$ and $y_n$, with WTs, ${\mathcal{W}}_{x_n}(s)$ and ${\mathcal{W}}_{y_n}(s)$, the wavelet coherence is defined as ${\mathcal{W}}_{x_n,y_n}={\mathcal{W}}_{x_n}{\mathcal{W}}_{y_n}^{*}$, where * denotes complex conjugation. In the same manner, the cross wavelet power will be ${|{\mathcal{W}}_{x_n,y_n}|}^2$. Torrence and Grinsted defined the wavelet coherence of two time series as^31,36

2.3.

Confidence Level

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.³⁵ Matteau-Pelletier studied noise in fNIRS and proved that it was a pink noise (a.k.a., $1/f$ noise), i.e., the power increased proportionally with a decrease in frequency.¹⁶ Kaulakys suggested that the noise can be analyzed by an autoregressive model.³⁸ Referring to the theory in the Fourier analysis, the normalized Fourier power spectrum is given by

2.4.

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.³⁹ 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 $p(n)$, then a local hemodynamic response can be thought of as the output of a linear time-invariant system.⁴⁰ 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 $p(n)$ with an impulse response, $h(n)$,⁴¹

This impulse response function in Eq. (7) consists of a linear combination of two Gamma functions.⁴⁰ Here, $\xi (n)$ stands for Gaussian noise.

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 $J$ channels that are used to collect the hemodynamic response of different regions of a subject’s brain, $x_j(n)$, $j=1,\dots ,J$, the WTs of the predicted response and the detected signals are denoted by ${\mathcal{W}}_{y_n}(n,s)$ and ${\mathcal{W}}_{x_{jn}}(n,s)$, respectively. Equation (3) can then be used to calculate the coherence between the fNIRS channel data and the predicted response. Since there are $J$ channels collecting hemodynamic data for a subject’s brain, we can obtain $J$ coherence distributions of channel data with the predicted response, $R_j^2(n,s)$, $j=1,\dots ,J$.

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 $J$ 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.

3.1.

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,⁴² and low frequency oscillations at 0.01 Hz were added into the predicted response.⁴³ 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.

Fig. 1

Predicted response in the simulation experiment. (a) is the HRF function and (b) is the task protocol.

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 ($t$, $x$ axis) and scale ($s$, $y$ 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.

Fig. 2

Wavelet transform of the channel 17 in (a) and the mask shown in (b). The upper left line plot in (a) presents the power spectral density of the response; the lower right line plot is the channel 17 data in the time domain; and the upper right color map is the scalogram of the wavelet transform of the channel.

The expected response was added into the measured resting-state data in channels 6, 17, and 19 with amplitudes of $2\times {10}^{-6}$, $8\times {10}^{-6}$, and $3\times {10}^{-6}$, 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.

Fig. 3

Data at channels 6, 7, 17, and 18 shown in (a), (b), (c), and (d), respectively.

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 $\pi$ phase lag between them.

Fig. 4

Wavelet coherence shown in (a) and (b) between the ideal response and the data from channels 17 and 18, respectively.

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.

Fig. 5

Normalized coherence of the channels (a) and the activation map in (b) and (c) by the method and GLM method, respectively.

3.2.

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).⁴⁴ 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.

Fig. 6

Experimental paradigm in a trial (a) and the probe positioned on the participant’s scalp (b).

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).

Fig. 7

Wavelet scalogram of an ideal response in (a) and the mask with significant power in (b). The wavelet transform spectrum in the time-frequency domain is shown in the upper right subfigure of (a). The scalogram in terms of the frequency is shown in the upper left subfigure. The bottom graph shows the predicted response.

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.

Fig. 8

Wavelet coherence between the predicted response with the data at channels 9 and 10 in (a) and (b), respectively.

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.

Fig. 9

Normalized coherence of the channels and the activation map.