Motion Artifacts Correction From EEG and fNIRS Signals Using Novel Multiresolution Analysis

Physiological signal measurement and processing are increasingly becoming popular in the ambulatory setting as the hospital-centric treatment is moving towards wearable and ubiquitous monitoring. Most of the physiological signals are highly susceptible to various types of noises, especially movement artifacts. The electroencephalogram (EEG) and functional near-infrared spectroscopy (fNIRS) signals are no exception to motion artifacts, which become prominent in the ambulatory setting. Since successful detection of various neurological disorders is greatly dependent upon clean EEG and fNIRS signals, it is a matter of utmost importance to remove motion artifacts from these two signal modalities using reliable and robust methods. This paper proposes three novel multiresolution analysis techniques: i) Variational mode decomposition (VMD), ii) VMD in combination with principal component analysis (VMD-PCA), and iii) VMD in combination with canonical correlation analysis (VMD-CCA), for motion artifact correction from single-channel EEG and fNIRS signals. The efficacy of these novel techniques is validated by computing the difference in the signal to noise ratio (<inline-formula> <tex-math notation="LaTeX">$\Delta SNR$ </tex-math></inline-formula>) and percentage reduction in motion artifacts (<inline-formula> <tex-math notation="LaTeX">$\eta$ </tex-math></inline-formula>). Among the three proposed novel methods, VMD-CCA decomposed with 15 intrinsic mode functions (IMFs) has shown the best denoising performance for EEG signals producing an average <inline-formula> <tex-math notation="LaTeX">$\Delta SNR$ </tex-math></inline-formula> and <inline-formula> <tex-math notation="LaTeX">$\eta $ </tex-math></inline-formula> values of 23.81 dB and 57.01%, respectively for all 23 EEG recordings. On the other hand, for the available 16 fNIRS recordings, VMD-CCA decomposed with 10 IMFs produced an average <inline-formula> <tex-math notation="LaTeX">$\Delta SNR$ </tex-math></inline-formula> and <inline-formula> <tex-math notation="LaTeX">$\eta $ </tex-math></inline-formula> values of 15.97 dB and 39.01%, respectively. The results reported using the proposed methods outperform most of the existing state-of-the-art techniques.


I. INTRODUCTION
Electroencephalography (EEG) measures the electrical activity of the human brain quantitatively originating from the The associate editor coordinating the review of this manuscript and approving it for publication was Humaira Nisar . firing of neurons [1]. Generally, such action is recorded using several electrodes which are positioned in different areas of the scalp [2]. The use of EEG, in the area of epileptic seizure detection, is the most stereotypical usage among its several utilizations [3], [4]. As epileptic activity is typically problematic to identify due to the unpredictability of the exact time of the occurrence, long-term monitoring is frequently being used to enhance the possibility of detection. This continuous monitoring for an extended period brings discomfort to the patients as they must stay stationary during the process. An in-home recording of the EEG and/or personal healthcare domain will elongate recording durations while ensuring more accuracy of the recording signal. Furthermore, such facilities can proliferate the user's level of comfort. Additional noteworthy utilization of EEG includes the measurement of drowsiness levels [5]- [8], emotion detection [9], cognitive workload [6], [10], and brain-computer interfaces (BCI) [11]- [16]. All these usages have potential applications in the personal healthcare domain. Recently, due to the inherent anti-spoofing capability of EEG signals, the implementation of biometric systems using EEG is being studied and has already shown promising outcomes [17].
The functional near-infrared spectroscopy (fNIRS), a noninvasive optical imaging technique, measures changes in hemoglobin (Hb) concentrations within the human brain [18]. Light of various wavelengths in the infrared band employed by fNIRS infiltrates the skull to measure the change in the concentration levels of oxygenated (oxy-Hb) and deoxygenated hemoglobin (deoxy-Hb) in the human brain by estimating the difference in the optical absorption [19]. The fNIRS signal distinguishes the activation changes in the cortex, allowing optical measurements to be used for imaging brain functions [20], [21]. Medical applications of fNIRS focus on the noninvasive measurement of the quantity and oxygen content of hemoglobin, cognitive tasks identification [22], [23], and BCI [24]- [26].
Both EEG and fNIRS, two very vital physiological signals, are prone to motion artifacts that occur due to the voluntary or involuntary movement of the test subject during signal acquisition. Constraining a test subject completely from physical movements, voluntary or involuntary, is always very difficult. Consequently, the EEG and fNIRS signals get corrupted to some extent by motion artifacts. In some cases, this contamination may become so prominent that the recorded signals may lose their usability unless the motion artifacts are removed or reduced by a significant amount [27].
Apart from motion artifacts, physiological signals suffer from other forms of artifacts as well. During simultaneous EEG-fMRI experiments, the EEG signal is usually contaminated by gradient artifacts (GA) and pulse artifacts (PA) [28]- [30]. Event-related fNIRS signals are normally contaminated by heartbeat, respiration, Mayer waves, etc., as well as extra-cortical physiological noises from the superficial layers [31].
Several approaches were investigated for the reduction of motion artifacts from the EEG [32], [33]. The performance comparison using several motion artifact removal techniques for the EEG signal was made by Sweeney et al. in [34]. Various multiresolution analysis methods, such as the discrete wavelet transform (DWT) [35], empirical mode decomposition (EMD) [36], ensemble empirical mode decomposition (EEMD) [37] were investigated. Authors of [34] also combined EMD with CCA (EMD-CCA), EMD with independent component analysis (EMD-ICA), EEMD with ICA (EEMD-ICA), and EEMD with CCA (EEMD-CCA). Maddirala and Shaik [38] used singular spectrum analysis (SSA) [39] to filter out the motion artifacts from the EEG signal. The wavelet-based techniques were proven to be an effective multiresolution approach for the analysis and decomposition of the EEG signal [40]. To filter out motion artifacts from the EEG signal, DWT along with the thresholding technique was utilized in [41]. Recently, to reduce motion artifacts from the EEG signals Gajbhiye et al. [42] employed wavelet-based transform along with the total variation (TV) and weighted TV (WTV) denoising techniques whereas in [43], wavelet domain optimized Savitzky-Golay filter has been incorporated.
In the last few decades, multiple motion artifacts removal techniques was developed [44]- [46] for the removal of motion artifacts from the fNIRS signal. Sweeney et al. [47] used adaptive filter, Kalman Filter, and EEMD-ICA. Scholkmann et al. [48] used the moving standard deviation and spline interpolation method whereas in [49] waveletbased method was employed. The authors of [34] used DWT, EMD, EEMD, EMD-ICA, EEMD-ICA, EMD-CCA, and EEMD-CCA. In [50], Barker et al. used an autoregressive model-based algorithm to reduce motion artifacts from the fNIRS signals. Kurtosis-based wavelet transform was proposed in [51] whereas Siddiquee et al. [52] utilized nine-degree of freedom (DoF) inertia measurement unit (IMU) data to estimate the movement artifacts in the fNIRS signals using autoregressive exogenous (ARX) input model. A hybrid algorithm was proposed in [53] to filter out the movement artifacts from fNIRS signals where both the spline interpolation method and Savitzky-Golay filtering were incorporated.
It is evident that the development of robust algorithms capable of reducing movement artifacts from EEG and fNIRS signals is of utmost importance; otherwise, the interpretation of the signals will be inaccurate. Also, the requirement of robust algorithms is required since EEG and fNIRS signals are highly non-stationary [34]. In this paper, three novel motion artifacts correction techniques have been proposed which are capable of removing motion artifacts from singlechannel EEG and fNIRS signals. The first novel method is a single-stage motion artifacts removal technique where variational mode decomposition (VMD), a robust multiresolution analysis, was employed. The rest two novel methods are VMD with PCA (VMD-PCA) and VMD with CCA (VMD-CCA). As the name suggests, both VMD-PCA and VMD-CCA are two-stage motion artifacts removal techniques. Since, both PCA and CCA algorithms can only process multi-channel signals, the VMD technique was used to generate multiple intrinsic mode functions (IMFs) and then fed to PCA and CCA algorithms for further processing. In all these three methods, single-channel EEG and fNIRS signals were decomposed into 5,10 and 15 IMFs separately and investigated resulting in nine different approaches. VOLUME 10, 2022 It should also be mentioned here that the VMD algorithm was used in removing motion artifacts from motion corrupted PPG signal [54] in our previous work which showed the potential of using VMD for these bio-signals. VMD in combination with PCA and CCA is also the novel contribution of this research work which reduces motion artifacts from EEG and fNIRS signals.
The remainder of this paper is organized as follows: Section II discusses the theoretical backgrounds of the different algorithms (VMD, VMD-PCA, VMD-CCA) investigated here, Section III provides brief information about the EEG and fNIRS benchmark dataset, and experimental methodology. Section IV provides the results of the artifact removal techniques and discusses the results. Finally, the paper is concluded in section V.

II. THEORETICAL BACKGROUND
In this section, the theoretical backgrounds of the VMD, PCA, CCA, VMD-PCA, and VMD-CCA are introduced.
A. VMD VMD [55] non-recursively decomposes time-series signal, X into k number of quasi-orthogonal IMFs which are compact around a center frequency ω k with limited bandwidth. The VMD algorithm, to assess the bandwidth of the time-series signal, is briefly described below: i. A unilateral frequency spectrum is generated for each mode x k using the Hilbert transform. The related analytic signal is calculated. ii. For each mode, the frequency spectrum of the mode is shifted to the baseband by mixing with an exponential which is tuned to the respective estimated center frequency. iii. Through the Gaussian smoothness of the demodulated signal, the bandwidth is estimated. The modes obtained by VMD show less instantaneous frequency fluctuations with better tone detection, tone separation, noise robustness as well as superior signal reconstruction capability in comparison to EMD and EEMD [55]. Let x 0 be the observed signal of the actual signal x contaminated by additive zero-mean Gaussian noise η and is given by: To solve this denoising problem, Tikhonov Regularization is used [32], [33], as follows: from which Euler-Lagrange equations could be obtained and solved in Fourier domain as: where,x (ω) → F {X (.)} (ω) → 1/ 2π R f (t)e jωt dt is the Fourier transform of the signal x(t). The solution corresponds to a convolution with Weiner filter, with α as the variance of the white noise and 1/ω 2 as lowpass power spectrum of the signal. VMD algorithm is the generalized form of Weiner filter with adaptive and multiple band methods and the band-limited IMFs are obtained by solving a constrained variational problem described by [55]: where x k for k = 1, 2, 3 . . . K are the band-limited IMFs having the center frequency ω k obtained by decomposition, with K defined a priori.
The constrain of Equation 4 is addressed by using a quadratic penalty term and Langrangian multiplier (λ). Hence, with δ t as Dirac distribution and ( * ) as convolution operator, the augmented Langrangian equation is given by [56]: where α is the balancing parameter for the data fidelity constraint. Alternate Direction Method of Multipliers (ADMM) [57] is utilized to estimate the saddle point of Equation 5 corresponding to the solution of Equation 4, with convergence tolerance of ε, where the convergence criterion is formulated by, The k th mode estimate is updated using Equation 7 and Equation 8 defined by [55]: PCA [58] is one of the most popular blind source separation (BSS) techniques that are being used for signal processing for the past several decades. PCA transforms the signal time-course among all N channels into a set of N uncorrelated variables using a linear orthogonal transformation. A simple way of computing PCA of a matrix A is to compute the 29762 VOLUME 10, 2022 eigenvalue decomposition of its covariance matrix. The principal components of the input matrix A are the eigenvectors associated with the largest eigenvalues. Since the eigenvalues in the diagonal matrix are sorted in decreasing order, the first d vectors in the matrix V (a matrix whose columns represents the corresponding eigenvectors) are the principal components of A : The use of VMD in combination with PCA for source separation from single-channel measurements is introduced in this study for the very first time to reduce motion artifacts from EEG and fNIRS modalities. As mentioned earlier, the VMD algorithm can be used to generate a multi-channel signal X, comprised of K number of IMFs, from a single channel recording x. This matrix X, can then be fed as the input to the PCA algorithm to produce principal components (as described in the previous subsection) matrix, Y. The columns of matrix Y which contribute to motion artifacts can then be selected as motion artifact components and set to zero. Finally, the artifact-free signal can be reconstructed by simply adding up the rest of the columns of matrix Y.

D. CCA
CCA [59] and Independent component analysis (ICA) [60] are the two most popular BSS methods for separating several contaminated signals. Assuming linear mixing, square mixing, and stationary mixing [61], the ICA technique computes an un-mixing matrix W which is used to find out the unknown independent componentsŜ.Ŝ = WX (9) where X is a matrix of the recorded multi-channel signals. CCA also estimates the unknown independent componentsŜ using Equation 9, but it is different in comparison with the ICA technique because CCA uses second-order statistics (SOS) to generate components whereas ICA uses higherorder statistics (HOS). Hence CCA is also less computationally complex in comparison to ICA. CCA solves the BSS problem by forcing the sources to be maximally autocorrelated and mutually uncorrelated [62]. Given an input signal x, let y be a linear combination of neighboring samples (i.e. y (t) = x (t − 1) + x(t + 1)) [63]). Consider the linear combinations of the components in x and y, called the canonical variates: CCA computes the weight matrices w x and w y which will maximize the correlation ρ between x and y [63]: where C xy is the between-sets covariance matrix and C xx and C yy are the nonsingular within-set covariance matrices. The maximum of ρ is calculated by setting the derivatives of Equation 12 to zero for w x and w y .
w x and w y can then be determined as the eigenvectors of the matrices C −1 xx C xy C −1 yy C T yx and C −1 yy C yx C −1 xx C T xy respectively and the corresponding eigenvalues ρ are the squared canonical correlations. The first pair of variates are the eigenvectors of w x and w y that correspond to the largest square correlation coefficient (or eigenvalue) ρ 2 . The following pairs of variates (w x , w y ) 2...m (with m recording sites) are then the remaining eigenvectors in decreasing order of correlation which are themselves uncorrelated with the previous eigenvectors. The CCA technique, therefore, creates a weight Matrix W x = [w x1 , w x2 , . . . , w xm ] that can be used to separate the recorded sources into the self-correlated and mutually uncorrelated sources.
Using the CCA algorithm, the un-mixing matrix is determined, and the underlying source signalsŜ can be estimated. The components seem to be artifacts that can be removed by setting the corresponding columns of theŜ matrix to zero. When the artifact-suppressed source signals are passed through the inverse of the mixing matrix W −1 , the resultant signal is the artifacts-free signal.

E. VMD-CCA
The single-channel signal x is converted into a multi-channel signal X using the VMD algorithm. This matrix X can then be fed as the input to the CCA algorithm to estimate the underlying true sourcesŜ (Equation 9). The individual sources determined to be artifacts are selected and the corresponding columns of the matrixŜ are set to zero. The source matrix is then passed through the inverse of the un-mixing matrix W −1 to return the multi-channel signalsX which are now, ideally, free of artifacts. The original single-channel signalx can be determined by simply adding the columns in the matrixX.

III. METHODS
This section discusses data description, pre-processing, study design, motion component identification, and evaluation metrics.

A. DATASET DESCRIPTION
In this research work, a dataset, publicly available in the PhysioNet [33], [34], [64] is used which contains ground truth and motion corrupted signals for both EEG and fNIRS modalities. The details of the data recording methodology for EEG and fNIRS modalities can be found in [47] where all the EEG and fNIRS recordings are of 9 minutes of duration.
The proposed denoising methods used a pair of EEG signals which are highly correlated as they were recorded from a very closely spaced location on the scalp. VOLUME 10, 2022 Twenty-three EEG recordings sampled at 2048 Hz, collected from six patients in four different sessions, are available in the publicly available database. Each recording contains two EEG signals: 1) motion corrupted EEG signal and 2) reference ''ground truth'' EEG signal. The average correlation coefficient is very high during the epochs when the motion artifacts are absent [33]. But, when the EEG signal has a motion artifact, the average correlation coefficient drops significantly. The overlaid reference EEG signal and motion corrupted EEG signal are depicted in Figure 1a. The reference EEG signal is utilized for the validation of the proposed denoising methods. As EEG signals can be divided into several sub-bands, namely delta (1-4 Hz), theta (4-8 Hz), alpha (8)(9)(10)(11)(12)(13), beta (13)(14)(15)(16)(17)(18)(19)(20)(21)(22)(23)(24)(25)(26)(27)(28)(29)(30), and gamma (30-80 Hz) [65], in this study, all the EEG recordings are downsampled from 2048 Hz to 256 Hz which ensures data fidelity for further applications while processing the signals becomes relatively convenient. The performance of the proposed techniques is evaluated for all 23 EEG recordings.
Unlike the EEG recordings, the fNIRS recordings contain two pairs of signals, recorded at 690 nm and 830 nm wavelengths. There were 16 fNIRS recordings (9 recordings at 830 nm wavelength and 7 recordings at 690 nm wavelength) in total from 10 test subjects at a sampling frequency of 25 Hz [34], [64]. The performance of the proposed denoising techniques is evaluated for all the 16 fNIRS recordings. The overlaid reference fNIRS signal and motion artifact contaminated fNIRS signal are depicted in Figure 1b.

B. SIGNAL PREPROCESSING 1) POWER LINE NOISE REMOVAL
Power line noise causes artifacts in physiological signals. The recorded EEG and fNIRS signals are no exception to this. To remove power line noise, a 3 rd order Butterworth notch filter with a center frequency of 50 Hz was chosen to remove 50 Hz and its subsequent harmonics as a pre-processing technique for all the EEG and fNIRS signals.

2) BASELINE DRIFT CORRECTION
It is observed that both the EEG and fNIRS signals suffered greatly from baseline drift i.e. unwanted amplitude shifts in the signal which would consequently give incorrect outcomes if not removed. To eliminate baseline drift from EEG and fNIRS recordings, a polynomial curve fitting algorithm (e.g., profit function of MATLAB) was utilized to approximate the baseline and then subtracted from the corresponding original signal.

C. STUDY DESIGN
The simulations of this study were carried out in a PC with Intel(R) Core (TM) i5-8250U CPU at 1.80GHz equipped with 8 GB RAM. In-house built MATLAB code was designed to pre-process the EEG and fNIRS data and three different multiresolution analysis (MRA) techniques: one stage (VMD) and two-stage (VMD-PCA and VMD-CCA) were deployed in ''MATLAB R2020a, The MathWorks, Inc''. Figure 2 outlines the motion artifacts removal framework which is proposed in this research work. An automated way for identifying motion corrupted components of the preprocessed signal is also discussed.
As mentioned earlier, VMD can generate k number of IMFs where k is user-defined. In this research work, three different values of k have been chosen namely 5, 10, and 15 separately. In the rest of the manuscript, the subscript VMD (5) would refer to the generation of 5 IMFs, while, the subscripts of VMD (10) and VMD (15) would refer to the generation of 10 and 15 IMFs, respectively. With the IMF's availability using the VMD technique, the artifact components can be selected and removed. All the remaining IMFs can then either be added up to reconstruct a cleaner signal or all the IMFs can be fed as inputs to PCA /CCA algorithms separately to generate the output PCA/CCA components. PCA technique needs the number of input channels to be at least two or greater. In this work, single-channel EEG and fNIRS signals are evaluated for the correction of motion artifacts. Hence, it is required to generate several IMFs which would be used as the inputs for PCA and CCA separately. VMD-PCA-based (VMD (5) -PCA, VMD (10) -PCA, VMD (15) -PCA) two-stage artifacts removal technique has been realized for three different values of k = 5, 10, and 15 for both single-channel EEG and fNIRS modalities. Also, VMD-CCA-based (VMD (5) -CCA, VMD (10) -CCA, VMD (15) -CCA) two-stage artifacts removal technique has been realized for three different values of k = 5, 10, and 15 for both single-channel EEG and fNIRS signals.

D. REMOVAL OF MOTION ARTIFACT COMPONENTS USING ''GROUND TRUTH'' METHOD
A frequently encountered issue in removing motion artifacts using the aforementioned artifact removal techniques is to reliably identify the motion corrupted components of the decomposed signal, remove those components, and reconstruct a cleaner signal using the remaining components. To evaluate the effectiveness of the proposed algorithms, the available reference ''ground truth'' signal of EEG and fNIRS modalities were utilized. If a component of the decomposed signal is removed and the signal is reconstructed using the remaining components, then the correlation coefficient between the newly reconstructed signal and the reference ''ground truth'' signal will only increase if that removed component suffers from motion artifacts. Using this simple yet effective idea, component/s of the decomposed signal suffering from motion artifacts were identified and removed to produce a cleaner signal which eventually ensured the optimal performance of each proposed technique during evaluation. Figure 3a shows an example motion corrupted EEG signal and Figure 3b represents the corresponding 10 bandlimited IMFs (BLIMFs) generated from the EEG signal using VMD (10) algorithm. Figure 4a represents an example motion corrupted EEG signal from which 5 BLIMFs are generated using VMD (5) and then these 5 IMFs are fed as 5 channel input to the PCA algorithm whereas Figure 4b represents the resultant PCA components. Figure 5a depicts an example motion corrupted EEG signal and Figure 5b represents resultant 5 CCA components where the input of the CCA method was 5 BLIMFs generated from the motion corrupted EEG signal using VMD (5) .
Similarly, Figure 6a, Figure 7a, and Figure 8a shows three separate motion corrupted fNIRS signals whereas Figure 6b, Figure 7b, and Figure 8b represent the BLIMFs generated from VMD (5) , 5 output components from the PCA algorithm, and 5 output components using the CCA algorithm, respectively for the corresponding motion corrupted fNIRS signals.

E. PERFORMANCE METRICS
Using the available reference ''ground truth'' signal for each modality, as described earlier, the efficacy and performance of each proposed artifact removal technique can be determined. Since the objective of each proposed technique is to reduce motion artifacts from the contaminated signal, computation of SNR and percentage reduction in motion artifacts can quantify the efficacy of that corresponding technique's ability in removing artifacts. Therefore, the difference in SNR before and after artifact removal ( SNR), and the improvement in correlation coefficient between signals, denoted by the percentage reduction in motion artifact η [34], are used  as performance metrics. For the calculation of SNR, the following formula was used which was given in [34]: where σ 2 x , σ 2 e before , and σ 2 e after represent the variance of the reference ''ground truth'' signal, motion corrupted signal, and cleaned signal, respectively.
The improvement in correlation coefficient between signals was used to estimate another performance metric, namely the percentage reduction in motion artifact η [34]:  where ρ before denotes the correlation coefficient between the reference ''ground truth'' and motion-corrupted signals and ρ after represents the correlation coefficient between the reference ''ground truth'' and the processed cleaner signal. ρ clean is the correlation between reference ''ground truth'' and motion corrupted signals over those epochs where motion artifact is absent. In this study, a modified version of Equation 15 has been used assuming ρ clean = 1 as in an ideal situation the correlation coefficient between the ground truth and the motion corrupted signal over the clean epochs would always be 1. Hence, the modified version of Equation 15 becomes:

IV. RESULTS AND DISCUSSION
The results obtained in this work, using the proposed novel motion artifact removal techniques are mentioned below where the performance metrics were calculated using Equations 14 and 16.

A. MOTION ARTIFACT CORRECTION FROM EEG
All the methods mentioned earlier were applied to all the 23 recordings of EEG. Figure 9 depicts an example EEG trial after the correction of the motion artifacts using VMD (5) , VMD (5) -PCA, VMD (5) -CCA.

1) VMD
When employing VMD (5) , VMD (10) , and VMD (15) , an average SNR of 25.91 dB, 24.58 dB, and 24.78 dB, respectively was found for all (23) recordings and 25.34 dB, 24.63 dB, and 24.28 dB, respectively was found over 21 recordings. The average percentage reduction in artifact for three separate VMD approaches was found to be 53

B. MOTION ARTIFACT CORRECTION FROM fNIRS
Again VMD, VMD-PCA, and VMD-CCA (with three different decomposition levels of IMFs i.e. 5, 10, and 15) techniques were applied and the performance matrices were calculated for all the 16 fNIRS recordings. Figure 10a, Figure 10b, and Figure 10c illustrate an example fNIRS trial after the removal of the motion artifact using VMD (10) , VMD (10) -PCA, VMD (10) -CCA, respectively.

1) VMD
For all trials, the average SNR improved to 15.63 dB, 15.39 dB, and 15.45 dB when VMD (5) , VMD (10) , and VMD (15) techniques were applied, respectively, and the average percent reduction in motion artifacts improved by 28 Table 1 summarizes the results obtained (average SNR and average percentage reduction in motion artifacts η) using the artifact removal techniques described in the paper i.e. VMD, VMD-PCA, and VMD-CCA for all the EEG (23) and fNIRS (16) recordings. The values inside brackets in Table 1 denote the corresponding VOLUME 10, 2022   (15) , VMD (15) -PCA, and VMD (15) -CCA for all 16 cleaned fNIRS signals. standard deviations. For both the EEG and fNIRS modalities, it is evident from the result presented in Table 1 that the two-stage artifacts removal technique i.e. VMD-PCA and VMD-CCA algorithms performed relatively better compared to the single-stage artifact removal technique i.e. VMD.  (10) , VMD (10) -PCA, and VMD (10) -CCA, and (c) VMD (15) , VMD (15) -PCA, and VMD (15) -CCA for 21 motion-corrected EEG signals. Figure 11a shows the box and whisker plot of the performance metrics ( SNR and percentage reduction in motion artifacts) using VMD (5), VMD (5) -PCA, and VMD (5) -CCA techniques for all the motion-corrected EEG signals. Figure 11b represents the box and whisker plot of the same performance metrics when VMD (10), VMD (10) -PCA, and VMD (10) -CCA algorithms were employed whereas Figure 11c depicts the performance metrics obtained using box and whisker plot from VMD (15), VMD (15) -PCA, and VMD (15) -CCA for all the EEG recordings.
It is evident from the result of Table 1 and Figure 11 that the EEG signal reconstructed from VMD (5) method has a greater SNR value (25.91 dB) and lesser η value (53.59%) compared to both VMD (10) and VMD (15) techniques. VMD (15) provided the largest η value (55.86%) while VMD (5) rendered the largest SNR (25.91 dB) among the three different single-stage motion artifact correction techniques (VMD (5), VMD (10), and VMD (15) ). The probable reason for the greater improvement in correlation is as VMD (15) generates 15 IMFs and VMD (5) produces 5 IMFs, there is a lesser probability that the last IMF which is considered as motion component for the VMD (15) has a lesser chance of getting mixed up with the signal components and movement artifacts components. In other words, VMD (15) generates finer decomposition of the single-channel EEG signal compared to VMD (5) .
When two-stage motion artifacts removal techniques are employed (VMD-PCA and VMD-CCA), the best average correlation improvement is produced by VMD (15) -CCA approach which is 57.01%. On the other hand, the best average SNR value (24.86 dB) is obtained from VMD (10) -CCA technique for EEG signals (Table 1, Figure 11b and Figure 11c). Since two-stage decomposition allows the low IMF number-based techniques to be further decomposed to extract noise components, VMD (10) -CCA performed well for the two-stage technique as well while the best correlation improvement is produced by VMD (15) -CCA approach. Figure 12a, Figure 12b, and Figure 12c represent the box and whisker plot of the SNR and percentage  Table 1 and Figure 12, it is observed that for fNIRS recordings, the best average SNR value is 15.97 dB and the greatest average percentage reduction in movement artifacts is 39.01% produced by VMD (10) -CCA. Also, twostage motion artifacts correction techniques (VMD-PCA and VMD-CCA) performed significantly better compared to the single-stage motion artifacts correction technique (VMD) in terms of the average percentage reduction in motion artifacts.
Authors of [38] found that brain activity was absent in the EEG trials of 12 and 15. Moreover, they found a poor correlation coefficient over the clean epochs of the recordings of trials 12 and 15 and carried out their investigation on the remaining 21 recordings of EEG. We have observed similar phenomena while investigating this study. In addition to this, similar phenomena were observed for the single-channel fNIRS signal of trial 8 (recorded at 830 nm wavelength). Therefore, the authors of this work are presenting a second table ( Table 2) that depicts the average SNR and average percent reduction in motion artifacts for the remaining 21 EEG signals (excluding trials 12 and 15) and 15 fNIRS signals (excluding trial 8 recorded at 830 nm wavelength) using VMD, VMD-PCA, and VMD-CCA based algorithms. Figure 13a, Figure 13b, and Figure 13c represent the box and whisker plot of the calculated performance parameters ( SNR and η) using Equations 14 and 16) when VMD (5) , VMD (5) -PCA, VMD (5) -CCA; VMD (10) , VMD (10) -PCA, and VMD (10) -CCA; VMD (15) , VMD (15) -PCA, and VMD (15) -CCA algorithms are employed, respectively for 21 motion corrupted EEG recordings. On the other hand, Figure 14a, Figure 14b, and Figure 14c represent the box and whisker plot of the same performance parameters which are obtained from VMD (5) , VMD (5) -PCA, VMD (5) -CCA; VMD (10) , VMD (10) -PCA, and VMD (10) -CCA; VMD (15) , VMD (15) -PCA, and VMD (15) -CCA techniques respectively for 15 motion artifact corrected fNIRS signals. Box and whisker plot of performance parameters using VMD (5) for 21 cleaned EEG trials. Here, for each trial, the first four IMFs were added to generate the clean EEG signals directly.
After excluding two recordings (trials 12 and 15) of EEG signals, the largest average SNR was found as 25.70 dB using VMD (5) -CCA and VMD (15) -CCA technique provided the highest average η which is 62.09%. For fNIRS signals, the highest average SNR (17.01 dB) was found when VMD (15) -CCA technique was incorporated. The largest η (39.48%) was generated by VMD (10) -CCA method.
As previously mentioned in the introduction section, it has been found that DWT, EMD, EEMD, EMD-ICA, EMD-CCA, EEMD-ICA, EEMD-CCA, SSA, Wavelet-based techniques along with approximation sub-band filtering, adaptive filtering (ARX model with exogenous input), etc. were mostly used as movement artifacts removal techniques. To obtain denoised signals from motion corrupted physiological signals using DWT-based techniques, choosing the appropriate wavelet is very vital and somewhat difficult as well. To date there is no concrete rule available for choosing the right wavelet for the specific physiological signal of interest, rather, in most cases, wavelets are chosen based on the morphology of the signal of interest. Hence, inappropriate selection of wavelets would lead to inefficient denoising. EMD based motion artifacts removal technique greatly suffers from the 'mode mixing' problem which would lead to an erroneous result. This mode mixing problem is solved in EEMD techniques. Although EEMD is free from the mode mixing problem, this algorithm requires an optimum parameter, the number of ensembles to be used. Selecting the number of ensembles is entirely done based on trial and error [34]. To use the SSA algorithm for the analysis of physiological signals, a prior declaration of the window length and the required number of reconstruction components is necessary, which makes SSA inefficient as well. In [42], the authors used DWT and approximation sub-band filtering using total variation (TV) and weighted TV. During the reconstruction of the signal, they have discarded the first three high-frequency detailed sub-band signals based on that those first three detailed sub-band signals contained no relevant information of EEG signal. Therefore, identifying the non-useful components while using DWT-based methods is critical to remove motion artifacts from EEG as well as fNIRS signal. Also, the value of the regularization factor to solve the optimization problem of TV and MTV techniques was chosen without giving any justification. In [52], the autoregressive exogenous input model (adaptive approach) was investigated to model the motion corrupted segments as output, and IMU data was used as exogenous input. The authors showed the efficacy of their proposed model only for 4 test subjects. One of the crucial parts presents in implementing this method is the proper synchronization of the fNIRS data and IMU data which is very challenging in a practical environment. Also, if the epoch duration of the motion artifacts is sufficiently large (more specifically, the sample size), then modeling the artifacts mathematically using the least square method would require higher-order polynomial models (stability issue), which makes this approach very difficult to implement in a real-life scenario.
The two-stage motion artifacts removal techniques (VMD-PCA and VMD-CCA) proposed in this paper will not be able to identify the motion corrupted PCA components/CCA components in the absence of ground truth signal. In the absence of ground truth signal, an alternative method can be realized proposed by Hassan et al. [66] where the authors used the autocorrelation function to identify the motion corrupted components. The automatic artifact component selection technique using the autocorrelation function proposed in [66] has not been presented in this work.
However, the proposed single-stage motion artifact removal technique (VMD) will give optimum results even in the absence of the ground truth signal. It has been visually observed that the maximum percentage of motion artifacts are present in the last IMF of VMD. To validate this statement, VMD (5) technique was again used for 21 noisy EEG (trials 12 and 15 were excluded) signals, where the first four IMFs were added up leaving the 5 th IMFs as it is. This technique produced an average SNR of 23.04 ± 6.22 dB and a 58.31% reduction in artifacts with a standard deviation of 25.16. Figure 15 depicts the box and whisker plot for the performance parameter obtained using this technique.
Applications where obtaining ground truth signal is difficult (EMG signal acquisition), could use this method for successful elimination of motion artifacts. On the other hand, in the physiological signal acquisition modalities, where obtaining the ground truth signal is relatively easier (ECG signal), making use of the two-stage motion artifacts removal technique (VMD-CCA) is suggested for obtaining greater accuracy and signal quality.
The results, as well as the box and whisker plots presented in this research work, are a clear indication of the efficacy of our proposed novel multiresolution techniques VOLUME 10, 2022 i.e. VMD, VMD-PCA, VMD-CCA for the successful removal of motion artifacts from EEG and fNIRS signals. Since all the techniques, proposed in this work, are reliable, robust, and computationally efficient, similar methodology with minimal alterations can be applied for the correction of motion artifacts from other physiological signals (EMG, ECG, PPG, etc.) as well. We strongly believe, the motionartifacts-free EEG and fNIRS signals obtained using our proposed models would facilitate in identifying neurological disorders by both medical doctors and machine learning based systems to a great extent. The only limitation of this study is that the two-stage motion artifact correction methods (VMD-PCA and VMD-CCA), proposed in this study, are relatively computationally expensive in comparison with the single-stage motion artifact correction method (VMD).

V. CONCLUSION
Three robust motion artifact removal techniques have been proposed in this paper, namely variational mode decomposition (VMD), variational mode decomposition in combination with principal component analysis (VMD-PCA), and variational mode decomposition in combination with canonical correlation analysis (VMDD-CCA) for EEG and fNIRS modalities. Further, the proposed algorithms are investigated by 9 different approaches. Both VMD-PCA and VMD-CCA techniques can be used on the single-channel recordings as the VMD algorithm can decompose a single-channel signal into a predefined number of IMFs that can be fed to the PCA/CCA algorithm. The performance parameters obtained from all these approaches are a clear indication of the efficacy of these algorithms. The novel VMD (15) -CCA and VMD (10) -CCA technique provided the best performance in terms of the percentage reduction in motion artifacts (62.09% and 39.48%) when analyzing the 21 EEG and 15 fNIRS recordings, respectively. On the other hand, the VMD (10) -CCA technique generated the highest average SNR (25.35 dB) for EEG signals, and VMD (15) -CCA produced the highest SNR (17.01 dB) for fNIRS signals. An alternative approach for removing motion artifacts from EEG signals has also been investigated that produced very good results. This claim is validated from the computation of the performance matrices also. In the future, deep learning models will be investigated for the automated detection and removal of motion artifacts in physiological signals (EEG, ECG, EMG, PPG, fNIRS, etc.). New methods based on the use of different multivariate signal processing approaches will be developed for the elimination of other artifacts from the EEG and fNIRS signals recorded using multiple electrodes.

ETHICAL STATEMENT
No ethical statement to be declared

ACKNOWLEDGMENT
The dataset used in this experiment is kindly shared in the PhysioNet database by Sweeney et al. [33], [34], [64]. The statements made herein are solely the responsibility of the authors.
SERKAN KIRANYAZ (Senior Member, IEEE) joined Qatar University, in 2015, as a Professor with the Department of Electrical Engineering. He was a Professor with the Department of Signal Processing, Tampere University of Technology, from 2010 to 2015. He has published two books, three book chapters, more than 60 journal articles in many IEEE TRANSACTIONS and other high-impact journals, and around 100 papers in international conferences. His principal research field is ''intelligent systems'' which comprises systems, devices, algorithms, software, and processes requiring intelligent processing of possibly ''big'' data and sensor signals together with artificial reasoning. He made significant contributions to bio-signal analysis, particularly ECG analysis and processing, classification and segmentation, signal processing and management, computer vision, scene analysis, automatic object segmentation, and machine learning strategies with applications to recognition, classification, detection and evolutionary machine learning, swarm intelligence, stochastic optimization, and deep learning. His team is among the pioneers in patient-specific ECG classification. His early work became one of the most-read articles in IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING and has currently more than 400 citations according to Google Scholar. His most recent work on real-time patient-specific ECG classification too became one of the most popular articles in the IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING. His team has won 2nd place in the PhysioNet Challenge Grand Competition (see http://www.cinc.org/physionet-cinc-challenge-awards/) where 48 international teams worldwide competed to achieve the most accurate systems for the classification of normal/abnormal heart sound recordings.
AMITH KHANDAKAR (Senior Member, IEEE) received the B.Sc. degree in electronics and telecommunication engineering from North South University, Bangladesh, and the master's degree in computing (networking concentration) from Qatar University, in 2014. He graduated as the Valedictorian (President Gold Medal Recipient) of North South University. After graduation, he was a Consultant in a reputed insurance company in Qatar and in a private company that is a sub-contractor to national telecom service provider in Qatar. He is a certified Project Management Professional and the Cisco Certified Network Administrator. He was a Teaching Assistant and a Laboratory Instructor for two years for courses, such as mobile and wireless communication systems, principle of digital communications, introduction to communication, calculus and analytical geometry, and Verilog HDL: modeling, simulation, and synthesis. Simultaneously, he was a Laboratory Instructor for the following courses: programming course ''C,'' Verilog HDL, and general physics course. He has been with Qatar University, since 2010. He is currently the General Secretary of the IEEE Qatar Section and also the Qatar University IEEE Student Branch Coordinator and an Adviser (Faculty Member).

MOHAMMED ALHATOU is working as a Senior
Consultant with the Department of Neurology, Al-Khor, Hamad General Hospital, from August 2015 to present. After completing his neurology training at the University of Virginia, he did a fellowship in clinical neurophysiology and electromyography at Emory University. He built the Neurophysiology Laboratory at the Neurology and Pain Clinic and directed the Clinical Neurophysiology Laboratory at the Regional Medical Center of Orangeburg and Calhoun Counties, from 2001 to 2015. He is an Assistant Professor at V-COM. He is board-certified in neurology, clinical neurophysiology, clinical electromyography sleep medicine, neuroimaging, and neuromuscular disorders. He participated as a PI in multiple clinical trials in neuropathy, epilepsy, and dementia supported by Pfizer, Ortho-McNeil, UCB, and King pharmaceutical. He has multiple peer-reviewed papers, journal articles, and abstracts. He has also served as the Head of the Credential Committee and a member of the Medical Executive and Peer Review Committee.
RUMANA HABIB received the M.B.B.S. and F.C.P.S. (neurology) degrees. She is a Neurology Specialist in Dhaka, Bangladesh. She is working as a Neurology Specialist Doctor at the BIRDEM General Hospital and Ibrahim Medical College, Dhaka. She has published several journal articles and abstracts in renowned peer-reviewed journals. Her research interests include autonomic, peripheral, proximal, focal neuropathy, and processing and interpreting 1D bio-signals. VOLUME 10, 2022