Accelerometer-based methods include adaptive filtering, active noise cancelation (ANC) (Kim et al., 2011), accelerometer-based motion artifact removal (ABAMAR) (Virtanen et al., 2011), acceleration-based movement artifact reduction algorithm (ABMARA) (Metz et al., 2015), multi-stage cascaded adaptive filtering (Islam et al., 2017), blind source separation, accelerometer-based artifact rejection, and detection (BLISSA2RD) (von Luhmann et al., 2019). The introduction of the accelerometer improves the feasibility of real-time rejection of MAs.

The Wiener filter approach is the first to remove motion artifacts without incorporating additional hardware devices (Izzetoglu et al., 2005; Orihuela-Espina et al., 2010; Li L. et al., 2021). The technique assumes that the measured fNIRS signals are a simple addition between the actual fNIRS signals, x(n), and motion artifacts, v(n). Moreover, it is assumed that x(n) and v(n) are stationary and uncorrelated.

Consequently, the Wiener filter, g(n), minimizes the mean square error between x(n) and x^(n), that is,

Therefore, the optimum filter can be obtained using the orthogonality principle and simplified using Eq. (19) as follows.

We obtain the Fourier transform of the Wiener filter by converting Eq. (21) into the frequency domain using the Fourier transform, which is as follows.

where p_x(w) and p_v(w) denote the power spectral densities (PSDs) of the actual fNIRS signals and motion artifacts.

In application, a prior experiment is required to determine the PSDs of g(n) by calculating the values of x(n) and v(n). Subsequently, we can determine g(n) on a time scale using the inverse Fourier transform and apply it to new experimental data.

The Wiener filter is the first attempt to remove motion artifacts without a reference signal from additional hardware devices, such as accelerometers. With the g(n) determined, the filter can be implemented for online applications. However, it requires prior knowledge of the PSDs of x(n) and v(n), which makes initial calibration more complex (a particularly designed paradigm is needed). The idea of building the filter model from a prior experiment inspired later research on motion artifact removal techniques.

The Kalman filter approach was also proposed based on the general idea of the Wiener filter but as a different model (Izzetoglu et al., 2010). The motion artifact-free fNIRS signal, x(n), was modeled using an autoregressive (AR) model. An AR model of order p can be written as

With the motion artifact-free data from the prior experiment, the parameters of the AR model, a_i, i = 1, …, p, can be determined using the Yule-Walker equations. Therefore, the process equation for the Kalman filter has the following form.

where ϕ(n) is composed of p motion artifact-free fNIRS signals, and ω_n denotes the zero-mean noises in the AR model with an error covariance matrix Q. Matrix A can be obtained using Eq. (23) as follows.

The measurement equation can be written as follows.

where ν_n denotes the measurement noise (such as the motion artifact) with an error covariance matrix R. z(n) denotes the motion artifact corrupted signal.

Subsequently, Eqs. (24) and (26) form the state-space model for the Kalman filter, see Table 4. The minus sign in the superscript denotes the prior estimate of a variable. Subsequently, the motion-free fNIRS signal can be obtained using the Kalman filter (Wan and Nelson, 2001). The Kalman filter method can be applied to any online application of optical intensities, optical densities, and concentration changes, and in the Kalman filter theory, both ω_n and ν_n are assumed to be zero-mean Gaussian white noise (Huang and Hong, 2021; Haghighi and Pishkenari, 2021; Li B. A. et al., 2021; Sun and Zhao, 2021; Yang et al., 2021; Lv et al., 2021; Tang et al., 2020; Pham et al., 2021). However, it is not the case for ν_n (motion artifacts do not observe the zero-mean Gaussian distribution). It may degrade the filter’s performance (Zhou et al., 2017). Moreover, matrices A and C were fixed once determined by the Yule-Walker method. Therefore, further development of the algorithm focuses on; (i) the compensation of the instrumental noises (Amian and Setarehdan, 2013) and (ii) adaptive adjustment of A and C, or using a more sophisticated and nonlinear model in place of Eqs. (24) and (26) (Dong and Jeong, 2018).

Additionally, the state-space model for the Kalman filter method is not unique. A typical example is the incorporation of the autoregressive iterative robust least-squares model (AR-IRLS), the general linear model (GLM), and two linear Kalman filters (Barker et al., 2016). A GLM and an AR-IRLS replaced the state-space model in the Kalman filters. The GLM was introduced to describe the dynamics of hemodynamic responses and physiological noises. The AR-IRLS was introduced as compensation for the MAs in the signals. A difference between the GLM-based method and the AR model-based method is that the GLM-based method requires the information of experimental paradigms, while the AR model-based method does not. Despite the significant adaptation ability of the Kalman filter regarding its state-space model, its applications are limited due to the effort required to set its initial parameters (e.g., the error covariance matrices for the state and the observation).

The spline interpolation method was first proposed by Scholkmann’s group and is referred to as the MARA (Scholkmann et al., 2010; Selb et al., 2015; Lee et al., 2021). The method made two fundamental assumptions: (i) The measured fNIRS signal is a linear addition of the motion artifacts and the motion-free fNIRS signal, and (ii) in the motion-corrupted segments in the signal, the motion artifact component dominates the measured fNIRS signal. Therefore, the proposed MARA comprised two parts: (i) Motion artifact detection and segmentation, and (ii) motion artifact removal. A flowchart of MARA is illustrated in Figure 4. The spline interpolation method encompasses six processing steps. First, the moving standard deviation (MSD) is calculated within a moving time window of W samples and stored as msd(n).

where N* is the set of natural numbers.

We can determine the start and end points of the motion artifacts and store the indices of the corresponding samples in a vector ξ(n) by comparing with the MSD [i.e., msd(n) in Eq. (28)] with a user-defined threshold value T. If the MSD is smaller than T, the corresponding msd(n) will be assigned as zero. The start points and endpoints of motion artifacts can be extracted by considering the first and last samples of non-zero values in msd(n). Next, let us suppose that there are L segments of motion artifacts and let the motion artifact segments, z_MA(n), and the non-corrupted segments, z_NC(n), be expressed, respectively, as

Using Eqs. (29) and (30), the measured fNIRS signals can be segmented into non-corrupted segments and motion artifact segments.

In the next part, the spline interpolation method corrects motion artifact segments. Because the motion artifact components dominate the MA segments, the spline interpolation fitting of z_MA(n) can be viewed as the motion artifact component. The difference between z_MA(n) and its spline interpolation fitting is stored as z_diff(n):

Because z_NC(n) and z_diff(n) may have different magnitude levels, the final step involves correcting the signal levels for the entire time series. Each segment is parallel-shifted according to the mean of the previous segment and that of the target segment. Two empirical constant thresholds, α = 3^–1 Hz^–1⋅fs and β = 2 Hz^–1⋅fs (where the variable fs denotes the sampling frequency), were chosen for comparison. The detailed shifting rules are listed in Table 1 in Scholkmann et al. (2010).

The spline-interpolation method is used widely for offline analysis in fNIRS studies. It has also been included in some open-source toolboxes, such as HomER2 and NIRSLAB (Balardin et al., 2017). Moreover, the method applies not only to the optical intensities and optical densities but also to concentration changes. However, both the segmentation procedure (the procedures in the blue box in Figure 4) and the parallel-shifting procedure (in the reconstruction process) increase the difficulty for online filtering applications. Moreover, the filter performance depends on artifact detection results (Brigadoi et al., 2014). Some variations of the spline interpolation method have also been proposed (Jahani et al., 2018; Zhou et al., 2021).

The wavelet-based method eradicates motion artifacts by removing the corresponding wavelet coefficients (Molavi and Dumont, 2012; Pinti et al., 2015), without requiring auxiliary devices. The method made the same assumption as the Wiener filter case, that is, z(n) = x(n) + v(n). Based on discrete wavelet transform (DWT), fNIRS signals can be expanded as follows.

where i denotes the dilation parameter, k indicates the translation parameter, and i₀ denotes the coarsest scale. The scaling function φ_i0k and the wavelet function ψ_ik are as follows.

The functions Φ(n) and Ψ(n) correspond to the mother scaling function and mother wavelet function, respectively. With the fast wavelet transform, the coefficients of the DWT of z(n) become

where N denotes the number of samples in the signal and N = 2^J. g(n) and h(n) denote the wavelet filter bank high- and low-pass filters, respectively. Therefore, Eqs. (35) and (36) can be written in the matrix form as follows.

where W = [B₁, B₂, …, B_J–1, A_i0]^T, B_i = [b_i0, …, b_i(2∧i–1)], A_i0 = [a_i00, …, a_{i0(2∧i0–1)}]. Ω is an N by N DWT matrix. For any element w_ik in the wavelet coefficient vector W, the probability of observing its value can be written as

where NCDF(⋅) is a normal cumulative distribution function. The estimated standard deviation of w_ik can be approximated empirically as

where W_i denotes B_i or A_i0 according to the value of i. Given a probability threshold α, if p_ik < α, we consider that the corresponding wavelet coefficient w_ik is dominated by artifacts and set to zero. The filtered signals can be obtained from the updated wavelet coefficient vector as

The hat over x(⋅) denotes the filtered fNIRS signal, and the hat over W denotes the updated wavelet coefficient vector after thresholding.

Similar to the spline interpolation method, the wavelet-based method can be applied to optical intensities, optical densities, and concentration changes. Initially, the wavelet-based method was first introduced as an offline artifact removal algorithm. It may be adjusted to be an online method using the windowing technique; however, its filtering performance will be degraded. Several variations of this method have also appeared in recent years (Chiarelli et al., 2015; Shukla et al., 2018; Perpetuini et al., 2021). Gao Y. et al. (2022) addressed that the tuning of the probability threshold is crucial, and Wei et al. (2018) designed a dual-threshold structure to improve the performance of the wavelet-based method.

Based on the negative correlation between measured oxyhemoglobin (HbO) and deoxyhemoglobin (HbR) concentrations, the correlation-based signal improvement (CBSI) method was proposed in 2010 (Cui et al., 2010; Haeussinger et al., 2014). The algorithm made three fundamental assumptions: (i) HbO concentrations are strictly negatively correlated with HbR concentrations (close to –1); (ii) motion artifacts have identical effects on HbO and HbR, subject to a constant factor; and (iii) motion artifacts are not correlated to the actual concentration changes. The measured concentration changes were modeled as follows.

where v denotes the motion artifact of the HbR signals, and χ is a constant. Based on the first assumption, we can reasonably obtain

β denotes the free ratio accounting for the magnitude difference between HbO and HbR. Equations (41) and (42) can be used to obtain the following:

Based on the third assumption, we set the correlation between v and x_HbO to zero, yielding

We assumed that χ = β. Then, Eq. (44) can be reduced to

where σ_HbO and σ_HbR denote the standard deviations of z_HbO and z_HbR, respectively. Therefore, the final CBSI filter is

Once χ is decided based on prior-experimental data, the CBSI filter can be easily implemented for online applications using Eq. (46). As mentioned in (Cui et al., 2010), CBSI can be applied after an exponential moving average (EMA) filter to facilitate real-time analysis. The behavior of an EMA filter is similar to an online version of a band-pass filter. The EMA is described as follows.

L(n) denotes the signal filtered by a large-window moving average (window size as W_L), and S(n) denotes the signal filtered by a moving average of window size W_S.

Owing to the fundamental assumptions of the filter, the CBSI filter can only be applied to concentration changes. However, the negative correlation assumption between the HbO and HbR concentrations may not always be true for experimental data. Typically, the CBSI method is used for comparison in the literature (Janani and Sasikala, 2017; Reyes et al., 2018; Sherafati et al., 2020).

The solution based on principal component analysis (PCA) was developed from a multivariate method to remove physiological noises. The PCA-based method, entitled targeted PCA (tPCA), is similar to spline interpolation, except for using PCA instead of spline interpolation to fit the noisy segments in the signals (Yucel et al., 2014a). In this subsection, we subsequently explain the mechanism of Zhang’s method and tPCA.

In accordance with the regular PCA method, the systemic spatial interference subspace (the actual hemodynamic responses) dominates in the baseline over the measured concentration changes and it spans during the stimulation periods for the first few components. The measured fNIRS signals are first converted to HbO and HbR concentration changes and are band-pass filtered. The filtered data from a chosen channel set are formatted into two N_ch × N matrices, Z_HbO and Z_HbR, where N_ch denotes the number of channels and N is the number of time samples. When segmented into two parts according to the baseline and stimulation periods, the two matrices turn out to be as follows.

For simplification, we omit the subscripts HbO and HbR, primarily because the filtering operations are identical for both chromophores. Accordingly, the spatial correlation matrices for the baseline signals are as follows.

The eigen-values and the eigenvectors can be obtained by the eigen-decomposition of spatial correlation matrices:

where the rows of U_base contains the eigenvectors, and Σ_base contains the eigenvalues in its diagonal. Based on the dominant interference subspace assumption, we assume that the first r spatial eigenvectors span in the stimulation period as U_base,r, an N_ch by r matrix padding by zeros. Accordingly, the actual hemodynamic responses will be the projection of the measured fNIRS signals onto the orthogonal subspace of the first r eigenvectors. The superscript ^⊥ denotes the orthogonal subspace. Then, the filtered fNIRS signals were

where X_stim denotes the filtered fNIRS signal and matrix Id is an N_ch by N_ch identity matrix. Finally, we can block the average X_stim around the stimulus to improve the signal-to-noise ratio (SNR).

To remove MAs, tPCA applies the multivariate PCA method on MA segments. The eigenvectors in U_base are ranked in decreasing order of the percent variance of eigenvalues in Σ_base. The number of spatial eigenvectors span, r, is the minimum integer obtaining not less than 97% percent variance. Then, the filtered fNIRS signals of the segments are stitched back into the original signals after shifting the average of adjacent data segments (identical to MARA). The processing procedures are repeated twice or three times until there is no further improvement in the results.

Targeted PCA applies to optical densities. Due to the iteration in MA identification, the method cannot apply to online filtering. tPCA can remove both spike-like MAs and baseline shifts. It can better suppress high-frequency MAs than MARA (Yucel et al., 2014a). However, during application, the performance of tPCA may degrade if the MAs are not present in multiple channels (Brigadoi et al., 2014). The method requires multiple channels, which is not an appropriate choice when the number of channels is limited.

The idea of the method is to subtract a best-fit linear combination of the noise references from each signal (Robertson et al., 2010). For a k-channel regression, the measured signal is modeled as follows.

where z ϵ ℝ^{1 × N}, a ϵ ℝ1 × k, H ϵ ℝ^{(k + 1) ×N}, x ϵ ℝ^{1 × N}. The row vector, a, is the weight for each reference channel in H. The regressor H is defined as follows.

where z_ri (i = 1, 2, …, k) is the measurement of the ith reference channel. It is recommended to take co-located channels as reference channels in priority. Therefore, a can be estimated using the least-squares method.

The motionless signal, x, can be obtained by subtracting aH from z.

The multi-channel regression method does not require an auxiliary device and applies to optical intensities, optical densities, and concentration changes. It was first designed as an offline filtering method, but it can apply to online applications simply by replacing the least-square estimation procedure with recursive least-square estimation. When different types of artifacts affect collocated channels differently, the multi-channel regression method may not remove MAs entirely.

Independent component analysis (ICA) can decompose a signal to the weighted sum of multiple independent sources statistically (Kohno et al., 2007; Santosa et al., 2013). Therefore, the ICA solutions belong to blind source separation methods (von Luhmann et al., 2019). The key to removing MAs using ICA is how to select the components relevant to motionless signals or MAs. Kohno et al. (2007) suggested identifying the components related to artifacts using the coefficient of spatial uniformity (CSU).

Additionally, Santosa et al. (2013) estimated motionless signals by weighted reconstructing from independent components according to the t-values between each component and the desired hemodynamic response. The desired hemodynamic response is the convolution between the canonical hemodynamic response and the experimental paradigm (a boxcar function).

The ICA applies to optical intensities, optical densities, and concentration changes. The input requirements may vary according to the reconstruction strategy. Some researchers concluded that the performance of ICA was not satisfactory in some applications because fNIRS signals contained both non-instantaneous and non-constant coupling, correlated noise, and source dependencies (von Luhmann et al., 2019).

The temporal derivative distribution repair (TDDR) is an online filtering method based on the temporal derivative of fNIRS signals. Three basic assumptions were made for this method: Non-motion fluctuations are assumed to be normally distributed; most fluctuations are independent of motion artifacts; and the derivatives of motion artifacts dominate in the derivatives of fNIRS signals when they are present (Fishburn et al., 2019). Therefore, the algorithm can be divided into five steps.

Step 1: For time instance n, the variation in the fNIRS signal can be calculated as follows.

where the D before z(n) denotes the temporal derivative of the measured signal at instance n.

Step 2: The observation weight w(n), corresponding to Dz(n), is initialized as one.

Step 3: The observation weight is adaptively updated using a robust estimator, M-estimator, and Tukey’s biweight function:

where σ_s(n) denotes the scaled deviation at instance n. An estimate of σ_s(n) can be obtained from the weighted mean of the fluctuations, the absolute residuals of the temporal derivative, and the median absolute residual scale for the normal distribution. Let μ denote the weight mean of fluctuations, res(n) denote the absolute residuals of temporal derivatives, and σ(n) represent the median absolute residuals scale. Subsequently, we can obtain the following:

Repeat (54)–(57) and consider σ_s(n) = σ(n) until μ converges.

Step 4: The filtered temporal derivative can be obtained using μ and w(n) obtained in Step 3.

Step 5: The final filtered signal is the integration of corrected temporal derivatives.

Importantly, the filter can be implemented offline and online. The only difference is that when implemented online, the sum of w(n) and w(n)Dz(n) should be stored in each loop to facilitate the calculation in Eq. (59). The TDDR algorithm can be used for concentration changes, optical intensities, and optical densities. The performance and computation time of the TDDR depend on the convergence criterion in Step 3. In particular, in Eq. (61), the iterative application of the median operator may result in a significant computation time. Furthermore, as the total number of data increases, the computation time increases according to regarding the complexity of the sorting algorithm in the median filter. Fishburn’s group claimed that 130 ms were required to process a dataset of 20 Hz (50 ms in period) (Fishburn et al., 2019). A long computation time may lead to a loss in data transfer or a delay in the data output. In addition, high-frequency noise may degrade the performance of TDDR.

The transient artifact reduction algorithm (TARA) is an offline artifact removal method. The algorithm attempts to estimate motion artifacts instead of filtering the original signal directly (Selesnick et al., 2014). Generally, motion artifacts can be divided into two types. The Type 1 artifacts are modeled as spikes, whereas Type 2 artifacts are modeled as additive step shapes. The corrupted fNIRS signals were modeled as a low-pass filter/compound noise denoising problem (LPF/CSD).

where ν₁ and ν₂ denote the two types of MAs, and ω denotes the Gaussian white noise. Besides, the desired uncorrupted fNIRS signal x is assumed to approximate an all-zero signal when filtered by an appropriately chosen high-pass filter, denoted by HF. Accordingly, when designed as a zero-phase recursive discrete-time filter, H can be factorized into matrices A and B such that

The first-order difference operator matrix D is defined as follows.

Inversely, the discrete-time integrator S is defined as the matrix satisfying DS = Id_{N × N}, where Id_{N × N} indicates an N by N identity matrix. Then, matrix B can be further factorized as

Subsequently, we can change the variables using the following equations.

Type 2 MAs are expressed in Eq. (69) because their first-order difference resembles Type 1 artifacts.

The two types of artifacts can be estimated using the following procedures when the factorization matrices A and B for the desired high-pass filter, two weight parameters λ₁ and λ₂, and two penalty functions r₀, r_1, and r₂ are available.

Step 1: Input fNIRS signals z, three regularization parameters λ₀, λ₁, and λ₂, penalty functions r₁ and r₂, and factorization matrices A and B for the desired high-pass filter (they can be calculated from a given degree, cutoff frequency, and length of z).

Step 2: The first additive terms for the estimator for both types of artifacts are as follows.

The vectors u₁ and u₂ will be initialized as zero vectors.

Step 3: Repeat the following calculation until u₁ and u₂ converge.

where Λ_i, (i = 0, 1, or 2) is diagonal, ψ_i = u_i/r_i’, [⋅]_i,j denotes the element located at the ith row and jth column of the matrix, and [⋅]_i denotes the ith element of the vector. The superscript prime notation, ′, is the derivative operator.

Step 4: The estimated artifacts can be obtained using Eqs. (68) and (69). Accordingly, the filtered fNIRS signals becomes the following.

The TARA method is applied to optical densities, optical intensities, and concentration changes. TARA exhibits a similar performance to the wavelet-based method but is better at removing Type 2 motion artifacts. TARA is limited to Type 1 and Type 2 artifacts but is not designed for oscillatory transients.

The dual-stage median filter (DSMF) method resembles an online method for motion artifact removal (Huang et al., 2022). Similar to the TARA, the design of a DSMF is based on Type 1 and Type 2 artifact categories.

The measured fNIRS signal, z, is first filtered using the first median filter with a window size of W₁.

Subsequently, the difference between z and x₁ is filtered using the second median filter with a window size of W₂.

The output of the second median filter conserves the variations in the uncorrupted signals; however, it is biased. Therefore, the final step involves compensation of the bias using the initial time point.

W₁ can be set as 8 s, while W₂ as 18 s.

The DSMF is applied to optical intensities, optical densities, and concentration changes. Unlike the TDDR method, although two median filters are used in the algorithm, they are called once within fixed window sizes for each estimate. Moreover, registers can store the sorted data sequence during real-time applications, thereby reducing computation time. Therefore, a long computational time problem can be solved relatively quickly. Similar to the TARA, the DSMF is not designed for oscillatory transients.

The convolution neural network (CNN) was adopted for MA removal in 2022 (Gao Y. et al., 2022; Kim et al., 2022). The denoising auto-encoder (DAE) model adopted a serial structure incorporating max-pooling and up-sampling layers (Gao Y. et al., 2022). The model enabled an offline MA removal by minimizing the linear combination of mean square error (L_mse), variance (L_var), standard deviation (L_std), and amplitude loss (L_amp).

where θ₁, θ₂, and θ₃ were set as 1, 1, and 10. The model parameters were estimated using simulated data generated from experimental data in the training process. The CNN structure is shown in Figure 1 of the reference (Gao Y. et al., 2022).

Kim et al. (2022) adopted the CNN in the U-net structure for offline MA removal. Different from the DAE model, the U-net structure accepted an input containing two feature channels.

where u(k) is the input of U-net of the kth package, J is the total number of packages, and dHR(k) is the desired hemodynamic response of the kth package. dHR(k) is the convolution between a canonical hemodynamic response function and a stimulus function.

The U-net solution adopted the mean square errors between the ground-truth motionless fNIRS signals and the estimated signals as its loss function in the training process. The configurations of the U-net structure are shown in Figure 1 of the reference (Kim et al., 2022).

The DAE model can be applied to optical intensities for offline filtering. It takes measured data in a full-time sequence, outputs the estimated motionless signals, and is trained using simulated data from an AR model. Therefore, the quality of the simulated data has a significant impact on the filtering performance of the model.

In contrast, the U-net structured CNN was applied to concentration changes. The training dataset was obtained through a semi-simulation process, and the size of each package was fixed as 1,024 samples. If the fNIRS data are obtained at a low sampling frequency, the package size is sufficient, but the filtering performance for high-sampling-frequency fNIRS data remains unknown. The matrix structure of the input data enables an easy adaptation to fNIRS measurement with multi-distance configurations. However, similar to the GLM-based method, the information of the stimulation paradigm is a prerequisite in the input. The inclusion of the stimulation paradigm is critical to obtain good performance of the method.

The methods presented above are the fundamental techniques available from the existing MA removal solutions. Among the new techniques under investigation, many solutions were developed from or in stack of one or several basic techniques. Table 5 demonstrates the relationship between other new solutions and these basic techniques. For a fair comparison among these techniques, systematic approaches are required and the results are presented.

Comparisons have been made among the trial rejection method, recursive least squares adaptive filtering method (Park et al., 2020), wavelet-based method, ICA, two- or multi-channel regression method, Kalman filter, MARA, PCA, moving average, band-pass filter, median filter, Savitzky–Golay filter, and CBSI for adults’ fNIRS data. Overall, the MA correction approach is better than the trial rejection approach (Brigadoi et al., 2014). MA correction can also help to improve classification accuracy (Janani and Sasikala, 2017). Wavelet-based method, ICA, multiple-channel regression, and MARA had a good performance in noise suppression. When processing children’s fNIRS signals, researchers compared PCA, MARA, wavelet-based method, moving average, CBSI, and tPCA. MARA, tPCA, and CBSI can retain a higher number of trials. The tPCA and MARA are more robust. Moving average, wavelet-based method, and tPCA have a better performance than other methods. Researchers also examined the filtering performance of infants’ fNIRS data. Trial rejection, wavelet-based method, tPCA, MARA, and different processing pipelines were compared in their studies. MA correction can retain many trials, but hemodynamic response functions were also suppressed. The optimal pre-processing pipeline depends on the dataset (Gemignani and Gervain, 2021). Similar to the conclusion for adults’ data, MA correction is better than trial rejection. The wavelet-based method alone or in a stack with MARA and the wavelet-based method (iqr = 0.5) had better performance in noise suppression. MARA in a stack with the wavelet-based method can recover a higher number of trials.

Reportedly, the wavelet-based filter has good performance when MAs are caused by sudden large movements (or spikes). Moreover, ICA and multi-channel regression also perform well in terms of noise suppression (Robertson et al., 2010). Cooper et al. (2012) reported that PCA, MARA, wavelet-based method, Kalman filter, and trial rejection method could reduce the MSE in corrupted fNIRS signals (Cooper et al., 2012). PCA, MARA, and wavelet-based methods exhibited much better performance than the other methods (Piper et al., 2014). In contrast, Brigadoi et al. (2014) claimed that the Kalman filter performed better using the area under the curve (AUC ratio) of the mean hemodynamic responses than PCA, and MARA, followed by the wavelet-based method and CBSI. They also claimed that the wavelet-based method exhibited a better performance in reducing the AUC for the first two seconds after stimulus onset (AUC_0–2).

Therefore, the evaluation of MA removal techniques depends on the evaluation metrics of the researchers. Due to this phenomenon, multiple metrics need to be considered when comparing various methods. Cooper et al. (2012) used MSE to assess how different solutions can suppress MAs and contrast-to-noise ratio (CNR) to evaluate how well the solutions can retain brain activation signals. Janani and Sasikala (2017) reported that a band-pass filter, median filter, and Savitzky-Golay filter could increase the CNR. Additionally, the performance of the MA solutions differs from each other when the MA types are different. Figure 5 showed an example of MA removal using different filters.

The optimal solution for different age groups can be different. Hu et al. (2015) tested the MA removal performance of six solutions using fNIRS data from children: wavelet-based method, PCA, MARA, moving average, CBSI, and the combination of moving average and wavelet-based methods. It was found that MARA successfully suppressed all types of MAs but had an unsatisfactory AUC ratio. Moving average and wavelet-based method turned out to be better solutions for children. Meanwhile, CBSI was efficient in removing spike-like and fast step-like MAs. Besides, Reyes et al. (2018) claimed that tPCA was the best solution for processing children’s fNIRS signals and the authors added that CBSI occasionally produced unstable hemodynamic responses. However, the authors also claimed that MARA was robust and performed well in both AUC and the average standard deviation of each trial-specific hemodynamic response. The contradiction between Reyes and Hu’s findings regarding the MARA originates from improper artifact correction (Hu et al., 2015). Another study by Di Lorenzo et al. (2019) used a combination of MARA and wavelet-based methods and found it to be more effective for infant data than using trial rejection, MARA, or wavelet alone. They further justified that using MA correction could retain many trials. Similar results were obtained by Behrendt et al. (2018). They claimed that trial rejection was not recommended for infants’ data. Table 6 summarizes the data type, age group, techniques in comparison, and main conclusions of all nine comparison papers.

According to Table 5, three types of data exist: (i) Simulated data, (ii) semi-simulated data, and (iii) experimental data. A semi-simulated dataset can be generated by adding desired concentration changes to MA-corrupted experimental data at resting state, or by adding uncorrupted experimental data at the resting state to the sum of desired concentration changes and designed artifact sequence. For both simulated and semi-simulated datasets, the ground-truth motionless signals are known. Nevertheless, this is not true for experimental datasets. Therefore, the input requirements of each metric decide the metric’s compatible data types. Table 7 summarizes the data type compatibility of different metrics.