Neural Network Classification of Brain Hemodynamic Responses from Four Mental Tasks

Abstract

We investigate subjects’ brain hemodynamic activities during mental tasks using a nearinfrared spectroscopy. A wavelet and neural network-based methodology is presented for recognition of brain hemodynamic responses. The recognition is performed by a single layer neural network classifier according to a backpropagation algorithm with two error minimizing techniques. The performance of the classifier varied depending on the neural network model, but the performance was usually at least 90%. The classifier usually converged faster and attained a somewhat greater level of performance when an input was presented with only relevant features. The overall classification rate was higher than 94%. The study demonstrates the accurate classifiablity of human brain hemodynamic useful in various brain studies.

Keywords:

• brain-computer interface
• functional near infrared spectroscopy
• mental task classification
• neural networks
• wavelet transforms

NOMENCLATURE

A	=	light attenuation measured in optical density
b	=	translation parameter of wavelet; bias parameter in a neural network
C	=	concentration of the absorbing compound (Oxy-Hb, deOxy-Hb or CytO) measured in μ molar
E	=	squared error function of a neural network
G	=	additive term reflecting the light scattering loss
I	=	transmitted light intensity measured by the fNIRS
k	=	localization index of wavelet transform, a
N c	=	number of correct classifications
N r	=	number of rejections
P e	=	probabilities of erroneous classification
R(.)	=	cross–correlation parameter
s(.)	=	noise free signal in x(n)
	=	universal threshold value applied to wavelet scales
w x (.)	=	wavelet transform of an observed signal x(n)
w s ()	=	noise free wavelet coefficients of an analyzed signal
x(.)	=	noisy model signal of hemodynamic response n = 1, N ∈ R
x .	=	Input–output vector of neural network classifier
α	=	specific extinction coefficient of the absorbing compound
a	=	wavelet scaling parameter
B	=	differential pathlength factor
d	=	interoptode distance defined in cm
g(.)	=	dditive noise term in x(n)
I o	=	medium light intensity incident given by the light of fNIRS
j	=	scale index of wavelet transform a
L	=	layers of neurons in neural network structure
N c	=	number of erroneous classification
ψ a, b (.)	=	wavelet function used for signal decomposition
P c	=	probabilities of correct classification
σ	=	noise level of a signal in wavelet domain
t	=	time sample
w x (.)	=	wavelet transform of an observed signal x(.)
w g (.)	=	noisy wavelet coefficients of an analyzed signal
ΔW	=	neural network weight parameter
	=	data set with n instances
y i	=	desired output vector in classification task
η	=	neural network learning rate parameter

1. INTRODUCTION

Recent advances in neuroimaging demonstrate a new way of accessing the brain's functional state by using a functional near infrared spectroscopy (fNIRS). It is an emerging noninvasive, safe and portable sensing modality for monitoring physiological changes that occur in the human brain (Izzetoglu et al. 2005). fNIRS presents the information about cortical hemodynamics and oxygenation status during functional activity through three parameters such as oxyhemoglobin, deoxydemoglobin, and total hemoglobin. The ratio of these concentrations is determined by a combination of oxygen consumption, supply of oxygenated arterial blood flow, and drainage of de-oxygenated venous blood (Elwell 1995). The concentrations of these parameters may provide an important information of subjects’ cognitive state for studying the brain anatomy and activity. One of the emerging application of fNIRS is in brain computer interface (BCI) technology (Matthews et al. 2008). For instance, in BCI the brain signals are studied to derive useful information about the users’ mental state to map it into an external action to control devices (Gerven et al. 2009). Many studies already propose the use of the fNIRS as a practical approach to realize the BCI (Coyle et al. 2007; Ayaz et al. 2007). For instance, recently, Sitaram et al. (2007) demonstrated the feasibility of decoding brain hemodynamics arising from right and left hand motor imagery. In addition, Utsugi et al. (2007) achieved real-time fNIRS-based toy train control by performing complex mental arithmetic tasks. Sassaroli et al. (2008) showed that simple k-means algorithm satisfies successful classification brain hemodynamics from five mental tasks. Truong et al. (2009) studied a neural network approach for the classification of brain hemodynamics with wavelet transform input features for BCI. Moreover, the brain hemodynamics and oxygenation have been studied to identify the brain areas for human speech recognition (Sato et al. 1999), and for brain disorders such as epilepsy, posttraumatic, and mood disorders (Watanabe et al. 2000; Matsuo et al. 2003).

These studies show that fNIRS is an interesting alternative approach to other widely used neuroimaging techniques in brain studies. Because of its simple and economical instrumentation and more efficient techniques for long-term recording, fNIRS also offers accurate and good signal-to-noise ratio with detailed information of neural activity in high spatial resolution. Examples of widely used neuroimaging techniques are brain electroencephalography (EEG) and magnetoencephalography (MEG) measured from scalpelectrodes or blood–oxygen-level dependent functional magnetic resonance imaging (fMRI) as well as positron emission tomography (PET). The techniques, such as MEG, fMRI, and PET, provide good spatial and temporal resolutions; however, they require heavy and expensive instrumentation under restricted conditions to utilize. In contrast, EEG devices are cheap and portabl. They also provide noninvasive recordings which are safe and easy to use in all environmental conditions. Coffey et al. (2010) provide a detailed table comparing the characteristics of fMRI, MEG, EEG, and fNIRS in terms of principle of operation (e.g., hemodynamics or electromagnetic), signal characteristics (e.g., temporal and spatial resolution), and portability and comfort. In their comparison table, Coffey et al. (2010) arrive at the conclusion that EEG and fNIRS are better BCI candidates for space applications. Therefore, other techniques are not suitable for evaluating the brain activity of subjects under normal working conditions. There are many EEG-based BCI studies that demonstrate its applicability in controlling devices (McFarland et al. 2006). Unfortunately, most independent EEG-based BCI systems' performance are limited to only a few commands due to the EEG's poor signal-to-noise ratio and spatial resolution.

In this study, we present our preliminary research results in classification of mental task relevant brain oxygenation changes with the aim to decode them for BCI. The methodology presented in this study is shown in Figure 1. It consists of the following three major steps. First, four healthy participants’ data are acquired while performing several mental tasks. Second, we apply continuous wavelet transforms (CWT) for preprocessing and feature extraction. In this stage, we test four wavelet functions and propose best basis function. Subsequently, we employ a soft thresholding technique to further reduce irrelevant (noisy) features of the signal. The important features are extracted based on wavelet statistics. Finally, a single-layer neural network classifier is trained and tested with extracted wavelet features of each mental task. We train and test various neural network models to obtain the best classifier by employing two error minimizing techniques of a backpropagation algorithm. The output of the classifier consists of four decoded mental tasks associated with brain oxygenation.

Figure 1 Flowchart of the methodology (color figure available online).

2. FUNCTIONAL NEAR-INFRARED SPECTROSCOPY

An fNIRS measures the physiological changes associated with brain activity by exploiting the optical properties of the brain tissue. Near infrared light in the range of 650–950 nm can pass through the skull and reach the cerebral cortex up to a depth of 3 cm (Elwell 1995). The light is attenuated due to a combination of absorption and scattering as it passes through the brain tissue. There are three types of varying chromophores in the brain that are present in variable concentrations: oxy-haemoglobin (Oxy-Hb), deoxy-haemoglobin (deOxy-Hb), and cytochrome oxidise (CytO). Any change in observed light attenuation is due to a change in the concentration of these chromophores. The quantitative assessment of the concentration of absorbing chromophores can be obtained using the modified Beer-Lambert law. To quantify cerebral oxygenation, the optical path length of the photons traveled through the tissue should be known and remain constant. The distance of emitter-detector optodes is known. The lights are scattered by tissue components and diffused through the brain. The modified Lambert-Beer Law in equation (1) allows conversion changes in absorbtion and attenuation in changes in concentration of the chromophores (Elwell 1995).

Here, A is the light attenuation measured in optical density; I ₀ is the medium light intensity incident given by the light of a fNIRS; and I is the transmitted light intensity measured by the fNIRS. An α is the specific extinction coefficientof the absorbing compound measured in μ molar⁻¹ cm⁻¹, and this is well known and shown in Figure 2 a. C is the concentration of the absorbing compound (Oxy-Hb, deOxy-Hb or CytO) in the solution measured in μ molar; d is the interoptode distance defined in cm; and B is the differential pathlength, factor. Because light scatters through tissue, d has to be multiplied by B to find the differential pathlength, which is the true optical distance. Here, G is the additive term reflecting scattering loss. Part of the light will be scattered so that the light measured is not only changed due to absorption but also because of scattering. As we cannot measure the scattering, G remains an unknown factor and no absolute values canbe measured.

Figure 2 The locations and settings of the emitter-receiver optodes of fNIRS system on the frontal cortex of a subject. (a) 19 channels locations according to 10-20 electrode placement system; (b) oxyHb and deOxyHb reconsctruction using Modified Beer Lambert law; (c) absorbtion spectrum of hemoglobin (color figure available online).

The only values that can be measured are differential values comparing one moment with another. This will give us the following equations.

Because the specific extinction coefficient is different for the three chromophores, the different instruments will use different laser diodes emitting different wavelengths (see Figure 2 c) to calculate the proportional concentration of the cerebral oxygenation. The amount of wavelength used will have an influence on the precision of the measurement.

3. EXPERIMENTAL DATA ACQUISITION

In our experiment, we used a multichannel optical brain-function imaging system FOIRE-3000 by Shimadzu Co. Ltd., Kyoto, Japan. The system uses safe near-infrared light to assess Oxy-Hb and deOxy-Hb contents on the surface of the brain. The system uses three lasers at wavelengths of 780 nm, 805 nm, and 830 nm. The sampling rate of the acquired data is 10 Hz. We measured the concentration levels of Oxy-Hb, deOxy-Hb, and total haemoglobin (total-Hb) from four healthy right-handed subjects (mean ages 22–29 years old). Nineteen channel recordings have been acquired from the frontal lobe of the prefrontal cortex according to the 10–20 electrode placement system. Figure 2 a shows the locations of the emitter-detector probes, where the red circled numbers represent the emitter and the blue circled numbers represent the detector probes. The corresponding fNIRS 19 channels are shown as numbers in white squares, and the interoptode distance is fixed to 30 mm. Figure 3 demonstrates our data acquisition protocol in pre-rest, and task, and post-rest paradigm. We set the onset and offset points of a specific mental task by referring to the time points.

Figure 3 Experimental setup for mental task recording; after the pre-task rest, an audio tone notifies the participant to prepare for the task and in 3 seconds mental tasks are provided (color figure available online).

Four mental tasks were prepared in our work: baseline, object rotation, math, and verbal fluency. The baseline state was recorded while a participant rested as a reference to compare with other mental tasks. In the object rotation task, several objects were presented to the participants for imaginary rotation. In addition, we provided tasks in the form of puzzles so that their brain activity could be easily assessed; During the math task, subjects were asked to multiply two random numbers as fast as possible; for instance, multiply 18 × 87. The math task was designed as such that it was difficult but could be accomplished within the given time segment. The verbal fluency task consisted of constructing a word with a randomly determined first letter. For instance, the operator would say the letter a [ei], and the participant would have to pronounce words starting with the same letter (e.g., apple, ant, etc.).

Before the data measurement session, the participants were asked to sit, relax, and rest in order to stabilize the blood flow in all 19 channels. We started the recordings after the haemodynamic concentration of the blood flow was normalized. Our data acquisition were arranged into sessions where in a single session, a participant performs the same mental task three times. After the data acquisition, we visually inspected for possible motion or experimental noises. We excluded all possible sessions where the data were suspect to noise (i.e., we, sorted our data from only three potential sessions), and there were no fixed number of sessions. The data was averaged at the end of each session from each three trials. Thus, we had the recordings of four mental tasks performed by four participants. In total, 48 data sets were collected, which is made of 12 data sets from each subject for four mental tasks.

In this study, we focused on the analysis of Oxy-Hb concentration levels because it was found to be constant in monitoring of neural activities. A study reports that the increase in the nervous activity results in an increase of consumption of local oxygen. As a result, oxygenated hemoglobin decreases and deoxygenated hemoglobin increases. Further, the blood vessel is enhanced to supply fresh blood. Then, the blood flow volume increases locally. So, oxy-Hb increases and deoxy-Hb decreases. These processes are reactions that occur within a few seconds (Obrig and Villringer 2003). The Blood-oxygen-level dependence (BOLD) signal is an index of the blood stream change in measurement by fMRI. There is a report that oxy-Hb strongly has the correlation with the BOLD signal (Strangman et al. 2002). To evaluate the functional brain state, we may mainly judge activation from the increase of oxy-Hb. Therefore, our focus in this study is in the analysis of only Oxy-Hb concentrations.

3.1. Signal Pre-Processing

One of the challenges in the accurate analysis of hemodynamic signals is the noise interference. The noise interferences in hemodynamic signals may arise from instrumentation, experimental setup, or physiological sources. We present a wavelet-based approach; in particular, focused on a soft-thresholding technique to separate a true hemodynamic response from background noise. The soft-thresholding technique was originally proposed by Donoho and Johnstone (1994; 1995), who studied wavelet denoising techniques. They proposed a remarkably good wavelet denoising method which has been shown to be a better denoising method than other classical methods, such as MiniMax, Hybrid method, and SURE method (Johnstone and Silverman 1997). Their technique is based on thresholding the wavelet coefficients of the original hemodynamic signal. Let s(n) be the noise-free signal and x(n) is a signal corrupted by non-stationary Gaussian noise g(n), which is defined as

where g(n) is a zero-mean Gaussian noise with time-varying variance σ and is independent of the noise-free signal s(n). It is s(n) that we wish to estimate from the noisy signal x(n) as accurately as possible. Let the wavelet transform of the observed signal be denoted by w _x (j, k) = w{x(t)}_, where j = 1, 2, …, N is a scale index and k ∈ Z is a localization index. It follows

where w _x (j, k) is the noisy wavelet coefficients, w _s (j, k) is the noise-free wavelet coefficients, and w _g (j, k) is the noise in the wavelet domain which is assumed to be a zero-mean Gaussian random variable. The continuous wavelet transform (CWT) of a signal x(t) is obtained by convolution of the signal time series x(n) with the mother wavelets ψ _{a, b} (t):

where the parameters a and b are the scale and translation, respectively (Mallat 1999). As a increases, the wavelet becomes more narrow and captures the low frequency contents of the signal, and by varying b the mother wavelet is displaced in time. Thus, CWT provides us with in depth analysis of the signal in multi-scale representation. We first perform dyadic CWT with different wavelet functions. Next, we apply the following soft thresholding technique to each scale to separate the neural signal from noise.

Here, t is the universal threshold estimation proposed by Donoho and Jonstone (1994; 1995) which is defined as,

where N is the total number of samples of t _{a, N} . Here, σ = median (w _x (j, k)/0.6745) is the noise level of a signal, which is computed at each scale. The soft thresholding method is applied to the individual scales a _j,k .

3.2. Mother Wavelet Selection

Notice that in equation (6), we have to define the type of the mother wavelet for a signal decomposition. All three wavelet analysis parameters have several or even an infinite number of options. The parameters are closely related and must be considered at the same time. The wavelet type should be set according to the signal being analyzed. The general selection of wavelets is based on the shape of the wavelet in time domain, its length (support), and smoothness. Equation (5) shows that the CWT is essentially a collection of dot product of a signal s(t), and the translated and scaled wavelet basis functions ψ(t) for all a and b, . For fixed scale a ₀ and translation b ₀ parameters the wavelet transform of the function f(t) is denoted by w(a ₀, b ₀), which represents its similarity index to the wavelet function. In other words, we can represent wavelet transform as a cross-correlation measure given as . The maximum cross-correlation will provide only intrinsic structures of the signal represented in a few coefficients. Thus, selecting a suitable wavelet function becomes an important condition for extractingtrue neural signals from background noise. Figure 4 shows the four candidate wavelets for hemodynamic signal analysis. These functions are Mexican hat (Mexh), Daubechies order 8 wavelet (Db8), and Biorthogonal wavelets (Bior1.3) and (Bior.15). Our criterion for choosing wavelet functions was based on the shape of wavelet and its similarity to the patterns of the original hemodynamic transients. Mexican hat wavelet function has a simple expression and it is piecewise smooth. It has good localization in both time and frequency domains with a minimum number of oscillations. The following Biorthogonal wavelets are an important class of functions. They provide generalization to orthogonality property that also offers more flexibility in the construction (Mallat 1999). Daubechies wavelets are also orthogonal functions with compact support of 2N-1, where N is an order of wavelet function.

Figure 4 Selected 4 wavelet functions for CWT decompositions; (a) mexh, (b) db8, (c) bior1.3, and (d) bior1.5.

3.3. Feature Selection

In general, feature selection aims to reduce the dimensionality of pattern for classification by selecting the most informative rather than irrelevant and/or redundant features. For our case, the important features of hemodynamic responses are located in low and high frequency regions of CWT multi-scale decompositions. Because CWT coefficients provide clear insight into the structure of a signal, we must extract a comparatively small number of features that are maximally informative from varying cognitive conditions across CWT decompositions. We have defined five potential features based on the statistics of the wavelet coefficients studied by Tzanetakis et al. (2001) and the highest amplitude value of the signal in wavelet domain. The following statistical features define the time-frequency distribution and the amplitude of the mental tasks.

	Mean of the absolute values of the coefficients in each scale.
	Average power of wavelet coefficients in each scale.
	Standard deviation of the coefficients in each scale.
	Ratio of absolute mean values of adjacent scales.
	Amplitude of the hemodynamic concentration in wavelet domain.

Here, features 1 and 2 represent the frequency distribution of the hemodynamic signal, and features 3 and 4 the amount of changes in frequency distribution. Finally, the amplitude level of the hemodynamic signal, was considered in addition to statistics of wavelet coefficients.

4. NEURAL NETWORK CLASSIFICATION

We perform a neural network-based classification of mental tasks, followed by wavelet pre-processing and feature selection. The formulation of the classification task is as follows. Suppose we are given a data set (feature set) with n instances, (x ₁, y ₁) ⁿ , …, (x _i , y _i ) ⁿ  ∈ X where X-is a nonempty set which contains the mental task features, and y _i  ∈ [0, 1] ^m is the desired output vector where m is set to the number of classes in classification task. We will consider a multilayer neural network with L = 1 layers of neurons and L + 1 layers of synaptic connections, as shown in Figure 1 c. Mathematically, the neural network is given as

where x ^l  ∈ R ∀ l = 0, …, L and W ^l is and n _l  × n _l − 1 matrix ∀ l = 0, …, L. Our task is to determine the weights W and biases b so that the actual output x ^L becomes closer to the desired output y. We apply back-propagation learning algorithm which consists of the following steps (Bishop 1997): 1) The input vector x ⁰ is transformed into the output vector x ^L , by evaluating the equation for l = 1 to L; 2) the difference between the target output d _i and the computed output x ^L is evaluated by an error function ; 3) the error signal at the output nodes is propagated backwards through the entire network; and 4) the synaptic weights and biases are updated using the results of the forward and backward passes by and . Notice that the last step can be performed by using variations of gradient descent algorithms; namely, the network parameters are continually changed in the direction of steepest descent with respect to error. In this study, we train a backpropagation neural network (BPNN) with two algorithms, such as gradient descent and scaled conjugate gradient algorithms. The gradient descent algorithm is one of the line search minimization procedures. This method smoothes the descent direction in the steepest decent method (Haykin 1999). The weights and biases are updated in the direction of the negative gradient of the performance function. The scaled conjugate gradient algorithm was developed by Moeller (1993) and used to avoid thetime consuming line search. It combines the Model-trust region approach with the conjugate gradient approach. Then, backpropagation is used to calculate the derivation of performance with respect to the weight and bias variables.

Training BPNN involves other parameters except weight and bias. Particularly, we defined a learning parameter which represents the percentage of the step taken towards minimum error (e.g., a too small learning rate produces a too long training time and vice versa). In addition, initial weights, number of neurons in a hidden layers, and the activation functions at each neuron were defined. The number of epochs needed to achieve a minimum error rate was also defined. Furthermore, we determined the optimal number of neurons in the hidden layer and number of epochs. The number of input nodes was determined by the size or dimension of input feature vector. The activation function used here was a sigmoid function given as . The details of BPNN parameters are given in Table 1.

Table 1. Parameters used in the training of a neural network classifier

Parameters	Properties
The number of neurons in the input layer	30
The number of hidden layers	1
The number of neurons in the hidden layer	10
The number of neurons in the ouput layer	4
Learning rate	0.01
Momentum parameter	0.003
The range of initial weight parameters	[−0.1, −0.0001] & [+0.1, +0.001]
Learning algorithm	GDA & SCG
Transfer function	Logarithmic sigmoid

The input feature vector consisted of five wavelet features of six channel, making a total of 30 input features. The target vector is specified for each of the network input features. For instance, for the math task, the target is specified as 0100. Next for the verbal fluency task the target is specified as 0010. Similarly, for each mental task the target is specified such that the particular task's target is classified as 1, and those of the other tasks are 0. There are four neurons in the output that corresponds to four mental tasks, as given in equation (10).

We perform n = 20–fold cross validation by dividing the input dataset into 20 disjoint approximately equal proportions. One proportion is used for testing, and the other n−1 proportions are used for training alternatively in the total n rounds ofclassifications. We evaluate the classifiers’ superiority and reproducibility on the basis of the probabilities of correct classification P _c and erroneous classification P _e , given as

where N _c is the number of correct classifications, N _c is the number of erroneous classifications, and N _r is the number of rejections. The validation was performed across the signal classes, which were indicated by the index of the highest output of the classifier. We repeated this validation several times, using a different train set each time untill all possible choices for the test set were eliminated.

5. EXPERIMENTAL RESULTS

At first, the brain activation regions with changes in the levels of oxy-Hb, deoxy-Hb, and total Hb were consistently observed at each channel in all subjects. It was noticed that the level of activation depends on the task difficulty. Higher activation responses are obtained when a subject considers a task difficult, and lower activation responses are obtained with easy tasks. For example, the math and verbal fluency tasks were considered difficult, while imaginary object rotations were considered easy. Figure 5a illustrates 19 channel hemodynamic concentrations during the math task. The pre-rest, task, and post-rest segments are shown divided by a vertical line; wherein, the task time is 10–25 time intervals in seconds. The relative changes in the oxy-Hb, deoxy-Hb, and total Hb concentrations can be observed at each corresponding channel. Significant changes in the oxy-Hb level during the mental task are seen in almost all channels. However, compared to other channels, 8, 9, and 17 show higher activation regions of the frontal cortex. Let's consider the other mental task during the mental rotation, as illustrated in Figure 5b . Compared to the math task, the hemodynamic concentration levels are small in some channels. As was mentioned, the difference in response and stimulus representation is because of the difference in task difficulty. Obviously, subjects do not tend to overload their minds to perform an easy task. Another difficult task could be the composition of a list of words beginning with a randomly assigned letter (especially for non-native speakers, as were our subjects). Regardless of change in the level of oxy-Hb concentration during each mental task, it is possible to extract task relevant stimulus features by using the wavelet transforms. Let's demonstrate some of the selected channels for multiscale wavelet decompositions.

Figure 5 Example of 19-channel hemodynamic concentrations during two mental tasks, a) Math task and b) object rotation task. Vertical lines represent [prerest – task – postrest] periods (color figure available online).

Consider a channel 17 of math task shown in the plane of Figure 6 a. During the task, the increase in oxy-Hb and total Hb is clearly seen; however, the signal tends to be noisy. Next, we obtained thresholded wavelet coefficients of the signal using a Mexican hat wavelet, as seen in Figure 6b . This figure demonstrates the time-scale representation of the signal. We may notice that this type of representation provides localized signal features across low and high frequency contents. Some features of the signal present at high frequency may not be detected at low frequencies. For instance, Figure 6b at time samples 380, shows most neural activity is located at a low frequency band. In contrast, around time samples 600 the neural activity is probably located at a higher frequency. Such ability of wavelet transform will enable us to combine important signal features from multi-scale representation for accurate classification. In principle, wavelet transform operates as a zooming function by computing the dot product between the signal with the scaled and shifted version of a wavelet function. With increase in a, the wavelet function becomes more narrow and captures the low frequency contents of the signal, and by changing b the mother wavelet is displaced in time. More intuitive visualization of wavelet transforms can be seen in Figure 7. The figure illustrates multi-scale decomposition of the same signal from Figure 6 by using four different wavelet functions. Scales (a2, a4) capture high frequency contents and scales (a10, a12) capture low frequency contents. We notice that wavelet representation of the signal varies depending on the type of wavelet function. As an example, Figure 7d shows that mexh wavelet mainly localizes low frequency signal contents while bior1.3 and bior1.5 capture high frequency signal characteristics in Figures 7a and 7b. With different wavelets we obtain different signal representation in wavelet domain. It is usually difficult to find a suitable wavelet for a particular signal analysis task because an unlimited number of wavelet functions exist in the wavelet library. We will further analyze the classification results of a neural network classifier based upon all these wavelet input features.

Figure 6 Analysis of hemodynamic concentration from math task, a) details of channel 17, oxy-Hb b) Time-scale wavelet decomposition of the signal using Mexican hat wavelet (color figure available online).

Figure 7 Comparative visualization of four selected wavelet functions in decomposition of oxy-Hb signal from figure 6; a) bior1.3 b) bior1.5, c) db8 and d) mexh wavelet functions (color figure available online).

We started neural network simulations with a 2-layer neural network model. Our neural network had single hidden-layer neurons with sigmoid transfer functions. The number of input nodes was fixed as equal to the number of input feature vectors, and the number of output vectors was fixed to 4, which equals the number of output classes. We repeated the network simulation to determine the parameters related to the optimal number of neurons in the hidden layer and training algorithm. The changes in the overall classification performance of neural networks were confirmed with varying numbers of hidden-layer neurons for a particular task. For instance, the results are shown in Table 2 and the classification rate changes with varying hidden neurons when the SCG algorithm is tested. Here, the input features consisted of wavelet coefficients using mexh function. The optimum number of neurons required in the hidden layer was found to be 10 and, hence, we selected the neural network configuration 30-10-4 (number of input nodes, neurons, and output vectors, respectively). The best performance was obtained for the training set; the validation test set with models whose hidden layer had two or more neurons. We tested the classifier with independent wavelet features of each wavelet function to investigate the classification accuracy.

Table 2. Performance of the various neural network architectures tested for the mental tasks

Model no	Network architecture	No. of. epochs	Training ime	Rate (%)
1	30-10-4	200	1.00	97.33
2	30-7-4	251	0.642	94.09
3	30-5-4	321	0.766	93.67
4	30-4-4	385	0.721	93
5	30-2-4	1202	1.903	85.1
6	30-1-4	5786	8.62	71.88

Figure 8 shows the comparative classification results obtained using the wavelet input features. The mexh wavelet offers better accuracy than the others, and bior1.3 is marginally better than db8. Therefore, the mexh and bior1.3 wavelets were proved to perform better and were selected for the application. In addition, we can see the good classification rates of even the other proposed wavelets; they achieved an accuracy of 90% and higher.

Figure 8 Mental task classification results of neural network with four different wavelet input features from five individual scales (color figure available online).

Following, both GDA and SCG algorithms are applied to obtain the classification performance of test set. Table 3 provides the classification rate of SCG and GDA algorithms with mexh wavelet input features. We can see that the SCG algorithm is a better compromise for this specific application. Further, overall classification rates of all subjects by using four wavelet function input features to a SCG based neural network classifier is illustrated in Figure 9. The classification rate varies depending on the wavelet features and the subject's mental workload over the same task. Besides, it should be noted that a neural network with Mexican hat wavelet transform features provides the most accurate classification accuracy than the other studied wavelets. The variation of the results from our method shown in Figure 9 may depend on at least two factors.

Figure 9 Overall classification results of 4 mental tasks from all subjects using the selected wavelet input features; a) bpnn with mexh input features, b) db8 input features, c) bior1.3 input features and d) bior1.5 input features. (color figure available online).

Table 3. Performance of SCG and GDA [30-10-4] neural network model in the classification of Mexh wavelet features from a test set

Method	Rest	Math	Verbal fluency	Rotation
SCG	91.2 (%)	92 (%)	94.2 (%)	89 (%)
GDA	85.5 (%)	90.4 (%)	87 (%)	86 (%)

One obvious one is the different brain activation levels when four subjects perform the same task, e.g., math. The other one remains in the classifier's performance. For instance, there's a possibility that subsets of two different input feature set may overlap; as a result, a classifier may fail in the decision making of these subsets properly. As mentioned earlier, the results obtained demonstrate that designing a task is also important to achieve a good neural activation. As an example, in the plane of Figure 9a math and verbal fluency tasks classified more accurately. Similarly, the other planes of the figure show the same outcome except some cases as in Figure 9c . When the bior1.5 function was used, it was found to be stable in the feature extraction of all subjects for one task which resulted in stable accuracy on average. However, the lowest classification rates goes to 82.3%, which is less compared to mexh features in Figures 9a or 9c, which is, about 84%. However, we can find some trade-offs by using these functions.

6. DISCUSSION

The biophysical characteristics of hemodynamic transients have been studied earlier by Elwell et al. (1999) and Obrig et al. (2000). Hemodynamic transients possess slow and smooth varying waves, which could somehow be characterized as wavelets. Coyle et al. (2004) explained the biophysical characteristics of the noisy hemodynamic response. Their assumption was based on the fact that the hemodynamic wave is quasi-periodic and can be represented as a sine wave. Further, hemodynamic changes are not stationary because it may actually be difficult to represent the signal with sine waves. In contrast wavelets provide a better localization of hemodynamic transients, which to our knowledge represent pure brain hemodynamic stimuli.

The use of wavelet transforms in the fNIRS signal analysis has been conducted earlier by Emir et al. 2003. They first applied wavelet denoising methods for fNIRS signals. Their empirical results indicate that wavelet denoising may fail in some circumstances compared to independent component analysis (ICA) methods. However, we have shown that with proper design of the wavelet algorithm it could be an effective filtering tool. Truong (2009) used discrete wavelet transforms for feature extraction. Their method used a similar approach, i.e., the use of wavelet coefficients as an input to a neural network classifier. The authors didn't consider the noisy coefficients of the signal in the wavelet domain instead relying, on a standard wavelet bank of filters. However, they have received relatively good accuracy in classification, even with raw wavelet coefficients as an input to a classifier. Compared to studies by Emir et al. (2003) and Truong et al. (2009), we use CWT and a denoising technique to eliminate the unwanted noisy components of the signal in wavelet domain. This actually increases the performance of the classifier. For example, classifiers may yield misleading results if adequately relevant features of the brain phenomena are not provided. In other words, irrelevant features may lead to garbage-in garbage-out types of results from the classifier with impractical solutions. One solution to this problem is to extract most of the relevant and clean features from the data; then, most types of classifiers will perform accurately. Aiming at efficient feature extraction and a classification framework for hemodynamic-based BCI, we have done our efforts on the following attributes. 1) We explore brain hemodynamic signals through overcomplete multi-scale decompositions or CWT. One of the advantages of CWT over discrete time wavelet transform (DWT) is its ability to compute the coarsest levels of the signal, which is crucial in the analysis of physiological signals where each signal component could be important information to retain. Even though DWT is computationaly efficient, there's a possibility to lose that important information while passing the signal through a number of filters. Besides, it also makes it difficult to look for the information in transform domain that makes use of the signal's approximations and details. In contrast, CWT transforms the signal into one time-scale domain, where all information can be searched by only measuring the correlation of the mother wavelet function and the analyzed signal. 2) We discuss the mother wavelet type that should be set according to the signal being analyzed. The general selection of wavelets is based on the shape of the wavelet in time domain, its length (support), and smoothness. Their common description can be found in wavelet literature. Four selected wavelet functions are based on the shape of the signal and smoothness. Our experimental results suggested that the Mexican hat function could match the signals’ intrinsic characteristics better than the other wavelet functions (e.g., Figure 8). 3) Further, we apply a universal thresholding technique to eliminate of noisy coefficients of the signal in the wavelet domain. We do not achieve inverse transform of the signal from the wavelet domain; instead, we try to get signal approximations on it. As a result, it makes our algorithm computational less. 4) Finally, for classification, we explore the optimal neural network structure for the given problem. Further, we demonstrate two neural network learning algorithms’ performance, GDA and SCG, trying to find a fast algorithm to generalize to the given problem. SCG performed fast in convergence, and provided better classification performance than the standard GDA algorithm.

ACKNOWLEDGEMENTS

This work was supported by the DGIST R&D Program of the Ministry of Education, Science and Technology of Korea (11-RS-01).