The Physyqx data, which was used in this study, is third-party data from the following published paper: Gupta R, Banville HJ, Falk TH. PhySyQX: A database for physiological evaluation of synthesised speech quality-of-experience. Proceedings of IEEE Workshop on Applications of Signal Processing to Audio and Acoustics (WASPAA'15). 2015;1--5. DOI: [10.1109/WASPAA.2015.7336888](https://doi.org/10.1109/WASPAA.2015.7336888). Researchers interested in accessing this data may contact Dr. Rishabh Gupta (<grishabhg@gmail.com>).

Introduction {#sec001}
============

Text-to-Speech (TTS) systems, which convert a written text into speech, and are becoming more widely implemented in mobile phones, car navigation systems, and other consumer electronics. Such systems play a critical role in many applications because speech is the most fundamental and easiest communication tool for human beings. Therefore, synthesized speeches must sound natural for good machine-to-human communications.

The research of TTS systems needs reasonable criteria that evaluate the qualities of synthesized speeches. Several current evaluation methods have their own advantages and disadvantages: (1) subjective ratings \[[@pone.0193521.ref001]--[@pone.0193521.ref003]\], (2) analyzing a speech signal itself \[[@pone.0193521.ref004]--[@pone.0193521.ref006]\], and (3) measuring the physiological responses of listeners to speech \[[@pone.0193521.ref007]--[@pone.0193521.ref014]\].

In the first approach, the two most common aspects for quality judgment are naturalness and intelligibility. Naturalness describes how close synthesized speech is to human speech, and intelligibility reflects how well the speech content can be heard. The former is usually measured by a mean opinion score (MOS) test \[[@pone.0193521.ref001]\], and the latter is gauged by semantically unpredictable sentences (SUS) \[[@pone.0193521.ref003]\]. In addition, valence and arousal are often used to evaluate the subjective impressions of speech \[[@pone.0193521.ref011], [@pone.0193521.ref013], [@pone.0193521.ref015], [@pone.0193521.ref016]\] and to model emotions \[[@pone.0193521.ref017]--[@pone.0193521.ref020]\]. Valence reflects a positive or a negative emotion. Arousal reflects the degree of intensity or activation. In a MOS test, subjects listen to speech and rate its relative perceived quality on some kind of a scale, for example, "excellent," "good," "fair," "poor," "bad." Then the scores are averaged across subjects. This is well established method for which references on how to perform it are available \[[@pone.0193521.ref002]\], making it the only standard way to evaluate the naturalness quality of synthesized speech. However, their appropriateness has not been fully proven because high inter- and intra-subject inconsistencies are often observed in the ratings, resulting in poor reproductivity \[[@pone.0193521.ref021]\].

In the second approach, speech quality is automatically evaluated at its signal level by software that inputs a speech file and outputs the estimated speech quality. Advantages of these methods include complete reproductivity and less time consumption after such software is developed. However, appropriateness is difficult to prove because the exact relationship between the acoustic features and the perceived quality of speech by a listener is not well understood \[[@pone.0193521.ref021]\]. In fact, speech quality must be evaluated not only physically but also psychologically because it is commonly defined as an assessment result within which a listener compares his/her perceptions with expectations \[[@pone.0193521.ref022], [@pone.0193521.ref023]\].

Last, quality estimation methods are emerging that measure the physiological responses of a listener \[[@pone.0193521.ref024]\]. Even though these methods have not been established yet, they are worth investigating because physiological signals can be recorded automatically and continuously to provide insight about listener's cognitive states without interruptions caused by directly asking him/her to answer questions. Among existing non-invasive physiological response measures, electroencephalography (EEG) has especially great potential to estimate a listener's perceived speech qualities for the following reasons. EEGs can be recorded at a higher temporal resolution, e.g., a millisecond range, than hemodynamic measures, including functional magnetic resonance imaging (fMRI) and functional near-infrared spectroscopy (fNIRS), both of which analyze the changes in blood flow that inherently take a few seconds until a brain response can be recorded. Temporal resolution is important to evaluate speech quality since the temporal structure of speech largely affects its perceived quality. In addition, an EEG recording equipment is relatively small, less expensive than other brain recording equipments, and can be even wireless, which allows us to use it in daily environments, whereas fMRI and magnetoencephalography (MEG) can only be used in experimental rooms because of the lack of portabity. Measuring physiological responses to speech in daily environments is critical because speech is everywhere. Despite the above advantages, the main disadvantage of physiological measures is the difficulty of data gathering. The amount of data that can be collected from a subject is limited for practical and ethical reasons. Conducting experiments is usually time-consuming and labor-intensive. In addition, physiological data are generally noisy and easily contaminated by artifacts. Furthermore, multi-channel EEG signals are usually highly correlated to each other, which makes the features extracted from them less informative compared to the height of their dimensions. These aspects of EEG (limited amount of data, noise, and high correlation and dimension) complicate training a predictive model with EEG data and require a sophisticated dimension reduction or regularization techniques \[[@pone.0193521.ref025]\].

Existing researches have analyzed EEG responses to speech stimuli using event-related potentials (ERP), which are time-locked responses to external or internal events in terms of a voltage change that are usually visualized and quantified after synchronous averaging of multiple epochs \[[@pone.0193521.ref007]--[@pone.0193521.ref009]\]. Due to its definition, measuring ERP need the instantaneous time-locking points at which an event occurs, complicating the use of ERP if stimuli onsets are gradual or unclear \[[@pone.0193521.ref026]\]. Therefore, ERP is not suitable for our purpose of the predicting perceived quality of speech whose length exceeds a second because it is usually unclear which time points affect a listener's perceived quality. Other research used power spectral density \[[@pone.0193521.ref014], [@pone.0193521.ref027]\] and their difference between EEG channels \[[@pone.0193521.ref011], [@pone.0193521.ref013]\] at multiple frequency bands. Neuroscience studies reported that EEG spectral changes in distinct regions and between hemispheres are related to emotions \[[@pone.0193521.ref028]--[@pone.0193521.ref031]\]. Other studies used EEG phase synchronization between EEG channel pairs and found a correlation to emotions \[[@pone.0193521.ref032], [@pone.0193521.ref033]\].

The purpose of this research is to predict the perceived qualities of synthesized speeches using only EEG. Interest is growing in the development of a machine learning algorithm that uses an input/output data structure as tensor formats \[[@pone.0193521.ref034]--[@pone.0193521.ref036]\]. Such tensor structured features were investigated in this study because EEG signals can have structures in time, frequency, space, experimental condition, and other modalities.

Materials {#sec002}
=========

We used the PhySyQX data set \[[@pone.0193521.ref010]\], which consists of speech files, their subjective rating scores from 21 subjects, and EEG signals from the same subjects recorded while they listened to the speech. The data recording protocol was approved by the INRS Research Ethics Office, and participants gave informed consent for their participation and to make their data anonymous and freely available online. The details of the data set and the experimental procedures are available in \[[@pone.0193521.ref010]\]. We obtained it by an e-mail request.

Speech stimuli {#sec003}
--------------

The speech stimuli presented to the subjects in the data set consist of speech collected from four humans and seven commercially available TTS systems. From each human and each TTS system, four English sentences were collected, whose durations ranged from 13 to 22 seconds. The 44 human and synthesized speeches were presented to each subject in random order.

Experimental procedure {#sec004}
----------------------

The experiment's timeline is shown in [Fig 1](#pone.0193521.g001){ref-type="fig"}. A 15-second rest period was provided before each stimulus presentation. It is followed by a subjective rating period during which the subjects evaluated the speech to which they had just listened. The subjective rating scales used in this study are shown in [Table 1](#pone.0193521.t001){ref-type="table"} and include overall impression (MOS), valence (VAL), and arousal (ARL). MOS was evaluated with a 5-scale rating and the others with a 9-scale using self-assessment manikin \[[@pone.0193521.ref037]\].

![Timeline of EEG and subjective evaluation data recording experiment.](pone.0193521.g001){#pone.0193521.g001}

10.1371/journal.pone.0193521.t001

###### Subjective rating scales.

![](pone.0193521.t001){#pone.0193521.t001g}

  Rating Scale         Abbreviation   Description
  -------------------- -------------- ------------------------------
  Overall Impression   MOS            1 (Bad) to 5 (Excellent)
  Valence              VAL            1 (Negative) to 9(Positive)
  Arousal              ARL            1 (Unexcited) to 9 (Excited)

### EEG recording and preprocess {#sec005}

EEG data were recorded throughout the experiment with 64 scalp channels. The sampling rate was 512 Hz, which was down-sampled to 256 Hz. All the channels were placed on scalp according to the modified 10/20 system \[[@pone.0193521.ref038]\]. Some channels were removed from further analysis because they were noisy. A band-pass filter was applied to all the data between 0.5--50 Hz and applied an independent component analysis based semi-artifact removal technique using the ADJUST toolbox \[[@pone.0193521.ref039]\]. After these preprocessing, the EEG signal of each subject was cut into 44 epochs corresponding to the stimuli listening periods.

Methods {#sec006}
=======

Feature extraction {#sec007}
------------------

All features were extracted at five frequency bands from a channel or a channel pair. The frequency bands include delta (*δ*: 1--4 Hz), theta (*θ*: 4--8 Hz), alpha (*α*: 8--12 Hz), beta (*β*: 12--30 Hz), and gamma (*γ*: 30--45 Hz). Let us denote the Fourier transformation at the frequency of *f*~*k*~ of the *n*-th trial recorded by the *m*-th channel by ***x***~*n*,*m*~(*f*~*k*~). An estimator of the power spectrum density and a phase spectrum denoted by *p*~*k*~ and *h*~*k*~ can be calculated using the periodogram method as follows: $$\begin{array}{r}
{p_{n,m}\left( f_{k} \right) = \frac{1}{T}\left| \mathbf{x}_{n,m}\left( f_{k} \right) \right|^{2}} \\
\end{array}$$ $$\begin{array}{r}
{h_{n,m}\left( f_{k} \right) = \text{angle}\left( \mathbf{x}_{n,m}\left( f_{k} \right) \right),} \\
\end{array}$$ where *T* is the number of time samples within a trial. Then, we averaged the power spectrum density over the frequency bins within the range of each frequency band to define channel-based features PSD~*n*~(*m*, *f*) as follows: $$\begin{array}{r}
{\text{PSD}_{n}\left( m,f \right) = \frac{1}{|D_{f}|}\sum\limits_{f_{k} \in D_{f}}p_{n,m}\left( f_{k} \right),} \\
\end{array}$$ where *D*~*f*~ is the index set of the frequency bins included in the range of the *f*-th frequency band and \|*D*~*f*~\| is the number of the elements in *D*~*f*~. The channel-pair-based features are also defined using the averaged power spectrum density and the phase spectrum as follows: $$\begin{array}{r}
{\text{PWD}_{n}\left( m_{1},m_{2},f \right) = \text{PSD}\left( m_{1},f \right) - \text{PSD}_{n}\left( m_{2},f \right)} \\
\end{array}$$ $$\begin{array}{r}
{\text{PHD}_{n}\left( m_{1},m_{2},f \right) = \frac{1}{|D_{f}|}\sum\limits_{f_{k} \in D_{f}}h_{n,m_{1}}\left( f_{k} \right) - h_{n,m_{2}}\left( f_{k} \right).} \\
\end{array}$$

If *M* EEG channels and *F* frequency bands are used (*F* = 5 in this study), *I* = *F*(*M*(*M* − 1) + *M*) features are calculated. The feature matrix X can be expressed as: $$\begin{array}{r}
{X = \begin{pmatrix}
{\mathbf{x}\left( 1 \right),\mathbf{x}\left( 2 \right),\ldots,\mathbf{x}\left( N \right)} \\
\end{pmatrix}^{\top} \in \mathbb{R}^{N \times I},} \\
\end{array}$$ where *N* is the number of training trials and ***x***(*n*) is a feature vector of the *n*-th trial and has all the features PSD~*n*~, PWD~*n*~, and PHD~*n*~.

To exploit structures of the features, we organized the features as a tensor $\mathcal{X} \in \mathbb{R}^{N \times M \times M \times F}$ as follows: $$\begin{array}{r}
{\mathcal{X}\left( n,m_{1},\mspace{360mu} m_{2},f \right) = \left\{ \begin{array}{ll}
{\text{PWD}_{n}\left( m_{1},m_{2},f \right)} & \left( m_{1} > m_{2} \right) \\
{\text{PHD}_{n}\left( m_{1},\mspace{360mu} m_{2},f \right)} & \left( m_{1} < m_{2} \right) \\
{\text{PSD}_{n}\left( m_{1},f \right)} & \left( m_{1} = m_{2} \right) \\
\end{array}\operatorname{}. \right.} \\
\end{array}$$

The feature matrix and tensor are depicted in [Fig 2](#pone.0193521.g002){ref-type="fig"}.

![Schematic image of vectorized and tensorial features.\
(a) Vectorized feature is included in a matrix whose first and second modes are trials and features. PWD, PHD, and PSD at each frequency band are lined as a vector in each trial. (b) Tensor structured dependent variable with four modes: trials, channel-1, channel-2, and frequency bands. Tensor elements with a larger channel-1 index than channel-2 are PWD, and a smaller channel-1 index than channel-2 are PHD, and identical channel indexes are PSD.](pone.0193521.g002){#pone.0193521.g002}

Regression analysis {#sec008}
-------------------

Higher order partial least square (HOPLS) \[[@pone.0193521.ref034]\] and standard partial least square (PLS) \[[@pone.0193521.ref040], [@pone.0193521.ref041]\] simultaneously perform dimension reduction and regression, which were used in this study. The former is a natural extension of the latter so that tensor-format features can be used.

Let us denote the response matrix by Y that has all the response variables of all training trials: $$Y = \left( {\mathbf{y}\left( 1 \right),\mathbf{y}\left( 2 \right),\ldots,\mathbf{y}\left( N \right)} \right)^{\top} \in \mathbb{R}^{N \times J},$$ where ***y***(*n*) is the *J* = 3 dimensional response vector of the *n*-th trial. All response variables were normalized to have zero mean and unit variance.

PLS performs a simultaneous decomposition of X and Y to find common latent variables $\mathbf{t}_{\mathbf{r}} \in \mathbb{R}^{N}$ as: $$X{= \sum\limits_{r = 1}^{R_{1}}\mathbf{t}_{r}\mathbf{p}_{r}^{\top} + E}$$ $$Y{= \sum\limits_{r = 1}^{R_{1}}\mathbf{t}_{r}\mathbf{q}_{r}^{\top} + F,}$$ where *E* and *F* are the residual matrices, and *R*~1~ is called the number of the components.

On the other hand, HOPLS can be similarly formulated as the problem to find latent variables as follows: $$\begin{array}{r}
{\mathcal{X} = \sum\limits_{r = 1}^{R_{2}}\mathcal{G}_{r} \times_{1}\mathbf{t}_{r} \times_{2}P_{r}^{(1)} \times_{3}P_{r}^{(2)} \times_{4}P_{r}^{(3)} + \mathcal{E}} \\
\end{array}$$ $$\begin{array}{r}
{\mathcal{G}_{r} \in \mathbb{R}^{1 \times M \times M \times F}, P_{r}^{(k)} \in \left\{ \begin{array}{ll}
\mathbb{R}^{M \times L_{k}} & \left( k = 1,2 \right) \\
\mathbb{R}^{F \times L_{k}} & \left( k = 3 \right) \\
\end{array}\operatorname{} \right.} \\
\end{array}$$ $$\begin{array}{r}
{Y = \sum\limits_{r = 1}^{R_{2}}\mathbf{t}_{r}\mathbf{q}_{r}^{\top} + V,} \\
\end{array}$$ where $\mathcal{G}_{r}$ is called the core tensor, $\mathcal{E}$ and V are the residuals, *R*~2~ is the number of the components, and ×~*k*~ denotes the *k*-mode product \[[@pone.0193521.ref042]\]. $P_{r}^{(n)}$ is called the loading matrix of the *r*-th component, and *L*~*k*~ is called the number of the *k*-mode loadings.

If data are plentiful, which is rare in EEG studies, the best approach for training and evaluating the performance of a predictive model is to randomly divide the dataset into three parts: training, validation, and test sets, which are respectively used to train a model, tune hyper-parameters or select a model, and evaluate the generalization error \[[@pone.0193521.ref043]\]. However, since the amount of data in this study is too small to exploit such an ideal protocol, we instead used leave-one-out cross-validation for each subject. The hyper-parameter *R*~1~ of PLS varied from 1 to 43, loadings of the channel-1 *L*~1~ and the channel-2 *L*~2~ ranged from 1 to 7. The loading gs of the frequency band *L*~3~ and the number of components *R*~2~ of HOPLS ranged from 1 to 5. The result of the models that achieved the best performance was reported in Results.

Evaluation metrics {#sec009}
------------------

Root mean squared error (RMSE) was used to quantify the predictability of the regression models for each subject, which are formulated as: $$\begin{array}{r}
{\text{RMSE} = \sqrt{\frac{1}{N}\sum_{i = 1}^{N}\left( {\hat{y}}_{i} - y_{i} \right)^{2}},} \\
\end{array}$$ where *N* is the number of test samples, ${\hat{y}}_{i}$ is the predicted value for the *i*-th test data, and *y*~*i*~ is the actual value.

Results {#sec010}
=======

[Table 2](#pone.0193521.t002){ref-type="table"} summarizes RMSE, and the numbers of latent factors identified by PLS and HOPLS, respectively. Predictions with tensorial features generally made smaller errors than the vectorized ones for all the rating scales. [Fig 3](#pone.0193521.g003){ref-type="fig"} reports the one hundred features that contributed to the prediction the most greatly, where feature contributions were calculated by taking the magnitude of the regression coefficients and PSD, PWD, and PHD are separately shown. PSD rarely appeared in the list of the top one hundred features list for all of the rating scales. Among the five frequency bands, the alpha band contributed the most to the MOS prediction, followed by the beta band. For the VAL and the ARL predictions, the beta band contributed the most, followed by the alpha band. The top ten channel pairs, which contributed the most to the MOS prediction extracted from subjects 1, 2, 3, and 4, are shown in [Fig 4](#pone.0193521.g004){ref-type="fig"}.

![Contribution of features.\
Feature contributions were calculated by taking magnitude of their regression coefficients. (Left) Numbers of PWD, PHD, and PSD among the one hundred features that most greatly contributed to the prediction of each rating scale among all features. (Right) Number of features of each frequency band among the one hundred features that most greatly contributed to each rating scale.](pone.0193521.g003){#pone.0193521.g003}

![Contribution of channel pairs.\
Top ten channel pairs that contributed the most to MOS predictions of subjects 1, 2, 3, and 4. This figure was made by modifying the original one, which is distributed under the public domain dedication \[[@pone.0193521.ref044]\].](pone.0193521.g004){#pone.0193521.g004}

10.1371/journal.pone.0193521.t002

###### Prediction results and the number of latent factors.

![](pone.0193521.t002){#pone.0193521.t002g}

  subject   RMSE (vector)   RMSE (tensor)   *R*~1~   *L*~1~   *L*~2~   *L*~3~   *R*~2~               
  --------- --------------- --------------- -------- -------- -------- -------- -------- --- --- --- ---
  1         1.091           1.090           1.086    0.961    0.957    0.962    1        3   4   5   3
  2         1.042           1.057           1.058    0.987    0.958    0.975    1        4   4   4   3
  3         0.997           1.018           1.020    0.907    0.936    0.931    2        7   3   4   2
  4         1.148           1.110           1.043    0.949    0.944    0.994    2        7   5   2   2
  5         1.187           1.225           1.229    1.082    1.177    1.150    1        6   7   4   5
  6         1.260           1.292           1.151    1.000    2.176    1.941    4        7   4   4   4
  7         0.979           0.981           1.007    0.869    0.935    1.167    1        2   3   5   5
  8         1.155           1.186           1.125    0.970    0.989    1.017    1        1   5   1   1
  9         1.215           1.221           1.160    0.996    0.994    1.023    2        7   7   1   2
  10        1.111           1.112           1.022    0.957    0.985    1.144    1        5   4   3   4
  11        1.125           1.243           1.0.19   1.013    1.047    0.920    3        7   7   1   1
  12        0.996           1.193           1.177    0.641    0.912    0.680    29       4   2   3   4
  13        1.258           1.227           1.234    1.051    1.050    1.035    3        4   2   5   4
  14        0.991           1.102           1.040    0.940    0.980    0.929    12       7   1   2   3
  15        1.022           0.989           0.969    0.965    0.927    0.934    1        7   3   2   5
  16        1.196           1.206           1.087    1.058    1.047    1.027    4        7   6   2   5
  17        1.021           1.083           1.087    0.884    0.886    0.924    3        4   7   4   3
  18        1.055           1.027           1.092    0.915    0.920    0.969    2        1   3   4   4
  19        1.130           1.142           1.126    1.021    1.081    1.020    1        5   1   5   4
  20        0.995           1.028           0.944    0.887    0.990    1.057    1        4   5   2   5
  21        1.121           1.157           1.103    0.900    0.969    0.997    1        1   1   4   2

Discussion {#sec011}
==========

Channel-pair-based features (PWD and PHD) contributed more to the predictions than channel-based ones (PSD), which agrees with a previous study \[[@pone.0193521.ref031]\] and suggests the importance of considering scalp EEG dynamics between brain regions, and that graph theory based features or functional connectivity analysis can be effective \[[@pone.0193521.ref045], [@pone.0193521.ref046]\]. The importance of spectral differences in caudality (DCAU) between the anterior and posterior \[[@pone.0193521.ref012], [@pone.0193521.ref047]\] or the front-posterior brain regions \[[@pone.0193521.ref031]\] as well as the lateral (left-right) spectral difference (DLAT) have been documented \[[@pone.0193521.ref028], [@pone.0193521.ref030]\]. In this study, both of DLAT and DCAU contributed to the predictions ([Fig 4](#pone.0193521.g004){ref-type="fig"}) although their effectiveness was dependent on the subjects.

Quality prediction models were independently trained for each subject in this study because emotion regulation is reportedly dependent on individuals \[[@pone.0193521.ref048]\]. The commonality of the channels/channel pairs, which greatly contributed to the predictions, was actually rather small ([Fig 4](#pone.0193521.g004){ref-type="fig"}). Therefore, creating subject-independent features is an interesting future work. However, note that the alpha and beta bands commonly contributed to the predictions, whereas the effective channels/channel pairs differed depending on the subjects. The alpha and beta bands contributed more largely to the predictions than the other frequency bands, which is in line with previous neurophysiological studies. The relationship between alpha band asymmetry and the withdrawal or disengagement from a stimulus or negative valence has been well documented in response to a variety of stimuli, including pictures \[[@pone.0193521.ref049], [@pone.0193521.ref050]\], music \[[@pone.0193521.ref031], [@pone.0193521.ref047], [@pone.0193521.ref051]\], movies \[[@pone.0193521.ref052]\], and speech \[[@pone.0193521.ref011], [@pone.0193521.ref013]\]. The beta band, which contributed the most to the ARL predictions, is reportedly associated with arousal and emotional experiences \[[@pone.0193521.ref053], [@pone.0193521.ref054]\].

Gupta et al. \[[@pone.0193521.ref013]\] predicted MOS values using the same data set that we used in this study. Their study used a simple linear regression model with not only EEG but also speech features. They reported the RMSE of their model was 0.117, which is much lower than our model, and suggests that speech features are much more informative than EEG features to predict subjective quality ratings.

Although we predicted the response values of MOS, VAL, and ARL, other perpetual dimensions were also proposed recently to model emotions or perceived quality-of-experiences \[[@pone.0193521.ref055], [@pone.0193521.ref056]\], which should be investigated in future research.

Neither previous work nor our current study advocate that physiological assessment methods of speech quality should replace subjective rating methods or signal analysis methods because, as stated in Introduction, each method has its own advantages and disadvantages and they can complement each other.

Several open questions remain. First, features were extracted and constructed as tensors as described in Feature Extraction and Regression Analysis, but other features and construction ways are also possible. For example, if time-frequency analysis is employed, times frames can be treated as one of the tensor modes. Second, this study analyzed the overall quality of each speech stimulus longer than ten seconds. However, parts of speech can affect much more largely its overall perceived quality. Therefore, analysis methods to specify such parts need to be studied.

Conclusion {#sec012}
==========

This study predicted the subjective quality ratings of synthesized speech solely based on EEG. We created vectorized and tensorial features for the regression that include channel-based and channel-pair-based features at multiple frequency bands. The experimental result showed that tensorial features more effectively predicted the subjective ratings than the other, and the trained predictive models were neurophysiologically plausible.

[^1]: **Competing Interests:**The authors have declared that no competing interests exist.
