{'fraction_non_alphanumeric': 0.03245743287491814, 'fraction_numerical': 0.010334806810740014, 'mean_word_length': 4.948265368228849, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Robust and Complex Approach of Pathological Speech Signal Analysis

This paper presents a study of the approaches in the state-of-the-art in the field of pathological speech signal analysis with a special focus on parametrization techniques. It provides a description of 92 speech features where some of them are already widely used in this field of science and some of them have not been tried yet (they come from different areas of speech signal processing like speech recognition or coding). As an original contribution, this work introduces 36 completely new pathological voice measures based on modulation spectra, inferior colliculus coefficients, bicepstrum, sample and approximate entropy and empirical mode decomposition. The significance of these features was tested on 3 (English, Spanish and Czech) pathological voice databases with respect to classification accuracy, sensitivity and specificity.


Introduction
Voice Pathology description and characterization has always demanded attention from the physiological and medical fields [1], as well as from the voice function point of view [2]. The term voice is described by [3] both in a broad and a narrow sense. In the broad sense voice may be taken as synonymous of speech, therefore terms as Voice over IP (VoIP) can be found in the literature and media with the meaning of speech data on internet. In the narrow sense voice refers to the vibration of the vocal folds. Speech sounds resulting from the interaction of this vibration with the Oro-Naso-Pharyngeal Tract (ONFT) are referred to, as voiced. Speech sounds produced by turbulent flow within the ONFT are termed voiceless. Phonation is the recommended term to refer to vocal fold vibration. A speaker showing an anomalous vocal fold vibration pattern is referred as dysphonic, and aphonic if there is no vocal fold vibration at all. If no anomalies are present the speaker is referred as normophonic.
Dysphonic voice is a perceptual and subjective term associated to Voice Pathology. Dysphonia is the perceptual quality of voice signaling that something wrong is happening in the phonation organs (mainly the larynx and its associated structures). Voice Pathologies and Dysphonic Voice are thus intrinsically related. The classification of Voice Pathologies or Disorders is described in [3] as tissue infection (e. g. laryngitis, bronchitis, croup,. . . ), systemic changes (e. g. dehydratation, pharmacological and drug effects, hormonal changes,. . . ), mechanical stress (e. g. vocal nodules, polyps, ulcers, granulomae, laryngocele, hemorrhage,. . . ), surface irritation (e. g. laryngitis, leukoplakia, gastroesophageal reflux,. . . ), tissue changes (e. g. laryngeal carcinoma, keratosis, papillomas, cysts,. . . ), neurological and muscular changes (e. g. bilateral and unilateral vocal fold paralysis, Parkinson's Disease (PD), Amyotrophic Lateral Sclerosis (ALS), myotonic dystrophy, Huntington's Chorea, myasthenia gravis,. . . ), and abnormal muscle patterns (e. g. conversion aphonia or dysphonia, spasmodic dysphonia, mutational dysphonia, ventricular phonation,. . . ). In [4] a description of vocal pathologies can be found. A last important group of neurological diseases which leave a correlate in voice and speech is that of cognitive origin, Alzheimer's Disease (AD) being the most relevant as polyps, nodules, edemae, cysts, etc., where full closure of the glottal gap is not granted by vocal fold adduction and produce acoustic correlates as open and close phase perturbations, recovery phase attenuation, harmonic display reduction, etc. Other pathologies, especially those of neurological origin produce deviations in biomechanical parameters, as vocal fold tension (hypo-and hyper-tonic) and tremor, which can also be manifested as modulations in amplitude or frequency, and as changes in the vocal fold tension. As many organic pathologies induce a hyper-tonic behavior, the main problem when dealing with correlates produced by neurological pathology is to differentiate their origin from that of organic origin. Contact defect pathologies can also show asymmetric vocal fold vibration, and this may also be the case in aging voice (presbyphonia). To establish differentiation criteria for using acoustic correlates when dealing with organic, neurologic or aging-induced perturbations, is a major issue in AVQA [21].
Another problem regarding AVQA is the lack of good reference baselines to establish the methodologies for voice pathology classification from acoustic analysis. Rigorous databases, acquired to represent each of the different pathologies under well-defined sample population size and recording conditions are scant.
Many times the problem comes from the simultaneous presence of different larynx pathologies in the same patient, either related or unrelated (e.g. it is common to find a counter-lateral lesion as a consequence of a unilateral polyp). The problem then is how to decide if a specific acoustic correlate is produced by one cause or another. Most of the times acoustic databases are produced by laryngological services as part of examination protocols [22,23], but this is not always the case, as many laryngologists prefer to depend on visual exploration, neglecting the possibilities offered by AVQA for different reasons [24]. Speech therapists rely more on acoustic exploration, but many times it is mainly restricted to measurements of vocal effort, long term frequency analysis, respiratory efficiency, or distortion measures. Therefore voice records from speakers ranked by etiology and severity index, using compatible standards (digitalization, channel, microphones) are scant. Most of the available databases are either incomplete, inconsistent (made up of recordings taken under different conditions) or deficient (many non-frequent pathologies are not well represented in sample size). Besides, there is a lack of good normative databases, as most of the records produced by medical services contain information from dysphonic voice, but normative speakers are absent or poorly represented (this is the case of the most widely used database [25]). This is a severe limitation for systematic AVQA. Another problem is cross-lingual representation. As far as sustained vowels are concerned, this would not be a problem, but it becomes a major obstacle when segmental parameters are involved, as in the use of passages (either read or spontaneous).
After all these considerations, it may be said that the aim of AVQA is to design the best methodology to use Voice Correlates in Voice Pathology detection, grading and classification. Several steps are to be assumed and methodologically formulated for such. Voice Pathologies and Correlates have to be well connected by adequate physiological and biomechanical modeling of larynx dynamics to understand how acoustic correlates and organic dysfunctions could be associated [26]. This requirement is of vital importance in meeting the challenge of pathology classification from acoustic correlates. On its turn, inverse methods to obtain robust and significant estimates of voice correlates are of great relevance [27]. These are to be combined with other parameterization methods based in direct time and frequency domain estimates to produce rich feature sets [28]. Finally, good Statistical Pattern Matching and Machine Learning methodologies have to be used to reduce redundancy and to make a selection by relevance within the feature set, to determine the best combinations for each specific task (detect, grade, classify), and procure specific results under quality criteria (sensitivity, specificity and accuracy) [29]. These methods have to be contrasted on generally-accepted benchmark databases [30].
Having in mind all these conditions the present study is aimed to describe a wide set of tests on different databases, using an exhaustive set of acoustic features, by means of well-known classifiers to estimate the performance of the methods and features regarding specificity, sensitivity and accuracy in dysphonic voice detection tasks. There are already many publications describing the methods of pathological speech detection. Usually the authors use just a limited set of speech features and most of the publications present the detection accuracies tested only on one monolingual database [31][32][33][34][35][36][37]. Probably the widest range of features describing different aspects of speech was tested in a work of Tsanas et al. [34]. Regarding the databases there are just a few publications considering 2 different data sets [38][39][40][41].
To sum up the state-of-the-art approaches in the field of pathological voice analysis there is still a lack of publications providing a complex overview of features quantifying pathological speech and providing strong conclusions supported by a robust testing. Therefore this work has 4 main goals: 1. According to complex parameterization and consequent robust testing identify features that have the largest discriminative power in the field of pathological speech analysis.
2. Design new features that can quantify hoarseness, breathiness and non-linearities in pathological speech signals.
3. Prove that the proposed large set parameterization approach can provide better classification results (with respect to classification accuracy, sensitivity and specificity) than those published in the field of pathological speech analysis by the other researchers. 4. Select a database that has high potential for the future, especially in terms of speech features design, tuning and testing.
The paper is organized as follows: the classically used features as well as the newly designed ones are discussed in section 2, section 3 describes the 3 databases, section 4 describes the testing procedure, section 5 is devoted to the experimental results and discussion. The conclusions are given in section 6.

Features
The aim of this work was to explore significance of the mostly used speech features when focusing on the ability of differentiation between healthy and pathological speech. The features that are usually known in different fields of speech signal processing (speech recognition, enhancement, denoising, identification etc.) and newly designed features originally introduced in this paper will be investigated as well.
Due to the limited size of this paper the features that were not originally introduced in this work will be mentioned without their deeper description by an algorithm. However each parameter will be accompanied by a reference where the reader can find further information.

Features Describing Phonation
Probably the most popular features describing pathological voice are fundamental frequency F 0 and parameters describing its variability in time (jitter): PPQ5 (five-point Pitch Perturbation Quotient), RAP (Relative Average Perturbation), jitt loc (average absolute difference between consecutive periods, divided by the average period), jitt abs (average absolute difference between consecutive periods), jitt ddp (average absolute difference between consecutive differences on neighbour glottal periods, divided by the average period) [39,[41][42][43]. These features are good candidates especially for quantification of voice tremor [44].
A disadvantage of previously mentioned measurements is that the values of the features are highly dependent on gender and variable acoustic environment. To overcome this disadvantage Little et al. proposed the PPE (Pitch Period Entropy) [45]. During the calculation of PPE, logarithmic semitone scale, inverse filtering and entropy estimation are incorporated. Another method which is close to jitter is a measure of standard deviation (std) of the time that vocal folds are apart (GQ open ) and in collisions respectively (GQ closed ) [34].
To effectively describe hypophonia or intensity perturbations, short-time energy E or pitch-level variations (shimmer) can be used [34,42,46]. In this work 6 kinds of shimmer will be considered (they are calculated similarly to jitter but the intensity is used): APQ3 (three-point Amplitude Perturbation Quotient), APQ5, APQ11, shimm loc , shimm ddp , shimm dB (average absolute base-10 logarithm of the difference between the amplitudes of consecutive periods) [43,47].
Another feature describing speech intensity is TKEO (Teager-Kaiser Energy Operator) [48]. The advantage over simple E is that it takes into account also signal frequency. It has been shown that speech contains dominant modulation frequencies in a range 2-20 Hz with maximum at approximately 4 Hz [49].
The 4 Hz modulation energy (ME) was selected as a feature which is related to a measure describing energy distribution in power spectra. Similarly MPSD (Median of Power Spectral Density), usually called median frequency, was selected [50]. The last feature in this category is LSTER (Low Short-Time Energy Ratio) [51]. This feature is usually used for differentiation between speech and music signals because speech exhibits higher variations in energy per 10-30 ms frames. However this feature was selected for analysis of pathological speech as well. It is considered that patients will reach higher values of this feature in the case of maintained vowels due to the inability to sustain the same amount of airflow during the whole phonation.

Features Describing Tongue Movement
Frequencies of first three formants F 1 , F 2 , F 3 and their bandwidths B 1 , B 2 , B 3 are related to volumes of vocal tract cavities. Especially the volume of throat and oral cavity is modified by the tongue position.
According to this fact it is possible to consider formants as a measure of tongue movement [52,53].

Features Describing Speech Quality
Signs of vocal fold dysfunctions are usually associated with breathiness or hoarseness [32]. From signal theory point of view these dysfunctions are characteristically decreasing voice quality, under a simplified consideration, by additive noise. This implies that methods based on speech quality measurements can suitably quantify vocal folds impairment and describe its progress. Besides breathiness and hoarseness these methods can be also used for analysis of hypernasality caused by an improper work of soft palate [54].
Probably the simplest quality measure feature is ZCR (Zero-Crossing Rate) and its modification HZCRR (High Zero-Crossing Rate Ratio) which takes into account a variation of ZCR in time [51]. Another simple measure is FLUF (Fraction of Locally Unvoiced Frames) which can describe an impossibility of carrying out periodical glottal closure [55].
The next three measures are based on the variation of spectrum values between adjacent frames. Specifically these are SF (Spectral Flux) [56], SDBM (Spectral Distance Based on Module) and SDBP (Spectral Distance Based on Phase) [55].
We have also used several features based on cepstral analysis. It was shown that cepstrum and its rahmonics correlate well with the perception of breathiness [57]. Moreover from the nature of cepstrum it can be said that it is a kind of periodicity measure and therefore it should also predict roughness. The most famous feature based on the real cepstrum is CPP (Cepstral Peak Prominence) originally introduced by Hillenbrand et al. [57]. Besides this feature, PECM (Pitch Energy Cepstral Measure) [55] and VR (Variation in Ratio between the second/first harmonic within the derived cepstral domain) were also used [55].
The last feature in this category is SSD (Segmental Signal-to-Dysperiodicity ratio). It was shown that this feature correlates strongly with the perceived degree of the speaker's hoarseness [38].

Segmental Features
Although some of the features mentioned in the other categories can be denoted as segmental as well (they are calculated from 10-30 ms segments), here the segmental features are considered as matrices (not only vectors) calculated from the whole signal. Probably the most popular segmental features in the field of speech signal analysis are MFCC (Mel Frequency Cepstral Coefficients) [33,36,42,54]. The advantage of these coefficients is that they can indirectly detect slight misplacements of articulators [34]. In this work 20 MFCC coefficients are extracted. The coefficient number zero is replaced by an estimate of log-energy.
Although MFCC are generally used, there is a lack of publications that compare these features to other segmental ones for the purpose of pathological speech analysis. Therefore we decided to include other segmental features. We extracted 20 mel frequency cepstral coefficients but in this case the bank of triangular filters was adjusted to the equal loudness curve [61]. We call these features MFCCE. The other two sets of features derived from MFCC are LFCC (Linear Frequency Cepstral Coefficients) and CMS (Cepstral Mean Subtraction coefficients). In the case of LFCC the bank of triangular filters is equidistantly spread in the frequency scale [62]. CMS is a kind of standardized z-score of MFCC (subtraction of mean and division by std over the time). [36,63]. These features can capture a class of source mechanism characteristics related to voice quality [33].

MSC (Modulation Spectra Coefficients) can provide information complementary to MFCC
In the next step features based on linear prediction were extracted. LPC (Linear Predictive Coefficients) [42], PLP (Perceptual Linear Predictive coefficients) [64], LPCC (Linear Predictive Cepstral Coefficients) [65], LPCT (Linear Predictive Cosine Transform coefficients) [42] and ACW (Adaptive Component Weighted coefficients) [65] were tested. The advantage of PLP over simple LPC or MFCC is that it also takes into account an adjustment to the equal loudness curve and intensity-loudness power law [64]. The advantage of LPCC and LPCT over "classic" LPC is that a transformation is used into the cepstral domain and thus the values do not correlate much. The advantage of ACW is that these coefficients are less sensitive to channel distortion [65].
The last segmental features used in this work are parameters that analyze amplitude modulations in voice using a biologically-inspired model of the inferior colliculus [66]. These features are called ICC (Inferior Colliculus Coefficients).
Segmental features are sometimes extended by their 1 st and 2 nd order regression coefficients (∆ and ∆∆ respectively). In this work used the ∆ coefficients.

Features Based on Bispectrum
Alonso et al. proved that a greater presence of quadratic coupling is observed in healthy voice when comparing it to the pathological one [55]. It is probably due to a fact that healthy voice is characterized by a vocal tract which is more non-linear than in the case of pathological voices. This quadratic coupling can be appropriately described by bispectrum and features derived from this 2D signal.

Features Based on Wavelet Decomposition
The wavelet transform is widely used especially in the field of coding and speech denoising. However its application can be found in the field of pathological speech analysis too [67]. Detail coefficients after the decomposition can be used to estimate the present noise and consequently it is possible to calculate SNR (Signal-to-Noise ratio). In fact this is just another method of voice quality measurement.
According to some measurements we have empirically selected 7 wavelets: 10 th , 15 th , 20 th -order Daubechies wavelets; 10 th , 15 th , 20 th -order symlet wavelets and 5 th -order coiflet wavelet. We will mark the features derived from the detail coefficients of the wavelet transform as SNRW(wvl ), where wvl corresponds to the specific wavelet (e. g. daub15).

Non-linear Dynamic Features
Tension on the vocal folds can significantly differ in the case of pathological speech. Voice becomes aperiodic, noisy-like, and it is very difficult to find any regularities in the signal. There is a frequent presence of sub-harmonics and chaos which can lead to a failure of conventional techniques of speech signal analysis. However it was shown that these kinds of signals can be sufficiently described by non-linear dynamical analysis [32,37,40,54].
The first representative in this category is CD (Correlation Dimension) which statistically measures attractor geometry in the phase space. This measure is related to a number of independent variables necessary for generating the attractor [32,45,47]. Another dimension measure FD (Fractal Dimension) is based on a number of basic building blocks that form a pattern [32,68]. We will also use complexity measures like ZL (Ziv-Lempel complexity) which quantify the regularity embedded in a time series [69].
Possible long-term dependencies in the speech signal will be described by HE (Hurst Exponent) [54].
Another set of measures is based on entropy. We will use SHE (Shannon Entropy), RE (second-order Rényi Entropy) [70], CE (Correlation Entropy) [40,70], RBE1 (first-order Rényi Block Entropy) [40], RBE2 (second-order Rényi Block Entropy) [40], AE (Approximate Entropy) [32,71,72], SE (Sample Entropy) [72] and RPDE (Recurrence Probability Density Entropy) [37]. Generally the entropy is a measure of uncertainty and it can be used to quantify the complexity of a system. Rényi entropies quantify the loss of information in time in a dynamic system [40], correlation entropy gives an indication of the predictability of the nonlinear time series [70] and RPDE represents the uncertainty in the measurement of the pitch period [37]. The only difference between AE and SE is that SE does not evaluate a comparison of embedding vectors with themselves.
Another measure we have considered in this work is FMMI (First Minimum of Mutual Information function) which was found by Henriquez et al. as a feature that better discriminates among the different voice qualities of the multiquality database [40]. To include also a measure of sensitivity to an initial condition, the LLE (Largest Lyapunov Exponent) was selected [32,54]. We also used detrended fluctuation analysis (DFA) to characterize the self-similarity of the graph of a signal from a stochastic process [34,37].
In this field NSE (Normalized Scaling Exponent) and FA (Fluctuation Amplitudes) were evaluated.

High-level Features
The features that are calculated for each speech segment separately form a vector or a matrix at the output of the parametrization process. This representation must be then transformed to a scalar value to be able to carry out the next processing like statistical analysis, classification etc. This is usually done by an extraction of some kind of statistics. These statistics are called high-level features while parameters extracted directly from the speech signal are called local features. If the local feature is represented by a matrix, then the high-level feature is calculated for each row separately (this monitors feature changes in time). We have extracted 60 high-level features: • max, min, position of max, position of min, relative pos. of max, relative pos. of min • range, relative range, interquartile range, rel. interquartile range, interdecile range, rel. interdecile range, interpercentile range, rel. interpercentile range, studentized range processing. However during the experiments these parameters will be extended by the other 36 features that are originally introduced in this work.

Features Based on Modulation Spectra
An initial step of modulation spectra calculation employs a short-time Fourier transform (STFT) of the discrete input speech signal s[n] with length N : where * corresponds to the average operator over m. Next the frequency analysis of subband envelopes is performed using the discrete Fourier transform (DFT): where p and l denote the acoustic and modulation frequency respectively. In the last step of modulation spectra extraction the second power of each subband is taken and normalized which partially suppresses the effect of training and testing conditions mismatch: The features we are proposing are based on a function ψ[l] which is extracted from the normalized modulation spectra according to: where l 5 Hz is a sample corresponding to 5 Hz modulation frequency (this limit was empirically found). MFP is defined as: where MFP is given in Hz. Finally RPHM can be calculated according to the following expression: where r is a linear regression line of ψ[l] for l = l 5 Hz , . . . , M − 1 and i is the index at which the function ψ[l] reaches maximum value.

Features Based on Inferior Colliculus Coefficients
Using inferior colliculus coefficients (ICC) it is possible to extract the frequency content of the modulation envelopes applied to different bands of an auditory stimuli [66]. Similarly to modulation spectra, the ICC extraction process employs at the beginning STFT and consequent calculation of power spectra |S[k, m]| 2 .
However in the next step instead of the bank of the triangular filters, the bank of P mel-spaced gammatone filters is applied. We used filters defined by the the impulse response: where f c is the center frequency in Hz, o the filter order (in our case 4) and f s the sampling frequency in Hz. The DFT on the subband envelopes was applied next, this time without prior to normalization: Finally the magnitude spectrum of each envelope |T [p, l]| is filtered by a bank of Q = 13 resonance filters with exponentially spaced frequencies from 12 to 107 Hz. In this work we used resonance filters defined by the following transfer function: At the output of this process a matrix Ξ[p, q] (q = 1, 2, . . . , Q) of ICC for the input signal s[n] is extracted.
The proposed features are based on a function ξ[p] which is extracted from Ξ[p, q] according to: Contrary to ψ[l] this function reflects the misplacement of articulators better than the instability of vocal fold vibrations. In Fig. 2  .
Similarly to RPHM the RPHIC is extracted using the regression line but in this case r[p] is calculated for p = 13, . . . , 20: heal. path.

Features Based on Bicepstrum
We will use a definition of the real bicepstrum c[n 1 , n 2 ] similar to the definition of the real cepstrum: where the bispectrum B[k 1 , k 2 ] is calculated as the DFT of triple correlation or circular triple correlation γ[n 1 , n 2 ] [75]. In our work we used the second approach based on γ[n 1 , n 2 ]: The first feature we are proposing in the present work is BCII (BiCepstral Index Interference): where c m [n 1 , n 2 ] is the bicepstrum calculated from the m th speech segment. The two features introduced next are HFEBC (High Frequency Energy of one-dimensional BiCepstral index) and LFEBC (Low Frequency Energy of one-dimensional BiCepstral index): where ρ[n] (for n = 0, 1, . . . , N − 1) is the one-dimensional bicepstral index and f max the maximum expected fundamental frequency (in our case f max = 350 Hz).
It is supposed that pathological voice contains much more white noise (symmetrically distributed) than healthy voice. This is especially due to incorrect glottal closure. This kind of noise disappears in the real cepstrum estimated by ρ[n] (see Fig. 3), therefore a difference between the real cepstrum c[n] (estimated using DFT) and the one-dimensional bicepstral index can estimate noise components of the analyzed signal.
According to this idea BCMII (BiCepstrum Module Interference Index) and BCPII (BiCepstrum Phase Interference Index) are proposed: where c m [n] and ρ m [n] are the m th frame real cepstrum and the one-dimensional bicepstral index respectively.
According to these equations we also propose the features LSBER (Low Spectra/Bispectra Energy Ratio) and − |c m [n 1 , n 2 ]|| , wherec m [n 1 , n 2 ] is the complex bicepstrum calculated without the absolute value in eq. (16). Similarly to eq. (29) and (30) features BMD (Bispectral Module Distance) and BPD (Bispectral Phase Distance) are also

Different Kernel Based Approximate and Sample Entropy
Approximate entropy (AE) is a measure of regularities inside the analyzed time series. The advantage of AE is that it can robustly estimate the system complexity using just a limited number of samples (100-5000) [76]. To define AE we need to firstly reconstruct a state space using Takens' embedding theorem: [77] x is an embedding vector, m is an embedding dimension and τ is a time delay (τ = 1 for AE).
Next, the regularity quantity of a particular pattern C[i, m, r] is defined as: Finally AE can be calculated according to expression: [32] A disadvantage of AE is its dependence on the signal length due to the self comparison of points in the attractor. This fact can be avoided using the sample entropy (SE) which does not evaluate the comparison of embedding vectors among themselves: [78] The usual values for embedding dimension are m = 1, 2. In this work we used m = 2 which provides more detailed reconstruction. Finally, a function must be defined κ (i, j, r) we call kernel function. Originally the AE or SE are based on: where Θ is the Heaviside function and r a radius in our case calculated according to: We will denote these entropies as AE(Heaviside) and SE(Heaviside) respectively. Orozco-Arroyave et al.
proposed AE(Gaussian) and SE(Gaussian) based on the Gaussian kernel: [78] κ In this work we propose AE and SE based on the other 6 kernels: exponential kernel Laplacian kernel, circular kernel, for x[i] − x[j] < r, zero otherwise and triangular kernel, where µ i is a parameter, or mean sequence value, calculated from the original i th IMF and I is the total number of the IMFs. In the case ofμ the i th IMF was logarithmized before the consequent parametrization.
To extract the sequence from IMF, Tsanas et al. used SEO (Squared Energy Operator) and TKEO. As a parameter they also used SHE.
We have extended this idea on these parameters: IMF-SNR RE (based on second-order Rényi Entropy), IMF-SNR ZCR (based on Zero-Crossing Rate) and IMF-NSR RE . In the case of IMF-SNR RE µ i is defined as: where p 2 x i j is the probability P IMF i = x i j and x i 1 , x i 2 , . . . , x i J are the possible values of the i th IMF. In the case of IMF-SNR ZCR µ i is calculated according to: where f i [n] is the i th IMF.
The time-varying high frequency components present in the 1 st IMF represent the noise part of the speech signal. We propose a new feature IMF-FD which is based on the fractal dimension calculated from the 1 st IMF. This complexity measure appropriately quantifies the amount of noise present in the signal and indirectly describes hoarseness or breathiness. The feature is defined as: IMF-FD = log 10 N log 10 N + log 10 Another feature based on the first intrinsic mode function is IMF-CPP (Cepstral Peak Prominence extracted from the 1 st IMF). It can be calculated as: 6. Pick the maximum from all the maxima in 5.

Databases
To provide robust results 3 (English, Spanish and Czech) databases have been used during the testing procedure. Each database represents a different language group (Germanic, Romanic and Slavic). This approach is advantageous from the cultural difference point of view. Speakers of different languages exhibit especially different prosodic characteristics. The aim of this work is to find features significant for the particular language, but we have focused on a selection of features that are language-independent as well.

MEEI Disordered Voice Database
The Massachusetts Eye and Ear Infirmary (MEEI) database [25] has been for many years used as are used.
Although MEEI consists of approximately 700 speakers in total, the age distributions between normal and pathological speakers are not matched. Therefore the amount of pathological speakers is usually limited to 173 according to criteria published by Parsa and Jamieson (we call this "limited version of MEEI database") [15]. The statistical characteristics of MEEI database used in this work can be found in Table 1.
Although this database is very popular and very often used, its size and data are considered insufficient.  Table 2 can be found in sec. 3.4. are sampled at f s = 25 kHz. The statistical characteristics of this database can be found in Table 3. Table 4 shows detection results (obtained on the PdA database) of two available works. As can be seen, the accuracies are not so high as in the case of MEEI database. Moreover, the PdA consists of more speakers than the limited version of MEEI. All these facts highlight the high potential of the PdA for future use.

PARCZ Database
The last database we have included in our test is the Czech Parkinsonian Speech Database (PARCZ) recorded at St. Anne's University Hospital in the Czech Republic. This database consists of 52 healthy   characteristics of the PARCZ database can be found in Table 5. As it can be seen, contrary to the MEEI or PdA, PARCZ is more focused on elder people.

The Rule of 30
To decide whether the corpus size is sufficient for the robust conclusions Doddington et al. introduced "the rule of 30" which comes directly from the binomial distribution, assuming independent trials [85].

Experiments
All the databases were resampled to f s = 16 kHz and in dependence on the next processing the data have been divided into 9 groups. Two approaches have been considered: gender-dependent and gender independent. Each database is randomly divided into 75 % and 25 % training and testing subsets respectively.
The classifier is evaluated consequently. This procedure (data split, classifier tuning and evaluation) is repeated 100 times. The resulting accuracy (ACC), sensitivity (SEN ) and specificity (SP E) are calculated according to: where T P (True Positive) and F P (False Positive) represent the number of correctly identified pathological speakers and number of speakers diagnosed as pathological, but being healthy. Similarly, T N (True Negative) and F N (False Negative) represent the total number of correctly identified healthy speakers, and pathological speakers evaluated as healthy controls.
Before classification the training data are z-score normalized. The testing data are normalized by subtracting the training set mean and dividing by the training set standard deviation for each feature.

Parameterization Programs
Several toolboxes and programs have been used for the purpose of feature extraction. Features based on the detrended fluctuation analysis have been calculated using FastDFA [86]. The glottal quotients (GQ open and GQ closed ) were extracted using the algorithm DYPSA implemented in VOICEBOX [87]. To estimate the largest Lyapunov exponent (LLE) TSTOOL has been used [88]. The software Praat has been used to estimate the fundamental frequency F 0 , all kinds of jitter and shimmer, harmonic-to-noise ratio (HNR) and  [91].
After the evaluation a list of the ten most significant features is drawn up. Finally the density estimation plots (computed using kernel density estimation with Gaussian kernels) of the most significant features in all scenarios are given.

Classification Methods
In our test we have used 2 classifiers: SVM (Support Vector Machine with a radial kernel) and RF (Random Forest). Regarding the SVM, the parameter kernel gamma γ and penalty parameter C were optimized using a grid search for possible values.

Results
A summary of pathological speech detection results expressed by accuracy, sensitivity and specificity can be found in Table 6. This table also provides some statistics related to the number of features selected in each scenario. In the case of MEEI database the accuracy, sensitivity and specificity were equal to 100.0 ± 0.0 % in all scenarios (considering both genders together and separately) when using an RF classifier. A discussion focused on the credibility of these results is given below. In comparison to RF the SVM classifier provided slightly worse results.
In the case of the PdA database the best results were found when classifying by RF too.  The ten most significant features selected by Mann-Whitney U test in all considered scenarios can be found in Tables 7 -15. The density estimation plots (computed using kernel density estimation with Gaussian kernels) of the most significant features in these scenarios can be seen in Fig. 4. The MEEI database itself introduces many issues related to the credibility of the results. The main problem is that both the healthy and disordered speech were recorded in a different way. First of all the healthy one was sampled at f s = 50 kHz and the disordered one at f s = 25 kHz. In our experiment all databases were resampled to f s = 16 kHz but some small differences will remain among signals. But probably the most serious problem is that the duration of a sustained vowel is 3 s and 1 s for healthy and disordered voice respectively. Therefore all high-level or local features that somehow reflect the signal length (entropy, range, std, duration, etc.) can provide very good discriminative results. This can also be a case of the relative interpercentile range mentioned in Table 7 and 9. Although it is a relative measure it is still dependent on the length of the input vector. We made a small experiment repetitively and randomly generating vectors with length N = 10 i for i = 1, 2, 3, 4 and with a normal distribution. After that we calculated the relative interpercentile range of these vectors and found out that the value of this measure is increasing with the decreasing vector length.
To avoid the problem of different vowel durations one can skip features dependent on this property but we would loose parameters that can be potentially very good candidates for pathological speech identification.
The other approach is to take only a one-second segment from the healthy speech but in our opinion this is not a good solution. The first second segment could be selected but then we are losing the sustained part of the signal (used for estimating jitter) and phonation trail. If we take only the sustained part, then we will loose information about phonation onset and offset. It has been shown that an analysis of these parts is useful for example for the dysarthric speech description [90,92]. Therefore there is no elegant solution that would not distort the results. Moreover, Malysla et al. mention that some speakers were recorded in different sites and over potentially different channels [66]. In conclusion, although the MEEI database is a very popular benchmark in the field of pathological speech analysis, the results obtained using this database should be taken very carefully. An introduction of a new English disordered voice database is thus highly important for the future evaluation of new state-of-the-art speech signal processing techniques.
In scenario P1 (PdA, males) all the most significant features are based on UCPP (Unsmooth Cepstral Peak Prominence). On the other hand, in the case of scenarios P2 (females) and P3 (both genders) all the      lower than in the case of MEEI or PdA databases. This is also closely related to the poor classification results mentioned in Table 6.
To summarize the discussion about the databases we can say that the MEEI database should no longer be used as a benchmark. The more challenging one is PARCZ database, however it has some disadvantages from the scientific point of view: it is in Czech language, which is not a widespread language in the world; it is focused only on PD people with hypokinetic dysarthria; and it is not really suitable for a binary classification (at least 4 classes should be considered). Probably the most suitable database for the evaluation of pathological speech identification methods is PdA. The classification accuracies are still challenging and they can be significantly improved. Another advantage of this database is that it is freely available for   when classifying the pathological speech generally (but still this should be proved by robust testing in future).  Our initial idea was also to try the inter-database classification (training the classifier on one database and testing it using the other one). However the classification results were very poor. It was caused mainly by these two facts: 1) When we trained the classifier on MEEI database, the features reflecting the signal length were selected. But these features were not so significant when testing them on PdA or PARCZ databases.
2) The PARCZ database is very different from MEEI or PdA, because it contains only one specific speech disorder (hypokinetic dysarthria), while the other two databases contain many different voice pathologies. In this work we have introduced 36 new speech features. Just one (IMF-CPP) has been mentioned among the 10 most significant ones (Table 11 and 12), however the rest of them are significant (p < 0.05) as well, see

Conclusions
This work provides an insight into the robust and complex approach of pathological speech analysis. To our best knowledge this is the first contribution providing a complex evaluation of feature significance from different fields of speech signal processing (e. g. speech analysis, recognition, coding, enhancement etc.). It is also the first contribution deriving conclusions according to robust tests where 3 (English, Spanish, Czech) databases were used. These languages belong to 3 different language groups (Germanic, Romanic, Slavic).
In general, the work has 4 goals, yet each of them has its conclusion.
1) According to complex parameterization and consequent robust testing identify features that have the largest discriminative power in the field of pathological speech analysis. Unfortunately most of the works published in the field of pathological speech analysis provides conclusions based on a limited set of parameterization methods. In other words there is still a lack of publications providing a complex overview of features quantifying pathological speech and providing strong conclusions supported by a robust testing.
Our work is unique in this way, because together it provides testing based on 128 local features. We made the wide overview of all these features used for pathological speech quantification so that the researchers can find out what exactly the specific feature quantifies and how it can be implemented. We used features describing phonation, tongue movement, speech quality (including non-linear dynamic features and features based on bispectrum/bicepstrum, empirical mode decomposition, wavelet decomposition) and segmental features.
Using the non-parametric Mann-Whitney U test we observed that among all parameterization techniques those based on segmental features provide the best classification results. To our best knowledge this is the first work that tested the significance of 11 segmental features. The largest discriminative power can be obtained thanks to especially segmental features like modulation spectra coefficients, adaptive component weighted coefficients and linear predictive cepstral coefficients. Although the segmental features are not very frequent when analysing pathological speech, their potential seems to be high. However, their disadvantage is that they are usually difficult to be interpreted clinically. This is probably the reason they are not frequently used. Although they can provide good classification results they don't say much about specific pathology or speech dysfunction.
2) Design new features that can quantify hoarseness, breathiness and non-linearities in pathological speech signals. Clinical signs of vocal fold dysfunctions are usually associated with breathiness or hoarseness.
Moreover voice can become aperiodic, noisy-like, and it is very difficult to find any regularities in the signal.
Sometimes there is a frequent presence of sub-harmonics and chaos, which can lead to a failure of conventional techniques of speech signal analysis and which requires new parameterization methods developed specifically for pathological speech description.
We introduced 36 new measures based on modulation spectra, inferior colliculus coefficients, bicepstrum, sample and approximate entropy and empirical mode decomposition. Features based on modulation spectra quantify instability of vocal fold vibrations and complements features based on inferior colliculus coefficients that reflect the misplacement of articulators. Due to incorrect glottal closure the pathological voice contains much more white noise that can be effectively quantified by proposed features derived from bicepstrum.
New features based on empirical mode decomposition are able to describe the noise component of analysed signal as well. Finally we proposed different kernel-based approximate and sample entropies to measure regularities inside the signal.
These novel features were statistically processed by the non-parametric Mann-Whitney U test. All of them have been identified as significant and they have passed the feature selection process in at least one database. Moreover, some of them were listed among top ten significant features selected in specific scenario (they outperformed the other conventional features). In other words they helped to improve classification results in terms of accuracy, sensitivity and specificity and due to their effective quantification abilities they have high impact on the future work.
3) Prove that the proposed large set parameterization approach can provide better classification results (with respect to classification accuracy, sensitivity and specificity) than those published in the field of pathological speech analysis by the other researchers. We tested the significance of all the mentioned parameters on 3 (English, Spanish, Czech) databases. In the case of the Massachusetts Eye and Ear Infirmary (MEEI) database we get accuracy, sensitivity and specificity equal to 100.0 ± 0.0 %, which are the best results that have been published in the frame of this database (Henriquez et al. reached 99.69 % [40]). However, we are very critical with these results. Therefore we discussed their trustability and the viability of using the MEEI database as a benchmark for pathological speech signal analysis.
The results obtained with the PdA database are more challenging. When we considered the singleclassifier approach, we reached the accuracy 82.1 ± 3.3 %, which is the best among the published numbers Regarding the last Czech Parkinsonian Speech Database (PARCZ), we obtained poor accuracy (67.9 ± 6.0 %), however, this is probably due to the fact that we used a binary classifier for a multi-class database.
We would like to do more experiments with this database in a near future splitting the data into healthy speech and mild, moderate and severe dysarthria.

4)
Select a database that has high potential for the future, especially in terms of speech features design, tuning and testing. Due to some issues related to vowels' length, different sampling frequencies, different recording conditions, etc., the MEEI database should no longer be used as a benchmark. Results obtained using this database are not trustable. The more challenging one is PARCZ database. Unfortunately it is focused only on parkinsonic people with hypokinetic dysarthria and it is not really suitable for a binary classification (at least 4 classes should be considered).
Therefore the most suitable database for the evaluation of pathological speech identification methods is PdA, where the classification accuracies are still challenging and they can be significantly improved. Moreover this database is freely available for research purposes.
There are many works that deal with the development of pathological voice identification methods and there is still a lot that can be improved in this field of science. However, the researchers should go further and focus not only on the pathological speech identification, but also on more sophisticated analysis that would be more helpful for doctors and that can make the treatment or diagnosis more effective. Probably the most important challenges to face in the next decade are: 1. Identification of particular voice pathology. Identification of pathological speech itself is not so interest-ing. There are many pathologies (adductor spasmodic dysphonia, erythema, hypokinetic dysarthria, etc.) and the issue is to classify them individually. This can be very problematic, therefore we propose to do some kind of clustering and split this big set into subsets, for example according to the way they are reflected in speech (problems with tongue movement, improper work of soft palate, disordered vocal folds, etc.).
2. Identification of voice pathology in its first stage or estimation of its progress. This would enable doctors to start the treatment very early and slow down the progress.
We are going to deal with these issues in future works. However this research is very dependent on good databases and especially in the case of voice pathology identification, as in its first stage, there is still a lack of suitable training data.