Real-Time Reduction of Task-Related Scalp-Hemodynamics Artifact in Functional Near-Infrared Spectroscopy with Sliding-Window Analysis

Functional near-infrared spectroscopy (fNIRS) is an effective non-invasive neuroimaging technique for measuring hemoglobin concentration in the cerebral cortex. Owing to the nature of fNIRS measurement principles, measured signals can be contaminated with task-related scalp blood flow (SBF), which is distributed over the whole head and masks true brain activity. Aiming for fNIRS-based real-time application, we proposed a real-time task-related SBF artifact reduction method. Using a principal component analysis, we estimated a global temporal pattern of SBF from few short-channels, then we applied a general linear model for removing it from long-channels that were possibly contaminated by SBF. Sliding-window analysis was applied for both signal steps for real-time processing. To assess the performance, a semi-real simulation was executed with measured short-channel signals in a motor-task experiment. Compared with conventional techniques with no elements of SBF, the proposed method showed significantly higher estimation performance for true brain activation under a task-related SBF artifact environment.


Introduction
Functional near-infrared spectroscopy (fNIRS) is an effective non-invasive neuroimaging technique that measures hemoglobin concentration changes as an indicator of activity in the cerebral cortex [1].Those changes are due to the relative changes of oxygenated and deoxygenated hemoglobin (∆Oxy-Hb and ∆Deoxy-Hb, respectively).The advantages of fNIRS include low cost, portability, safety, and fewer constraints for participants, which facilitates its application to a wide variety of fields such as neuro-rehabilitation and brain computer interfaces [2][3][4][5][6].In such applications, real-time signal processing is necessary to simultaneously monitor true brain activation and exclude any negative effects of artifacts in order to give effective feedback stimuli or to control some devices.
However, scalp blood flow (SBF), which exists on the near-infrared light signal path, contaminates task-related cerebral blood flow (CBF).Activation of the autonomic nervous system, which is easily affected by experimental tasks, leads to physiological changes (i.e., cardiac, respiratory, Mayer waves and low frequency) that affect the task-related evocation of SBF as a specific global artifact of fNIRS (scalp-hemodynamics artifact).Even though this negative-effect artifact shows task-related changes, it does not reflect brain activation [7][8][9][10].Furthermore, it is distributed over the whole head, thus masking true cerebral activity (for review, [1]).Because SBF has similar temporal profiles and is highly correlated with CBF, it is difficult to remove this effect with conventional methods such as filtering and block averaging.Therefore, the reduction of SBF global interference from fNIRS signals remains as a challenging task.
To prevent this effect, previous studies have employed shorter-distance channels to only measure the temporal pattern of SBF [11][12][13].Particularly, Sato et al. proposed ShortPCA GLM (expressed as off-line ShortPCA GLM in this paper) to solve the problem with a few 15-mm wide multi-short channels [14].They extracted a global temporal pattern of SBF from short-channels, which can be explained by the first principal component of principal component analysis (PCA), and removed SBF signals from long-channels using a general linear model (GLM).Even though their method did not need the same number of short-channels as long-channels, it performed better at rejecting the global SBF artifact with long-distance channels, which were possibly contaminated with SBF when compared with conventional techniques.However, their method is still impeded by an off-line restriction.
Because the aims of real-time signal processing of fNIRS-based applications are to estimate true brain activation and reduce artifacts from raw signals, we proposed a real-time ShortPCA GLM (rt-ShortPCA GLM) with sliding-window analysis (SWA), which is sometimes used for parameter dynamics estimation in functional magnetic resonance imaging (fMRI) [15,16] and fNIRS analysis [5] or removing motion artifact from fNIRS signals [17].This simple technique may retain the advantages of off-line ShortPCA GLM, i.e., simplicity and low physical constraint.For the performance assessment, we executed a semi-real simulation with artificial fNIRS signals constructed with a mathematical model of CBF, measured real SBF and white noise, and then compared it with conventional GLM methods and off-line.We also investigated the effect of sliding-window length on estimation accuracy.

Sliding-Window Analysis
For the overall proposed method, SWA is important for real-time estimation.SWA sets a certain sample range sliding-window, and the most recent data are used for computing statistics [16].When a new sample is measured, the sliding-window slides to the next sample, which includes the newest and discards the oldest sample (Figure 1).As the sample window is maintained at the same size, it has the advantage of constant dynamics estimation for the same parameters.In this study, we set a 350-sample wide sliding window and it covered at least a single-trial [5,16].The sliding-window was applied for each measuring time (i.e., excluding only the oldest sample).We also tested different window length simulations to assess the effect of length (ranging from 50 to 2077 samples with 50-samples step).Diagram of the sliding-window analysis (SWA).Samples placed in the sliding-window, including the newest sample, are applied signal processing or used for analysis.A sliding-window will move to the next when the time instant k is increased.

Real-Time Signal Processing
Figure 2 shows a flow diagram of our present method (rt-ShortPCA GLM based on the method of off-line ShortPCA GLM).The key aspect of ShortPCA GLM is the design matrix, which contains a global temporal pattern of SBF explained by the first principal component of PCA.However, due to real-time restrictions, we cannot define this element before measuring samples.For this reason, we applied SWA in each signal step individually.
When a new raw sample was measured, at first, only L-sample-length-windowed long-and short-channel samples had a 2nd-order Butterworth bandpass filtering with cut-off frequency of 0.01-0.3Hz applied (these values are similar to those used in previous studies; for review, see [1]), then the newest sample was stored as filtered data (Figure 2I-II).Before storing, z-scoring (mean µ and standard deviation σ during all filtered samples were 0 and 1, respectively) was applied only to the newest samples of short-channels (Figure 2II), in order to prevent extracted principal components for the next PCA step from being biased to a specific channel signal because of motion artifacts, signal-to-noise ratios, or different path lengths [14,18].Furthermore, to estimate a global temporal pattern of SBF (global scalp-hemodynamics model; GSHM), we applied PCA to the windowed and preprocessed samples and extracted principal components.The newest sample of the first principal component was then stored to the design matrix of GLM (Figure 2III-IV).
For the final step, we employed GLM with SWA to estimate the newest sample of CBF by removing the weighted GSHM.The GLM expresses response variable y Long (time sample N T × number of channels N CH ) in terms of a linear combination of explanatory variables (design matrix) X (N T × number of elements N): where β is an unknown weight matrix (N × N CH ) as relative contributions of the assumed elements in long-distance channels.ε is a normally distributed model error term with a mean of zero and variance of σ 2 (i.e., ε ∼ N(0, σ 2 )).In this study, we constructed a design matrix with five elements similar to the off-line ShortPCA GLM: cerebral-hemodynamic model (CHM) with its temporal and dispersion derivatives (CHM (1) ,CHM (2) ), constants (Const.)and a global temporal pattern of SBF (GSHM) estimated in real-time PCA.Here, the CHM was generated by convolving the mathematical hemodynamic response function [19] and experimental task function (box-car function).All elements without the GSHM were normalized by scaling the maximum peak amplitudes to one.Note that all elements without the GSHM were able to prepare before real-time processing, thus we only had to update the GSHM for each measuring time.The optimized estimation of weight vector at time k, βk , was calculated by ordinary least squares (Figure 2IV) as follows: Here, N T was restricted to L-samples by SWA so that each element of the matrices in Equations ( 2) and (3) had the samples ranging from k-L+1 to k (see Figure 1 again).Owing to the changing data samples in the sliding-window and the transition of k, we could estimate weight β dynamics corresponding to the amplitudes at time k of the assumed elements in the design matrix.After the weight estimation, we can remove global scalp-hemodynamics at time k from each long channel by subtracting the real-time estimated GSHM multiplied by the βGSHM,k for the corresponding channel.

Artificial Functional Near-Infrared Spectroscopy (fNIRS) Signals
Some studies have assumed that the global interference of SBF can be expressed as periodic physiological changes with variable amplitudes and frequencies of sinusoidal waves without considering task-related changes [2,[20][21][22].In our assumption, the task-related fluctuation of SBF, which is easily affected by experimental tasks and individual differences, is completely contradictory.Accordingly, a task-related SBF formulation is unknown.
To acquire task-related SBF data for the performance assessment, we executed a block-designed experiment [14] and synthesized SBF data using measured short-channel data.The experiment was conducted according to the guidelines outlined within the Declaration of Helsinki and were approved by the institutional ethical committee of Nagaoka University of Technology.All subjects provided written informed consent prior to participation in the experiment, which consisted of 6 trial loops constructed as rest-task-rest (15 s each).Thirteen subjects (22-34 years old, one female) were asked to perform a motor task (finger tapping or ball grasping) with their right hand in the task section.During the experiment, 18 short-channels were placed on the motor-related areas and measured with a 100 ms sampling period at wavelengths of 780, 805, and 830 nm by a multichannel optical imaging system (FOIRE-3000, Shimadzu Corp., Kyoto, Japan).The modified Beer-Lambert Law was applied to quantify ∆Oxy-Hb and ∆Deoxy-Hb, using the least-squares method [23,24].Furthermore, these quantities had baseline correction by subtracting mean of samples in pre-rest section of each channel, a 2nd-order Butterworth bandpass filter with a 0.01-0.3Hz cut-off frequency, Gaussian spatial interpolation, and piecewise cubic Hermite interpolation.All of the processing was applied off-line.Since ∆Oxy-Hb reflects task-related brain activation more than ∆Deoxy-Hb, the estimation improvement by removing the global SBF artifact of ∆Oxy-Hb is a more important challenge [20,25].Moreover, due to the difficulty of measuring ∆Deoxy-Hb with our wavelength settings and its low signal-to-noise ratios [5,18], we only used ∆Oxy-Hb to generate artificial fNIRS signals in this study.Note that the unit of ∆Oxy-Hb is nanomole (nM).
Artificial fNIRS signals were generated by assuming that 43 long-channels and 4 short-channels were placed over the motor-related areas in both hemispheres (bilateral primary sensorimotor cortices, dorsal premotor cortices, and supplementary motor areas) with a sampling period 130 ms.Each channel signal, y, consisted of the linear combination of 3 elements: the mathematical model of CBF [19], measured SBF, and white noise as follows: where a i and y denotes an amplitude adjustment parameter and time-series vectors for each signal, respectively.y CBF had common temporal pattern across eleven brain-activation channels and randomized amplitudes during trials with normal distributions, N(300, 100 2 ).Furthermore, we applied spatial randomization selected from N(300, 100 2 ) across whole channels to the peak amplitudes of trial-averaged signals.Therefore, the amplitude ratios of channels were always the same.
For simplification, the randomized CBF templates were common across the thirteen subject samples.Signal processing was applied to the measured SBF signal, y SBF , as mentioned above.Then, the peak of block-averaged SBF was selected to be equal to CBF (i.e., N(300, 100 2 )).The standard deviation of white noise, ε wn , was randomly selected from the normal distribution N(400, 180 2 ).The randomized parameters were based on those in previous studies, which used executed performance assessment with artificial data [20,22].Finally, fNIRS signals were synthesized with those elements using Equation (4) in the following three conditions: (a) equal amount of CBF and SBF (a 1 :a 2 = 1:1), (b) CBF varied from zero to five times larger than SBF (i.e., a 1 /a 2 = 1 ∼ 5 at intervals of 0.5), and (c) the opposite of (b) (i.e., a 2 /a 1 = 1 ∼ 5 at intervals of 0.5).The results of condition (b) and (c) were used to investigate the robustness of our proposed method under several CBF and SBF amplitude-ratios because the actual ratio of the fNIRS signal is unknown.Additionally, in the (b) and (c) conditions, the magnitude of the white noise components was fixed as described above (for details see the experiment, simulation, and acquired data from the study detailed in [14] including its supplementary data).Note that the range of amplitudes, trial-amplitude randomization, and their units were changed from the original.

Performance Evaluation
To investigate the superiority of our proposed method, we prepared two comparative conventional design matrices used in fMRI and fNIRS studies that did not have the following SBF elements: simple (CHM + Const.) and standard (CHM + CHM (1) + CHM (2) + Const.)[5,16,[26][27][28][29].Both of the comparative design matrices were also applied to SWA for real-time CBF estimation.After pseudo real-time processing, for performance quantification we calculated the root-mean-squared error (RMSE) for the true CBF, y CHM , and estimated CBF, ŷCHM , and adjusted R 2 as an estimation error and GLM fitting score, respectively.Both of the performance indicators were obtained by averaging across the selected long-channels (all: 43 long-channels; activated: 11 long-channels assuming CBFs; non-activated: 32 long-channels assuming no CBFs) and subjects.To assess the statistical superiority of our method, we applied two-tailed paired t-tests to the performance indicators between rt-ShortPCA GLM and each of the other two design matrices (with Bonferroni correction, n = 2, significance level α = 0.05/2).Note that we excluded the first trial data from statistical analysis due to potential instability of how the first-trial-window sized samples were stored.Furthermore, a performance comparison of on-line (by applying sliding-window) and off-line (by batch analysis) GLM analyses was also carried out.In off-line GLM, the time series of task-related CHM and its derivatives were devided by task to represent each individual task.By this transformation, the number of elements in the off-line design matrix was 20, including 3 elements × 6 repetitions of the task.Evaluation statistics were performed using two-tailed paired t-tests (α = 0.05).

Performance Comparison
Figure 3A,B shows the statistical comparison of RMSE and adjusted R 2 between the proposed method and comparative design matrices that have no elements to explain SBF with a 350-sample-wide SWA.Those indicators expressed in the performance of our proposed method was significantly higher than the conventional design matrices: lower estimation error (for vs. simple, paired t = 3.41, p < 0.003; for vs. standard, paired t = 3.64, p < 0.002) and higher GLM fitting (for vs. simple, paired t = 17.96, p < 0.0001; for vs. standard, paired t = 12.16, p < 0.0001) in task-related SBF situations.In contrast, Figure 3C, comparing with off-line ShortPCA GLM, shows that the RMSE of rt-ShortPCA GLM (on-line) was significantly higher (paired t = 5.72, p < 0.0001) even though there was no significant difference in adjusted R 2 between on-line and off-line (paired t = 1.71, p < 0.1).

Performance of Several Window Length
Figure 4 shows RMSE, obtained by different-window-length SWA ranging from 50 to 2077 samples at 50 sample intervals and averaging across all long channels.Estimation errors, which were estimated by SWA with shorter than 350-sample closing to the single-trial sample length (i.e., L < 350 samples), were substantially increased.In contrast, RMSE decreased gradually when we employed longer-sample-wide sliding-window than single-trial-samples window (i.e., L ≥ 350 samples).We also checked that the means of RMSE across activated-or non-activated-channels showed similar results shown in Figure 4.

Performance under Several Cerebral-and Scalp-Hemodynamics Amplitude Ratios
The actual ratio between CBF and SBF is still unknown.To investigate the robustness of the estimation accuracy of our proposed method under several amplitude ratio situations, we executed simulations under the two-conditions, (b) and (c), described in Section 3.1.The result obtained by applying real-time GLM with 350-sample-length SWA is shown in Figure 5.Under condition (c), the larger amplitude of with a SBF than of CBF, the proposed method reduces estimation error compared with the conventional design matrices.By contrast, the RMSE averaged across activated-channels in condition (b) was more sharply increased along with the increase of CBF-amplitude ratios.As a result, RMSEs averaged across all-channels of rt-ShortPCA and standard but not simple was increased.  .RMSE in several amplitude ratios between cerebral blood flow (CBF) and scalp blood flow (SBF).RMSE of each case was averaged across three-channel sets: all (43 long-channels), activated (11 long-channels containing evoked CBFs) and non-activated (32 long-channels containing no CBF).

Advantages of Proposed Method
Our proposed method, using a semi-real simulation with the assumption that fNIRS signals were possibly contaminated with global task-related SBF, achieved a significantly higher estimation accuracy than conventional design matrices of GLM, which have been used in several studies [5,16,[26][27][28][29].However, brain activity estimations with real fNIRS data are more challenging.To investigate the performance of the present method in real data, we partially tested the data obtained by the real experiment which was executed with the same experimental settings as the simulation (see details in Experiment 2 [14]).Figure 6 shows an example of the estimated cerebral-hemodynamics from Ch 21 and 27, placed on the left (ipsilateral) right (contralateral) primary sensorimotor cortices, respectively.The estimated waveform of rt-ShortPCA (represented by On-Line ShortPCA in the figure ) shows the most similarity with the off-line.In contrast, the waveforms of the other design matrices (i.e., simple and standard) show particularly larger amplitudes of cerebral-hemodynamics in the 4th-6th task sections.Sato et al. reported that the off-line ShortPCA GLM achieved a performance very close to the fMRI estimation of the identical real data [14].Thus, by assuming that the estimated wave by off-line ShortPCA GLM is true, the proposed method reduces the negative effect of task-related scalp-hemodynamics and achieved better estimation accuracy than the conventional design matrices in real time.This result partly indicates that the present method possibly effects for scalp-hemodynamics artifact reduction in real-time.To evaluate the quantitative performance, further studies should evaluate the time course of t-maps of cerebral activity and compare these with fMRI data for more practical use.Applying SWA for real-time artifact reduction improves estimation accuracy without reducing the benefit of off-line ShortPCA GLM (i.e., simplicity and few physical constraints).Hence, this method is possibly useful for real-time applications.Previous studies have shown that fNIRS may be helpful, for both diagnosis and therapy, in epilepsy, dementia and stroke [30][31][32][33].In clinical practice, the measurement area should be as broad as possible, in order to detect all potential brain activation differences in patients compared with healthy subjects, and stress should be minimized as much as possible.On-and off-line ShortPCA GLM can provide such broad measurement areas with few physical constraints because they need a limited number of short channels to extract the global temporal pattern of SBF.Furthermore, with real-time monitoring of brain activation, diagnosis times should be reduced, which may cause participants to experience less stress.Diagnosis quality should also be improved by real-time artifact reduction.The proposed method may also be utilized in therapeutic situations, especially with neurofeedback training.Kober et al. reported that the effectiveness of neurofeedback training can be reduced if feedback signals are not related to subject-dependent true brain activation [4].Thus, false positive brain activation, affected by task-related scalp-hemodynamics, probably reduces the quality of these training effects.For such cases, our proposed method would be powerful and suitable for supplying high quality feedback signals to a subject.Additionally, we can expect advantages for high-precision brain computer interfaces using task-related cerebral activity as the input signal.It should also be noted that the proposed method can be easily implemented in a conventional fNIRS system because the algorithms (PCA, GLM and SWA) are simple and the number of short channels required is small.In summary, the method should be very practical for real-time applications.

Effect of Sample Length of Sliding-Window
Previous SWA-based GLM often used single-trial or longer length sliding-windows to avoid lacking samples of the mathematical model of CBF in a design matrix corresponding to the task section [5,16].However, the influence of changing the fixed length of a sliding-window on brain activity estimation accuracy is still unknown.The result in Section 4.2 gives us some possible answers for the question.Interestingly, even though the sample size was restricted by SWA, the length of the sliding-window was not significant for improving estimation accuracy.For example, with a sliding-window set at a two-trials-sample length, the weight for CBF, which represents the amplitude of cerebral-hemodynamics at time k, would be estimated as a summarized value across two-trials of CBF though the amplitudes between current and previous trial are different.This leads to lower tracking capabilities for the signal fluctuations during trials and retains the estimation errors, even though a sufficient length to explain task-related CBF evocation was employed.In contrast, with a shorter than 350-sample-window length, the RMSEs dramatically increased because the CHM element lacked some samples by sliding-window cannot explain whole of task-related evocation of CBF.Thus, our results indicate that by using a one-trial sample length (350 samples in this simulation setting), SBF artifacts can be sufficiently removed.Additionally, sample length has an effect on computational complexity.A calculation complexity using ordinary least squares is O(N 3 T ) depended on time samples in the sliding-window [34].This means we should avoid using longer-sample length sliding-windows than necessary for fast computation in real-time.It should be noted that the averaged process time for each time sample using a 350-sample-wide sliding-window of our proposed method was approximately 26 ms, sufficiently shorter than the sampling period of the simulation setting (130 ms) when we processed it using the MATLAB software (version 9.2.0.556344 (R2017a), Mathworks Inc., Natick, Massachusetts, United States) on a general desktop personal computer (Intel R Core TM i7-4770K, 3.50 GHz, Intel Corp., Santa Clara, CA, USA).

Limitations
Unfortunately, there are still limitations corresponding to SWA achieving good estimation performance.One limitation arose when the ratio of CBF amplitude was more than two times larger than that of SBF, obtained in Section 4.3.Even though our proposed method is based on off-line ShortPCA GLM, whose estimation is robust even under larger CBF ratio conditions (see further details in the simulation section of the supplementary files in [14]), the result we obtained was not robust.It was possibly due to the standard design matrix which has a similar increase in activated channels with rt-ShortPCA GLM.Both of these design matrices but not the simple design matrix contained temporal and dispersion derivatives (CHM (1) and CHM (2) , respectively).Those elements represent individual differences of width and latency [19], even though those are regardless of our simulation settings, no differences during the subject CBFs.Conversely, those elements will possibly work under subject-dependent differences.However, an applicable method for both conditions is desired, thus further studies should enhance the present performance in such situations.
Next, fixed-length SWA-based least squares showed slow convergence for the weight estimation [34].In particular, it arose owing to the difference in large amplitudes between trials.At the start of the current trial, samples in the sliding-window were mostly occupied by data from the previous trial.In such a case, β was affected by previous CBF activities, leading to a false estimation.β would gradually converge to the true value with the transition of k.To overcome this weakness, an enhanced, faster convergence method, for example, by variable forgetting factors or variable-length sliding window analysis, should be utilized [34][35][36][37].
Finally, the proposed method does not account for the serial correlation that exists in the residual error term of the GLM equation, Equation (1).Due to the high sample rate of fNIRS, the signals have serially correlated errors which may lead to an inflated t-value calculated by sample-wise t-statistics.In this way, the brain activation may be overestimated.To overcome this negative effect, further studies should implement correction methods, for example, precoloring or an autoregressive model in real-time [38][39][40].

Conclusions
Due to the global task-related SBF artifact, which contaminates fNIRS signals, it is difficult to estimate true brain activation accurately.Our study aimed to remove this artifact in real-time with rt-ShortPCA GLM using SWA.From the results of a semi-real simulation using measured short-channel signals as SBFs, our proposed method achieved better CBF estimation accuracy than that of conventional design matrices with no elements of SBF by removing global SBF patterns, while retaining the advantages (simplicity and low physical constraint) of the original off-line ShortPCA GLM.

Figure 1 .
Figure 1.Diagram of the sliding-window analysis (SWA).Samples placed in the sliding-window, including the newest sample, are applied signal processing or used for analysis.A sliding-window will move to the next when the time instant k is increased.

Figure 3 .
Figure 3. Comparisons of performance indicators: root-mean-squared error (RMSE) and adjusted R 2 .(A) RMSE and (B) adjusted R 2 of rt-ShortPCA GLM were compared with on-line methods and with (C) off-line.Open circles show means of each subject across all channels.Filled circles and error bars represent the mean and standard deviation across subjects, respectively.

Figure 4 .
Figure 4. RMSE estimated by different-window-length SWA ranging from 50 to 2077 samples at 50 sample intervals.Colored circles at 350 samples represent the same results as in Section 4.1.

Figure 5
Figure5.RMSE in several amplitude ratios between cerebral blood flow (CBF) and scalp blood flow (SBF).RMSE of each case was averaged across three-channel sets: all (43 long-channels), activated (11 long-channels containing evoked CBFs) and non-activated (32 long-channels containing no CBF).

Figure 6 .
Figure 6.Typical temporal estimated waveforms of cerebral hemodynamics in ipsilateral and contralateral motor area symmetrically located by rt-ShortPCA GLM for real functional near infrared spectroscopy (fNIRS) data of a representative subject during a left-handed finger tapping task.The first trial samples from 0 to 45 s are not shown because of unstable estimated waves.