1.

Introduction

Stress is known as one of the main factors threatening human health. It has inevitable consequences on performance in daily life and in the workplace, especially where risky situations are present.^1,2 People with low psychosocial resources are more susceptible to illness and mood disturbance in stressful situations, even if little stress exists in their lives.³ Long-term exposure to stress causes a variety of health problems, such as heart disease, obesity, diabetes, stroke, and depression.⁴ Therefore, detecting and managing stress in its early stages is vital to limiting its damaging consequences on health.

Stress detection is useful not only for research and empirical studies but also for nonclinical applications. For example, such systems can be employed by people who are in a stressful working environment in order to evaluate occupational stress. In this regard, psychological games have been designed in the field of biofeedback, which put the user in a tense and stressful virtual environment while the physiological signals are recorded simultaneously, in order to recognize online the mental stress and feedback the related information to the user.⁵

Appropriate standardized protocols for studies of stress are required to induce stress in a reliable and credible way. They have been categorized into three groups: physical, psychological, and mixed stressors, which have different physiological effects on the human body and are used in the various applications and studies. Physical stressors alter systemic circulation and make direct biochemical changes, while psychological stressors induce behavioral, physiological, and biochemical alterations in the brain. Mixed stressors include the other two stressors’ alteration and therefore make changes in both the systemic circulation and brain regions.⁶ Montreal Imaging Stress Task (MIST) is one of the popular psychological stressor techniques employed in stress assessment studies.⁷ In the MIST, mental arithmetic calculations are performed during a stressful condition (in a limited time) while participants’ responses are assessed.⁸

The autonomic nervous system (ANS) is responsible for regulating involuntary body functions such as heartbeat, blood flow, and breathing. The ANS is divided into two branches: the sympathetic nervous system (SNS) and the parasympathetic nervous system.⁹ Stress activates the SNS and hence raises the cortisol hormone in the adrenal cortex. As a result of cortisol release in the body, function, and structure of the brain are affected. Therefore, stress can be assessed by brain signals including electroencephalography (EEG) and functional near-infrared spectroscopy (fNIRS).^10,11 In the past decade, fNIRS has been used as a measure to quantify the mental stress and workload,^12,13 showing increased concentration of oxyhemoglobin (${\mathrm{HbO}}_2$) and decreased concentration of deoxyhemoglobin (HHb) in the prefrontal cortex (PFC) during a mental stress. In addition, fNIRS has been used in comparison and in combination with EEG signal for stress assessment.^4,14

Under stress, the ANS affects the cardiac function and will act to raise the heart rate (HR) and respiration activity.¹⁵ Therefore, scanning cardiac activity is an informative measurement in order to evaluate the state of ANS. Heart rate variability (HRV) is a measure of the variations in the time intervals between the heartbeats, which is controlled by the ANS. The gold standard of HR tracking is to analyze the interbeat intervals (IBIs) determined using an electrocardiogram device.¹⁶ Alternatively, HR signal can be measured based on blood flow pulsation derived from photoplethysmography (PPG) or fNIRS signals.^17,18 Many studies have been performed in order to assess the stress by measuring the HR or HRV.^1,19–23 In recent years, HR or HRV measurement has been employed for stress assessment in combination with the measurements from other physiological signals including galvanic skin response, skin temperature, pupil diameter, electrooculogram, and EEG.^24–29

NIR devices measure the hemodynamic changes, which are caused not only by cognitive activities of the brain but also by the so-called physiological artifacts such as HR, breathing, and Mayer signal, and the nonphysiological noises such as motion artifact and electrical noises. While the physiological artifacts and other noises can be eliminated using different algorithms including wavelet analysis, auto-regression models, conventional or adaptive filtering, Kalman filter, adaptive Viener filter, and independent component analysis (ICA),^30–32 the physiological artifact produced by the heart beat can be extracted from the fNIRS signals and effectively used for the HR analysis as another source of information. Therefore, by means of only one portable device, it is possible to extract two different signals, which contain important functional information from both brain and heart. Since 2011, several studies have been performed in order to compute the HR signal from the fNIRS signals. Trajkovic et al. used empirical mode decomposition and parameter estimation of a model for almost periodic signals algorithms to estimate normal-to-normal (NN) intervals from the fNIRS signal.¹⁸ Although the accuracy of their proposed method was relatively acceptable and appropriate, due to the nature of their algorithm, it cannot be applied in the real-time applications. Perdue et al.³³ used the Nakajima et al.³⁴ band-pass filtering method, which had been applied to PPG signal, to extract HR from the fNIRS signal in infants.

In previous studies, although they could derive HR from the fNIRS signal with an acceptable accuracy, due to the fact that each algorithm has been applied in diverse conditions and to different signals, the presented results are not comparable. Therefore, it is necessary to apply all methods to one sample dataset for a better comparative study.

Efficient derivation of the HR from the fNIRS signal gives an opportunity to use the HR measurement as a heart response in combination with the fNIRS signal, which modulates the brain activity. Since stress can functionally affect the brain and heart, it is possible to assess the amount of stress by analyzing the brain and heart measurements, which can be derived by using only the fNIRS signals.

In this study, combining and improving peak detection algorithms including AMPD,³⁵ S function,³⁶ and M2D,³⁷ a new real-time algorithm was introduced to extract the heartbeat signal from the fNIRS signal. First, this algorithm was applied to the simulated fNIRS and ECG signals to identify each heartbeat, which was then used for IBIs calculation. Next, this algorithm was applied to the in vivo signals recorded in the laboratory environment. The HR derived from the fNIRS signal (EHR) was then compared to the reference HR derived from the ECG signal (RHR) in both simulation and real data.

Several feature conditioning and classification algorithms, including principal component analysis (PCA),³⁸ ICA,³⁹ support vector machine (SVM),⁴⁰ and multilayer perceptron (MLP),⁴¹ were developed and used for stress level classification. In this study, a modified version of the MIST⁸ was designed and used to induce a higher level of mental stress.

As it was expected, the proposed combination of the peak detection algorithms computes the EHR derived from both simulated and in vivo measured fNIRS signals more accurately rather than employing the algorithms individually. In addition, on the classification part, the EHR could outperform the RHR and fNIRS measurements in the stress assessment.

The rest of the paper is organized as follows. In Sec. 2, the materials and methods are introduced. In Sec. 3, after illustrating the results of HR derivation in both simulated and in vivo measured fNIRS signals, the result of stress classification is presented. The discussion and conclusion are drawn in Secs. 4 and 5, respectively.

2.

Materials and Methods

2.2.

Data Acquisition Protocol

23 fNIRS channels were used to monitor the activity of the brain in the PFC region. Figure 1 shows the locations of the transmitters, receivers, and hence the arrangement of channels. The OxyMon NIRS system (Artinis Medical Systems, The Netherlands), which was equipped with light sources at wavelength 845 and 762 nm, and a sampling rate of 10 Hz, was employed. The distance between each of the transmitters and receivers was set to 3 cm.

Fig. 1

The location of the fNIRS detectors and sources, and the layout of the 23 fNIRS channels. The source–detector distance was considered 3 cm.

ECG signals were recorded using three Ag/AgCl electrodes placed on the sites RA, LA, and LL (Fig. 2) with reference to the electrode located on the right earlobe in a monopolar configuration and with a sampling rate of 512 Hz. Meanwhile, Nasion was selected as the ground for the ECG and NIR devices. In this study, the ECG signal recorded from the LL electrode is processed while the other two signals were recorded to be used in case the LL signal is noisy and inappropriate.

Fig. 2

Three-lead ECG placement, in which three signals were recorded with reference located on the right earlobe.

2.3.

Experiment Procedure

The HR of the participants was measured before entering the laboratory environment. In order to reduce the adverse impact of anxiety caused by the laboratory and medical devices on the experiment, the signal recording process was started in that environment when the HR of the participants reaches its normal range determined before entering the laboratory.

Signals were registered at the Iranian National Brain Mapping Lab. The laboratory environment was silent, ventilated, and free from any environmental stress. Participants were asked to sit on a comfortable chair, which was located $\sim 1.5\mathrm{m}$ far from the monitor and avoid shaking during the experiment as much as possible (Fig. 3). The experiment has an ethics code of IR.IUMS.REC.1396.810194120 from Iranian University of Medical Sciences.

Fig. 3

The experiment procedure for participant performing the task.

2.4.

Task Sequence

In this study, MIST^4,8 is designed to induce mental stress. Several modifications are performed in order to induce a higher level of stress during the task. This task includes four levels, as follows:

• 1. In the first level, a brief explanation is given to the participant to get familiar with the task sequence.

• 2. Then, in the second level (training level), the participants perform a task including simple mathematical calculations (summation, subtraction, multiplication, or division of two or three operands) for 1 min in order to accustom to the task. The answers are integer numbers between zero and nine. Operands and operators related to each question are randomly selected with a uniform distribution from sets [0-9999] and {+, −, *, /}, respectively. At this level, the participants can estimate the response time to the questions. Since this level is only for training, no signal is taken from the individuals.

• 3. The third level (control level) is similar to the previous one. Wrong answers are not shown to the participants and there is also no time limitation. However, there are two main differences. First, the fNIRS and ECG signals are recorded simultaneously and second, the average of the answering time to each question is also calculated.

• 4. In the fourth level (stress level), a time limitation is added to the third level. The time considered to solve each question is 90% of the mean time calculated in the control level. Moreover, while incorrect answer or the completion of the time to answer each question is marked with “incorrect” and “time’s up” messages, they are also emphasized by changing the background color to red. At this level, the arrangement of the buttons 3 to 9 is changed in order to induce a higher level of stress. Changing the location of the buttons will add more shock and stress to the participants because during the previous two levels the user has been accustomed to the sequential order of the buttons.

Time scheduling in the data collection protocol for each participant is as follows. The total time of the protocol is about 11 min. The first 1 min is aimed to train the participants in order to get familiar with the task. Then there are two 5 min period namely “control” and “stress” levels, each performed under different conditions. Each level lasts 5 min, including 20 s “rest” and 30 s “calculation” which are irritated five times. During the cycles of calculation in the control level, to reduce the stress, there is no stop at the end of the 30 s period till he/she responds the final question. In addition, in both control and stress levels going from the rest part to the calculation is controlled by the participant (clicking on the keyboard or mouse) so not to add excess stress. The block diagram of the training, control, and stress levels is shown in Fig. 4.

Fig. 4

The block diagram of the MIST, including training, control, and stress levels. Control and stress levels have five two-part sections, which consist of 20s rest and 30s mathematical calculation under the presented circumstances.

Before starting the task, the participants were encouraged to achieve more points in the final level (stress level) and they were asked to do their best to solve the questions. In addition, the percentages of the correct answers in the control and stress levels were calculated to obtain the participants’ engagement at each level. According to Ref. 4, the less percentage of correct answers, the less stress is transferred to the subject. For this reason, signals related to those participants whose correct answers in the control and stress levels are less than 90% and 30%, respectively, are removed from the dataset.

In order to synchronize the task with the NIR and ECG devices, parallel communication interface is used. The starting and ending of the calculation parts in the control and stress levels are marketed by sending trigger pulses.

2.5.

Signal Processing

The method considered to achieve the objectives of this study consists of two main parts: HR derivation and stress classification. These two parts are presented in the block diagram shown in Fig. 5.

Fig. 5

Schematic showing the signal processing methods for stress assessment consists of two main parts: HR derivation and stress classification.

2.5.1.

HR derivation

The process of HR derivation includes three subsections as follows: preprocessing, derivation process, and correction process. At the end of this section, fNIRS signal simulation method is explained.

Preprocessing

Two preprocessing steps are employed in order to eliminate the artifacts and maintain the heart beat component in the fNIRS signal. First, the algorithm proposed in Ref. 37 is employed to eliminate the sharp changes of the baseline due to such factors as motion artifact or undesired light intensity variations of the sources. Then, a band-pass filter consisting of two zero-phase Butterworth filters, a 1.9-Hz low-pass filter and a 1-Hz high-pass filter is used to remove the unwanted components including Mayer waves with a frequency component around 0.1 Hz, respiratory fluctuations, and other frequency components unrelated to the heart beat fluctuations of the fNIRS signal.

In order to compute HR from the ECG signal (RHR), it is necessary to first eliminate the artifacts, noises, and unrelated frequency components. It was observed that a zero-phase low-pass Butterworth filter with a cutoff frequency of 0.44 Hz can effectively identify the $R$-peak in the ECG signal recorded from the LL electrode reference to the electrode placed on the earlobe.

Derivation process

Detection of the peak points of the signal is the most important part of the HR derivation process. After the extraction of the heart beat signal, HR is measured by calculating the time interval between the two adjacent peaks. There are real-time algorithms to detect the peaks of a signal. In this study, four algorithms including AMPD, S1 and S5 functions, and M2D are employed to detect the peak points in the fNIRS signal.

Two reforms are made to the S1 and S5 functions, and M2D algorithms in this study in order to improve their performance as follows:

• 1. In addition to the calculations performed in these algorithms, first-difference of the signal is also obtained. In this correction process, a point of a signal will be determined as a peak point, when it satisfies two conditions as follows: First, the value of its objective function is larger than the predetermined threshold, and second, the value of the first-difference in that point is also 0.1 times greater than the maximum first-difference.

• 2. Some of the points detected by the algorithms might be relative maximum and nonpeak. A reform is done in order to exclude them from the detected peak points as follows: after detecting point $n1$ as a peak, 0.25-s searching before the time $n1$ is done (which is equal to $F_s/4$ samples before the time $n1$, where $F_s$ is the sampling frequency). If there is another peak, which has been previously detected in the 0.25-s period, the point $n1$ is considered as a relative maximum and is removed from the detected peak points (otherwise $n1$ will be remarked as a peak).

In order to combine the performance of the algorithms, the average of the IBIs measured by the mentioned four techniques is computed based on a weighted averaging procedure named weighted mean algorithm. The weights are selected based on the simulation results. The simulation process is explained in Sec. 2.5.1.4. According to Eq. (1), the weight of each algorithm is calculated by inverting the normalized BAR (Bland Altman ratio)⁴² computed when it is applied to the simulated signal

Eq. (1)

$$W_j=\frac{1}{({\mathrm{BAR}}_j/{\sum }_{i=1}^4{\mathrm{BAR}}_i)};j=1,\dots ,4.$$

The recommended weights for the AMPD, S5 function, S1 function, and M2D are 5.04, 4.39, 3.78, and 3.24, respectively, which have been determined based on the performed simulation in this study.

Correction process

Time intervals calculated between the two adjacent peaks might be apart from their correct values due to false positive and false negative errors of the peak detection. In order to reduce these errors, a correction process is needed to specify these cases and fix the time intervals. In this study, a correction process is introduced to eliminate the outliers of the IBIs based on a window method. In addition to total standard deviation (SD) of the IBIs calculated from the start of the recording to the current time, two other metrics including the mean (m) and standard deviation (sd) of the IBIs are calculated in the recent 4 s interval. These three parameters are used to correct the IBI which is outside the permissible range. The permissible range in this correction process is based on the study.³⁴ The sd and SD parameters represent the local and total signal variations, respectively. Since the motion artifact causes sudden local oscillations in the signal, it is locally considered in the configuration of the permissible range by using the sd calculated during the 4-s time intervals. If the fluctuations of the IBIs during the 4-s time interval are high, it is better to reduce the permissible range in order to make further corrections. Therefore, sd must be considered in the denominator when calculating the permissible range. Moreover, in order to consider the natural oscillations of the HR signal in the correction process, SD is included as a multiplicand in the permissible range. The HR correction process is described step by step as follows.

• 1. The current IBI is examined.

• 2. The mean (m) and standard deviation (sd) of IBIs in the recent 4 s interval are calculated.

• 3. The standard deviation of IBIs from the start to the present time is computed (SD).

• 4. Parameter $k$ is defined as follows: $k=0.8*15/\mathrm{s}\mathrm{d}$.

• 5. Finally, if the IBI is outside the range of $m-k*\mathrm{SD}<\mathrm{IBI}<m+k*\mathrm{SD}$, it will be considered as outlier and the corrected value will be obtained by a cubic interpolation of the closest neighboring IBIs which are within the permissible range.

In the descriptions aforementioned, the constants 0.8 and 15 used in calculating the parameter $k$ have been experimentally obtained. If there is no permissible IBI in the 4-s neighborhood of the IBI which is outside the permissible range, the fluctuations of HR signal in that range are normal, not caused by the instantaneous noises, and hence the HR will not be corrected.

In this study, the oxyhemoglobin signals are employed to compute EHR using the proposed derivation and correction process. Although each participant has only one reference HR signal derived from his/her ECG signal using the AMPD algorithm, 23 different EHR signals are derived from 23 fNIRS channels (each single HR signal is obtained from one channel).

Simulation

In addition to the in vivo measurement in this study, simulation is added to compare and analyze the applied HR derivation algorithms. The fNIRS signal simulation method performed in this paper is based on the simulation procedure proposed in Ref. 43. Several changes are made to generate a simulated fNIRS signal, which is more realistic and closer to the real ones. In addition, the measurements proposed in article⁴⁴ are used to initialize the required parameters for simulation. In this simulation, instead of instantaneous changes in the frequency and amplitude of the sinusoidal signals, each $N$ periods of the sinusoidal component, frequency and amplitude are changed ($N$ is chosen proportional to the signal length). As a consequence, the high-frequency artifact, which is due to the instantaneous variations of the frequency and amplitude will be eliminated. Moreover, when changing the frequency of the signal, the sinusoidal phase of the signal is reorganized to remove the abrupt changes in the simulated signal.

2.5.2.

Stress classification

In this regard, the stress classification is separately performed using fNIRS, EHR, and RHR signals recorded during the MIST. Initially, the features which are described in the following are extracted. Then, they are pruned by applying the PCA and ICA feature conditioning algorithms, and features with less correlation and dependence are obtained. Finally, the dataset is categorized using SVM and MLP classifiers.

Feature selection is considered as one of the most important steps in the classification process. So that the proper extraction which means considering comprehensive features of the signal, even if simple classifiers are used, can greatly increase the classification accuracy. For this reason, in addition to defining some features in the time domain, it is better to extract some other features in the frequency domain of the signal. One of the advantages of feature selection in the frequency domain is that the variations and fluctuations of the signal can be revealed more completely with fewer features than time domain.

In this study, nine different features, including four features in the time domain ($m_t$, $s{d}_t$, $S_t$, $K_t$) and five features in the frequency domain ($m_f$, $s{d}_f$, $S_f$, $K_f$, LF/HF) are extracted from the HR measurements. In this regard, $m_t$, and $s{d}_t$, which are representing AVNN (average of NN intervals) and SDNN (standard deviation of NN intervals), respectively, are time-domain HR metrics.³⁷ In order to more efficiently demonstrate the distribution of the HR data in the feature extraction process, in addition to the first two statistical moments of the HR data, the third and fourth moments are also calculated ($S_t$, $K_t$). In the frequency domain, recent studies have reported that LF/HF, which means the ratio of the low frequency components (LF) to the high frequency components (HF) of the HR measurements, does not accurately measure the autonomic balance.^45,46 However, it is employed in this study as a feature since it has been able to provide a better characterization of HR derived from fNIRS signal compared to the individually used LF and HF metrics under stressful condition.⁴⁷ The other features in the frequency domain of the HR measurements are extracted as follows: the Fourier transform of the data is considered as the probability density function (PDF) of random variable $f$ (frequency), and therefore from the first to fourth statistical moments of random variable $f$ are obtained.

Regarding the fNIRS signal, the presented four features in the time domain, which are the first four statistical moments of the signal, are also extracted. These features characterize the time-history data and have been employed in brain–computer interface (BCI) studies to represent the fNIRS signal.⁴⁸ In addition, the other four statistical moments in the frequency domain ($m_f$, $s{d}_f$, $S_f$, and $K_f$) are extracted in order to efficiently reveal the fluctuations of the fNIRS signal in the extracted features. The selected features are presented in Table 1. The descriptions of the statistical moments are according to Ref. 49.

Table 1

Features that are extracted from the time and frequency domains of HR and fNIRS signals.

Feature	Description
$m_t$	$\frac{1}{N}\sum _{i=1}^N{x}_i$
$s{d}_t$	$\sqrt{\frac{1}{N}\sum _{i=1}^N{(x_i-m_t)}^2}$
$S_t$	$\frac{1}{N*s{d}_t^3}\sum _{i=1}^N{(x_i-m_t)}^3$
$K_t$	$\frac{1}{N*s{d}_t^4}\sum _{i=1}^N{(x_i-m_t)}^4$
$m_f$	$\sum _{i=1}^N{X}^{\prime }(f_i)*f_i$
$s{d}_f$	$\sqrt{\sum _{i=1}^N{(f_i-m_f)}^2*X^{\prime }(f_i)}$
$S_f$	$\frac{1}{s{d}_f^3}\sum _{i=1}^N{(f_i-m_f)}^3*X^{\prime }(f_i)$
$K_f$	$\frac{1}{s{d}_f^4}\sum _{i=1}^N{(f_i-m_f)}^4*X^{\prime }(f_i)$
LF/HF	Ratio of low-frequency component (0.04 to 0.15 Hz) to high-frequency component (0.15 to 0.4 Hz) of HR measurement

Where $x_i;i=1,\dots ,N$ is the desired signal (fNIRS or HR) in the time domain which has been windowed with a particular length ($N$). The signal $X^{\prime }(f_i)$; $i=1,\dots ,N$ is the normalized Fourier transform amplitude of the windowed signal, which $f_i$ is the digitized frequency in the range of [0,Fs/2] (Fs is the sampling frequency). Signal $X^{\prime }(f)$ is obtained according to Eq. (2). Where $|X(f_i)|$ is the amplitude of the signal Fourier transform in the frequency of $f_i$,

Eq. (2)

$$X^{\prime }(f_i)=\frac{|X(f_i)|}{{\sum }_j|X(f_j)|}.$$

As mentioned in Sec. 2.4, control and stress levels have five two-part periods: 20-s rest, and 30-s calculation under the presented circumstances. In this study, same as classification procedure conducted by Al-Shargie et al.,⁴ a classification of the 30-s parts between control level and their corresponding at the stress level is desired. To form the dataset, 10 30-s parts are separated from each signal. Then each part is divided into two 15-s halves, and the mentioned features are extracted from each of these 15-s halves. This process is applied to both oxyhemoglobin and deoxyhemoglobin signals of fNIRS channels, and hence 10 data (five control data and five stress data) with 32 features are extracted from each fNIRS channel. Regarding the HR measurements, 10 data with 18 features are constructed from each measurement.

3.

Results

3.1.

Results of HR Derivation Algorithms

3.1.1.

Results of applying HR derivation algorithms to simulated fNIRS signal

The comparison of an example-based simulated and in vivo measured fNIRS signal is shown in Fig. 6. It is observed that the baseline and high-frequency oscillations of both signals are similar in terms of structure.

Fig. 6

An example of (a) in vivo measurement and (b) simulated fNIRS signal.

The EHR signal computed by applying the weighted mean algorithm is shown in Fig. 7 compared with the reference HR (RHR). A quite good agreement can be seen between the simulated RHR and EHR signals.

Fig. 7

(Blue) EHR and (red) RHR signals derived from the simulated fNIRS signal.

To quantitatively compare the algorithms used for derivation of the HR from the fNIRS signal, following metrics were computed: mean error (mean difference between EHR and RHR), SD error (standard deviation of error), RMS error (root mean square of error), BAR (the ratio of half the range of limits of agreement to the mean of the pairwise measurements), maximum error, mean, and standard deviation of the EHR and RHR, and Spearman linear correlation between EHR and RHR. These metrics have been introduced as the measures to verify the agreement between HR measurements derived from PPG and ECG.⁴²

The weights of the proposed weighted mean algorithm were obtained using the BAR calculated for each of the S1 function, S5 function, AMPD, and M2D algorithms. Comparing the values reported in Table 2, the best algorithm in each metric is identified. The weighted mean algorithm in five metrics and AMPD algorithm in two metrics are better than the other algorithms.

Table 2

Quantitative comparison between the algorithms applied to the simulated fNIRS signal.

Quantity	S1 function	S5 function	AMPD	M2D	Weighted mean
Mean error (s)	0.026	0.020	0.014	0.052	0.022
SD error (s)	0.040	0.035	0.030	0.047	0.023
RMS error (s)	0.048	0.040	0.033	0.070	0.032
BAR (%)	7.206	6.234	5.393	8.445	4.210
Maximum error	0.150	0.153	0.158	0.222	0.136
Mean of EHR and RHR (s)	1.088, 1.089	1.088, 1.089	1.090, 1.089	1.090, 1.089	1.089, 1.089
SD of EHR and RHR (s)	0.099, 0.106	0.096, 0.106	0.098, 0.106	0.105, 0.106	0.093, 1.062
Linear correlation (%)	88.823	92.117	95.801	76.594	95.684

Note: The best results for each metric is written in bold.

3.1.2.

Results of applying HR derivation algorithms to in vivo measured fNIRS signal

In this section, the results of applying the HR derivation algorithms to the signals recorded during the MIST are evaluated. Figure 8 shows the average of the EHR signal computed from one of the participants’ signals using the weighted mean algorithm, comparing with the RHR signal derived from the ECG signal using the AMPD algorithm. The linear correlation between the two signals is 98.48%.

Fig. 8

An example of (blue) the averaged EHR and (red) RHR signals computed from the fourth participant.

The statistical metrics proposed in the previous section were obtained for EHR derived from the in vivo measurements in order to quantitatively compare the algorithms. In Table 3, the mean and standard deviation of the metrics calculated for each algorithm are summarized. It is observed that the weighted mean algorithm in seven metrics and the M2D algorithm in one metric are better than other algorithms.

Table 3

Quantitative comparison between the algorithms applied to the in vivo fNIRS signal.

Quantity	S1 function	S5 function	AMPD	M2D	Weighted mean
Mean error (s)	0.028 ± 0.004	0.028 ± 0.003	0.029 ± 0.004	0.030 ± 0.004	0.018 ± 0.004
SD error (s)	0.026 ± 0.006	0.026 ± 0.005	0.026 ± 0.006	0.029 ± 0.008	0.018 ± 0.005
RMS error (s)	0.038 ± 0.007	0.039 ± 0.005	0.039 ± 0.007	0.042 ± 0.008	0.026 ± 0.006
BAR (%)	6.925 ± 1.293	7.056 ± 1.067	6.961 ± 1.214	7.989 ± 2.451	4.934 ± 1.841
Maximum error (s)	0.292 ± 0.119	0.299 ± 0.128	0.260 ± 0.105	0.347 ± 0.173	0.190 ± 0.078
Mean of EHR and RHR (s)	0.739, 0.738	0.739, 0.738	0.739, 0.738	0.739, 0.738	0.737, 0.738
SD of EHR and RHR (s)	0.074, 0.067	0.074, 0.067	0.073, 0.067	0.073, 0.067	0.065, 0.067
Linear correlation (%)	86.225 ± 0.054	85.922 ± 0.054	85.051 ± 0.047	83.892 ± 0.047	92.254 ± 0.034

Note: The best results for each metric is written in bold.

The results show that the weighted mean algorithm performs a better derivation of EHR in comparison to the individual algorithms, producing a higher percentage of Spearman linear correlation of $92.254\pm 0.034\%$ between the EHR and their corresponding RHR. Therefore, it is more suitable for real-time applications where the motion artifact is not considerable. The weighted mean algorithm is employed to compute HR for the aim of stress classification.

3.2.

Results of Stress Classification

An example of filtered fNIRS signals averaged between channels is shown in Fig. 9. The applied filter is the band-pass filter introduced in Sec. 2.5.1.1. The first five yellow parts and the second five purple parts are related to the control and stress levels, respectively. It is observed that there is a greater signal variation in the stress level.

Fig. 9

An example of averaged fNIRS signals recorded from the fifth participant during the MIST. The first five yellow parts are corresponding to the control level and the second five purple parts are related to the stress level. The black, red, and blue curves are related to the total hemoglobin, oxyhemoglobin, and deoxyhemoglobin signals, respectively.

The average of the EHR signal computed from to the same participant and his RHR signal are shown in Fig. 10. The Spearman linear correlation between the two signals is 97.98%. As the previous figure, higher changes related to the stress level in the HR trend are observed in this figure.

Fig. 10

An example of (blue) averaged EHR signal derived from the fifth participant in comparison with his (red) corresponding RHR signal during the MIST. The first five yellow parts are corresponding to the control level and the second five purple parts are related to the stress level.

SVM (with RBF kernel function) and MLP (with conjugate gradient optimization) classifiers and PCA and ICA feature conditioning algorithms were used to classify the features extracted from each of the fNIRS, EHR, and RHR measurements. For this purpose, 75% of the data was considered as train data, 15% as the validity data, and the remaining 15% as the test data. Since the train, validity, and test data were randomly selected with a uniform distribution, fNIRS and HR were categorized separately 20 times in order to reduce the effect of random selection and the variability of classification accuracy. The mean and standard deviation of the classification accuracy obtained by each classifier are reported in Table 4.

Table 4

Results of stress classification using different classifiers constructed by SVM/MLP classifier, and PCA/ICA feature conditioning under employing fNIRS, EHR, and RHR features.

Classifier	SVM			MLP
Data	fNIRS	EHR	RHR	fNIRS	EHR	RHR
Accuracy (%) without feature conditioning	78.8 ± 2.4	92.6 ± 1.1	62.2 ± 10.1	72.4 ± 4.1	93.9 ± 1.3	50.0 ± 13.4
Accuracy (%) using PCA	77.1 ± 1.9	94.6± 0.9	56.1 ± 9.7	68.5 ± 4.8	93.7 ± 1.6	52.9 ± 11.9
Accuracy (%) using ICA	77.2 ± 3.1	94.4 ± 0.8	54.3 ± 9.6	71.1 ± 3.3	94.3 ± 2.2	51.7 ± 12.8

Note: The best accuracy for each feature type and each classifier is written in bold.

The EHR features were classified using the SVM and MLP classifiers with the accuracies of $92.6\pm 1.1\%$ and $93.9\pm 1.3\%$, respectively, which are more than the accuracy obtained by each of the fNIRS and RHR features (Table 4). Applying the PCA algorithm to the features extracted from the EHR, the stress classification accuracy using the SVM and MLP classifiers increased by 2% and decreased by 0.2%, respectively. The reason for reducing the classification accuracy by MLP after applying the PCA is that the MLP classifier itself has a nonlinear PCA in the first layer, and therefore it does not need the PCA feature conditioning. Applying the ICA algorithm (the number of independent components was obtained using the PCA algorithm) to the features of EHR, stress classification accuracy using the SVM and MLP classifiers was increased by 1.8% and 0.4%, respectively. In total, the greatest stress classification accuracy during the MIST was obtained using the EHR features, which was classified by SVM classifier and PCA feature conditioning with the accuracy of $94.4\pm 0.8\%$. Considering both the fNIRS and EHR features, the average classification accuracy by SVM and MLP classifiers were calculated 91.3% and 88.1%, respectively, and then having applied the PCA algorithm to these combined features, the average accuracy reduced by 0.4% and increased by 1.4%.

4.

Discussion

This study investigated if the HR derived from the fNIRS signal has the same efficiency as the HR derived from the ECG signal in the mental stress assessment. The mental states of 10 healthy male subjects during the MIST were assessed by EHR, RHR, and fNIRS measurements. ECG and fNIRS signals were recorded simultaneously in order to measure the heart and brain responses. Real-time peak detection algorithms were employed including AMPD, M2D, S1 function, and S5 function in order to compute the HR from the fNIRS signal. The weighted mean algorithm was proposed to combine the performance of these four algorithms. The time and frequency domain features were extracted from the fNIRS, EHR, and RHR. They were then classified using different classifiers and the results of which are then compared to each other.

In the HR derivation from fNIRS signal, it was shown that the weighted mean algorithm is more accurate compared to the other four algorithms when applied to both simulation and in vivo data. In simulation, the percentage of BAR calculated for AMPD and weighted mean algorithms were 5.393% and 4.210%, respectively. In the case of real data, the percentages of BAR calculated for AMPD and weighted mean algorithms were $6.961\pm 1.214$ and $4.934\pm 1.841$, respectively.

Correlation coefficients such as Pearson, Spearman, and intraclass either measure the linear correlation or detect the agreement only according to the rank order between two variables. However, they are unable to detect additive or multiplicative bias. BAR is defined as an efficient way to quantify the agreement between two measures. Agreements were ranked as good when the percentage of BAR is lower than 10%.⁴² In the present study, in addition to the Spearman correlation coefficient, BAR has been used to verify the agreement between EHR and RHR. The result of evaluation metrics used in this study show that a high agreement ($r>0.9$, $\mathrm{BAR}<5\%$) was quantified between the EHR and RHR, which is in accordance with the previous studies.^18,34,47

Our results show that the EHR provides more effective means for stress classification during the MIST compared to the RHR and fNIRS in the sense of providing higher classification accuracy. Using the features extracted from EHR, stress was classified employing the combined SVM classifier and PCA feature conditioning with an accuracy of $94.4\pm 0.8\%$, which was 15.6% and 32.2% greater than the accuracy obtained using features extracted from fNIRS and RHR, respectively. Our stress classification results are in accordance with those reported in Ref. 47, in which EHR was compared to the HR derived from the PPG signal under a physical stressor. However, in the present study, a psychological stressor was induced, because it has become much more prominent and pervasive with advances in modern technologies and lives.⁵⁰ ECG signal was employed in this study to measure the reference HR, which is the gold standard for investigating the HR.¹⁶ In addition, a quantitative procedure was considered to assess the mental stress using different classifiers.

The stress classification accuracy acquired by fNIRS features in the present study is 5.35% lower than what was reported in Ref. 4. In addition to the difference in the experiment procedure, this lower accuracy might be resulted from our simple preprocessing method for fNIRS signal. Fusing the features extracted from the EHR and fNIRS using joint analysis techniques such as joint independent component analysis⁵¹ and parallel independent component analysis⁵² could be included in the future analysis.

As an interesting point, a noticeable difference was observed comparing the correlation values computed between the EHR and the RHR measurements for different fNIRS channels located at different positions over the PFC. The HR derived from the lateral positions of PFC showed more correlation with the corresponding reference HR. It can be concluded that blood fluctuations in those regions are more correlated to the heart beat. In contrast, the difference between the EHR and RHR measurements over fNIRS channels in some regions can be due to other physiological causes. We will further study this subject in order to find these differences.

The compatibility of the HR derived from fNIRS signal with the HR computed from ECG signal in the present study is in accordance with the results reported in Refs. 18 and 34. Derivation of the HR from fNIRS signal gives an opportunity to use the HR measurement as a heart response besides or in combination with the fNIRS signal which modulates the brain activity by means of only one portable device. In addition to the stress assessment, the EHR and fNIRS measurements can be employed for applications in which both the brain and heart are affected, including detection of the mental workload,⁵³ BCI,⁵⁴ biofeedback,^55,56 and diagnosis of some diseases,⁵⁷ especially in children and infants because of eliminating the need for additional HR sensors.^58,59

5.

Conclusion

fNIRS and ECG were recorded simultaneously in order to assess the mental stress induced by the MIST. The present study indicates that the EHR outperforms the RHR and fNIRS measurements in the stress assessment. In particular, our results showed that the stress was classified using the EHR features with the accuracy which was 15.6% and 32.2% greater than those obtained by the individual measurements of fNIRS and RHR, respectively. This suggests that the EHR can be employed for stress assessment as a compatible heart response besides or in combination with the fNIRS signal by means of only one portable device.