# Cortical Tasks-Based Optimal Filter Selection: An fNIRS Study

Authors: Rayyan Azam Khan; Noman Naseer; Sajid Saleem; Nauman Khalid Qureshi; Farzan Majeed Noori; Muhammad Jawad Khan
Journal: Journal of Healthcare Engineering
 (2020)

---

## Abstract

Functional near-infrared spectroscopy (fNIRS) is one of the latest noninvasive brain function measuring technique that has been used for the purpose of brain-computer interfacing (BCI). In this paper, we compare and analyze the effect of six most commonly used filtering techniques (i.e., Gaussian, Butterworth, Kalman, hemodynamic response filter (hrf), Wiener, and finite impulse response) on classification accuracies of fNIRS-BCI. To conclude with the best optimal filter for a specific cortical task owing to a specific cortical region, we divided our experimental tasks according to the three main cortical regions: prefrontal, motor, and visual cortex. Three different experiments were performed for prefrontal and motor execution tasks while one for visual stimuli. The tasks performed for prefrontal include rest (R) vs mental arithmetic (MA), R vs object rotation (OB), and OB vs MA. Similarly, for motor execution, R vs left finger tapping (LFT), R vs right finger tapping (RFT), and LFT vs RFT. Likewise, for the visual cortex, R vs visual stimuli (VS) task. These experiments were performed for ten trials with five subjects. For consistency among extracted data, six statistical features were evaluated using oxygenated hemoglobin, namely, slope, mean, peak, kurtosis, skewness, and variance. Combination of these six features was used to classify data by the nonlinear support vector machine (SVM). The classification accuracies obtained from SVM by using hrf and Gaussian were significantly higher for R vs MA, R vs OB, R vs RFT, and R vs VS and Wiener filter for OB vs MA. Similarly, for R vs LFT and LFT vs RFT, hrf was found to be significantp<0.05. These results show the feasibility of using hrf for effective removal of noises from fNIRS data.


---

## Body

## 1. Introduction

Brain-computer interface (BCI) also known as human-machine interface (HMI) or brain-machine interface (BMI) provides a communication mean between the user and external devices through a combination of hardware and software systems [1–3]. These systems are trained to generate control commands based on a specific set of patterns of brain signals [4].Brain signal acquisition is categorized between invasive and noninvasive techniques. However, due to surgical risks and limited access to the cortical region, noninvasive techniques are common in practice [5]. Noninvasive modalities include functional magnetic resonance interference (fMRI), functional near-infrared spectroscopy (fNIRS) [6], and electroencephalography (EEG) [4]. fNIRS is a comparatively new modality that has better spatial resolution and low artifacts, cost, and portability [4, 7]. So far, promising results have been shown by fNIRS-BCI [8–10]. Acquired brain signals for a specific task may contain noises that can contaminate signals and can effect informative data. These noises are categorized between physiological noise, experimental noise, and instrumental noise [8, 11]. In fNIRS, the experimental noises are removed prior to the change of the raw signal to its magnitude through the modified Beer–Lambert law [8, 12, 13]. Noises produced due to hardware or by surrounding are known to be instrumental noises. These noises usually have high frequency that can be removed using low-pass filter; furthermore, keeping isolation from external sources such as light can reduce such type of noises. Experimental noises include motion artifacts such as head motion during signal acquisition that can cause the dislocation of optodes from the assigned position, thus generates a spike-like noise due to change in light intensity. Various studies [13–15] have utilized commonly developed filtering techniques randomly for noise removal. However, obtained signals can be corrupted by different kinds of noises that could affect further analysis. Noises can be physiologically produced due to Mayer waves (∼0.1 Hz), respiration (0.2∼0.5 Hz), and heartbeat (1∼1.5 Hz), mainly due to fluctuations of blood pressure [8, 16–19]. These noises can be removed using adaptive or bandpass filtering [20, 21]. After preprocessing, useful information is extracted from the filtered data afterward classified using different classifiers mainly named as linear discriminant analysis (LDA), support vector machine (SVM), quadratic discriminant analysis (QDA), and naïve bayes (NB) [22] to generate control commands, hence completing the loop for BCI. Previous studies [13, 14, 23] show that an appropriate filter for correctness of data is the key to achieve more accurate results.In this study, we hypothesized to find an optimal filter for commonly used cortical tasks, owing to a particular cortical region. Hence, we compared six commonly used filters to remove previously discussed noises. These filters include discrete Kalman [24], time-varying Wiener [25], 4th order Butterworth, hemodynamic response filter (hrf), Gaussian [26], and window-based finite impulse response (FIR) [27]. For the said purpose, cortical data were acquired from the three main regions of the brain, namely, prefrontal (PFC), motor (MC), and visual cortex (VC). Since the data acquired from PFC relate to thinking tasks [28–31], hence arithmetic and object rotation tasks were performed for this cortical region [32–34]. Similarly, tasks related to movement of limbs or fingers is related to the motor cortex [35]; therefore, the finger tapping tasks were performed for desired data acquisition [36–38]. Likewise, flickering of checker box was performed for visual cortex data [39]. Keeping in view the target of an optimal filter for a specific cortical region, a previous study [22] reported the different combinations of statistical features. Therefore, for consistency of extracted data, statistical features were kept the same for all experimental tasks, hence making combinations of six features, namely, signal mean (SM), signal slope (SS), signal peak (SP), signal skewness (SK), signal kurtosis (KR), and signal variance (SV). For classification, a number of studies [22, 40, 41] reported nonlinear SVM classifiers for comparatively better accuracies, hence all experimental tasks were classified using a nonlinear SVM classifier. Therefore, the main contribution of this work is (1) to analyze the effect of the six most commonly used filtering techniques and (2) to propose the optimal one among the most frequently discussed noise removal techniques. For the aforementioned experiment, we select three main cortical regions: prefrontal, motor, and visual cortex with seven various paradigms. The canonical hemodynamic response filter (hrf) [26] performed overall best among the opted techniques. On the basis of these systematic and explicit analyses, it can be seen that selection of an optimal filter has a significant role in enhancing accuracies. Hence, these observations can serve as a standard guide for others to test the effect of noise correction algorithms for fNIRS experiments, and therefore can select a significant methodology.


## 2. Materials and Methods

### 2.1. Experimental Setup

To acquire experimental data, seven paradigms owing to three main cortical regions were designed and explicitly performed using a dynamic near-infrared optical tomography (DYNOT-232; NIRx Medical Technologies, NY, USA) device at Pusan National University. It operates on two wavelengths that are 760 and 830 nm where the signal acquisition sampling frequency was 1.81 Hz. Five healthy subjects with normal or corrected-to-normal vision took part in the experiment with a verbal consent before experimentation. All subjects were right-handed with an age range of 26 ± 3. Right-handed subjects were selected to minimize hemodynamic response variation due to hemispheric-dominance difference. The experimental participants had no history of alcoholism, psychiatric, neurological, and visual disorder, cardiovascular and respiratory disease, mental illness, or any motor disability. Moreover, three hours before the commencement of the study, participants were asked to refrain from caffeinated drinks. As discussed in the literature [28, 30, 31, 35, 37, 38], for thinking-related task, signals were acquired from PFC, similarly for motor execution tasks from the primary motor cortex (PMC), and visual stimuli task from VC. Performed experiments were according to the latest Declaration of Helsinki.


### 2.2. Experimental Paradigms

In accordance with the literature [35, 39], subjects were seated on a comfortable chair and were asked to take rest with restricted movements as they can, so that the hemodynamic response activation owing to previous activities can be avoided. Hence, each paradigm related to PFC, MC, and VC starts with a rest of 20 s period to set up the baseline conditions. As the literature [42–45] show, (10∼12)  s task is adequate to acquire hemodynamic response of brain activity, hence 20 s initial rest was followed by 10 s task, and this was followed in turn by another 20 s rest period permitting signals to return to their baseline values before the start of the next trial in paradigm (a). For paradigm (b), 20 s rest after 10 s task 1 was again followed by 10 s task 2. The 20 s rest between two 10 s tasks was added to differentiate two classes through the baseline value. Figure 1 depicts paradigm (a) and paradigm (b). For optimal filter selection, tasks were selected concerning specific cortical regions such that for PFC, three different experiments were performed, rest (R) vs mental arithmetic (MA), R vs object rotation (OB), and OB vs MA task. Similarly, for MC, R vs left finger tapping (LFT), R vs right finger tapping (RFT), and LFT vs RFT. Likewise, for VC, R vs visual stimuli (VS) task.Figure 1
Experimental paradigms for cortical tasks. (a) Task vs rest. (b) Task 1 vs task 2.
(a)(b)


#### 2.2.1. Prefrontal Tasks

(1) Mental Arithmetic. According to previous studies [44, 46, 47], for R vs MA task, subjects were asked to perform a series of arithmetic calculations for 10 s based on the pseudorandom order, such that to subtract the random two-digit number (between 10 and 20) successively from the previous result of the three-digit number subtraction appearing on the screen (e.g., 400-11, 389-17, and 372-14). Afterwards, the screen was turned black so that subjects does not go beyond the 10 s task.(2) Object Rotation. For the object rotation task, subjects were asked to imagine a cube rotating for the 10 s task while the “object rotation” word appeared for 10 s on the screen at about 2 m distance [48–50].(3) Mental Arithmetic vs Object Rotation. In this protocol, subjects were asked to perform the aforementioned MA task versus OB task in between 20 s of rest to distinguish two tasks. The experimental paradigm is depicted in Figure 1(b).


#### 2.2.2. Motor Tasks

(1) Finger Tapping. According to the literature [8, 17, 19, 23, 31, 37, 38, 51], subjects were asked to tap the self-paced index finger of one hand for 10 s afterwards the 20 s rest task was performed allowing signals to return to their reference values. Also, repetition for 10 times was performed as depicted in Figure 1(a). Similar trials were performed on the other hand, while for the LFT vs RFT task, 20 s rest was performed in between two tasks for the restoration of the signal to the baseline level as shown in Figure 1(b).


#### 2.2.3. Visual Task

(1) Checker Box Flickering. In this experiment [4, 18, 39], a screen was placed in front of the subjects at a distance of approximately 2 m, and also subjects were requested to avoid eye blinking during the experiment. The 10 s task of checker box flickering at 4 Hz was performed followed by 20 s rest of the black screen. The sound was also generated during the transition between rest and task. The paradigm followed for visual stimuli is shown in Figure 1(a).


### 2.3. Experimental Setup

Since the mental imagery task activates the PFC [34, 43], a total 11 of near-infrared (NI) light optodes were placed on PFC, 3 of which were detectors and 8 were the source in accordance with the literature [46, 47]. Similarly, for the motor execution task, the primary motor cortex (PMC) is activated [37, 43], hence 15 optodes were placed on PMC out of which 8 were the source and 7 were detectors. To extract data for LFT, optodes were placed on the right hemisphere, while for RFT on the left hemisphere [43]. Similarly, for the visual stimuli task data acquisition from the visual cortex [4, 18, 39], eleven optodes were placed having eight sources and three detectors. The distance between the source and the detector was 3 cm. Optode placement with channel configuration for MC, VC, and PFC is shown in Figure 2.Figure 2
Optode placement configuration in accordance with the 10–20 international system. (a) Motor cortex. (b) Visual cortex. (c) Prefrontal cortex.
(a)(b)(c)


#### 2.3.1. Signal Acquisition

In accordance with the literature [6, 8, 29], the raw optical density signal is converted to oxyhemoglobin (ΔcHbOt) and deoxyhemoglobin (ΔcHbRt) concentration using the modified Beer–Lambert law (MBLL) as decribed in the following equation:(1)ΔcHbOtΔcHbRt=βHbOλ1βHbRλ1βHbOλ2βHbRλ2−1Δψt,λ1Δψt,λ2d∗l,where βHbXλ is the HbX extinction coefficient in µM−1cm−1, d is the differential path length factor for the curved path in (mm), l is the detector and emitter distance in (mm), and ΔψHbXt is the absorbance difference of the light emitter wavelength of λi.


### 2.4. Signal Processing

The acquired raw signals of the brain contain various noises that can be categorized into physiological, experimental, and instrumental noise [8, 52]. In fNIRS, the instrumental and experimental noises are removed prior to change of the raw signal to its magnitude through MBLL [8, 12, 14].


#### 2.4.1. Instrumental Noises

Noises produced due to hardware or by surrounding are categorized under instrumental noises. These noises usually have high frequency that can be removed using a low-pass filter; furthermore, keeping isolation from external sources such as light can reduce such type of noises [8].


#### 2.4.2. Experimental Noises

These errors contain motion artifacts due to unintentional body movements like head motion during signal acquisition. It can cause the dislocation of optodes from the assigned positions that result in a spike-like noise due to change in light intensity. Filters like Kalman and Wiener can be used to remove such type of noises [12–15].


#### 2.4.3. Physiological Noises

These noises are produced due to Mayer waves (∼0.1 Hz), respiration (∼0.5 Hz), and heartbeat (1∼1.5 Hz) that occur due to fluctuations of blood pressure [16–19]. These noises can be removed using adaptive or bandpass filtering [8, 20, 38].


### 2.5. Data Analysis

NIRS-SPM is a toolbox designed for fNIRS data analysis. For signal processing, it provides common filtering techniques, namely, Butterworth, Gaussian, and hrf [26]. The comparative data analysis was performed by implementing Gaussian and hrf filtering using the NIRS-SPM toolbox, while other techniques on MATLAB® 2017b. The generalized mathematical models with necessary details are as follows.


#### 2.5.1. Gaussian Filtering

A Gaussian filter is used in various forms depending upon the nature of the signal. Generally, a Gaussian filter is based on a Gaussian function which defines the probability distribution of noise or data. It can also be used as a smoothing operator. A Gaussian kernel is used for smoothing the signal in which each value is replaced with the weighted average of itself and its neighboring values [14, 21, 26]. A simple representation of the 2D Gaussian filter can be defined in the following equation as(2)Gx,y=12πσ2e−x2+y2/2σ2,where x and y are the distance from the origin in horizontal and vertical axis and σ is the known standard deviation of the distribution.


#### 2.5.2. Hemodynamic Response Filter (hrf)

The hrf is based on the canonical representation of the hemodynamic response functions (HRF) and is used for the temporal smoothing of the fNIRS time series signal. In NIRS-SPM, the given functional data were smoothen using the least square estimate with ideal HRF. The hrf and Gaussian filter model details are in accordance with the literature [26].


#### 2.5.3. Butterworth Filter

Butterworth filter is a model-based bandpass filter which performs on frequency attenuation using high and low-pass filter. The filtered value not only depends on the weighted average of the unfiltered time series, but also recursively on the previous values of the filtered time series. This filter aims to have a flat frequency response in the desired pass band [26]. The 4th order Butterworth filter with a band pass of (0.01∼0.1) Hz was applied by MATLAB® build in a library for the desired experiments.


#### 2.5.4. Finite Impulse Response Filter

FIR filter is designed by finding the coefficients and filter order so that it performs a cross-correlation between the input signal and the time reversed impulse response; therefore, by sampling the pulse shape, coefficients of the filter are designed [14]. FIR filter of order N can be defined as the following equation:(3)xn=∑i=0Naiyn−i,where yn is the input signal, xn is the output signal, N is the filter order, and bi is the value of impulse response at ith instance. Here, the 4th order FIR filter with a low-pass band of 0.1 Hz was utilized. The coefficients were estimated using the least square estimate. However, a time-varying Wiener filter, based on the short-time Fourier series, was implemented as in [25, 53].


#### 2.5.5. Kalman Filter

Kalman filters the input signal containing statistical and other noises by linear quadratic estimation. The estimated unknown variables based on the Bayesian inference and joint probability distribution tend to be more accurate [13, 14]. Its simple model can be seen in the following equation:(4)xk=Fkxk−1+Bkuk+wk,where Fk is the state transition model applied to previous state xk−1, Bk is the control input model applied to the control vector uk, and wk is the process noise. A discrete model of Kalman was implemented in accordance with [24].


### 2.6. Feature Selection

For the consistency between extracted data across all paradigms, six statistical features (SV, KR, SS, SM, SK, and SP) were used to extract information across data [28, 43, 46, 54, 55]. SM is calculated in the following equation as(5)SM=1n∑X=1nAx,where Ax is the input signal such that ΔcHbOt and n is the total number of observations. Signal variance is calculated in equation (6) as(6)σ2=1n−1∑x=1n−1Zx−μ2,where Zx is the input signal, μ is the mean found from equation (6), n is the number of samples, and σ is the standard deviation. For KR, calculation was made by the following equation:(7)kurtX=EX−μσ4,where X is the input signal and E is the expected value of X.SK is the asymmetry of values relative to normal distribution around the mean, hence calculated in the following equation:(8)skewX=EX−μσ3.MATLAB®polyfit function fits the line to all input data points, therefore used to calculate SS. Similarly, max function was used to calculate SP. Statistical features were rescaled between 0 and 1 using the following equation:(9)Z′=Z−minZmaxZ−minZ,where Z′ is the rescaled feature and Z refers to the original feature values. The scatter plot across all six statistical features for the OB vs MA task of subject 1 is shown in Figure 3.Figure 3
Extracted feature scatter plot for OB vs MA.


### 2.7. Support Vector Machine

Statistical significance of accuracy is analyzed for selection of an optimal filter for a specific cortical region; therefore, for higher classification performance, nonlinear SVM is used [35, 38]. It can rescale high-dimensional data and can control errors explicitly by maximizing the margins between two or more classes thus creating hyperplanes named as support vectors [40, 41, 56].The cost function that is to be maximized for the SVM classifier gives a correlation between training data and hyperplane as defined in equations (10) and (11), respectively,(10)Minimize12z2+k∑x=1Nεx,(11)ProvidedyxZTξx+a3≥1−εx,εx≥0,where z,ξxεR2,bεR1,z2=zTz, k is the positive regularization parameter, εx is the measure of the training error, and yx is the class label for the nth sample. Here, the third-degree polynomial kernel function with k=0.5 and 10-fold cross-validation was applied for the estimation of classification accuracies [53].


## 2.1. Experimental Setup

To acquire experimental data, seven paradigms owing to three main cortical regions were designed and explicitly performed using a dynamic near-infrared optical tomography (DYNOT-232; NIRx Medical Technologies, NY, USA) device at Pusan National University. It operates on two wavelengths that are 760 and 830 nm where the signal acquisition sampling frequency was 1.81 Hz. Five healthy subjects with normal or corrected-to-normal vision took part in the experiment with a verbal consent before experimentation. All subjects were right-handed with an age range of 26 ± 3. Right-handed subjects were selected to minimize hemodynamic response variation due to hemispheric-dominance difference. The experimental participants had no history of alcoholism, psychiatric, neurological, and visual disorder, cardiovascular and respiratory disease, mental illness, or any motor disability. Moreover, three hours before the commencement of the study, participants were asked to refrain from caffeinated drinks. As discussed in the literature [28, 30, 31, 35, 37, 38], for thinking-related task, signals were acquired from PFC, similarly for motor execution tasks from the primary motor cortex (PMC), and visual stimuli task from VC. Performed experiments were according to the latest Declaration of Helsinki.


## 2.2. Experimental Paradigms

In accordance with the literature [35, 39], subjects were seated on a comfortable chair and were asked to take rest with restricted movements as they can, so that the hemodynamic response activation owing to previous activities can be avoided. Hence, each paradigm related to PFC, MC, and VC starts with a rest of 20 s period to set up the baseline conditions. As the literature [42–45] show, (10∼12)  s task is adequate to acquire hemodynamic response of brain activity, hence 20 s initial rest was followed by 10 s task, and this was followed in turn by another 20 s rest period permitting signals to return to their baseline values before the start of the next trial in paradigm (a). For paradigm (b), 20 s rest after 10 s task 1 was again followed by 10 s task 2. The 20 s rest between two 10 s tasks was added to differentiate two classes through the baseline value. Figure 1 depicts paradigm (a) and paradigm (b). For optimal filter selection, tasks were selected concerning specific cortical regions such that for PFC, three different experiments were performed, rest (R) vs mental arithmetic (MA), R vs object rotation (OB), and OB vs MA task. Similarly, for MC, R vs left finger tapping (LFT), R vs right finger tapping (RFT), and LFT vs RFT. Likewise, for VC, R vs visual stimuli (VS) task.Figure 1
Experimental paradigms for cortical tasks. (a) Task vs rest. (b) Task 1 vs task 2.
(a)(b)


### 2.2.1. Prefrontal Tasks

(1) Mental Arithmetic. According to previous studies [44, 46, 47], for R vs MA task, subjects were asked to perform a series of arithmetic calculations for 10 s based on the pseudorandom order, such that to subtract the random two-digit number (between 10 and 20) successively from the previous result of the three-digit number subtraction appearing on the screen (e.g., 400-11, 389-17, and 372-14). Afterwards, the screen was turned black so that subjects does not go beyond the 10 s task.(2) Object Rotation. For the object rotation task, subjects were asked to imagine a cube rotating for the 10 s task while the “object rotation” word appeared for 10 s on the screen at about 2 m distance [48–50].(3) Mental Arithmetic vs Object Rotation. In this protocol, subjects were asked to perform the aforementioned MA task versus OB task in between 20 s of rest to distinguish two tasks. The experimental paradigm is depicted in Figure 1(b).


### 2.2.2. Motor Tasks

(1) Finger Tapping. According to the literature [8, 17, 19, 23, 31, 37, 38, 51], subjects were asked to tap the self-paced index finger of one hand for 10 s afterwards the 20 s rest task was performed allowing signals to return to their reference values. Also, repetition for 10 times was performed as depicted in Figure 1(a). Similar trials were performed on the other hand, while for the LFT vs RFT task, 20 s rest was performed in between two tasks for the restoration of the signal to the baseline level as shown in Figure 1(b).


### 2.2.3. Visual Task

(1) Checker Box Flickering. In this experiment [4, 18, 39], a screen was placed in front of the subjects at a distance of approximately 2 m, and also subjects were requested to avoid eye blinking during the experiment. The 10 s task of checker box flickering at 4 Hz was performed followed by 20 s rest of the black screen. The sound was also generated during the transition between rest and task. The paradigm followed for visual stimuli is shown in Figure 1(a).


## 2.2.1. Prefrontal Tasks

(1) Mental Arithmetic. According to previous studies [44, 46, 47], for R vs MA task, subjects were asked to perform a series of arithmetic calculations for 10 s based on the pseudorandom order, such that to subtract the random two-digit number (between 10 and 20) successively from the previous result of the three-digit number subtraction appearing on the screen (e.g., 400-11, 389-17, and 372-14). Afterwards, the screen was turned black so that subjects does not go beyond the 10 s task.(2) Object Rotation. For the object rotation task, subjects were asked to imagine a cube rotating for the 10 s task while the “object rotation” word appeared for 10 s on the screen at about 2 m distance [48–50].(3) Mental Arithmetic vs Object Rotation. In this protocol, subjects were asked to perform the aforementioned MA task versus OB task in between 20 s of rest to distinguish two tasks. The experimental paradigm is depicted in Figure 1(b).


## 2.2.2. Motor Tasks

(1) Finger Tapping. According to the literature [8, 17, 19, 23, 31, 37, 38, 51], subjects were asked to tap the self-paced index finger of one hand for 10 s afterwards the 20 s rest task was performed allowing signals to return to their reference values. Also, repetition for 10 times was performed as depicted in Figure 1(a). Similar trials were performed on the other hand, while for the LFT vs RFT task, 20 s rest was performed in between two tasks for the restoration of the signal to the baseline level as shown in Figure 1(b).


## 2.2.3. Visual Task

(1) Checker Box Flickering. In this experiment [4, 18, 39], a screen was placed in front of the subjects at a distance of approximately 2 m, and also subjects were requested to avoid eye blinking during the experiment. The 10 s task of checker box flickering at 4 Hz was performed followed by 20 s rest of the black screen. The sound was also generated during the transition between rest and task. The paradigm followed for visual stimuli is shown in Figure 1(a).


## 2.3. Experimental Setup

Since the mental imagery task activates the PFC [34, 43], a total 11 of near-infrared (NI) light optodes were placed on PFC, 3 of which were detectors and 8 were the source in accordance with the literature [46, 47]. Similarly, for the motor execution task, the primary motor cortex (PMC) is activated [37, 43], hence 15 optodes were placed on PMC out of which 8 were the source and 7 were detectors. To extract data for LFT, optodes were placed on the right hemisphere, while for RFT on the left hemisphere [43]. Similarly, for the visual stimuli task data acquisition from the visual cortex [4, 18, 39], eleven optodes were placed having eight sources and three detectors. The distance between the source and the detector was 3 cm. Optode placement with channel configuration for MC, VC, and PFC is shown in Figure 2.Figure 2
Optode placement configuration in accordance with the 10–20 international system. (a) Motor cortex. (b) Visual cortex. (c) Prefrontal cortex.
(a)(b)(c)


### 2.3.1. Signal Acquisition

In accordance with the literature [6, 8, 29], the raw optical density signal is converted to oxyhemoglobin (ΔcHbOt) and deoxyhemoglobin (ΔcHbRt) concentration using the modified Beer–Lambert law (MBLL) as decribed in the following equation:(1)ΔcHbOtΔcHbRt=βHbOλ1βHbRλ1βHbOλ2βHbRλ2−1Δψt,λ1Δψt,λ2d∗l,where βHbXλ is the HbX extinction coefficient in µM−1cm−1, d is the differential path length factor for the curved path in (mm), l is the detector and emitter distance in (mm), and ΔψHbXt is the absorbance difference of the light emitter wavelength of λi.


## 2.3.1. Signal Acquisition

In accordance with the literature [6, 8, 29], the raw optical density signal is converted to oxyhemoglobin (ΔcHbOt) and deoxyhemoglobin (ΔcHbRt) concentration using the modified Beer–Lambert law (MBLL) as decribed in the following equation:(1)ΔcHbOtΔcHbRt=βHbOλ1βHbRλ1βHbOλ2βHbRλ2−1Δψt,λ1Δψt,λ2d∗l,where βHbXλ is the HbX extinction coefficient in µM−1cm−1, d is the differential path length factor for the curved path in (mm), l is the detector and emitter distance in (mm), and ΔψHbXt is the absorbance difference of the light emitter wavelength of λi.


## 2.4. Signal Processing

The acquired raw signals of the brain contain various noises that can be categorized into physiological, experimental, and instrumental noise [8, 52]. In fNIRS, the instrumental and experimental noises are removed prior to change of the raw signal to its magnitude through MBLL [8, 12, 14].


### 2.4.1. Instrumental Noises

Noises produced due to hardware or by surrounding are categorized under instrumental noises. These noises usually have high frequency that can be removed using a low-pass filter; furthermore, keeping isolation from external sources such as light can reduce such type of noises [8].


### 2.4.2. Experimental Noises

These errors contain motion artifacts due to unintentional body movements like head motion during signal acquisition. It can cause the dislocation of optodes from the assigned positions that result in a spike-like noise due to change in light intensity. Filters like Kalman and Wiener can be used to remove such type of noises [12–15].


### 2.4.3. Physiological Noises

These noises are produced due to Mayer waves (∼0.1 Hz), respiration (∼0.5 Hz), and heartbeat (1∼1.5 Hz) that occur due to fluctuations of blood pressure [16–19]. These noises can be removed using adaptive or bandpass filtering [8, 20, 38].


## 2.4.1. Instrumental Noises

Noises produced due to hardware or by surrounding are categorized under instrumental noises. These noises usually have high frequency that can be removed using a low-pass filter; furthermore, keeping isolation from external sources such as light can reduce such type of noises [8].


## 2.4.2. Experimental Noises

These errors contain motion artifacts due to unintentional body movements like head motion during signal acquisition. It can cause the dislocation of optodes from the assigned positions that result in a spike-like noise due to change in light intensity. Filters like Kalman and Wiener can be used to remove such type of noises [12–15].


## 2.4.3. Physiological Noises

These noises are produced due to Mayer waves (∼0.1 Hz), respiration (∼0.5 Hz), and heartbeat (1∼1.5 Hz) that occur due to fluctuations of blood pressure [16–19]. These noises can be removed using adaptive or bandpass filtering [8, 20, 38].


## 2.5. Data Analysis

NIRS-SPM is a toolbox designed for fNIRS data analysis. For signal processing, it provides common filtering techniques, namely, Butterworth, Gaussian, and hrf [26]. The comparative data analysis was performed by implementing Gaussian and hrf filtering using the NIRS-SPM toolbox, while other techniques on MATLAB® 2017b. The generalized mathematical models with necessary details are as follows.


### 2.5.1. Gaussian Filtering

A Gaussian filter is used in various forms depending upon the nature of the signal. Generally, a Gaussian filter is based on a Gaussian function which defines the probability distribution of noise or data. It can also be used as a smoothing operator. A Gaussian kernel is used for smoothing the signal in which each value is replaced with the weighted average of itself and its neighboring values [14, 21, 26]. A simple representation of the 2D Gaussian filter can be defined in the following equation as(2)Gx,y=12πσ2e−x2+y2/2σ2,where x and y are the distance from the origin in horizontal and vertical axis and σ is the known standard deviation of the distribution.


### 2.5.2. Hemodynamic Response Filter (hrf)

The hrf is based on the canonical representation of the hemodynamic response functions (HRF) and is used for the temporal smoothing of the fNIRS time series signal. In NIRS-SPM, the given functional data were smoothen using the least square estimate with ideal HRF. The hrf and Gaussian filter model details are in accordance with the literature [26].


### 2.5.3. Butterworth Filter

Butterworth filter is a model-based bandpass filter which performs on frequency attenuation using high and low-pass filter. The filtered value not only depends on the weighted average of the unfiltered time series, but also recursively on the previous values of the filtered time series. This filter aims to have a flat frequency response in the desired pass band [26]. The 4th order Butterworth filter with a band pass of (0.01∼0.1) Hz was applied by MATLAB® build in a library for the desired experiments.


### 2.5.4. Finite Impulse Response Filter

FIR filter is designed by finding the coefficients and filter order so that it performs a cross-correlation between the input signal and the time reversed impulse response; therefore, by sampling the pulse shape, coefficients of the filter are designed [14]. FIR filter of order N can be defined as the following equation:(3)xn=∑i=0Naiyn−i,where yn is the input signal, xn is the output signal, N is the filter order, and bi is the value of impulse response at ith instance. Here, the 4th order FIR filter with a low-pass band of 0.1 Hz was utilized. The coefficients were estimated using the least square estimate. However, a time-varying Wiener filter, based on the short-time Fourier series, was implemented as in [25, 53].


### 2.5.5. Kalman Filter

Kalman filters the input signal containing statistical and other noises by linear quadratic estimation. The estimated unknown variables based on the Bayesian inference and joint probability distribution tend to be more accurate [13, 14]. Its simple model can be seen in the following equation:(4)xk=Fkxk−1+Bkuk+wk,where Fk is the state transition model applied to previous state xk−1, Bk is the control input model applied to the control vector uk, and wk is the process noise. A discrete model of Kalman was implemented in accordance with [24].


## 2.5.1. Gaussian Filtering

A Gaussian filter is used in various forms depending upon the nature of the signal. Generally, a Gaussian filter is based on a Gaussian function which defines the probability distribution of noise or data. It can also be used as a smoothing operator. A Gaussian kernel is used for smoothing the signal in which each value is replaced with the weighted average of itself and its neighboring values [14, 21, 26]. A simple representation of the 2D Gaussian filter can be defined in the following equation as(2)Gx,y=12πσ2e−x2+y2/2σ2,where x and y are the distance from the origin in horizontal and vertical axis and σ is the known standard deviation of the distribution.


## 2.5.2. Hemodynamic Response Filter (hrf)

The hrf is based on the canonical representation of the hemodynamic response functions (HRF) and is used for the temporal smoothing of the fNIRS time series signal. In NIRS-SPM, the given functional data were smoothen using the least square estimate with ideal HRF. The hrf and Gaussian filter model details are in accordance with the literature [26].


## 2.5.3. Butterworth Filter

Butterworth filter is a model-based bandpass filter which performs on frequency attenuation using high and low-pass filter. The filtered value not only depends on the weighted average of the unfiltered time series, but also recursively on the previous values of the filtered time series. This filter aims to have a flat frequency response in the desired pass band [26]. The 4th order Butterworth filter with a band pass of (0.01∼0.1) Hz was applied by MATLAB® build in a library for the desired experiments.


## 2.5.4. Finite Impulse Response Filter

FIR filter is designed by finding the coefficients and filter order so that it performs a cross-correlation between the input signal and the time reversed impulse response; therefore, by sampling the pulse shape, coefficients of the filter are designed [14]. FIR filter of order N can be defined as the following equation:(3)xn=∑i=0Naiyn−i,where yn is the input signal, xn is the output signal, N is the filter order, and bi is the value of impulse response at ith instance. Here, the 4th order FIR filter with a low-pass band of 0.1 Hz was utilized. The coefficients were estimated using the least square estimate. However, a time-varying Wiener filter, based on the short-time Fourier series, was implemented as in [25, 53].


## 2.5.5. Kalman Filter

Kalman filters the input signal containing statistical and other noises by linear quadratic estimation. The estimated unknown variables based on the Bayesian inference and joint probability distribution tend to be more accurate [13, 14]. Its simple model can be seen in the following equation:(4)xk=Fkxk−1+Bkuk+wk,where Fk is the state transition model applied to previous state xk−1, Bk is the control input model applied to the control vector uk, and wk is the process noise. A discrete model of Kalman was implemented in accordance with [24].


## 2.6. Feature Selection

For the consistency between extracted data across all paradigms, six statistical features (SV, KR, SS, SM, SK, and SP) were used to extract information across data [28, 43, 46, 54, 55]. SM is calculated in the following equation as(5)SM=1n∑X=1nAx,where Ax is the input signal such that ΔcHbOt and n is the total number of observations. Signal variance is calculated in equation (6) as(6)σ2=1n−1∑x=1n−1Zx−μ2,where Zx is the input signal, μ is the mean found from equation (6), n is the number of samples, and σ is the standard deviation. For KR, calculation was made by the following equation:(7)kurtX=EX−μσ4,where X is the input signal and E is the expected value of X.SK is the asymmetry of values relative to normal distribution around the mean, hence calculated in the following equation:(8)skewX=EX−μσ3.MATLAB®polyfit function fits the line to all input data points, therefore used to calculate SS. Similarly, max function was used to calculate SP. Statistical features were rescaled between 0 and 1 using the following equation:(9)Z′=Z−minZmaxZ−minZ,where Z′ is the rescaled feature and Z refers to the original feature values. The scatter plot across all six statistical features for the OB vs MA task of subject 1 is shown in Figure 3.Figure 3
Extracted feature scatter plot for OB vs MA.


## 2.7. Support Vector Machine

Statistical significance of accuracy is analyzed for selection of an optimal filter for a specific cortical region; therefore, for higher classification performance, nonlinear SVM is used [35, 38]. It can rescale high-dimensional data and can control errors explicitly by maximizing the margins between two or more classes thus creating hyperplanes named as support vectors [40, 41, 56].The cost function that is to be maximized for the SVM classifier gives a correlation between training data and hyperplane as defined in equations (10) and (11), respectively,(10)Minimize12z2+k∑x=1Nεx,(11)ProvidedyxZTξx+a3≥1−εx,εx≥0,where z,ξxεR2,bεR1,z2=zTz, k is the positive regularization parameter, εx is the measure of the training error, and yx is the class label for the nth sample. Here, the third-degree polynomial kernel function with k=0.5 and 10-fold cross-validation was applied for the estimation of classification accuracies [53].


## 3. Results and Discussion

### 3.1. Results

In this study, an optimal filter was chosen based on cortical tasks. Activities were categorized based on three main regions such as PFC, MC, and VC.ΔcHbOt signals were filtered using six filters. Figures 4–6 show the averaged ΔcHbOt filtered signals across trials of prefrontal, motor execution and, visual stimuli tasks. The three different signals show various paradigms for each brain region, while the horizontal axis is aligned with one complete event. The pictorial analysis of filtered responses shows that using the 4th order Butterworth changes data form altogether while discrete Kalman and Wiener filtered the signals, but there remain some noises in the output response, whereas, canonical hrf and Gaussian give much smoother response as compared to any other technique. Moreover, previous studies have commonly utilized Butterworth for signal processing, while the obtained visual does not show any significant improvement in comparison. Likewise, the Gaussian filter is also utilized in various studies that assume the acquired data with some normal distribution; therefore, it considers Gaussian function which defines the probability distribution of noise or data. It can also be used as a smoothing operator. A Gaussian kernel is used for smoothing the signal in which each value is replaced with the weighted average of itself and its neighbouring values. However, the hrf filter considers ideal hemodynamic response as a reference for smoothing the data that are closely related to chi-squared distribution. Due to this fact, hrf outperforms any other considered technique.Figure 4
Filtered averaged HbO of PFC tasks.Figure 5
Filtered averaged HbO of MC tasks.Figure 6
Filtered averaged HbO of VS task.For consistency among extracted information, classification accuracies were obtained across six statistical features. The accuracies obtained across these features are shown in Table1. The two-tailed independent t-test was performed to check statistical significance with Holm–Bonferroni for multiple comparison of filters. Results obtained in Table 2 show significance using hrf or the Gaussian filter for R vs MA and R vs OB while for the OB vs MA task, no significant filter was found; however, the time-varing Wiener filter outperformed others. Similarly, in motor execution tasks, hrf was found to be significant for R vs LFT, likewise for R vs RFT hrf and Gaussian. Moreover, for LFT vs RFT, the hrf filter shows better performance. In the VS task, hrf and Gaussian were significant. Moreover, we plot the mean accuracies of filters with respect to cortical regions as depicted in Figure 7. The Gaussian and canonical hrf performed consistently better across all three cortical regions. However, comparatively better accuracies were obtained using hrf only. The statistical results validate our hypothesis of optimal filter selection by assuring that hrf generally has better performance for fNIRS-based studies, which is in accordance with a previous study [53].Table 1
Accuracies across six features for multiple cortical paradigms.
Prefrontal tasksR vs MASubjectsButterFIRGaussianhrfKalmanWienerS157.0362.2281.4882.9682.9680S257.7759.2571.8583.768.14869.62S358.5154.0788.1490.3773.3382.22S464.4471.8582.2282.9672.5974.81S555.5569.6285.1882.2261.4874.81R vs OBS163.768.148086.6671.4875.18S254.8168.8878.518068.8873.33S36065.9282.2289.6289.2274.81S463.771.1197.0397.0390.3794.81S5606082.9677.0368.1465.18OB vs MAS15064.4468.8875.5585.5590S253.3366.6672.2272.2284.4486.66S354.4468.8885.558076.6674.44S456.6671.1181.1183.3385.5586.66S557.7775.5582.2282.2285.5588.88Motor execution tasksR vs LFTS165.1857.0377.0381.4857.0366.66S26057.7778.5185.1845.1850.37S355.5549.6267.472.5959.2559.25S453.3356.2968.1473.336061.48S545.1857.7776.2976.2945.1845.18R vs RFTS161.4843.778.5189.6266.6665.92S259.2558.5168.8879.2559.2559.25S357.7765.9268.8862.9646.6660.74S454.8151.8566.6674.076047.4S562.9658.5171.8577.0355.7756.29LFT vs RFTS154.4446.6687.7792.2264.4460S255.5561.1164.4476.6653.3351.11S358.8858.8864.4468.8863.3345.55S454.4457.7778.888048.8852.22S56062.22708055.5556.66Visual stimuli taskR vs VSS153.3360.7471.1165.1848.8851.85S254.8166.6674.0775.5559.2562.96S359.2564.4476.2982.2271.1169.62S457.7747.468.1473.3349.6257.77S568.8862.9673.3378.5156.2961.48Table 2
Statistical significance of filters across multiple cortical paradigms.
Filtersp valuesPrefrontal tasksR vs MAButter vs hrf, FIR, Gaussian, Kalman, Wiener0.089, 0.012, 0.012, 0.012, 0.012hrf vs Gaussian, Wiener, FIR, Kalman0.089, 0.012, 0.012, 0.022Gaussian vs FIR, Kalman, Wiener0.012, 0.089, 0.012Kalman vs Wiener, FIR0.089, 0.089Wiener vs FIR0.012R vs OBButter vs FIR, Gaussian, Kalman, Wiener, hrf0.022, 0.012, 0.012, 0.012, 0.012hrf vs Kalman, Gaussian, Wiener, FIR0.012, 0.158, 0.012, 0.012Gaussian vs FIR, Kalman, Wiener0.012, 0.089, 0.012Kalman vs Wiener, FIR0.327, 0.022Wiener vs FIR0.012OB vs MAButter vs FIR, Gaussian, Kalman, Wiener, hrf0.012, 0.012, 0.012, 0.012, 0.012hrf vs Kalman, Gaussian, Wiener, FIR0.089, 0.281, 0.089, 0.012Gaussian vs FIR, Kalman, Wiener0.012, 0.089, 0.089Kalman vs Wiener, FIR0.089, 0.012Wiener vs FIR0.012Motor execution tasksR vs LFTButter vs FIR, Gaussian, Kalman, Wiener, hrf0.65, 0.025, 1, 0.158, 0.025hrf vs Kalman, Gaussian, Wiener, FIR0.012, 0.022, 0.012, 0.025Gaussian vs FIR, Kalman, Wiener0.012, 0.012, 0.012Kalman vs Wiener, FIR0.041, 1Wiener vs FIR0.32R vs RFTButter vs FIR, Gaussian, Kalman, Wiener, hrf0.089, 0.012, 0.5, 0.5, 0.012hrf vs Kalman, Gaussian, Wiener, FIR0.012, 0.089, 0.012, 0.089Gaussian vs FIR, Kalman, Wiener0.012, 0.012, 0.012Kalman vs Wiener, FIR0.5, 0.32Wiener vs FIR0.32LFT vs RFTButter vs FIR, Gaussian, Kalman, Wiener, hrf0.158, 0.012, 0.32, 0.089, 0.012hrf vs Kalman, Gaussian, Wiener, FIR0.012, 0.012, 0.012, 0.012Gaussian vs FIR, Kalman, Wiener0.012, 0.012, 0.012Kalman vs Wiener, FIR0.32, 0.32Wiener vs FIR0.089Visual stimuliR vs VSButter vs FIR, Gaussian, Kalman, Wiener, hrf0.32, 0.012, 0.32, 0.5, 0.012hrf vs Kalman, Gaussian, Wiener, FIR0.012, 0.089, 0.012, 0.012Gaussian vs FIR, Kalman, Wiener0.012, 0.012, 0.012Kalman vs Wiener, FIR0.089, 0.32Wiener vs FIR0.32Figure 7
Filters’ mean accuracies.


### 3.2. Discussion

Despite of the fact that fNIRS offers portability, low cost, and ease of equipment setup, there still remains a challenge of removing noises like systematic physiological (Mayer waves, muscle activity, blood pressure, and respiration and heart rate) and artifacts [17, 57], as fNIRS signals are highly contaminated by measurement noises and physiology-based systemic interference [58]. Several studies have proposed methodologies that can remove noises robustly. In [59], temporal filtering using a low-pass filter with 0.6 Hz cutoff frequency and canonical hemodynamic response function with 4 s full width at half maximum was applied. In [60], exponential moving average and Chebyshev filter were used to remove artifacts from the fNIRS data. Similarly, in [61], only low-pass filter with a 0.14 Hz cutoff was applied to remove physiological noises from fNIRS signals. Aqil et al. [62] used recursive least square estimation for online imaging. This adaptive approach provided a spatial filtering with low and high pass, detrending, and baseline correction. Similarly, Seo et al. [63] assessed the utility of NIRS in removal of physiological noises in fMRI data by reducing variance of the residual error in the baseline fMRI signal through the NIRS signal in the model. Similarly, in [64], adaptive filtration with the affine projection algorithm was used to accelerate convergence with colored noise, but it increases computational cost. In [65], a bandpass filter based on the 5th order Butterworth filter was used to filter motor execution signals based on EEG. In literature [52], fNIRS-based walking signals and walking signals while talking were acquired. These signals were preprocessed using a low-pass with finite impulse response filter, while the talking task results in low-amplitude artifacts with a similar frequency of hemodynamic response, hence it may affect cortical activity [57]. In [66], multiple filters were applied for EEG- and ECG-based signal preprocessing such as for removal of motion artifacts, a median filter with 5-point, for systemic component removal, a second-order Chebychev type-2 filter, high-pass filter for signal drift removal, 1st order low-pass Butterworth filter for blood pressure signal, and the 4th order Butterworth bandpass filter for ECG were applied, while using too many filters for every aspect may result in cost of computation. In [8], previously mentioned aspects in Section 2.4 were briefly discussed. To remove such noise methods like MBLL, eigen-based vector approach [19], a bandpass filter was introduced in a MATLAB-based graphical user interface-based program, HomER. However, no such statistical significance comparison was seen for generic tasks. In literature [14, 15], cancellation of motion artifacts using Wiener and discrete Kalman filter was discussed; meanwhile, t-tests performed in comparison showed better performance with Kalman. Similarly, in [12, 13], four techniques, namely, Kalman, principal component analysis, wavelet analysis, and spline interpolation, were compared to remove NIRS data artifacts. Results showed that spine interpolation and wavelet analysis were significant for such noise removals. In [67], systemic noise was removed using wavelet minimum description length detrending approach and artifact using moving standard deviation and spline interpolation. Eggenberger et al. [68] removed movement artifacts by using visual inspection, and Mayer waves were avoided using averaged blocks. In [69], physiological noise was removed with a low-pass filter, and sliding window was applied for motion artifact rejection. Holtzer et al. [70] combined independent component analysis and principal component analysis to remove noise and signal drifts. In [70–77], the signal was low-pass filtered with a cutoff frequency at 0.14 Hz, and in [78], with a cutoff frequency of 0.2 Hz and in [78], a low-pass filter set at 0.67 Hz, meanwhile a moving average filter with a width of 4 s was used to smoothen the signal. In [79], Gaussian smoothing with a full width at half max of 2 s was applied, while motion artifacts were removed using the wavelet minimum description length detrending algorithm. In [80], a bandpass filter (0.01–1.25 Hz) was applied to the signal. Likewise, in [81], a bandpass filter (0.01–0.2 Hz) was applied while motion artifacts were removed through principal component analysis and spike rejection. Similarly, in [82], data were bandpass filtered (0.01 Hz to 0.14 Hz), while the wavelet filter and correlation-based signal improvement were applied to remove motion artifacts. In [83], data were filtered with a 0.01 Hz high-pass filter and a 5.0 s moving average filter and also principal component analysis were applied to reduce physiological noise. Metzger et al. [84] removed channels with large movement or technical artifacts. Smaller artifacts were corrected with the correlation-based signal improvement method while a low-pass filter with a 5.0 s moving average filter was applied. In [85–87], signals were analyzed with SPM99 (Statistical Parametric Mapping software; Wellcome Department of Cognitive Neurology, London, UK). In [88], moving standard deviation-based artifact removal (moving artifact reduction algorithm: MARA) with a threshold of 0.45 for HbO and 0.18 for HbR was applied and signals were linearly detrended and low-pass filtered at 0.1 Hz.Keeping in view the discussions in the literature, various methodologies were adopted for artifact removal while for systematic physiological noises, mostly high and low-pass filter were applied. Therefore, to make use of the best optimal filter for a specific cortical task, we performed an explicit and systematic analysis to find the filter based on statistical significance for mostly used cortical tasks in fNIRS study applications. The result shows that hrf outperforms the discussed techniques due to the fact that it considers ideal hemodynamic response distribution for smoothing the data. Keeping this in view, our future work involves designing of an adaptive filter technique, considering different distributions of data and noise. Therefore, opting more accurate separation of data.


## 3.1. Results

In this study, an optimal filter was chosen based on cortical tasks. Activities were categorized based on three main regions such as PFC, MC, and VC.ΔcHbOt signals were filtered using six filters. Figures 4–6 show the averaged ΔcHbOt filtered signals across trials of prefrontal, motor execution and, visual stimuli tasks. The three different signals show various paradigms for each brain region, while the horizontal axis is aligned with one complete event. The pictorial analysis of filtered responses shows that using the 4th order Butterworth changes data form altogether while discrete Kalman and Wiener filtered the signals, but there remain some noises in the output response, whereas, canonical hrf and Gaussian give much smoother response as compared to any other technique. Moreover, previous studies have commonly utilized Butterworth for signal processing, while the obtained visual does not show any significant improvement in comparison. Likewise, the Gaussian filter is also utilized in various studies that assume the acquired data with some normal distribution; therefore, it considers Gaussian function which defines the probability distribution of noise or data. It can also be used as a smoothing operator. A Gaussian kernel is used for smoothing the signal in which each value is replaced with the weighted average of itself and its neighbouring values. However, the hrf filter considers ideal hemodynamic response as a reference for smoothing the data that are closely related to chi-squared distribution. Due to this fact, hrf outperforms any other considered technique.Figure 4
Filtered averaged HbO of PFC tasks.Figure 5
Filtered averaged HbO of MC tasks.Figure 6
Filtered averaged HbO of VS task.For consistency among extracted information, classification accuracies were obtained across six statistical features. The accuracies obtained across these features are shown in Table1. The two-tailed independent t-test was performed to check statistical significance with Holm–Bonferroni for multiple comparison of filters. Results obtained in Table 2 show significance using hrf or the Gaussian filter for R vs MA and R vs OB while for the OB vs MA task, no significant filter was found; however, the time-varing Wiener filter outperformed others. Similarly, in motor execution tasks, hrf was found to be significant for R vs LFT, likewise for R vs RFT hrf and Gaussian. Moreover, for LFT vs RFT, the hrf filter shows better performance. In the VS task, hrf and Gaussian were significant. Moreover, we plot the mean accuracies of filters with respect to cortical regions as depicted in Figure 7. The Gaussian and canonical hrf performed consistently better across all three cortical regions. However, comparatively better accuracies were obtained using hrf only. The statistical results validate our hypothesis of optimal filter selection by assuring that hrf generally has better performance for fNIRS-based studies, which is in accordance with a previous study [53].Table 1
Accuracies across six features for multiple cortical paradigms.
Prefrontal tasksR vs MASubjectsButterFIRGaussianhrfKalmanWienerS157.0362.2281.4882.9682.9680S257.7759.2571.8583.768.14869.62S358.5154.0788.1490.3773.3382.22S464.4471.8582.2282.9672.5974.81S555.5569.6285.1882.2261.4874.81R vs OBS163.768.148086.6671.4875.18S254.8168.8878.518068.8873.33S36065.9282.2289.6289.2274.81S463.771.1197.0397.0390.3794.81S5606082.9677.0368.1465.18OB vs MAS15064.4468.8875.5585.5590S253.3366.6672.2272.2284.4486.66S354.4468.8885.558076.6674.44S456.6671.1181.1183.3385.5586.66S557.7775.5582.2282.2285.5588.88Motor execution tasksR vs LFTS165.1857.0377.0381.4857.0366.66S26057.7778.5185.1845.1850.37S355.5549.6267.472.5959.2559.25S453.3356.2968.1473.336061.48S545.1857.7776.2976.2945.1845.18R vs RFTS161.4843.778.5189.6266.6665.92S259.2558.5168.8879.2559.2559.25S357.7765.9268.8862.9646.6660.74S454.8151.8566.6674.076047.4S562.9658.5171.8577.0355.7756.29LFT vs RFTS154.4446.6687.7792.2264.4460S255.5561.1164.4476.6653.3351.11S358.8858.8864.4468.8863.3345.55S454.4457.7778.888048.8852.22S56062.22708055.5556.66Visual stimuli taskR vs VSS153.3360.7471.1165.1848.8851.85S254.8166.6674.0775.5559.2562.96S359.2564.4476.2982.2271.1169.62S457.7747.468.1473.3349.6257.77S568.8862.9673.3378.5156.2961.48Table 2
Statistical significance of filters across multiple cortical paradigms.
Filtersp valuesPrefrontal tasksR vs MAButter vs hrf, FIR, Gaussian, Kalman, Wiener0.089, 0.012, 0.012, 0.012, 0.012hrf vs Gaussian, Wiener, FIR, Kalman0.089, 0.012, 0.012, 0.022Gaussian vs FIR, Kalman, Wiener0.012, 0.089, 0.012Kalman vs Wiener, FIR0.089, 0.089Wiener vs FIR0.012R vs OBButter vs FIR, Gaussian, Kalman, Wiener, hrf0.022, 0.012, 0.012, 0.012, 0.012hrf vs Kalman, Gaussian, Wiener, FIR0.012, 0.158, 0.012, 0.012Gaussian vs FIR, Kalman, Wiener0.012, 0.089, 0.012Kalman vs Wiener, FIR0.327, 0.022Wiener vs FIR0.012OB vs MAButter vs FIR, Gaussian, Kalman, Wiener, hrf0.012, 0.012, 0.012, 0.012, 0.012hrf vs Kalman, Gaussian, Wiener, FIR0.089, 0.281, 0.089, 0.012Gaussian vs FIR, Kalman, Wiener0.012, 0.089, 0.089Kalman vs Wiener, FIR0.089, 0.012Wiener vs FIR0.012Motor execution tasksR vs LFTButter vs FIR, Gaussian, Kalman, Wiener, hrf0.65, 0.025, 1, 0.158, 0.025hrf vs Kalman, Gaussian, Wiener, FIR0.012, 0.022, 0.012, 0.025Gaussian vs FIR, Kalman, Wiener0.012, 0.012, 0.012Kalman vs Wiener, FIR0.041, 1Wiener vs FIR0.32R vs RFTButter vs FIR, Gaussian, Kalman, Wiener, hrf0.089, 0.012, 0.5, 0.5, 0.012hrf vs Kalman, Gaussian, Wiener, FIR0.012, 0.089, 0.012, 0.089Gaussian vs FIR, Kalman, Wiener0.012, 0.012, 0.012Kalman vs Wiener, FIR0.5, 0.32Wiener vs FIR0.32LFT vs RFTButter vs FIR, Gaussian, Kalman, Wiener, hrf0.158, 0.012, 0.32, 0.089, 0.012hrf vs Kalman, Gaussian, Wiener, FIR0.012, 0.012, 0.012, 0.012Gaussian vs FIR, Kalman, Wiener0.012, 0.012, 0.012Kalman vs Wiener, FIR0.32, 0.32Wiener vs FIR0.089Visual stimuliR vs VSButter vs FIR, Gaussian, Kalman, Wiener, hrf0.32, 0.012, 0.32, 0.5, 0.012hrf vs Kalman, Gaussian, Wiener, FIR0.012, 0.089, 0.012, 0.012Gaussian vs FIR, Kalman, Wiener0.012, 0.012, 0.012Kalman vs Wiener, FIR0.089, 0.32Wiener vs FIR0.32Figure 7
Filters’ mean accuracies.


## 3.2. Discussion

Despite of the fact that fNIRS offers portability, low cost, and ease of equipment setup, there still remains a challenge of removing noises like systematic physiological (Mayer waves, muscle activity, blood pressure, and respiration and heart rate) and artifacts [17, 57], as fNIRS signals are highly contaminated by measurement noises and physiology-based systemic interference [58]. Several studies have proposed methodologies that can remove noises robustly. In [59], temporal filtering using a low-pass filter with 0.6 Hz cutoff frequency and canonical hemodynamic response function with 4 s full width at half maximum was applied. In [60], exponential moving average and Chebyshev filter were used to remove artifacts from the fNIRS data. Similarly, in [61], only low-pass filter with a 0.14 Hz cutoff was applied to remove physiological noises from fNIRS signals. Aqil et al. [62] used recursive least square estimation for online imaging. This adaptive approach provided a spatial filtering with low and high pass, detrending, and baseline correction. Similarly, Seo et al. [63] assessed the utility of NIRS in removal of physiological noises in fMRI data by reducing variance of the residual error in the baseline fMRI signal through the NIRS signal in the model. Similarly, in [64], adaptive filtration with the affine projection algorithm was used to accelerate convergence with colored noise, but it increases computational cost. In [65], a bandpass filter based on the 5th order Butterworth filter was used to filter motor execution signals based on EEG. In literature [52], fNIRS-based walking signals and walking signals while talking were acquired. These signals were preprocessed using a low-pass with finite impulse response filter, while the talking task results in low-amplitude artifacts with a similar frequency of hemodynamic response, hence it may affect cortical activity [57]. In [66], multiple filters were applied for EEG- and ECG-based signal preprocessing such as for removal of motion artifacts, a median filter with 5-point, for systemic component removal, a second-order Chebychev type-2 filter, high-pass filter for signal drift removal, 1st order low-pass Butterworth filter for blood pressure signal, and the 4th order Butterworth bandpass filter for ECG were applied, while using too many filters for every aspect may result in cost of computation. In [8], previously mentioned aspects in Section 2.4 were briefly discussed. To remove such noise methods like MBLL, eigen-based vector approach [19], a bandpass filter was introduced in a MATLAB-based graphical user interface-based program, HomER. However, no such statistical significance comparison was seen for generic tasks. In literature [14, 15], cancellation of motion artifacts using Wiener and discrete Kalman filter was discussed; meanwhile, t-tests performed in comparison showed better performance with Kalman. Similarly, in [12, 13], four techniques, namely, Kalman, principal component analysis, wavelet analysis, and spline interpolation, were compared to remove NIRS data artifacts. Results showed that spine interpolation and wavelet analysis were significant for such noise removals. In [67], systemic noise was removed using wavelet minimum description length detrending approach and artifact using moving standard deviation and spline interpolation. Eggenberger et al. [68] removed movement artifacts by using visual inspection, and Mayer waves were avoided using averaged blocks. In [69], physiological noise was removed with a low-pass filter, and sliding window was applied for motion artifact rejection. Holtzer et al. [70] combined independent component analysis and principal component analysis to remove noise and signal drifts. In [70–77], the signal was low-pass filtered with a cutoff frequency at 0.14 Hz, and in [78], with a cutoff frequency of 0.2 Hz and in [78], a low-pass filter set at 0.67 Hz, meanwhile a moving average filter with a width of 4 s was used to smoothen the signal. In [79], Gaussian smoothing with a full width at half max of 2 s was applied, while motion artifacts were removed using the wavelet minimum description length detrending algorithm. In [80], a bandpass filter (0.01–1.25 Hz) was applied to the signal. Likewise, in [81], a bandpass filter (0.01–0.2 Hz) was applied while motion artifacts were removed through principal component analysis and spike rejection. Similarly, in [82], data were bandpass filtered (0.01 Hz to 0.14 Hz), while the wavelet filter and correlation-based signal improvement were applied to remove motion artifacts. In [83], data were filtered with a 0.01 Hz high-pass filter and a 5.0 s moving average filter and also principal component analysis were applied to reduce physiological noise. Metzger et al. [84] removed channels with large movement or technical artifacts. Smaller artifacts were corrected with the correlation-based signal improvement method while a low-pass filter with a 5.0 s moving average filter was applied. In [85–87], signals were analyzed with SPM99 (Statistical Parametric Mapping software; Wellcome Department of Cognitive Neurology, London, UK). In [88], moving standard deviation-based artifact removal (moving artifact reduction algorithm: MARA) with a threshold of 0.45 for HbO and 0.18 for HbR was applied and signals were linearly detrended and low-pass filtered at 0.1 Hz.Keeping in view the discussions in the literature, various methodologies were adopted for artifact removal while for systematic physiological noises, mostly high and low-pass filter were applied. Therefore, to make use of the best optimal filter for a specific cortical task, we performed an explicit and systematic analysis to find the filter based on statistical significance for mostly used cortical tasks in fNIRS study applications. The result shows that hrf outperforms the discussed techniques due to the fact that it considers ideal hemodynamic response distribution for smoothing the data. Keeping this in view, our future work involves designing of an adaptive filter technique, considering different distributions of data and noise. Therefore, opting more accurate separation of data.


## 4. Conclusion

To the best of our knowledge, this is the first study to propose filter selection for commonly used cortical tasks of functional near-infrared spectroscopy (fNIRS) contaminated with artifacts and systematic physiological noises. Six filters, namely, Gaussian, hemodynamic response filter (hrf), Butterworth, time-varying Wiener, discrete Kalman, and window-based finite impulse response, were tested. The results obtained have validated the overall statistical significance of the hrf for prefrontal, motor, and visual cortex tasks. Furthermore, signals acquired from different cortical regions may contain different types of noises. Hence, the prime goal of this study was to suggest an optimal filter for a specific task, owing to the specific cortical region so that further studies can achieve maximum accuracies by reliably improving the recovered hemodynamic response function. Outcome of this study shows that there is a significant impact of filter selection on the accuracy of the classified data, therefore facilitating users to avoid analysis of complex signal techniques by themselves.



---

*Source: 9152369-2020-05-22.xml*