1.

Introduction

Functional near-infrared spectroscopy (fNIRS) is an emerging noninvasive brain imaging technique that uses light in the near-infrared (NIR) range to measure local changes in the cerebral concentration of oxygenated hemoglobin ($\mathrm{\Delta }[{\mathrm{HbO}}_2]$) and deoxygenated hemoglobin ($\mathrm{\Delta }[\mathrm{HbR}]$) associated with the underlying brain activities.^1–4 Compared to functional magnetic resonance imaging (fMRI), which is only sensitive to $\mathrm{\Delta }[\mathrm{HbR}]$, fNIRS is able to detect both $\mathrm{\Delta }[{\mathrm{HbO}}_2]$ and $\mathrm{\Delta }[\mathrm{HbR}]$, providing additional information related to brain activity. In addition, compared to fMRI, fNIRS is relatively compact, inexpensive, and is a less restraining imaging system, allowing for brain studies that are conducted in naturalistic settings. Because of these advantages, fNIRS has been widely used in a variety of neuroscience studies including spatiotemporal mapping of brain activities,^5–11 investigating functional connectivity of brain networks,^12–16 and brain computer interface applications.^17–21

The averaging operation is performed at different stages of analysis (e.g., across trials, blocks, subjects, and channels) in fNIRS brain imaging studies, with the objective of enhancing the signal strength associated with task-induced brain activities, and reducing noise and randomness. For example, in task-based fNIRS studies, the averaging operation is used at the early stages of analysis. Task-based experimental paradigms are categorized as “block design” and “event-related design.” In block design experiments, to increase the detection power for estimating the location of task-induced active regions, the general idea is to present multiple trials of the same type to the subject within each block, and repeat the experiment across multiple blocks. Blocks of different experimental conditions are often interleaved, and recorded fNIRS signals across blocks of similar conditions are averaged through the conventional point-by-point averaging technique.^{10,11,22–26} In event-related design studies, brain activities associated with individual trials are recorded, allowing for estimating the brain’s hemodynamic response related to the stimulus. This hemodynamic response can be obtained through averaging recorded activities in response to several discrete events of the same type. In both categories, to estimate the location of brain regions associated with the task, the stimulus or the event of interest, various statistical tests (e.g., student’s $t$-test) are used to evaluate the statistical significance of features (e.g., amplitude) of the averaged signals.²⁷

As the averaging operation is conducted at the early stages of the analysis,^{10,11,23,28–30} inaccuracies in the averaged signal could lead to type I (incorrectly detecting a region as active) or type II (incorrectly detecting a region as inactive) errors in the statistical analysis, resulting in inaccurate conclusions about the functionality of the brain. As stated before, to perform the averaging operation, typically a conventional point-by-point averaging technique is used in fNIRS studies. Previous work, however, has shown that there exists variability (e.g., latency differences) in the brain response to trials of the same type (e.g., trial-to-trial variability) due to for example differences or delays in neural responses, or individual’s performance.^31–34 Furthermore, in patient populations (e.g., patients with autistic spectrum disorders [ASD]), several studies have reported variable latencies in their responses to stimuli.^35,36 While preserving information related to variability would be important in identifying parameters related to behavioral variability (e.g., understanding neural mechanisms related to variability in response time),^31,34,37 in several neuroimaging studies (e.g., those interested in functional specificity) conclusions are made based on an “aggregate-then-analyze” approach. In these type of studies, using conventional point-by-point averaging approach, which does not take temporal variation of signals into account, can result in inaccuracies (e.g., blurring the peaks and valleys) in the averaged signal. Examples of three scenarios are conceptually visualized in Fig. 1. In the scenario shown in Fig. 1(a) (left column), no temporal variation exists in the alignment of the three signals to be averaged. In the scenario shown in Fig. 1(b) (middle column), a “linear” temporal variation exists among the three signals. In the scenario shown in Fig. 1(c) (right column), a “nonlinear” temporal variation exists among the three signals. The averaged signals obtained through a point-by-point averaging technique for each scenario are shown in the last row. It can be seen that, compared to Fig. 1(a), the presence of linear temporal variation for signals in Fig. 1(b) has resulted in attenuation of the amplitude of the averaged signal. The presence of nonlinear temporal alignment for signals in Fig. 1(c) has caused nonlinear distortion in the averaged signal obtained through conventional point-by-point technique. In reality, such inaccuracies in the averaged signal can lead to misunderstandings about the brain function.³⁸

Fig. 1

Visual illustration of using point-by-point conventional averaging technique for three scenarios: (a) when three signals are temporally aligned, (b) when there exists linear temporal variation in the alignment of three signals, and (c) when there exists nonlinear temporal variation in the alignment of three signals.

To address the problem of temporal variation across fNIRS signals, in this paper, a dynamic time warping (DTW)-based averaging technique is presented for fNIRS-recorded time series. DTW algorithm was originally introduced in the field of speech processing³⁹ and has been widely used for measuring the similarity of two time series in various fields of research, such as biometric, data mining, gene expression analysis, human motion recognition, and EEG signal analysis.^32,40–44 Various forms of DTW-based averaging approaches have also been proposed for applications, such as EEG³² and satellite image time series.⁴⁵

To extensively examine the performance of the proposed DTW-based averaging technique in fNIRS-based detection studies, experiments and simulation are performed. For the experimental part, both block design and event-related design experiments are considered. For the simulation study, a dataset is generated based on the standard model of hemodynamic response function (HRF), and receiver operating characteristic (ROC) curves are used to compare the performance of conventional- and DTW-based averaging techniques. It is shown from both experimental and simulation studies that the DTW-based averaging technique outperforms the conventional averaging technique in terms of accurately estimating the location of task-induced brain regions. While for the proof of concept, in this paper, we focus on the problem of localizing task-induced active brain regions, the presented averaging technique can be used in other averaging stages of fNIRS neuroimaging studies (e.g., computing grand averages).

The remainder of this paper is organized as follows: in Sec. 2, the DTW-based averaging technique is described. In Sec. 3, the experimental and simulation studies are presented and results are discussed. Finally, the paper is concluded in Sec. 4.

2.

DTW-Based Averaging

Let ${\mathbf{b}}_{\mathbf{k}}=[b_k(1),\cdots ,b_k(N)]$ represent the $k$’th time series of a group of $K$ time series, each with $N$ time points. In fNIRS neuroimaging experiments, each of these time series corresponds to the signal from a channel associated with a given block (in a block design experiment, $K$ being total number of blocks) or associated with a given trial (in an event-related design experiment, $K$ being total number of events). The aim is to obtain the averaged representation of these $K$ time series using DTW.⁴⁵

To obtain the DTW-based averaged representation of $K$ signals, first, the best alignment between each signal and a “reference” signal $\mathbf{c}=[c(1),\cdots ,c(N)]$ is found. The reference signal could be, for example, the conventional averaged signal of all $K$ time series. To optimally align signal ${\mathbf{b}}_{\mathbf{k}}$ and the reference signal $\mathbf{c}$, a cost matrix ${\mathbf{D}}_{\mathbf{k}}$ needs to be determined. ${\mathbf{D}}_{\mathbf{k}}$ is a $N\times N$ matrix in which its elements are obtained through a cost function representing the discrepancy between the $i$’th and $j$’th samples ($i=1,\cdots N$ and $j=1,\cdots N$) of signals $\mathbf{c}$ and ${\mathbf{b}}_{\mathbf{k}}$. Measures of the Euclidean distance, or the square of difference between normalized samples, can be used as the cost function.³²

Next, from the cost matrix, an optimal alignment path, ${\mathbf{W}}_{\mathbf{k}}^{\mathbf{opt}}={[w_1,\cdots ,w_L]}^{\mathrm{T}}$, $N\le L\le 2N-1$, where $w_l=[i(l),j(l)]$, $1\le l\le L$, $1\le i(l)$, $j(l)\le N$, must be determined so that the overall similarity between the two signals is maximized. The optimal alignment path ${\mathbf{W}}_{\mathbf{k}}^{\mathbf{opt}}$ shows how the mapping between the indices of the two time series $\mathbf{c}$ and ${\mathbf{b}}_k$ must be made to achieve the best alignment. For example, if ${\mathbf{W}}_{\mathbf{k}}^{\mathbf{opt}}={[(1,1),(2,3),(2,4),\dots ,(N,N)]}^{\mathrm{T}}$, then the sample $c(1)$ is aligned with $b_k(1)$, $c(2)$ is aligned with $b_k(3)$ and $b_k(4)$, etc. To obtain ${\mathbf{W}}_{\mathbf{k}}^{\mathbf{opt}}$ the solution to the following optimization problem³²

Once ${\mathbf{W}}_{\mathbf{k}}^{\mathbf{opt}}$ is obtained, a new $N$-points time series, ${\mathbf{b}}_{\mathbf{k}(\mathbf{\text{aligned}})}=[b_{k(\text{aligned})}(1),\cdots ,b_{k(\text{aligned})}(N)]$, is formed as follows

As an example, if

This procedure will be performed for all $K$ signals. Once all signals are aligned with the reference signal $\mathbf{c}$ and their aligned representations are determined, the DTW-based averaged signal is obtained as

Figure 2(a) shows an example of aligning signal ${\mathbf{b}}_1$ with a reference signal $\mathbf{c}$. The color-coded cost matrix and the obtained optimal warping path (in blue) are shown in Fig. 2(b). It can be observed that when a sample in ${\mathbf{b}}_1$ is aligned with several samples in $\mathbf{c}$, the warping path has a vertical direction, whereas when a sample in $\mathbf{c}$ is aligned with several points in ${\mathbf{b}}_1$, the warping path follows a horizontal direction. Note that the optimal warping path is mostly along the antidiagonal elements of the cost matrix, illustrating that the two signals experience temporally variable latencies. Table 1 summarizes the steps involved in the DTW-based averaging procedure using a reference signal.

Fig. 2

Illustration of the alignment process for an exemplary signal ${\mathbf{b}}_{\mathbf{1}}$ and a reference signal $\mathbf{c}$ through DTW: (a) time series ${\mathbf{b}}_{\mathbf{1}}$ and $\mathbf{c}$ where their aligned points are connected, and (b) color-coded cost matrix, representing the distance between the two time series. The warping path is shown in blue. Plots are created using R programming package.⁴⁶

Table 1

Summary of steps involved in DTW-based averaging procedure using a reference signal.

It should be noted that the DTW-based averaging technique can be realized in various forms.⁴⁵ For example, in Ref. 32, instead of using a reference signal, the average is obtained sequentially in a pairwise manner. As such, errors at early stages of computation could propagate throughout the averaging process, resulting in loss of information.⁴⁵ To address the problems associated with pairwise averaging, a global averaging strategy was introduced in Ref. 45, in which the averaged signal is obtained considering all signals and is updated through an iterative process. The method presented here also computes the average by considering all signals and therefore, will not be sensitive to ordering effects.

3.

Evaluation

The performance of the proposed averaging framework for analyzing fNIRS-recorded time series is evaluated through both experimental and simulated data. On the experimental end, two experiments, one block design and one event-related design, are considered. The recorded time series are analyzed through both conventional and DTW-based averaging techniques. As it will be shown, compared to conventional point-by-point averaging, when DTW-averaging is used, an increase in the detection power for the block-design task and a decrease in the false positive rate for the event-related task are observed. For the simulation study, data using a widely-used equation of the hemodynamic response is generated, and ROC curves for the conventional and the DTW-based averaging approaches are obtained and compared. Furthermore, the impact of the choice of reference signal on the performance of DTW-based averaging technique is investigated.

3.1.

Experimental Studies

Two experiments, one block design (experiment I), and one event-related design (experiment II) are performed. We first describe the experimental setup for each case, and then present and discuss the results.

3.1.1.

Experimental Setup

Experiment I: block design paradigm

Five healthy right-handed volunteers (one female) participated in experiment I. All volunteers gave their informed consent approved by the Rutgers University Institutional Review Board prior to the experiment.³⁸ The paradigm for this experiment was the 2-back task [see Fig. 3(a)], which has been widely used in determining brain activities related to working memory.^25,26,47 Three blocks were presented to participants. Each block lasted for 32 s and consisted of 16 letters (computerized in a pseudorandom order with four being target stimuli) with intertrial interval (ITI) of 2 s. The order of presentation of the target stimulus was different across blocks. Participants were asked to respond by pressing the right button on the mouse when the presented letter matched the letter shown two stimuli back.

Fig. 3

(a) Experimental paradigm, and (b) optode setup (red: light emitter, green: detector, and blue: fNIRS channel), for the 2-back task (experiment I).

It is worth noting that the block design experiments are performed with the goal of continuously engaging the brain in the cognitive process of interest over the duration of the block, to increase the detection power for localizing the related activities in the brain. Here, the cognitive process of interest is “working memory.” The process of memorizing and recovering letters is carried out on a continuous basis over the duration of the block for each letter (target or nontarget). For letters that satisfy the 2-back condition, other additional processes and functions such as pressing the button (activating the somatosensory region) will become involved. Through the averaging process across blocks, the signal-to-noise ratio (SNR) related to these additional functions would be small because of the variations in the timing of the target stimulus from one block to other. But information related to “working memory,” which is the main objective of the $n$-back task, should be preserved over the duration of the block.

fNIRS data were collected using a 52-channel [Hitachi ETG-4000, 17 sources (lasers, 695 and 830 nm) and 16 detectors] at a sampling rate of 10 Hz. The source–detector separation was 3 cm. The headband optode holder was placed on the forehead of volunteers to cover the prefrontal cortex, as shown in Fig. 3(b).

Experiment II: event-related paradigm

Six right-handed healthy volunteers (all males) participated after giving their informed consent. An event-related modified visual oddball task^44,48 consisting of three graphical stimuli (plus, square, and circle), presented in random order, was used as the paradigm [see Fig. 4(a)]. A total of 220 stimuli [30 target stimuli (“plus”), and 190 nontarget stimuli] were presented. Each stimulus was presented for 50 ms, with ITI of 10 to 12 s. To minimize the periodic systemic effects,⁴⁹ the ITI was randomized to prevent the subjects from predicting the presentation time. Participants were asked to press the left button of the mouse once the target stimulus was shown on the screen.

Fig. 4

(a) Experimental paradigm, and (b) optode setup (circle: fNIRS channels in prefrontal cortex, triangle: fNIRS channels in visual cortex), for the modified visual oddball task (experiment II). Locations of optodes are visualized using MATLAB, BrainNet Viewer.⁵⁰

fNIRS data were measured by an NIRx System [NIRx NIRScout, 16 sources (LEDs, 760 and 830 nm) and 16 detectors] at the sampling rate of 10.42 Hz. Measured signals from source–detector pairs with a separation distance of 3 cm were considered, resulting in a total of 38 channels. Optodes were placed over the prefrontal and visual regions of the cortex, as shown in Fig. 4(b).

3.1.2.

Preprocessing

Recorded signals were visually inspected. Bad channels (for example due to loose contact to skin) were excluded from further analysis. In experiment II, the response to the first trial for the majority of participants included traces of the subject’s movement, and so it was excluded from the analysis, for all subjects. Signals were then segmented. For experiment I, signals were segmented by blocks. For each block, the segmentation time window began at 5 s prior to the onset of the first stimulus in the block and ended 15 s after the end of the block. For experiment II, signals associated with “target” trials and “nontarget” trials were segmented separately. Trials associated with missed response (i.e., target trials without subject’s response) were excluded from further analysis. The segmentation window began at 1.9 s prior to the onset of the stimulus and had a duration of 13.9 s.

Next, signals related to $\mathrm{\Delta }[{\mathrm{HbO}}_2]$ and $\mathrm{\Delta }[\mathrm{HbR}]$ with respect to the baseline were extracted using the modified Beer–Lambert law,⁵¹ based on the following equations

Eq. (7)

$$\begin{gathered} \ln \left(\frac{I_{{\lambda }_1}}{I_{\text{baseline},{\lambda }_1}}\right)=-({\epsilon }_{{\mathrm{HbO}}_2,{\lambda }_1}·\mathrm{\Delta }[{\mathrm{HbO}}_2]+{\epsilon }_{\mathrm{HbR},{\lambda }_1}·\mathrm{\Delta }[\mathrm{HbR}])·{\mathrm{DPF}}_{{\lambda }_1}·x; \\ \ln \left(\frac{I_{{\lambda }_2}}{I_{\text{baseline},{\lambda }_2}}\right)=-({\epsilon }_{{\mathrm{HbO}}_2,{\lambda }_2}·\mathrm{\Delta }[{\mathrm{HbO}}_2]+{\epsilon }_{\mathrm{HbR},{\lambda }_2}·\mathrm{\Delta }[\mathrm{HbR}])·DP{F}_{{\lambda }_2}·x. \end{gathered}$$

In Eq. (7), $I_{{\lambda }_i}$ and $I_{\text{baseline},{\lambda }_i}$ ($i=1$, 2) are the optical intensities measured at the detector location at wavelength ${\lambda }_i$, during task and during prestimuli baseline period, respectively, $x$ is the distance between the light source and the light detector, $DP{F}_{{\lambda }_i}$ is the differential pathlength factor, and ${\epsilon }_{{\mathrm{HbO}}_2,{\lambda }_i}$ and ${\epsilon }_{\mathrm{HbR},{\lambda }_i}$ are the extinction coefficient of ${\mathrm{HbO}}_2$ and HbR at wavelength ${\lambda }_i$, respectively.

Band-pass filters (0.01 to 0.2 Hz) were then implemented to remove artifacts and low-frequency drift. Furthermore, following a procedure in Ref. 52, signals showing sudden fast changes were detected, and excluded from further analysis.

3.1.3.

Results—Experiment I

For each subject and each channel, averaged signals, using both the conventional-based and the DTW-based averaging techniques, were first obtained. Next, for each averaged signal, the time point at which the signal reaches its maximum value following the onset of stimuli was identified. A temporal window with a duration of 2.1 s⁵³ around this time point was considered, and the mean of the averaged signal within this window was calculated. This number, to which we refer to as activation index (AI), was used in subsequent statistical tests to determine if the region associated with the channel was active in response to the external stimuli.

To determine active regions, one-sample $t$-test (with $p<0.05$) was performed on AIs (obtained through each averaging technique), with the null hypothesis being the region that is not active.⁵³ Figure 5(a) shows the result of the statistical test. As can be seen five channels (channels 1, 4, 9, 14, and 19) were identified as being significantly active in response to the task when DTW-based averaging technique is used for obtaining the averaged signal. Table 2 summarizes the result of the statistical test. These channels were not detected as active channels when conventional-based averaging technique is used. The activation pattern obtained through the DTW-based averaging technique appears to be more consistent with the results reported in previous fMRI studies where activation in bilateral prefrontal cortex was observed in response to working memory tasks.^54,55

Fig. 5

(a) Statistical activation map for 2-back task based on the averaged signal obtained from the conventional-based averaging technique and the DTW-based averaging technique. Blue-colored circles indicate channels with a significant increase in $\mathrm{\Delta }[{\mathrm{HbO}}_2]$ identified based on averaged signal obtained through both averaging techniques. Red-colored circles indicate channels with a significant increase in $\mathrm{\Delta }[{\mathrm{HbO}}_2]$ that were only detected from the averaged signal obtained through the DTW-based averaging. The statistical significance level is $p<0.05$. (b) Statistical map illustrating the difference between the DTW-based and conventional-based averaged signals.

Table 2

Mean and standard deviation (SD) of AIs across subjects (units in μM) as well as results of t-test for channels showing statistically more significant activation when DTW-based averaging technique is used as compared to when conventional point-by-point averaging technique is used (df=4).

Channel	Conventional averaging technique			DTW-based averaging technique
	Mean (SD)	p-value	t-value	Mean (SD)	p-value	t-value
1	0.65 (0.59)	0.068	2.48	0.84 (0.67)	0.048	2.82
4	0.45 (0.40)	0.068	2.49	0.61 (0.41)	0.035	3.14
9	0.50 (0.48)	0.081	2.32	0.73 (0.49)	0.029	3.34
14	0.44 (0.42)	0.082	2.31	0.61 (0.49)	0.048	2.83
19	0.98 (0.95)	0.082	2.31	1.30 (0.97)	0.040	3.00

While Fig. 5(a) shows the statistical activation map by employing a fixed threshold ($p<0.05$), it would be informative to also investigate the statistical significance of the difference of the outcomes of the two averaging approaches. To achieve this, for each channel and each subject, we subtracted the conventional-based averaged signals from their DTW-based averaged counterparts and conducted statistical test on the difference signals, with the null hypothesis that there are no significant differences. Figure 5(b) presents the $t$-map obtained from the one-sample right-tailed student’s $t$-test. It is shown that for almost all regions, the obtained $t$-values are significantly larger than zero, indicating that there are statistically significant differences between the averaged signals obtained from the two techniques. Next, we used the metric of contrast-to-noise ratio (CNR)^56,57 to quantify the SNR for the averaged signals obtained from each averaging technique. We considered both $\mathrm{\Delta }[\mathrm{HbR}]$ and $\mathrm{\Delta }[{\mathrm{HbO}}_2]$ signals. Only channels that were identified as active through both averaging techniques [shown in blue in Fig. 5(a)] were considered in this analysis. Denoting $\text{mean}(\mathbf{DUR})$ and $\text{var}(\mathbf{DUR})$ as the mean and variance of the signal amplitude during 5 to 15 s after the onset of the first stimulus of the block, and $\text{mean}(\mathbf{ITI})$ and $\text{var}(\mathbf{ITI})$ as the mean and variance of the signal amplitude corresponding to 10 to 15 s after the presentation of the last stimulus of the block, the CNR is computed as

Eq. (8)

$$\mathrm{CNR}=\frac{|\text{mean}(\mathbf{DUR})-\text{mean}(\mathbf{ITI})|}{\sqrt{\mathrm{var}(\mathbf{DUR})+\mathrm{var}(\mathbf{ITI})}}.$$

The results are shown in Fig. 6. It is shown that for $\mathrm{\Delta }[{\mathrm{HbO}}_2]$, mean CNR across subjects is significantly higher when DTW-based averaging technique is used, as compared to when conventional averaging technique is used ($p<0.05$, $df=15$). Since $\mathrm{\Delta }[\mathrm{HbR}]$ signals are generally weaker than $\mathrm{\Delta }[{\mathrm{HbO}}_2]$ signals, we observe a less significant difference between the CNR in $\mathrm{\Delta }[\mathrm{HbR}]$ signals obtained from the two approaches, though the mean CNR obtained through the DTW-based averaging technique is still higher than that obtained through the conventional-based approach. These results show that the DTW-based averaging technique has increased the SNR in the averaged signals and can enhance the detection power in studies that aim to identify active brain regions associated with external stimuli.

Fig. 6

Comparison of mean CNR values in averaged $\mathrm{\Delta }[{\mathrm{HbO}}_2]$ and $\mathrm{\Delta }[\mathrm{HbR}]$ signals, for experiment I, obtained through conventional- and DTW-based averaging techniques, respectively.

3.1.4.

Results—Experiment II

Experiment II is an event-related experiment. The histogram of response time for the “target” stimuli for each subject is shown in Fig. 7. Variability in response time across trials is observed for all individuals. This variability can also be an indication of the presence of temporal variation in the timing of the occurrence of task-induced events in the recorded fNIRS signals, and if not addressed, would affect the accuracy of the averaged signal.

Fig. 7

Histogram of the response time for “target” stimuli for each subject in experiment II.

For each subject, $\mathrm{\Delta }[{\mathrm{HbO}}_2]$ signals from each channel were averaged for each condition (target and nontarget), using the conventional and DTW-based averaging techniques, separately. Figure 8 shows examples of $\mathrm{\Delta }[{\mathrm{HbO}}_2]$ signals and their averaged signals obtained through both averaging techniques, under the “target” condition (shown in left) and “nontarget” condition (shown in right) for channel 1, for a representative subject. It is observed that in both cases, the peak amplitude of the conventional-based averaged signal is lower than that of the DTW-based averaged signal. Furthermore, comparing the two conditions, the peak amplitude of the conventional-based averaged signal for the nontarget condition is lower than the peak amplitude of the conventional-based averaged signal for the target condition, which may indicate that the brain region under channel 1, has been active for the target response (and hence sensitive to the oddball effect). However, caution should be taken when making an inference from this observation to avoid false alarm. Indeed, when the averaging process is conducted using DTW-based averaging technique, the difference in the peak amplitude of the two conditions is minimized, and such conclusion cannot be made.

Fig. 8

Exemplary recorded $\mathrm{\Delta }[{\mathrm{HbO}}_2]$ signals from one channel (channel 1) under “target” condition (shown in left) and “nontarget” condition (shown in right) in experiment II. Each trace represents a signal associated with a trial. The averaged signals obtained through conventional (shown in blue) and DTW-based (shown in red) averaging techniques are also shown. The vertical bar represents the timing of the onset of the stimulus.

Next, for each subject and each channel, the AIs under the target and nontarget conditions, denoted as ${\mathrm{AI}}_{ta}$ and ${\mathrm{AI}}_{nt}$, respectively, were estimated using the window length of 21 samples (2.02 s). The difference in AIs across two conditions, defined as ${\mathrm{AI}}_{\mathrm{diff}}={\mathrm{AI}}_{ta}-{\mathrm{AI}}_{nt}$, was then computed. For each channel, ${\mathrm{AI}}_{\mathrm{diff}}$ obtained from all subjects, were pooled and tested using one-sample right-tailed student’s $t$-test, to determine whether the channel is active in response to the oddball effect, with the null hypothesis of $\mathrm{mean}({\mathrm{AI}}_{\mathrm{diff}})=0$. With a threshold of $p<0.05$, the identified active channels, using the conventional-based averaging technique and the DTW-based averaging technique are shown in Figs. 9(a) and 9(b), respectively. One can observe that, compared to the patterns obtained through the conventional averaging, the active regions identified through the DTW-based averaging are located mostly in the prefrontal and occipital cortices, which appear to be consistent with previous fMRI studies.⁵⁸

Fig. 9

Location of channels (shown in red) that are specifically sensitive to the target condition compared to the nontarget condition (sensitive to the oddball effect) in experiment II, using the conventional-based averaging [shown in (a)] and the DTW-based averaging [shown in (b)], respectively.

3.2.

Simulation Studies

Given that in simulations the “ground truth” in terms of the location of brain activation is known, simulations are performed to compare the performance of the conventional averaging and the DTW-based averaging techniques. Furthermore, we also investigate the impact of choice of reference signal on the performance of DTW-based averaging procedure.

3.2.1.

Simulation platform

We considered a scenario similar to experiment II where there are 50 fNIRS channels, and the task is an oddball task, with 20 target trials, and 150 nontarget trials. The designed “ground truth” was that 10 channels are sensitive to the oddball effect ($\text{target}>\text{nontarget}$).

Simulated signals were generated based on the equation that is typically used to model the HRF⁵⁹

Eq. (9)

$$\mathrm{HRF}(t)=A_1{\left(\frac{t-d}{{\tau }_1}\right)}^{{\delta }_1}{e}^{[-({\delta }_1/{\tau }_1)(t-{\tau }_1)]}-A_2{\left(\frac{t-d}{{\tau }_2}\right)}^{{\delta }_2}{e}^{[-({\delta }_2/{\tau }_2)(t-{\tau }_2)]},$$

where parameters $A_1$, $A_2$, ${\tau }_1$, and ${\tau }_2$ determine the amplitude of the peak and the undershoot, $d$ represents the time delay, and ${\delta }_1$ and ${\delta }_2$ form the general shape of the peak and the undershoot. Through fitting this model to a typical measured oxygenation response from fNIRS recording, we used $A_1=1$, $A_2=0.4$, ${\delta }_1=10$, and ${\delta }_2=20$ for signals under target and nontarget conditions. The amplitude of the signals for the target condition was set to be 3% larger than that of the non-target condition.

To simulate the nonlinear variability in the latency, parameters ${\tau }_1$, ${\tau }_2$, and $d$ were designated as normally distributed random variables. Their variance was set to 20, 25, and 10, respectively. Using this model, 20 target trials and 150 nontarget trials were created 40 times (corresponding to 40 participants). Additive white Gaussian noise was added to each of the simulated signals such that SNR equals to 10 dB.

3.2.2.

Performance Comparison of DTW-Based Averaging with Conventional Averaging

The analysis procedure for the simulated dataset followed the one we used for experiment II. First, the signals associated with each channel, each condition, and each subject were averaged separately using the conventional and DTW-based averaging techniques. Next, ${\mathrm{AI}}_{\mathrm{diff}}$ values were computed for each channel and each subject. One-sample right-tailed student’s $t$-tests were conducted across subjects for each channel, with the null hypothesis being $\text{mean}({\mathrm{AI}}_{\mathrm{diff}})=0$. As the “ground truth” is known, we expected that for the 10 “true” active channels, the estimated ${\mathrm{AI}}_{\mathrm{diff}}$ be significantly larger than zero.

The ROC curves were computed for each of the averaging techniques to evaluate their performance. Here, the ROC curve illustrates the fraction of detected active channels and the associated false positive rate, when the threshold (given as $t$-values) varies. Figure 10(a) shows the resultant ROC curves for the conventional (blue dashed line) and DTW-based (solid red line) averaging techniques using the simulation dataset. It is clearly shown that the DTW-based averaging technique outperforms the conventional one.

Fig. 10

ROC curves for the conventional and the DTW-based averaging techniques. For the right figure, a standard HRF is used as the reference signal for the DTW-based averaging technique.

Furthermore, to examine the effects of the choice of the reference signal on the outcome of the DTW-based averaging technique, we repeated the simulation study where instead of the point-by-point averaged signal, we used the standard hemodynamic response as the reference signal. The result is shown in Fig. 10(b). It illustrates that the performance of the DTW-based averaging, with the standard HRF as reference, is still better than that of the conventional-based averaging. This result shows that regardless of the choice of the reference signal, the DTW-based averaging technique outperforms the conventional-based averaging technique in fNIRS-based detection studies.

4.

Discussions and Conclusion

The averaging operation is performed at different stages of a wide range of fNIRS brain imaging studies (e.g., across trials, blocks, and subjects), with the objective of enhancing the signal strength associated with task-induced activities, and reducing noise and randomness. The averaging approach commonly used for fNIRS signals is the point-by-point averaging technique. As shown in this paper, due to the existence of variable latencies across fNIRS signals, the use of conventional point-by-point averaging technique can lead to inaccuracies in the averaged signal, and, consequently, incorrect conclusions about the functionality of the brain. To address this problem, a DTW-based averaging technique for fNIRS signals was presented. The technique optimally aligns each fNIRS signal with a reference signal such that their similarity is maximized. Once all fNIRS signals are optimally aligned, the averaged signal is computed.

To compare the performance of the DTW-based and conventional point-by-point averaging techniques for fNIRS signals, we focused on the problem of localizing task-induced active regions in the brain. Results from both block design and event-related design experiments showed that the location of task-induced active regions estimated based on the DTW-based averaged signals are better aligned with the results reported from prior fMRI studies. Furthermore, since in the presented algorithm fNIRS signals are individually aligned with a reference signal, we also investigated the question of whether choosing a different signal as the reference would impact the performance of the algorithm. With a reference signal modeled based on a standard HRF, the DTW-based averaging technique continued to show improved performance compared to the conventional point-by-point averaging technique, in localizing active brain regions. The results also show that the performance of the DTW-based averaging can be further improved with a proper choice of reference signal. In our simulation setup, all generated fNIRS signals were created based on the HRF model (with randomized parameters and added noise), hence the choice of the standard HRF as the reference signal resulted in an improved performance, compare to the scenario where point-by-point averaged signal was chosen as the reference signal.

Improvement in the accuracy of the average of fNIRS signals is expected to introduce a significant impact in various fNIRS-based neuroscience and clinical research studies. For example, in multimodal EEG-fNIRS experiments,^23,60,61 the location of active regions estimated from fNIRS signals, has been used as constraint for the EEG source localization problem.²³ Errors in localizing active brain regions due to inaccuracies in the average signal can negatively influence the outcome of the EEG source localization problem, and therefore, the method presented here can be employed to avoid such errors. While for the proof of concept, in this paper, we focused on the problem of localizing task-induced active brain regions, DTW-based averaging framework can be employed in other steps of the analysis of fNIRS signals, to avoid loss of information. For example, several studies report grand averages of the hemodynamic response across subjects.^8,62–66 We suggest DTW-based averaging be used in these averaging steps, instead of the commonly used point-by-point technique, since variations in latencies for signals recorded from different individuals are inevitable. This issue is of particular importance when the study focuses on patient population, as several investigations have confirmed the existence of variable latencies in responses of patients (e.g., those with ASD) to stimuli.^35,36

In conclusion, our results suggest that the conventional point-by-point averaging technique, commonly used in fNIRS brain imaging studies, can result in inaccurate conclusions about the brain function, and therefore, use of advanced averaging techniques such as DTW-based averaging that consider the temporal variations in the alignment of recorded fNIRS signals is highly recommended. Simulation results showed that the performance of the DTW-based averaging technique can be further improved based on the choice of reference signal. Future work involves investigating this issue in depth and developing quantitative statistical measures and iterative procedures to obtain an optimum choice for a reference signal that maximizes the accuracy of the averaged signal. We also plan to implement a MATLAB toolbox with capabilities of computing the DTW-based averaged of fNIRS signals and make it available to the broader fNIRS neuroimaging scientific community.