2.1.

Participants and Experimental Protocols

Eight healthy participants participated in the study (seven males and one female, age: 26.8 ± 1.6 years). None of the participants had histories of neurological or psychiatric issues. All the eight participants performed the experiments related to the four-class classification while only three participants among them were able to finish the experiments to determine the task duration. The experiments were approved by the Institutional Review Board of Korea University (KU-IRB-10-34-A-2). Written consent was obtained from all participants. Considering quaternary classification, it is necessary to obtain the hemodynamic responses evoked by four kinds of distinct behaviors. Since lifting up and down an arm or leg are more active movements than finger tapping, we adopted the active movements in our experiment. Figure 1 indicates the graphical motions of the active movements. Hemodynamic responses evoked by the motor tasks were collected from each participant while performing overt motor execution. Figure 2 shows the experimental design for motor execution. The task consisted of 40 cycles of 6, 8, 10, or 15-s motor task and rest. The durations of the task period ($T_t$) and rest period ($R_t$) were kept the same. The participant sat on a chair in front of a screen displaying instructions and the progress of the experiment. In the prewaiting period, a red light was displayed. At the beginning of the task period, a green light was turned on. During this period, self-paced motor execution was performed. In the rest period, the color of the light was changed from green to red. The participant was asked not to think or move while staying in the rest state. The green and red lights were turned on and off repeatedly in the task and rest periods. In the postwaiting period, a yellow light was displayed and a red light was turned on when all the periods were finished.

Fig. 1

Graphical motions of (a) arm lifting and (b) knee extension. Participant sits on a chair and lifts up and down his/her arm or leg.

Fig. 2

Experimental design of motor tasks for the measurement of hemodynamic responses using functional near-infrared spectroscopy (fNIRS). Each experiment was composed of a prewaiting time ($T_t+R_t$ seconds), repetitions of rest periods for $R_t$ seconds and task periods for $T_t$ seconds, and postwaiting time ($T_t+R_t$ seconds). Hemodynamic responses evoked by four kinds of motor tasks (left/right arm lifting and knee extension) were measured separately.

A single experimental run was composed of a prewaiting period, 40 repetitions of rest and task periods, followed by a postwaiting period. Each participant first completed 20 repetitions of the rest and task periods in the first half. After a short rest, the participant completed the same number of repetitions again in the last half. Each half experimental run (20 repetitions) lasted 264, 352, 440, and 660 s for 6, 8, 10, and 15 s per one task period ($T_t$), respectively. Each participant conducted a single task repeatedly in each experimental run. Hence, four sets of experimental runs (i.e., for left/right arm lifting and knee extension) were completed for each participant. The tasks were performed in the order of left arm lifting, right arm lifting, left knee extension, and right knee extension. Between each experimental run, the participants took a rest for a minute. Each participant was allowed to do some stretching and speak lightly, but was not allowed to leave the seat during the rest period.

2.2.

System Setup and Data Processing

All windows of the experiment room were covered by blackout curtains, and lights were turned off to keep out ambient light. An eight-channel custom-built continuous-wave fNIRS system was used to measure the hemodynamic responses. The fNIRS system uses LEDs operated at two different wavelengths of 670 nm (L6112-01, Hamamatsu photonics, Hamamatsu, Shizuoka, Japan) and 890 nm (L2656-03, Hamamatsu) which are similar to the wavelength of LEDs used by Coyle et al., because the fNIRS system used in this investigation was implemented based on their system.¹¹ The LEDs were placed in direct contact with the scalp. The lights were modulated with 4- and 5-kHz sinusoidal waves. The detectors used in the system were APDs (C5460-01, Hamamatsu), which collected the lights via 3-mm diameter multimode optical fibers. The detected light signals were first demodulated by a lock-in amplifier implemented by LabVIEW software (National Instruments, Austin, Texas). The data were recorded by a 16-bit data acquisition card (PXI-8106, National Instruments) and downsampled with a sampling frequency of 12.5 Hz. The whole system architecture is schematically illustrated in Fig. 3.

Fig. 3

Block diagram of an fNIRS-based brain–computer interface (BCI) system. 670- and 890-nm LEDs are modulated at 4 and 5 kHz, respectively. Avalanche photodiodes collect the lights, which are demodulated by lock-in amplifiers. The data are recorded by a 16-bit data acquisition card. Normalization, band-pass filtering, and Savitzky–Golay smoothing are applied for signal processing.

Figure 4 shows the arrangement of the fNIRS optodes. An fNIRS cap was made of elastic Lycra and shaped very similarly to a conventional EEG cap. The receiver optodes were fixed at a distance of 30 mm from the source optodes. The optodes were arranged on both sides of a hemisphere around the motor cortex areas around the $C_3$, $C_4$, and $C_z$ positions of the international 10 to 20 systems. The source and detector optode pairs consist of eight channel arrangements. The concentration changes of oxy-hemoglobin were calculated from the changes in raw optical intensity data received using the MBLL given by³¹

Fig. 4

Optode arrangements for fNIRS. Yellow and green circles represent the source (S1 and S2) and detector optodes (D1 to D7), respectively. The source and detector are spaced 30 mm apart. The optodes were placed around the motor cortex area according to the international 10 to 20 systems.

The calculated concentration changes of oxy-hemoglobin of the task and rest periods, $\mathbf{X}$, were normalized by subtracting the mean $E(\mathbf{X})$ of the amplitude of the signal in the whole periods of each channel, and dividing by the standard deviation $\sigma$ of the amplitude of the signal in the whole periods of each channel [i.e., $\{\mathbf{X}-E(\mathbf{X})\}/\sigma$]. We used a zero-phase distortion filter (third-order Chebyshev I with a ripple factor of 0.5 dB) with a passband of 0.03 to 0.07 Hz to remove the baseline drift and physiological noise because the fundamental frequency of the evoked hemodynamic response was expected to be in the passband. If we increase the bandwidth of the BPF, the evoked hemodynamic responses can be distorted by low frequency oscillation signals.^23,32 The correlation-based signal improvement method was used for minimizing the effect of head movements,¹⁴ and a third-order Savitzky–Golay smoothing filter was used to determine the smoothed values for each point.¹⁵ Figure 5 shows the measured hemodynamic responses ($\mathrm{\Delta }{C}_{\mathrm{HbO}}$) from all the participants evoked by left arm lifting. Participant 1 result shows the clear change of the hemodynamic responses because the participant has shaved his/her hair so as not to weaken the intensity of the light. Since participants 2 and 3 have thin hair, we could obtain the relatively distinct evoked hemodynamic responses. Figure 6 shows the time-frequency analysis results of $\mathrm{\Delta }{C}_{\mathrm{HbO}}$ measured from participant 1 using a continuous-time wavelet transform (Morlet wavelet). The fundamental frequency components of $\mathrm{\Delta }{C}_{\mathrm{HbO}}$ are explicitly presented.

Fig. 5

Hemodynamic responses ($\mathrm{\Delta }{C}_{\mathrm{HbO}}$) measured from eight participants evoked by left arm lifting ($T_t=R_t=8\mathrm{s}$). (a)–(h) $\mathrm{\Delta }{C}_{\mathrm{HbO}}$ from participants 1 to 8, respectively. Shaded regions represent task durations. Twenty motor task executions are repeated. Yellow, magenta, cyan, red, green, blue, and black lines indicate the hemodynamic responses measured at D1, D2, D3, D4, D5, D6, and D7 positions, respectively.

Fig. 6

Time-frequency analysis of the hemodynamics response ($\mathrm{\Delta }{C}_{\mathrm{HbO}}$ at D2 of participant 1) using a continuous wavelet transform (Morlet). Low (around 0.1 Hz) and very low (around 0.04 Hz) frequency oscillations are much weaker than the evoked hemodynamic response.

2.3.

Feature Space

Among (1) mean, (2) max, (3) slope, (4) delay, (5) var, and (6) median, several features were chosen individually to obtain the highest classification accuracy. Individually chosen features are presented in Table 1. We set the calculation window with a length of $T_t$ seconds to calculate the features after 8 s from the onset of a motor task. For example, when $T_t=15\mathrm{s}$, the features were calculated for 15 s after 8 s from the onset of a motor task. Forty-two features were calculated for each participant [$6\text{features}\times 7\text{detectors}\times 1\text{signal}$ ($\mathrm{\Delta }{C}_{\text{HbO}}$)]. To the best of our knowledge, this is the first investigation using the max, slope, and median in the four-class classification of fNIRS signals. Figure 7 shows the grand average of the spatial hemodynamic responses of all the participants evoked by the four kinds of motor tasks. The contralateral hemispheres were more activated by the motor tasks. Since the knee and leg areas of Homunculus are near $C_z$, the knee extension activated more around $C_z$ (D4 position) than the arm lifting.³³ From this physiological reason, we can classify the hemodynamic responses evoked by arm lifting and knee extension using the $\mathrm{\Delta }{C}_{\text{HbO}}$ at D4.

Table 1

Chosen features of fNIRS signals for classification.

Fig. 7

Grand average of the spatial hemodynamic responses ($\mathrm{\Delta }{C}_{\text{HbO}}$) of all the participants evoked by motor task of (a) left arm lifting, (b) right arm lifting, (c) left knee extension, and (d) right knee extension. (a)–(d) The relatively large $\mathrm{\Delta }{C}_{\text{HbO}}$ on the contralateral hemisphere channels. (c) and (d) The large $\mathrm{\Delta }{C}_{\text{HbO}}$ compared with (a) and (b). Knee extension evokes larger $\mathrm{\Delta }{C}_{\text{HbO}}$ than arm lifting at D4.

2.4.

Dimension Reduction Method

A high-dimensional feature set is difficult to handle, since high dimensionality may increase the noise in the data set, which may degrade the performance of the classifier. This is known as “the curse of dimensionality.”³⁴ To prevent this phenomenon, principal component analysis (PCA) and linear discriminant analysis (LDA) were used to reduce the dimension of each channel data.³⁵ PCA converts a set of the correlated data into a set of the uncorrelated components. LDA maximizes the ratio of between-class variance to within-class variance in the data set.³⁶ After that, we considered every possible dimensionality of the feature set, and the dimensionality of the feature set which shows the highest classification accuracy was found and then used. Dimension reduction was performed using the entire measured data. Figure 8 shows the signal processing flow of fNIRS signals.

Fig. 8

Signal processing flow of fNIRS signals. Dimensionalities of the data used in each step are represented. $N_{\text{data}}$, $M_{\mathrm{ch}}$, $N_{\text{smpl}}$, $K_{\text{feat}}$, $S_{\mathrm{set}}$, and $P_{\text{dim}}$ are the number of data points, channels, samples, features, and dimensions, respectively. $k_{\text{feat}}$ is the number of features selected for classification ($1\le k_{\text{feat}}\le K_{\text{feat}}$).

2.5.

Classification and Outlier Elimination

The NB classification based on the Bayesian theorem was used. NB assumes that all the features are independent. Based on the Bayes’ theorem, the sample data were classified under the maximum likelihood rule (ML rule). Outliers were eliminated to remove the unreliable data. Ten out of 40 sample data were selected as outliers. Outliers were the farthest data from the average $\mathrm{\Delta }{C}_{\text{HbO}}$ of the measured data. After excluding the outliers, 10 runs of the threefold cross-validation method were used to validate the results using the training and validation (test) subsets. Thirty samples were randomly divided into three subsets. Each subset consisted of 10 samples. Two subsets were designated as the training subsets and the remaining subset was used for validation of the trained classifier. In the training step, NB estimates the parameters of the probability distribution of the samples by using the training samples in the training subsets. In the classification step, using the test samples in the validation subsets, NB calculates the posterior probability of the sample and determines the class of the test samples in accordance with the maximum posterior probability (i.e., he ML rule).³⁷ Because NB only requires a small amount of training data to estimate the parameters necessary for classification, NB is extremely fast compared to other more sophisticated methods.³⁸

2.6.

Data Transfer Rate

The ITR was calculated in bits/trial as follows:²⁹

3.1.

Task Duration

Naseer and Hong³⁹ considered various task durations with constant rest durations based on a binary classification. We estimate the proper task duration to increase BPM while maintaining a quaternary classification accuracy of larger than 70%. Figure 9 shows the average quaternary classification accuracy of four different tasks for participants 1–3 as a function of task duration. The highest classification accuracies were observed at different task durations. For participant 1, a maximum classification accuracy of 91.5% was observed at a task duration of 8 s. At a task duration of 6 s, the classification accuracy was estimated to be below 70% for all the participants according to the threshold of the sufficient classification accuracy for BCI.^40–42 Among the task durations that achieve classification accuracy above 70%, $T_t$ values of 8 s each can be chosen as the task duration.

Fig. 9

Average quaternary classification accuracies as a function of the task duration. Square, circle, and triangle symbols indicate the classification accuracies of participants 1, 2, and 3, respectively. Error bar indicates the standard deviation. At a window length of 6 s, the classification accuracies of all the participants were below 70%.

3.2.

Multiclass Classification

The classification performance of the binary, ternary, and quaternary classifications was investigated at the $T_t=8\mathrm{s}$. For binary classification, the six combinations (${{}_4\mathrm{C}}_2=6$) of two data classes among four data classes (LA: left arm lifting, RA: right arm lifting, LK: left knee extension, and RK: right knee extension) were compared [i.e., LA versus RA), (LA versus LK), (LA versus RK), (RA versus LK), (RA versus RK), and (LK versus RK)]. Four combinations (${{}_4\mathrm{C}}_3=4$) of three data classes [i.e., LA versus RA versus LK), (LA versus RA versus RK), (LA versus LK versus RK), and (RA versus LK versus RK)], and one combination (${{}_4\mathrm{C}}_4=6$) of four data classes [i.e., LA versus RA versus LK versus RK)] were considered for ternary and quaternary classifications, respectively.

Figure 10 shows the average classification accuracy of four different classes of the data for eight participants according to the number of classes of the data. The classification accuracy decreases with increases in the number of data classes. The binary classification shows the highest classification accuracy for each participant, ranging from 87.1% to 95.5%, and the quaternary classification shows the lowest classification accuracy for each participant, from 73.8% to 91.5%. Table 2 shows the average classification accuracies according to the number of data classes. The features set of each participant which showed the highest classification accuracy was selected. Because 10 runs of cross-validation were conducted, the results are presented as the $\text{means}\pm \text{standard deviation}$. The highest classification accuracy was obtained using PCA for participants 1, 3, 6, 7, and 8, because the LDA could not be applied to the high-dimensional feature vector because of the singular scatter matrix problem.⁴³ However, the classification accuracy obtained by LDA and applied to a low-dimensional feature vector was comparable to that obtained by PCA.

Fig. 10

Average classification accuracies of eight participants according to the number of data classes ($N=2$, 3, and 4) for $T_t=R_t=8\mathrm{s}$. Error bar indicates the standard deviation. The highest and lowest classification accuracies are observed when the number of classes of the data is 2 ($N=2$) and 4 ($N=4$), respectively.

Table 2

Average classification accuracies according to the number of data classes for binary, ternary, and quaternary classifications for Tt=Rt=8  s. Principal component analysis (PCA) and linear discriminant analysis (LDA) are applied, and the reduced dimensions by PCA or LDA are represented in the brackets. Among mean (1), max (2), slope (3), delay (4), var (5), and median (6), the feature sets are determined for obtaining the highest classification accuracy. Because 10 runs of cross-validation were conducted, the classification accuracies are presented in the form of mean±standard deviation (SD).

3.3.

Data Transfer Rate

The ITR for all the participants with the number of data classes ranging from two to four are presented as a function of classification accuracy in Fig. 11. The maximum and minimum classification accuracies were achieved when the number of data classes was two ($N=2$) and four ($N=4$), respectively. For all the participants, the ITR increases as the number of classes increases. The maximum ITR values were 1.44, 1.08, 1.22, 1.08, 0.75, 0.80, 1.08, and $1.06\text{bits}/\text{trial}$ for participants 1 through 8, respectively. At $T_t=8\mathrm{s}$, the available BPMs were 2.81 to $5.40\text{bits}/\min$ [0.75 to $1.44\text{bits}/\text{trial}\times 3.75\text{trials}/\min$ ($0.27\min /\text{trial}$)].

Fig. 11

Achievable information transfer rates (ITRs) as a function of the different numbers of data classes for all the participants ($N=2$, 3, and 4). Dash-dotted, dotted, and solid lines are the theoretical values of the ITR using BCI when the number of data classes is 2, 3, and 4, respectively.

4.1.

Comparison to Previous Studies

In many studies, classification has been performed using finger tapping as the motor imagery task. For an EEG-based BCI, quaternary classification using motor imagery tasks has been investigated, while binary or ternary classifications have been considered for fNIRS-based BCI.^29,30 Although motor imagery is more suitable for BCI, in this study, overt motor execution was performed for motor tasks to obtain distinct hemodynamic responses, because overt execution leads to more distinct activated signals. In spite of the quaternary classification, the classification accuracies achieved in this study were comparable to those of the binary classification in other studies.^17,24,26

4.2.

Experimental Consistency

Since the classification results can potentially be affected by the global changes in the signals by the contact condition of the optodes, ambient light, participant fatigue, attention, and physiological change, a custom-built head cap was used to maintain the optode contact conditions and to totally prevent the effects of ambient light on the experimental results. Physiological changes of the participants were filtered out at the preprocessing step. Because the experiments lasted for a long time, an exhausted participant could not participate in some experiments, as shown in Fig. 9.

4.3.

Data Transfer Rate

An experiment was designed to investigate four-choice fNIRS-based BCI, since the achievable BPM using a binary system is just 1 to $2\text{bits}/\min$. The suggested quaternary system could achieve a BPM of up to $5.33\text{bits}/\min$ with lower accuracy as a trade-off. Although Blankertz et al. achieved $35\text{bits}/\min$ using a 128-channel EEG system,⁴ it is impossible for an fNIRS-based BCI system to achieve a comparable ITR because of the inherent hemodynamic delay. A task duration and rest duration time longer than 16 s ($T_t+R_t$ is enough to obtain sufficient accuracy based on the results. Therefore, $7.5\text{bits}/\min$ [$2\text{bits}/\mathrm{trial}$ (theoretically maximum B) $\times 3.75\text{trials}/\min$] may be the maximum BPM using a four-choice fNIRS-based BCI system. This result shows the limitation of the fNIRS-based BCI using the quaternary system.