Activation Modeling and Classification of Voluntary and Imagery Movements From the Prefrontal fNIRS Signals

The trends in movement-related functional activity measurement for brain-computer interface (BCI) are mostly associated with the central lobe of the brain. This consideration may be a faulty approach for the paralyzed patient. This limitation demands an alternative approach for movement-related BCI. For the first time, we propose the prefrontal hemodynamics for implementing movement-related BCI. This paper aims to model the activation pattern and the classification performances of the prefrontal hemodynamics regarding the movement-related events. Utilizing functional near-infrared spectroscopy (fNIRS) the changes in the concentration of the oxidized hemoglobin and deoxidized hemoglobin regarding voluntary and imagery movements are acquired. With necessary preprocessing, the fNIRS signals are statistically analyzed to localize the most significant activated regions regarding the applied stimuli. The experiment shows that movement-related events have a strong correlation with the prefrontal hemodynamics. The patterns of the movement-related hemodynamics are modeled by polynomial regression and used to classify the voluntary and imagery events based on the maximum similarity approach. The resulting classification accuracies are found promising that proves the effectiveness of the prefrontal fNIRS signal to be effective in movement-related brain functionality analysis.


I. INTRODUCTION
The neural activations regarding the movement-related events are one of the most important research areas for implementing the practical brain-computer interface (BCI). The neural activations can be measured from the brain based on two fundamental methodsi) based on electrical activity estimation from the scalp of brain and ii) based on the hemodynamics estimation. Depending on the electrical activities, the brain function measuring The associate editor coordinating the review of this manuscript and approving it for publication was Kemal Polat . modality is electroencephalography (EEG) and magnetoencephalography (MEG). EEG has a major limitation regarding its spatial resolution (5 to 9 cm) [1]. MEG is not widely used for functional brain imaging due to its high degree of noise sensitivity [2].
The hemodynamics means the changes in the concentration of oxidized hemoglobin (HbO 2 ) and deoxidized hemoglobin (dHb) in the total blood volume of the brain. Based on the hemodynamics, functional magnetic resonance imaging (fMRI) is a gold modality in the field of functional brain imaging. The fMRI provides an excellent spatial resolution (3∼6mm) but its temporal resolution is poor VOLUME 8, 2020 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ (3∼5sec) [3]. In the field of BCI, it has a very slight scope of operation because of its very high cost, motion sensitivity, and heavyweight. On this contrary, Functional Near-Infrared Spectroscopy (fNIRS) is another technique for functional brain imaging that provides very good spatial resolution (∼1-1.5cm), moderate temporal resolution (up to 100Hz), portability to use, high value of signal to noise ratio (SNR), and less motion artifact compared to the others [4]- [6]. Furthermore, fNIRS is not as physically confining as fMRI and it allows more movements compared to fMRI during imaging. Recent publications [7], [8] demonstrate that fNIRS is comparable to fMRI and this method is reliable and valid for cortical activations measurement. As a result, fNIRS is getting the most attention for the recent researches in the field of functional neuroimaging and BCI. The fNIRS is a noninvasive optical method that utilizes the trans-illumination technique with NIR (850∼1000nm) [4] light to measure the change in the concentration of HbO 2 , dHb, and total blood volume. The brain activation level and corresponding functional region can be evaluated by these parameters. A number of research works have been performed by fNIRS in different valid fields of research and applications like mental workload assessment [9], [10], BCI [11], [12], human body monitoring for diagnosis and treatment [13], [14], cognitive skill [15], [16], emotions assessment [17], [18], and many more. In addition, fNIRS proves itself as a suitable modality for the brain-based experiment to develop the mechanisms for neurorehabilitation [19]. The motor activity (both voluntary and imagery) related to brain functionality plays an important role in BCI. Around one billion people in the world's population experience some form of disability [20]. If it could be possible to measure the proper activation pattern of motor action (in the case of voluntary and imagery) from the brain, BCI would be implemented to decline their physical limitations through rehabilitation.
Many researchers investigated the brain functionality related to voluntary and imagery movements by fNIRS. The research works [21]- [24] investigated the properties of voluntary and imagery movements through fNIRS modality and reported the activation strength of motor execution and imagination tasks. These studies suggested that a number of areas in motor cortex become active due to the voluntary and imagery movements and showed an increased amount of HbO 2 concentration. In [25], fNIRS based neuronal activations of imagery movements of the left and right-hand wrists were classified by linear discriminant analysis method. Three classes of brain responses are classified in [26] where two classes were motor imagery tasks. Voluntary movementrelated tasks of left and right hand were classified in [27] from the fNIRS signals of optimally selected channels. Multiple motor imagery tasks (movements of the left hand, right hand, left foot, and right foot) are classified in [28].
Since the most predictable region related to motor execution and imagery tasks is the motor cortex (central portion of the brain) [29], all of the aforesaid research works [21]- [28] conducted the investigations on the activations of the motor cortex. There is a problem of considering central lobe of the brain while the persons are paralyzed or previously disabled because these patients have the inactiveness of the central part of their brains. This issue can be a major obstacle to designing motor imagery based BCI for paralyzed patients. Though motor areas of the paralyzed persons are damaged, in most of cases their frontal areas of the brains remain sound. However, it is reported in [30]- [32] that the motor action planning occurs in the frontal lobe. In addition, some works [23], [33]- [36] reported that there is an indirect relationship of motor actions (especially the motor imagery) with the frontal lobe. Based on aforementioned idea, a simple survey was performed to classify the left and right hands movement (voluntary movement) from the prefrontal fNIRS data and on average we have found 79% classifying accuracy from five participants [37]. In addition, our statistical research work [38] has established the relation between the motor area and prefrontal area of the human brain in case of motor action. These works indicate that the hemodynamic information of movements can be measured from prefrontal cortex by fNIRS. Most of the previous research works [21]- [26], [37] regarding the movement-related stimuli by fNIRS are based on the classification of two class movements. The works also did not model the activation pattern of the different voluntary and imagery movements with their localizations. Therefore, multiple motor activity patterns of voluntary and imagery movements have not been analyzed by fNIRS from the prefrontal hemodynamics. So far our knowledge, the localization of the multiple movement-related stimuli and their classification accuracies has not been performed considering the prefrontal hemodynamics. As a result, there arises a scope to investigate the neural activations of the prefrontal cortex with respect to the stimuli of the voluntary and imagery movements by fNIRS.
In this work, four different movements of the hands and feet have been considered as both voluntary and imagery manner. According to these stimuli, fNIRS data were recorded from the prefrontal cortex of several young and healthy male subjects. With proper signal pre-processing, data were separated based on their class. Then, analysis of variance (ANOVA) was used to localize the most significantly activated area regarding the stimuli. The ANOVA results were verified by effect size (ES) estimation.
The temporal pattern of the change in the concentration of both HbO 2 and dHb signals regarding the most significant areas were modeled by polynomial regression. Utilizing the model activation patterns, a simple statistical classifier has been proposed based on the spatiotemporal maximum similarity concept. This proposed algorithm has been utilized to classify the movement-related fNIRS signals. For classification of the fNIRS signal, two, four, and six class problems were considered. The overall idea of the proposed research framework has been presented in Figure 1.
We have demonstrated the conventional time-domain feature extraction and classification strategy to classify the signals by four different classification algorithms: Artificial Neural Network (ANN), Support Vector Machines (SVM), Linear Discriminant Analysis (LDA), and k-Nearest Neighbor (kNN). We found that the classification accuracies of the proposed and conventional methods are almost similar. Eventually, this work contributes -• To reveal the most significant activation regions of the prefrontal cortex for the different voluntary and imagery movements, • To model the hemodynamic activation pattern based on the change of HbO 2 and dHb concentration using polynomial regression, • To propose a classification algorithm using the proposed hemodynamic activation model, and • Comparing the classification accuracy of the proposed model with the accuracies of the conventional classifiers The rest of the paper is organized as follows: The materials and applied mathematical methods are discussed in Section II. In Section III, the results are presented with discussions. We have concluded our research findings briefly in Section IV.

A. DATA ACQUISITION PROTOCOL
In this work, every participant performed eight different movement-related tasks (four voluntaries + four imageries). Since most significant voluntary movements of the human beings are performed by hands and feet, the movements by hands and feet (by means of voluntary and imagery) were considered for the neural stimuli. The data acquisition protocol was checked and permitted by the ''Data Acquiring Ethics Evaluation Committee (DAEEC)'' of Khulna University of Engineering & Technology (KUET). The subjects were verbally informed about the protocol of the data acquisition and according to the protocol all the subjects lifted their left hand, right hand, left foot, and right foot, sequentially. These four movements were performed in two kinds: i) voluntary and ii) imagery. The time scheduling for the data acquisition protocols is given in Figure 2. In one session this unit protocol was performed four times by a participant. After every session, each participant took rest at least five minutes. Eventually, every participant performed 20 trials for each movement-related task. A MATLAB based graphical protocol (as Figure 3) was designed for this research work that helped the participants by providing the instructions to make the data acquisition procedure easy with proper scheduling. The code of the designed graphical protocol is freely available in [39]. This program blinked according to the scheduled tasks and instructed the participant to follow and perform the tasks. In this graphical program, five different tasks were arranged and those were movements of the left hand, right hand, left foot, right foot, and rest. Finally, eight different tasks have been considered for analysis those are voluntary left hand (LH), right hand (RH), left foot (LF), right foot (RF) and imagery left hand (iLH), right hand (iRH), left foot (iLF), and right foot (iRF)).

B. DATA ACQUISITION
Thirty-five male subjects (age range=20 to 25) participated in the data acquisition procedure. Among them the recorded signals of four participants were excluded due to their poor signal quality. All participants were tested and found right-handed depending on the recommendation of Edinburg Handedness Inventory [40]. No participants had any history of psychiatric, neurological or visual disorder. In addition, no participant had pain in both hands and feet. VOLUME 8, 2020 FIGURE 2. Time schedule of data acquisition protocol for each participant regarding both the voluntary and imagery movements. This is a unit task performing schedule that was repeated four times in each session to complete 20 individual trials of every task.   Also, their verbal consents were taken before the data acquisition related to this research work as the rule of the university. All data acquisition procedures were completed in the Neuroimaging Laboratory of the Biomedical Engineering Department of KUET obeying the declaration of Helsinki [41]. For this work, a 16 channel continuous-wave fNIR system (model: Biopac 1200 fNIR imager) was used. By the system, hemodynamic signals from the prefrontal cortex were acquired from all the participants. The optode band for data acquisition for this fNIR system contains four NIR light sources and 16 detectors. The physical configurations of the optode band on prefrontal cortex are given in Figure 4. The data sampling rate was 2 Hz and COBI (cognitive optical brain imaging) software [42] was used for data acquisition. The total hardware configuration associated with the data acquisition is presented in Figure 5.
In this fNIRS system, two different sources of NIR light (λ 1 =730nm & λ 2 =850nm) were used those are almost transparent to skin, bone, and brain tissue [43]. The chromophores (HbO 2 and dHb) of the blood absorb NIR with dissimilar absorption coefficients. Therefore, from the scattered backlight after absorption, the amount of existed chromophores was measured modified Beer-Lambert law (mBLL), i = i 0 10 −ζ λ where, i = intensity of scatteredback light, i 0 = intensity incident light and ζ λ is the optical density which depends on the wavelength [44]. In addition, the concentrations of dHb and HbO 2 were calculated as [45], [46], In (i) ζ λ 1 and ζ λ 2 are optical density for two corresponding NIR wavelength, λ 1 and λ 2 . Here, α is the excitation coefficient of HbO 2 and dHb in µMole-1mm-1, d PL is the unit-less differential path length factor, and D EM is the distance between emitter and detector. The unit of D EM is mm. Besides, C is the concentration of HbO 2 and dHb in µMole.

C. PREPROCESSING
At first all the fNIRS signals were separated according to the tasks and thereafter signals of the same tasks for all the participants were arranged in individual arrays. Since we used 16 channel fNIRS system, the 16 channel fNIRS data carry the information of 16 different spatial positions of the prefrontal cortex. The positions and corresponding channels are shown in Figure 3. Interestingly, all these channels are not linearly uncorrelated. To check the linearly correlated and uncorrelated channels, principal component analysis (PCA) was applied and found that the signals are of four dimensions. The linearly uncorrelated profiles of the channels have been shown in Figure 6 where the Eigen-space is formed by the three principal components. In this consequence, the prefrontal cortex has been divided into four major region of interest ( It was found that the 4 signals corresponding to four channels in every defined region show a very strong correlation (0.95 < r < 1). Therefore, these four signals can be averaged to reduce the total dimension i × 16 to i × 4. Here, i indicate the sample number. If, a signal, X (N ) i×16 is of i × 16(here N is all samples of that signal), then the procedure to reduce its dimensions according to the previous description as, This selection of ROI is very important to localize the neural activities [48] and to reduce the feature dimension, which helps to achieve the high classification accuracy.
Savitzky-Golay (SG) filter with frame length 21 and order 3 was used to filter the noisy fNIRS signal. For smoothing the fNIR signals, we have used SG filter because of its special benefit i.e., SG filter is better than FIR filter for removing high-frequency contents from signals. Because, in the case of FIR filter, the Euclidian distance between original and filtered signal is more than SG filter [49]. Furthermore, employing SG filter needs no delay correction as an FIR filter. For each trial of the fNIRS data were corrected by subtracting the baseline from the original signal. Baseline was calculated from the average of the first 3 seconds of the task. This consideration [23], [50] ensures the initial signal points regarding each trial to remain at the zero level or the value closed to zero.

D. MODELING
This research work scopes to model the activation pattern with its localization. Consequently, modeling is done by two steps: 1) Statistically finding the localization of the activation, and 2) Activation pattern modeling by polynomial regression. The methods regarding the activation localization and modeling are presented as follows.

1) LOCALIZATION BY STATISTICAL ANALYSIS
One-way repeated ANOVA was used to find the significant ROI for each activity, separately. The signal mean was taken from three sample window (0-10, 10-20, and 20-30 samples) for one-way repeated measures of ANOVA. Student's t-distribution statistics were also used to find the significant difference between the hemispheres regarding any voluntary or imagery event. Since the p values of ANOVA is not enough VOLUME 8, 2020 to find the significance level for a big dataset [51], we also calculated ES of the event to localize the hemodynamic activation. Statistically, we can find the difference in activation level using ES calculation. ES is a simple statistical method of measuring the difference between groups of mental activation due to external excitation and resting-state activation. The widely used method of ES estimation is Cohen's d method. Suppose that, x 1 and x 2 are the mean values of the entries of group 1 and 2, respectively, then, the ES can be estimated between x 1 and x 2 by the following method [52], Here, N 1 and N 2 are the numbers of entries in each vector and S 1 and S 2 are the standard deviation of the dataset, respectively. The ROI's satisfied the significance by both ANOVA and ES test were considered for most significant activated ROI for a specific event.

2) ACTIVATION MODELING
After the confirmation of the statistically discovered activation localization, the fNIRS signals of all trials from all subjects were averaged to fit as model activation patterns regarding the event by polynomial regression. To fit a time series of data by polynomial fitting or regression, we usually contemplate a polynomial equation as an estimation function and assume that the estimation function, E(x) is of k th degree polynomial that can be presented as, Therefore, the difference between the actual value, Y and the estimated value resulting from the proposed model estimating For attaining the finest fitted estimated model equation, it is the foremost target to minimize the value of the residual. It is reported in [53] that hemodynamic activations can be modeled by polynomial fitting with the value of k = 5. In this work, k = 5 has been taken to get the minimum value of residual and for the coefficient, a i at the minimum error condition, the partial derivative of R 2 is zero. To achieve the regression with k th polynomial we get, The previous relation can be represented as, Here, X is a Vandermonde type matrix. This can be solved as, The order of the polynomial equation was estimated from the error performance of the fitted results. A satisfactory level of error was taken as the threshold for different activation modeling.

E. CLASSIFICATION
For classification, the voluntary and imagery tasks from the fNIRS signal, conventional feature extraction and classification technique can be applied. Moreover, in this work, we have proposed a statistical way to classify the fNIRS signal utilizing the proposed activation model of the HbO 2 and dHb concentration. Although the proposed classification method is not too powerful to the conventional classifiers, it is very simple to implement and significant in the results. Here, both the proposed method and the conventional classification method have been applied to compare the results.

1) PROPOSED CLASSIFICATION METHOD
Suppose, an fNIRS signal is to be tested whether it is the signal of LH, RH, iLH, iRH, iLF, or iRF. The testing signal can be presented as S k = [S LL S LM S RM S RL ] that contains 4 columns of data of 4 different ROI's of the prefrontal cortex and every column of signal contain N number of the data sample. The value of k = 1, 2, 3, and 4 indicates the positions (LL, LM, RM, and RL, respectively). Again, we have six model equations of six different stimuli (voluntary and imagery). The temporal activation model can be represented as, M c τ where the notation τ = 1, 2, 3, 4, 5, and 6 (LH, RH, iLH, iRH, iLF, and iRF, respectively) represents the index of the stimuli and c = (1, 2) = (HbO 2 , dHb) represents the concentration information. Therefore, M 2 3 means the activation model of iLH in dHb concentration. To measure the maximum similarity index of the testing signal S k with the models M c τ , the following correlation can be calculated.
Therefore, testing with six different models we get two final correlation coefficient matrices, r 1 τ,k and r 2 τ,k for the HbO 2 and dHb respectively, where, r c τ,k = [r c 1,k r c 2,k r c 3,k r c 4,k r c 5,k r c 6,k ]. The size of r c τ,k is 1 × 24. According to the proposed methodology, the maximum similarity pattern will show the minimum error. Therefore, the maximum value in the elements of the correlation coefficient matrix r c τ,k indicates the signal class. So, the index of the maximum valued element, I dx of the correlation coefficient matrix r c τ,k can be found by applying the following argument, If the model activation patterns are loaded in the algorithm of data testing in the proposed sequences, the proposed model suggests that the minimum valued index of a specific task will follow a fixed pattern. The aforementioned methodology has been briefly presented with a flow diagram in Figure 7. The  index for the value of τ will be 3, 6, 11, 14, 20, and 21 for all six classes, respectively. This same design can be regarded as two to six classes. For that proposition, the value of τ is to set according to the number of the classes. In this work, classifications of the two, four, and six class approaches were conducted.

2) CONVENTIONAL CLASSIFICATION METHOD
In the conventional classification approaches, feature extraction is a primary step to classify the data set. Although there are many techniques to extract features from the signals like independent component analysis, principal component analysis, common spatial pattern etc., these methods are often used in the feature extraction of the complex signals like EEG, MEG, EOG, etc. for its time and frequency domain properties. Since fNIRS signals exhibit simple time-domain characteristics, most of the classification techniques need only time domain features like mean, slope, skewness, median, maximum, minimum, etc. [37], [54]- [56] to classify the signals. Mean, slope, skewness, median, maximum, and minimum was considered as features. For classification purpose, the signals are generally separated into two parts: i) for training purpose and ii) for testing purpose. For a 5-fold cross-validation technique, 4/5 portion of the data is used to train the classifier to make a model and the rest of the data are used to predict the performance of the model. The training and testing files are separated with 5-different combinations from the original data set. Then the overall classification accuracy is found by averaging the results of the 5 different classification accuracies. This training and the testing procedures are shown in Figure 8 where there are two phases: training phase Figure 8(a) and testing phase Figure 8(b).
Although ''No Free Lunch'' theorem claims that no classification mechanism is entirely superior to the other [57], the most commonly used classifiers in fNIRS based event classification was utilized in this work. The review work [58] found that four classifiers are repeatedly used by the researchers of fNIRS based BCI and those are LDA, SVM, kNN, and ANN. These four classifiers were used to classify the events. Though the detail mathematical explanations of these machine learning based algorithms are out of the scope of this paper, a short note on the working method of each individual classification method is given here.
LDA: LDA is widely used classification technique in fNIRS signal classification due to its low computational complexity and high speed [50], [58]. To classify or separate the two or more than two-class data, LDA employs discriminant hyper-plane(s). Since the main mechanism of the LDA is dimension reduction, it chooses the hyperplane(s) by minimizing the ratio of within-class variance and maximizing the ratio of between-class variance (i.e., Fisher's criterion [58]) assuming the data of the classes are Gaussian distributed with equal covariance. Based on Fisher's criterion, the effective projection matrix P is calculated in LDA as [50], [58], In (13), τ b and τ w stand for the between-class scatter matrix and the within-class scatter matrix, respectively. Besides, det(•) represents the determinant of the matrix. We can define τ b and τ w by the following relations [50], [59].
Here, x's are the samples of the feature vector of a class, µ i andX represent the sample mean of class i and the grand mean of the total samples of m classes, respectively. The number of total samples is represented by ν where, ν i represents the number of samples in class, i. The solutions of (14) and (15) can be found considering them as an eigenvalue problem that leads to finding the optimum values of the projection matrix, P. In Matlab 2017a, fitclda() was used to construct the LDA based predictive model which was further utilized with 5-fold cross-validation to check the classification performance. SVM: SVM is an extensively employed classifier for its high prediction accuracy in high-dimensional features [60]- [62]. SVM can be used as a linear or nonlinear method. The main mechanism followed by SVM is to generate the hyper-planes that maximize the margin among the classes. The nearest points of the hyper-planes are called support vectors. The discriminating hyper-plane in a 2D feature space can be formularized as, In (10), x, d ∈ 2 and c ∈ 1 . To find the optimal results of d * maximizes the distance between the hyper-plane and the support vectors. This maximization procedure is obtained by minimizing the following cost function (17) cogitating the VOLUME 8, 2020 restrictions given in (18) [60]- [64].
In (17) d T d = d 2 and is a regularization parameter that can be chosen by the users based on the penalty factor of classification errors. Besides, ξ n represents the measurement of error during the training period, Z represents the number of samples those are misclassified, and y n represent the class label for the n th sample (for a binary classification problem, it is +1 and −1). In this work, we used the Matlab toolbox as one versus all approach of SVM. The SVM structure was defined with polynomial kernel function with default order value 3. kNN: Although kNN is known as a lazy nonparametric classifier, till now it is chosen by the researchers of various fields because of its simplicity. This method does not need any explicit training phase to generalize the training feature vectors. Therefore, the training phase is precisely fast. During the training period, kNN actually keeps all the training features with their labels for the testing phase. The kNN algorithm finds the points from the training data those are nearest to be considered for the selection of the class of a new testing observation. To take decision on the nearest points, there are several distance calculating formulas like Euclidean, Minkowoski, Cityblock, Manhattan, Mahalanobis, Cosine, Chebychev, etc. Therefore, we find there key steps of the kNN approach: i) a set of training feature vectors with label information ii) a distance metric to measure the distance between objects, and iii) the number of the nearest neighbors, k. Suppose, we have a training set, (T (ϕ, y) ∈ T ) that contains the feature vectors, ϕ with their labels, y and a test object t = (ϕ , y ) where ϕ is the feature vectors of the test object and y is its class. Now, the kNN algorithm measures the distance between (ϕ , y ) ∈ t and the training objects (ϕ, y) ∈ T to estimate the nearest neighbor list,(ϕ i , y i ) ∈ T t . From the list of the nearest neighbors, the class of the object will be decided by the following majority voting condition [65], where ν is a class label. On the other hand, y i is the class label for the i th nearest neighbors and I ( * ) is a function that indicates the value 1 for the true argument and 0 otherwise. In our proposed kNN based predictive model, the distance calculation was performed by the Euclidean method with k = 3 for two-class and k = 5for four and six-class classification. ANN: ANN is a complex but very efficient classifier. This algorithm was also used in this research work to find the highest classification accuracy of fNIRS based BCI system [37], [66]. ANN has the quality to mimic the comportment of the human brain. For the feedforward networks, commonly multilayer perceptors consist of three type layers: input, output, and hidden layers. The objective of the input layer is to buffer the distribution of the input signals x n (n = 1, 2, 3...) towards the hidden layer neurons. Each hidden layer neuron adds the input signals (x n ) after weighing the input signals by the strengthsW jn from the input layer and calculate the output, Y where Y is a function of their summations [37].
Here j is neuron numbers, W jn is the weight of a connection between n and j according to their relation, W jn = ηδ j x n . Here, η is the rate of learning parameter and the factor, δ j depends on the condition whether j is an input or hidden neuron. The adjustments of the weights are generally estimated by the back-propagation algorithm. Let, V i be the prediction for j th observation in an ANN system and is a function of the network weights vector W = (W 1 , W 2 , ...). Therefore, e, the total prediction error will also be a function of W as, [66]. For every weightW i , according to the gradient descent algorithm the updating formula is considered as Here, α is the learning parameter and the value of α lies between 0 and 1. In this work, we used two hidden layers for two class classifications and four hidden layers for the four and six class classification problems. In every case, the classification accuracy was calculated with a 5-fold crossvalidation technique.

III. RESULTS AND DISCUSSIONS
The activations of the four voluntary and imagery movements have been presented in Figure 9 and Figure 10, respectively. This is a result of an arbitrarily selected participant. The results have been presented after separating the data based on the ROI's, correcting baseline, and filtering. In the figures, both HbO 2 and dHb activation patterns have been presented. From the graphical depiction of the neural activations, the most significant activated areas and the activation patterns can be assessed. From the hemodynamic responses (HbO 2 and dHb) of four voluntary and four imagery tasks have been given regarding four ROI's (LL, LM, RM, and RL).
From the results of Figure 9, we get voluntary hand movements to create significant activation in the prefrontal cortex. The hemodynamics due to voluntary movements of the left and right hand show contralateral activation in the opposite hemispheres. Due to the left-hand movement, the activation is noticeable in the RM region of the prefrontal cortex and oppositely the LM region is activated due to the right-hand movement. The other regions show the random activation which does not indicate the clear hemodynamic activation. On the other hand, the voluntary movements of the left and right feet also show random activations as given  in Figure 9(c) and Figure 9(d). From the results of imagery movements given in Figure 10, we find that the hemodynamic activations due to the imagery hand movements show similarity with the voluntary movements but the activation strength is lower than that of the voluntary movements. In addition, the imagery feet movements show significant activation in the lateral area of the prefrontal cortex. The imagery left and right foot activates the right and left hemisphere of the prefrontal cortex, respectively where the activations of the RL and LL regions are most noticeable as given in Figure 10(c) and Figure 10(d). To examine the significant neuro-activation considering the data of the total population VOLUME 8, 2020  involved in this research, statistical analysis, ANOVA was performed. One way repeated measures (three levels: 0-5, 5-10, and 10-15 sec) ANOVA was performed on the fNIRS data of the tasks (left-hand, right-hand, left-foot, and right-foot movement as voluntary and imagination manner) to reveal the significant activation localization of the region of interests (LL, LM, RM, and RL). The ANOVA was conducted on the mean value of HbO 2 and dHb concentration. From the results of ANOVA, the following significant hypothesis has been found:   showed a significant difference (t = 13.4297, p < 0.001 for movement execution and t = 4.2874, p < 0.001 for imagery movement) with the right hemisphere (RM+RL). Moreover, the activated region (LM) due to the task (movement of the right-hand as voluntary and imagery manner) showed In case of voluntary and imagery movements of LH and RH some other regions showed significant activation (p < 0.01) based on the results of repeated ANOVA (see Table 1) but for large observations the p-value is not enough [51] for taking a statistical decision. Therefore, we also considered the value of ES to confirm the activation strength of the concerned ROI. The regions showed the significant ES have been considered as the significantly activated regions for the corresponding task (see the Table 1).

LF & RF:
The data regarding the voluntary movements of the left and right foot showed no significant activation either in HbO 2 or dHb concentration with insignificant ES. Since the activation region of the lower body part is situated in the deep brain, this inactiveness may occur. This result suggested us to advise the participants to concentrate deeply to imagine the feet movement during data acquisition so that there a significant cognitive load may occur in the prefrontal cortex.
In case of left and right foot imagery movement, several ROI's were found as significant according to the results of ANOVA. The LF imagery movement showed the significant activations in RL (F(2,90) 2 and F(2,90)= 12.87, p < 0.001 for dHb). Although three different regions were found as activated according to the ANOVA results, two of them showed significant ES values (LL: 1.1872 and RL: 0.6232) and between them, only the LL showed the highest (around two times than RL) ES value. Therefore, LL is the most significant region of activation due to the imagined movement of RF. All the results regarding ANOVA and the ES have been tabulated in Table 1. From the results of ANOVA and ES, the most significant areas for different tasks were selected. The concentration changing pattern of HbO 2 and dHb regarding the most activated area were averaged from the data of all the trials of all participants. The activation patterns of HbO 2 and dHb of the corresponding stimuli have been presented with their error disparity by Figure 11.
Moreover, to show the activation pattern due to the voluntary and imagery movements in the total area of the prefrontal cortex the functional neuroimages are given in Table 2. These functional neuroimages have been prepared based on the grand average value of the total population dividing five equal windows from the beginning of the task to the end of the trial. The concentration changes of the HbO 2 were considered here to construct the neuroimages where the activation level was registered on the MRI brain images. The total area covered with the activation level of the 16 channels was done with 20 points B-Spline interpolation technique using the fNIRSoft [67]. The color bar is placed just aside the neuro images to understand the activation levels.
The mean value of HbO 2 concentration due to voluntary hand movements (both hands) are more than that of the imagery movements of hands. In addition, the activation level of imagery feet movements is comparatively lowest than the voluntary and imagery movements of the hand. These results reveal that the voluntary and imagery movements of upper limbs show more activation than that of lower limbs. To observe the difference in the concentration of HbO 2 and dHb corresponding to different stimuli the mean concentrations are as of Figure 12.   (21) and (22), as shown at the bottom of the page, respectively.
It can be hypothesized that the most activated region for a typical stimulus shows the activation pattern closely correlated with its model equation. Therefore, the maximum temporal similarity pattern with its proper spatial region according to the proposed activation model gives us the class of the signal. Implementing the proposed methodology of the signal classification technique, subject dependent task classification was conducted based on the two, four, and six class perspectives. The subject dependent classification accuracies of the proposed method have been given in Figure 13, Figure 14, andFigure 15. Here, the results of two, four, and six class movement-related tasks are given. Moreover, the data classes have been oriented as two classes with iLH and iRH, four classes with iLH, iRH, iLF, and iRF, and six classes with LH, RH, iLH, iRH, iLF, and iRF. From the results we found that utilizing the model activation pattern of the proposed work, the average classification accuracies are 75.16±7.12 (2-class), 57.58±6.69 (4-class), and 51.11±7.74 (6-class). Since the spatiotemporal activation pattern of the LH & iLH and RH & iRH are similar (see Figure 11 and Figure 12), six class classification accuracies are slightly inferior. Further processing is necessary to improve the classification accuracies of six class problems. On the other hand, the average classification accuracies of the two and four class aspects are quite convincing.
We also classified the fNIRS data of voluntary and imagery movements considering the time domain features (mean, slope, variance, and maximum). The HbO 2 and dHb signals   (22) of proposed ROI's were considered for these time domain feature extraction. The signal window for feature extraction was considered as 5-15s since this window is the mostly activated period [68]. The features of two, four, and six class data were used to train and test the LDA, kNN, SVM, and ANN classifiers. The classification accuracies of the two, four, and six class data by the four different classifiers are given in Figure 13, Figure 14, and Figure 15, respectively along with the results of the proposed classification method. The classification accuracies have been presented as the subject dependent approach.
The average classification accuracies of the conventional classifiers along with the proposed activation model are given in Figure 16. From the results, we can see that the classification accuracies found by ANN and SVM are better than the kNN and LDA. In the case of ANN and SVM, the nonlinear kernel function has been used and the accuracies are quite good for four and six class problems. On the other hand, the proposed classification method shows significant classification accuracies. Although it provides slightly lower classification accuracies than that of the ANN and SVM, the proposed method has no training stage like the conventional classifiers (ANN, SVM, kNN, and LDA). Therefore, by utilizing the proposed mechanism it is easier to implement to find the class of a signal from its activated region (spatial information) and the temporal pattern of the activation.

IV. CONCLUSION
This paper report the localization of activations in the prefrontal cortex due to voluntary and imagery movements and modeled the activations of different movement-related events with polynomial regression. This work opens the doorway to measure the voluntary and imagery movements from the prefrontal hemodynamics by fNIRS modality. From the overall results of the proposed research work, we reveal that the upper limb movements (LH and RH) by voluntary and imagery movements are easily measurable from the prefrontal cortex due to their activation strength than the lower limb (iLF and iRF) activation. The voluntary movements of the lower limb (LF and RF) did not create activation significantly in the prefrontal cortex. This may be happened due to being the functional area of lower lamb activities situated in the deep brain. The hemodynamic activations based on the concentration change in HbO 2 and dHb regarding the significant activities have been also modeled by the polynomial regression. These model equations are useful to assess the prefrontal hemodynamic activation pattern of the proposed movement-related events. In addition, this work proposes a classification technique utilizing the proposed activation models to classify the fNIRS data. We classified the two, four and six class fNIRS data including the voluntary and imagery tasks. The resulting classification accuracy of the proposed approach has been found convincing. The same signals were also classified by the conventional classifiers from the temporal features of the signals. From the comparisons of the classification accuracies by the conventional and the proposed approach, it has been found that the proposed method provides accuracies slightly lower than the ANN and SVM but it provides a better result than that of the kNN and LDA. Another benefit of the application of the proposed method is its non-necessity of the training phase. The matrices of the polynomial coefficients (see (17) and (18)) regarding the propose activation models can be used as the initial marker to compare a signal to find its class. So, it is easier to implement the proposed classification methods for the voluntary and imagery movement-related task classification from the prefrontal hemodynamics. Since this work has found the very convincing classification accuracy from the prefrontal hemodynamics regarding voluntary and imagery movements, this result will play a great role to implement the BCI for the paralyzed people as well as for those people has a major injury in their central lobe of the brain. Therefore, the proposed research work contributes significantly to model and classify the voluntary and imagery movements from the prefrontal hemodynamics.
As all the participants of this study were normal, the actual scenarios of the physically challenged or paralyzed persons were not taken into this research which can be considered as the limitation of this study. Nonetheless, this limitation may be a future direction of research to find the actual scenario of motor imagery activation in the prefrontal cortex regarding the physically challenged persons. In our future work, aforesaid undone work will be performed to present a comparative scenario of voluntary and imagery movement-related activations of the prefrontal cortex concerning both the normal and physically challenged persons.