Exploiting Nonlinear Recurrence and Fractal Scaling Properties for Voice Disorder Detection
Max A Little, Patrick E McSharry, Stephen J Roberts, Declan AE Costello, Irene M Moroz
aa r X i v : . [ n li n . C G ] J u l Exploiting Nonlinear Recurrence and Fractal Scaling Prop-erties for Voice Disorder Detection
Max A Little ∗ , , , Patrick E McSharry , Stephen J Roberts , Declan AE Costello and Irene MMoroz Systems Analysis, Modelling and Prediction Group,Department of Engineering Science, University of Oxford, Parks Road, OxfordOX1 3PJ, UK Pattern Analysis Research Group,Department of Engineering Science, University of Oxford, Parks Road, Oxford OX1 3PJ, UK Applied Dynamical Systems Research Group,Oxford Centre for Industrial and Applied Mathematics, MathematicsInstitute,University of Oxford, Oxford OX1 3JP, UK Milton Keynes General Hospital, Standing Way, Eaglestone, Milton Keynes, Bucks MK6 5LD, UKEmail: Max A Little ∗ - [email protected]; Patrick E McSharry - [email protected]; Stephen J Roberts - [email protected];Declan AE Costello - [email protected]; Irene M Moroz - [email protected]; ∗ Corresponding author
Abstract
Background:
Voice disorders affect patients profoundly, and acoustic tools can potentially measure voicefunction objectively. Disordered sustained vowels exhibit wide-ranging phenomena, from nearly periodic to highlycomplex, aperiodic vibrations, and increased “breathiness”. Modelling and surrogate data studies have shownsignificant nonlinear and non-Gaussian random properties in these sounds. Nonetheless, existing tools are limitedto analysing voices displaying near periodicity, and do not account for this inherent biophysical nonlinearity andnon-Gaussian randomness, often using linear signal processing methods insensitive to these properties. They donot directly measure the two main biophysical symptoms of disorder: complex nonlinear aperiodicity, andturbulent, aeroacoustic, non-Gaussian randomness. Often these tools cannot be applied to more severedisordered voices, limiting their clinical usefulness.
Methods:
This paper introduces two new tools to speech analysis: recurrence and fractal scaling, whichovercome the range limitations of existing tools by addressing directly these two symptoms of disorder, togetherreproducing a “hoarseness” diagram. A simple bootstrapped classifier then uses these two features to distinguishnormal from disordered voices. esults: On a large database of subjects with a wide variety of voice disorders, these new techniques candistinguish normal from disordered cases, using quadratic discriminant analysis, to overall correct classificationperformance of . ± . . The true positive classification performance is . ± . , and the true negativeperformance is . ± . ( confidence). This is shown to outperform all combinations of the mostpopular classical tools. Conclusions:
Given the very large number of arbitrary parameters and computational complexity of existingtechniques, these new techniques are far simpler and yet achieve clinically useful classification performance usingonly a basic classification technique. They do so by exploiting the inherent nonlinearity and turbulentrandomness in disordered voice signals. They are widely applicable to the whole range of disordered voicephenomena by design. These new measures could therefore be used for a variety of practical clinical purposes.
Background
Voice disorders arise due to physiological disease or psychological disorder, accident, misuse of the voice, orsurgery affecting the vocal folds and have a profound impact on the lives of patients. This effect is evenmore extreme when the patients are professional voice users, such as singers, actors, radio and televisionpresenters, for example. Commonly used by speech clinicians, such as surgeons and speech therapists, areacoustic tools, recording changes in acoustic pressure at the lips or inside the vocal tract. These tools [1],amongst others, can provide potentially objective measures of voice function. Although acousticexamination is only one tool in the complete assessment of voice function, such objective measurement hasmany practical uses in clinical settings, augmenting the subjective judgement of voice function byclinicians. These measures find uses, for example, in the evaluation of surgical procedures, therapy,differential diagnosis and screening [1, 2], and often augment subjective voice quality measurements, forexample the GRB (Grade, Roughness and Breathiness) scale. [3] These objective measures can be used toportray a “hoarseness” diagram for clinical applications [4], and there also exists a variety of techniques forautomatically screening for voice disorders using these measures [5–7].Phenomenologically, normal and disordered sustained vowel speech sounds exhibit a large range ofbehaviour. This includes nearly periodic or regular vibration, aperiodic or irregular vibration and sounds2ith no apparent vibration at all. All can be accompanied by varying degrees of noise which can bedescribed as “breathiness”. Voice disorders therefore commonly exhibit two characteristic phenomena:increased vibrational aperiodicity and increased breathiness compared to normal voices [4].In order to better characterise the vibrational aperiodicity aspects of voice disorders, Titze [8] introduced atypology (extended with subtypes [1]). Type I sounds are those that are nearly periodic: coming close toperfect periodicity. Type II are those that show qualitative dynamical changes and/or modulatingfrequencies or subharmonics. The third class, Type III are those sounds that appear to have no periodicpattern at all. They have no single, obvious or dominant period and can described as aperiodic . Normalvoices can usually be classed as Type I and sometimes Type II, whereas voice disorders commonly lead toall three types of sounds. This is because voice disorders often cause the breakdown of stable periodicity invoice production. The breathiness aspect of disordered voices is often described as dominatinghigh-frequency noise. Although this original typology covered sounds of only apparently deterministicorigin, a very large proportion of voice signals seen in the clinic are so noisy as to be better consideredstochastic rather than deterministic [2]. Methods that are based upon theories of purely deterministicnonlinear dynamical systems, although they can be appropriate for sounds of deterministic origin coveredby the original typology, cannot in principle be applied to such noise-dominated sounds, that is, to soundsthat would be better modelled as stochastic processes rather than deterministic. This makes it impossibleto characterise the full range of signals encountered in the clinic. For these reasons, in this paper, when werefer to Type III sounds we include random, noise-like sounds (which, in keeping with original typology,have no apparent periodic structure, by virtue of their randomness).There exists a very large number of approaches to the acoustic measurement of voice function. The mostpopular of these are the perturbation measures jitter and shimmer and variants, and harmonics-to-noiseratios (HNR) [1, 4]. The effect of a variety of voice disorders on these measures has been tested under bothexperimental and clinical conditions [4, 9], showing that different measures respond to different disorders indifferent ways [10]. For example, in some disorders, jitter will increase with severity of the disorder, and forother disorders jitter is unaffected. Thus, although these measures can have value in measuring certainlimited aspects of voice disorders such as speech pitch variability, there is no simple relationship betweenthe extent or severity of voice disorder and these measures [4, 10]. Therefore, they cannot be used todirectly quantify the two main biophysical symptoms of voice disorders: increasingly severe aperiodicityand breath noise, a quantification required to differentiate normal from disordered voices.Another limitation of existing measures is that they are only properly applicable when near periodicity3olds: in Titze’s typology only Type I sounds satisfy this property [1]. The behaviour of the algorithms forother sound types is not known theoretically and limited only to experimental results and informalarguments [4]. The source of this limitation is that they make extensive use of extraction of the pitchperiod , or fundamental frequency (defined as the inverse of the pitch period) from the acoustic signal [1].Popular pitch period extraction techniques include zero-crossing detection, peak-picking and waveformmatching [1]. The concept of pitch period is only valid for Type I sounds, therefore, application of thesemethods based upon periodicity analysis, to any other type of sound is problematic [6]. Type II and IIIhave therefore received much less attention in the literature [8], such that there exist few methods forcharacterising these types, despite the fact that they exist in great abundance in clinical settings. Thisprecludes the proper use of these tools on a large number of disordered voice cases, limiting the reliabilityof the analysis, since in fact some algorithms will not produce any results at all for Type II and III sounds,even though they contain valuable clinical information [2]. Another reason for the limitations of thesemethods is that they are based upon classical linear signal processing methods that are insensitive to theinherent biophysical nonlinearity and non-Gaussianity in speech [1]. These linear methods includeautocorrelation, the discrete Fourier transform, linear prediction analysis and cepstral processing [11].Since standardised, reliable and reproducible results from acoustic measures of voice function are requiredfor clinical applications, these limitations of perturbation methods are problematic in clinical practice [2].It is clear that there is a clinical need for reliable tools that can characterise all types of disordered voicesounds for a variety of clinical applications, regardless of whether they satisfy the requirements of nearperiodicity, or contain significant nonlinearity, randomness or non-Gaussianity [8].Furthermore, analysis techniques are complicated by the use of any arbitrary algorithmic parameters whosechoice affects the analysis method – changing these parameters can change the analysis results. The valuesof such parameters must be chosen in order to apply the algorithms, and to be principled, it is better,when making that choice, to have a theoretical framework to apply which offers specific guidance on thatchoice. Often however, no such general theory exists, and therefore these values must be chosen byexperimental and empirical evaluation alone [12]. There exists the danger then that these parameters arehighly “tuned” to the particular data set used in any one study, limiting the reproducibility of the analysison different data sets or clinical settings. It is necessary therefore to reduce the number arbitraryparameters to improve the reproducibility of these measures.To address these limitations, empirical investigations and theoretical modelling studies have beenconducted which have lead to the suggestion that nonlinear dynamical systems theory is a candidate for a4nified mathematical framework modelling the dynamics seen in all types of disordered vowel speechsounds [1, 8, 13]. The motivation for the introduction of a more general model than the classical linearmodel is the principle of parsimony: the more general model explains more phenomena (more types ofspeech) with a smaller number of assumptions than the classical linear model [12, 14].These suggestions have lead to growing interest in applying tools from nonlinear time series analysis tospeech signals in order to attempt to characterise and exploit these nonlinear phenomena [1, 13]. Fornormal Type I speech sounds, fractal dimensions, Lyapunov exponents and bispectral methods have allbeen applied, also finding evidence to support the existence of nonlinearity [15, 16]. Extracting dynamicalstructure using local linear predictors, neural networks and regularised radial basis functions have all beenused, with varying reported degrees of success. Algorithms for calculating the correlation dimension havebeen applied, which was successful in separating normal from disordered subjects [17]. Correlationdimension and second-order dynamical entropy measures showed statistically significant changes before andafter surgical intervention for vocal fold polyps [18], and Lyapunov exponents for disordered voices werefound to consistently higher than those for healthy voices [19]. It was also found that jitter and shimmermeasurements were less reliable than correlation dimension analysis on Type I and unable to characteriseType II and (non-random) Type III sounds [20]. Mixed results were found for fractal dimension analysis ofsounds from patients with neurological disorders, for both acoustic and electroglottographic signals [21].Instantaneous nonlinear AM and FM formant modulations were shown effective at detecting muscletension dysphonias [22]. For the automated acoustic screening of voice disorders, higher-order statisticslead to improved normal/disordered classification performance when combined with several standardperturbation methods [7].In order to categorise individual voice signals into classes from the original typology (excluding severelyturbulent, noise-like sounds), recent studies have found that by applying correlation dimensionmeasurements to signals of these types, it was possible, over a sample of 122 stable vowel phonations, todetect a statistically significant difference between the three different classes [23, 24]. This provides furtherevidence in favour of the appropriateness of nonlinear signal analysis techniques for the analysis ofdisordered voices.These studies show that nonlinear time series methods can be valuable tools for the analysis of voicedisorders, in that they can analyse a much broader range of speech sounds than perturbation measures,and in some cases are found to be more reliable under conditions of high noise. However, very noisy,breathy signals have so far received little attention from nonlinear time series analysis approaches, despite5hese promising successes. Common approaches such as correlation dimension, Lyapunov exponentcalculation and predictor construction require that the scaling properties of the embedded attractor are notdestroyed by noise, and for thus very noisy, breathy signals, there is the possibility that such nonlinearapproaches may fail. There are also numerical, theoretical and algorithmic problems associated with thecalculation of nonlinear measures such as Lyapunov exponents or correlation dimensions for real speechsignals, casting doubt over the reliability of such tools [21, 25–27]. For example, correlation dimensionanalysis shows high sensitivity to the variance of signals in general, and it is therefore necessary to checkthat changes in correlation dimension are not due simply to changes in variance [28]. Therefore, as withclassical perturbation methods, current nonlinear approaches cannot yet directly measure the two mostimportant biophysical aspects of voice disorder.A limitation of deterministic nonlinear time series analysis techniques for random, very noisy signals is thatthe implicit, deterministic, nonlinear dynamical model, which is ordinarily assumed to represent thenonlinear oscillations of the vocal folds [29] is no longer appropriate. This is because randomness is aninherent part of the biophysics of speech production [30, 31]. Thus there is a need to expand the nonlineardynamical systems framework to include a random component, such that random voice signals and breathnoise can also be fully characterised within the same, unified framework.This paper therefore introduces a new, framework model of speech production that splits the dynamics intoboth deterministic nonlinear and stochastic components. The output of this model can then be analysedusing methods that are able to characterise both nonlinearity and randomness. The deterministiccomponent of the model can exhibit both periodic and aperiodic dynamics. It is proposed to characterisethis component using recurrence analysis . The stochastic components can exhibit statistical timedependence or autocorrelation, which can be analysed effectively using fractal scaling analysis . This paperreports the replication of the “hoarseness” diagram [4] illustrating the extent of voice disorder, anddemonstrates, using a simple quadratic classifier, how these new measures may be used to screen normalfrom disordered voices from a large, widely-used database of patients. It also demonstrates that these newmeasures achieve superior classification performance overall when compared to existing, classicalperturbation measures, and the derived irregularity and noise measures of Michaelis [4].
Methods
In this section we first discuss in detail the evidence that supports the development of a newstochastic/deterministic model of speech production. 6he classical linear theory of speech production brings together the well-developed subjects of classicallinear signal processing and linear acoustics (where any nonlinearities in the constitutive relationshipsbetween dynamic variables of the fluid are considered small enough to be negligible) to process and analysespeech time series [11]. The biophysical, acoustic assumption is that the vocal tract can be modelled as alinear resonator driven by the vibrating vocal folds that are the source of excitation energy [11]. However,extensive modelling studies, experimental investigations and analysis of voice signals have shown thatdynamical nonlinearity and randomness are factors that should not be ignored in modelling speechproduction.For the vocal folds, there are two basic, relevant model components to consider. The first is the vocal foldtissue, and the second is the air flowing through that structure. The vocal folds, during phonation, act as avibrating valve, disrupting the constant airflow coming from the lungs and forming it into regular puffs ofair. In general the governing equations are those of fluid dynamics coupled with the elastodynamics of adeformable solid. In one approach to solving the problem, the airflow is modelled as a modifiedquasi-one-dimensional Euler system which is coupled to the vocal fold flow rate, and the vocal folds aremodelled by a lumped two mass system [32]. A widely used and simplified model is the two-mass model inIshizaka [33], further simplified in asymmetric [29, 34] and symmetric versions [35]. These modelsdemonstrate a plausible physical mechanism for phonation as nonlinear oscillation: dynamical forcing fromthe lungs supplies the energy needed to overcome dissipation in the vocal fold tissue and vocal tract air.The vocal folds themselves are modelled as elastic tissue with nonlinear stress-strain relationship.Classical nonlinear vocal fold models coupled with linear acoustic vocal tract resonators appear to accountfor the major part of the mechanisms of audible speech, but from these mechanisms an importantcomponent is missing: that of aspiration , or “breath” noise [36]. Such noise is produced when the air isforced through a narrow constriction at sufficiently high speeds that “turbulent” airflow is generated,which in turn produces noise-like pressure fluctuations [37]. Aspiration noise is an unavoidable, involuntaryconsequence of airflow from the lungs being forced through the vocal organs, and can be clearly heard invowel phonation [38]. Certain voice pathologies are accompanied by a significant increase in suchaspiration noise, which is perceived as increased “breathiness” in speech [4]. This noise is therefore animportant part of sound generation in speech.One significant deficiency in the classical linear acoustic models is due to the assumptions about fluid flowupon which their construction is based [39]. These classical linear models make very many simplifyingassumptions about the airflow in the vocal organs, for example, that the acoustic limit [40] holds in which7he fluid is nearly in a state of uniform motion. Similarly, the
Bernoulli’s equation assuming energyconservation along streamlines, upon which many classical vocal fold models rely, applies only if the fluid isassumed inviscid and irrotational [41, 42]. The important point for this study is that these classicalassumptions forbid the development of complicated, “turbulent” fluid flow motion, in which the flowfollows convoluted paths of rapidly varying velocity, with eddies and other irregularities at all spatialscales [43]. This breakdown of laminar flow occurs at high
Reynolds number , and for the relevant lengthscales of a few centimetres in the vocal tract and for subsonic air flow speeds typical of speech [38], thisnumber is very large (of order 10 ), indicating that airflow in the vocal tract can be expected to beturbulent. Under certain assumptions, turbulent structures, and vortices in particular (fluid particles thathave rotational motion), can be shown to be a source of aeroacoustic sound [37]. Turbulence is a highlycomplex dynamical phenomenon and any point measurement such as acoustic pressure taken in the fluidwill lose most of the information about the dynamics of the fluid. Consequently, even if turbulence has apurely deterministic origin, it is reasonable to model any one single dynamical variable measured at a pointin space as a random process [43].There are two broad classes of physically plausible mathematical models of the effects of aeroacoustic noisegeneration in speech. The first involves solving numerically the full partial differential equations of gasdynamics (e.g. the Navier-Stokes equations), and the second uses the theory of vortex sound [37]. Forexample, the study of Zhao [44] focused on the production of aspiration noise generated by vortex sheddingat the top of the vocal folds, simulated over a full vocal fold cycle. This study demonstrates that thecomputed sound radiation due to vorticity contains significant high frequency fluctuations when the foldsare fully open and beginning to close. On the basis of these results, it can be expected that if the folds donot close completely during a cycle (which is observed in cases of more “breathy” speech), the amplitude ofhigh frequency noise will increase.The second class of models, which makes use of Lighthill’s acoustic analogy [37], are based around thetheory of vortex sound generated in a cylindrical duct, [37] where, essentially, each vortex shed at anupstream constriction acts as a source term for the acoustic wave equation in the duct, as the vortex isconvected along with the steady part of the airflow. The resulting source term depends upon not only theattributes of the vortex itself, but also upon the motion of the vortex through the streamlines of theflow [31, 37]. One model that uses this approach involves the numerical simulation of two overallcomponents: the mean steady flow field and the acoustic wave propagation in the vocal tract [38]. Vorticesare assumed to be shed at random intervals at constrictions at particular locations in the vocal tract, for8xample, at the vocal folds or between the roof of the mouth and the tongue. Each vortex contributes tothe acoustic source term at each spatial grid point. Here, an important observation is that simulatedpressure signals from this model are stochastic processes [45], i.e. a sequence of random variables. It is alsonoticeable from the spectra of simulated pressure signals that although the signals are stochastic, theyexhibit significant non-zero autocorrelation since the spectral magnitudes are not entirely constant, leadingto statistical self-similarity in these signals [15, 43].Other potential sources of biophysical fluctuation include pulsatile blood flow, muscle twitch and othervariability and randomness in the neurological and biomechanical systems of the larynx [30].The typical nonlinear dynamical behaviour of models of the vocal folds, such as period-doubling(subharmonics), bifurcations [29], and transitions to irregular vibration [35] have been observed inexperiments with excised larynxes, a finding that helps to support the modelling hypothesis that speech isan inherently nonlinear phenomena [29, 46]. Similarly, models of turbulent, random aeroacoustic aspirationnoise have been validated against real turbulent airflow induced sound generated in acoustic ductexperiments [38]. Such studies show that the models of vortex shedding at random intervals are plausibleaccounts for the dynamical origins of breath noise in phonation.Complementing modelling and experimental studies, the final source of evidence for nonlinearity andrandomness in speech signals comes from studies of voice pressure signals. Using surrogate data analysis, ithas been shown that nonlinearity and/or non-Gaussianity is an important feature of Type I sounds [47–49].Nonlinear bifurcations have been observed in excised larynx experiments [46], and period-doublingbifurcations and chaos-like features have been observed in signals from patients with organic and functionaldysphonias [13]. Aspiration noise has been observed and measured as a source of randomness in voicedphonation, in both normal [50] and disordered speech sounds [1, 4]. The fractal, self-similarity properties ofaspiration noise as a turbulent sound source have also been observed and quantified in normal [15] anddisordered speech [27].Taken as a whole, these modelling, simulation, validation and signal analysis studies suggest that duringvoiced phonation there will be a combination of both deterministic and stochastic elements, thedeterministic component attributable to the nonlinear movements of the vocal fold tissue and bulk of theair in the vocal tract, and the stochastic component the high frequency aeroacoustic pressure fluctuationscaused by vortex shedding at the top of the vocal folds, whose frequency and intensity is modulated by thebulk air movement (and other sources of biophysical fluctuation). During voiced phonation, thedeterministic oscillations will dominate in amplitude over the noise component which will show high9requency fluctuations around this oscillation. During unvoiced or breathy pathological phonation, theturbulent noise component will dominate over any deterministic motion.In order to capture all these effects in one unified model, we introduce a continuous time, two componentdynamical model that is taken to generate the measured speech signal. The state of the system at time t ∈ R is represented by the vector u ( t ) of size d . Then the equation of motion that generates the speechsignal is the following vector stochastic differential equation, commonly known as a Langevin equation [25]. ˙u ( t ) = f ( u ( t )) + ε ( t ) , (1)where ε ( t ) is a vector of stochastic forcing terms. It is not necessary to assume that these fluctuations areindependent and identically distributed (i.i.d.). The function f : R d → R d is unknown and represents thedeterministic part of the dynamics. Given an initial condition vector u (0) then a solution that satisfiesequation (1) is called a trajectory . Ensembles of trajectories can be shown to satisfy the properties of astochastic process with finite memory (a higher-order Markov chain) [25].Of importance to both deterministic and stochastic systems is the notion of recurrence in state space.Recurrent trajectories are those that return to a given region of state space [51]. Recurrence time statisticsprovide useful information about the properties of both purely deterministic dynamical systems andstochastic systems [51]. Recurrence analysis forms the basis of the method of recurrence plots in nonlineartime series analysis [25].In the context of the model (1), state-space recurrence is defined as: u ( t ) ⊂ B ( u ( t + δt ) , r ) , (2)where B ( u , r ) is a closed ball of small radius r > u in state-space, and u ( t ) B ( u ( t + s ) , r ) for 0 < s < δt . Each different t ∈ R may be associated with a different δt , called the recurrence time . We will define aperiodic as recurrent in the sense of (2) but not periodic . Periodicity isrecovered from the definition of recurrence in the special case when r = 0 and δt is the same for all t , sothat the system vector u ( t ) assumes the same value after a time interval of δt : u ( t ) = u ( t + δt ) , (3)for all t ∈ R . Then δt is the period of the system. Therefore, although these concepts of periodicity andaperiodicity are mutually exclusive, both are special cases of the more general concept of recurrence. Therequirement of periodicity is central to many linear signal processing methods (the basis of the Fourier10ransform, for example), but aperiodicity is a common feature of many voice disorders. It can be seentherefore that in order to characterise the inherent irregularity of disordered speech sounds, we requiremore general processing techniques that can directly account for such departures from strict periodicity.Using this analysis, nearly periodic speech sounds of Type I can be described as recurrent for some small r >
0, with δt nearly the same for each t . Type II sounds are more irregular than Type I, and for the same r , the δt will assume a wider range of values than for Type I. Similarly, Type III sounds which are highlyirregular and aperiodic, will have a large range of values of δt again for the same r .Similarly, of importance in the analysis of stochastic processes is scaling analysis [25]. Stochastic timeseries in which the individual measurements in time are not entirely independent of earlier time instantsare often self-similar , that is, when a sub-section of the time series is scaled up by a certain factor, it hasgeometrical, approximate or statistical similarity to the whole time series [25]. Such self-similarity is aproperty of fractals . [43] The tools of dimension measurement and scaling analysis may be used tocharacterise the self-similarity in signals such as speech. As discussed above, theoretical models ofaeroacoustic turbulent sound generation in speech predict randomly-timed impulses convolved with anindividual impulse response for each vortex that induces autocorrelated random noise sequences [31], sothat turbulent speech signals at one time instant are not independent of earlier time instants.It has also been found experimentally that changes in the statistical time dependence properties ofturbulent noise in speech, as measured by a particular fractal dimension of the graph of the speech signal,are capable of distinguishing classes of phonemes from each other [15]. As introduced above, disorderedvoices are often accompanied by increased “breathiness”, due in part to the inability of the vocal folds toclose properly, so that air escapes through the partial constriction of the vocal folds creating increasedturbulence in the airflow [31]. Thus scaling analysis (and more general graph dimension measures) could beuseful for characterising vocal fold disorders.Recent experiments have shown that the use of recurrence analysis coupled with scaling analysis candistinguish healthy from disordered voices on a large database of recordings with high accuracy [27]. Thesetechniques are computationally simple and furthermore substantially reduce the number of arbitraryalgorithmic parameters required, compared to existing classical measures, thus leading to increasedreproducibility and reliability. 11 ime-Delay State-Space Recurrence Analysis In order to devise a practical method for applying the concept of recurrence defined earlier and measuringthe extent of aperiodicity in a speech signal, we can make use of time-delay embedding [25]. Measurementsof the output of the system (1) are assumed to constitute the acoustic signal, s n : s n = h ( u ( n ∆ t )) , (4)where the measurement function h : R d → R projects the state vector u ( t ) onto the discrete-time signal attime instances n ∆ t where ∆ t is the sampling time, and n ∈ Z is the time index.From these sampled measurements, we then construct of an m -dimensional time delay embedding vector : s n = (cid:2) s n , s n − τ , . . . s n − ( m − τ (cid:3) T . (5)Here τ is the embedding time delay and m is the embedding dimension .We will make use of the approach introduced in Ragwitz [52] to optimise the embedding parameters m and τ such that the recurrence analysis produces results that are close as possible to known, analyticallyderived results upon calibration with known signals. (We use this as an alternative to common techniquesfor finding embedding parameters, such as false-nearest neighbours, which explicitly require purelydeterministic dynamics [25, 52]). Note that, under this approach, for very noisy signals, we will not alwaysresolve all signals without self-intersections. However, in the context of this study, achieving a completelynon-intersecting embedding is not necessary. For very high-dimensional deterministic or stochastic systems,any reconstruction with self-intersections due to insufficiently high embedding dimension can be consideredas a different stochastic system in the reconstructed state-space. We can then analyze the stochasticrecurrence of this reconstructed system. This recurrence, albeit different from the recurrence properties inthe original system, is often sufficient to characterize the noisy end of the scale of periodicity andaperiodicity.Figure 1 shows the signals s n for one normal and one disordered voice example (Kay Elemetrics DisorderedVoice Database). Figure 2 shows the result of applying the above embedding procedure for the samespeech signals.In order to investigate practically the recurrence time statistics of the speech signal, we can make use ofthe method of close returns [53], originally designed for studying recurrence in deterministic, chaoticdynamics. Here we adopt this method to study stochastic dynamics as well as deterministic dynamics. Inthis method, a small, closed ball B ( s n , r ) of radius r > s n .12hen the trajectory is followed forward in time s n +1 , s n +2 . . . until it has left this ball, i.e. until k s n − s n + j k > r for some j >
0, where k . k is the Euclidean distance. Subsequently, the time n at whichthe trajectory first returns to this same ball is recorded (i.e. the first time when k s n − s n k ≤ r ), and thedifference of these two times is the (discrete) recurrence time T = n − n . This procedure is repeated forall the embedded data points s n , forming a histogram of recurrence times R ( T ). This histogram isnormalised to give the recurrence time probability density : P ( T ) = R ( T ) P T max i =1 R ( i ) , (6)where T max is the maximum recurrence time found in the embedded state space. The choice of r is criticalto capture the properties of interest to this study. For example, if the trajectory is a nearly periodic (TypeI), we require that r is large enough to capture all the recurrences, but not too large to find recurrencesthat are due to spurious intersections of B ( u , r ) with other parts of the trajectory, violating the conditionsfor proper recurrence. The appropriate choice of embedding delay τ has a role to play: selecting τ toosmall means that any trajectory lies close to a diagonal in the reconstructed state-space, potentiallycausing spurious recurrences. Thus τ must be chosen large enough to avoid spurious recurrences. Similarly,if τ is too large then the points in the embedded state-space fill a large cloud where recurrences will bedifficult to find without using an inappropriately large value of r . There will be an optimum value of τ which traditionally is set with reference to autocorrelation or time-delayed mutual information estimates,for more details see [25].In order to understand the behaviour of this algorithm (partly for optimising the embedding parameters),we consider two extreme forms that the density (6) may assume. The first is the ideal limiting case in whichthe recurrence distance r tends to zero for a periodic trajectory. The recurrence time probability density is: P ( T ) = (cid:26) T = k , (7)where k is the period of the trajectory. In the second extreme case we consider a purely random, uniformi.i.d. signal which is normalised to the range [ − , P ( T ) ≈ T max . (8)Proofs for the results in equations (7) and (8) are given in the Appendix.We can then optimise m , τ and r such that for a synthetic signal of perfect periodicity, P ( T ) is determinedusing the close returns method such that it is as close as possible to the theoretical expression (7). This13ptimisation is carried out by a straightforward systematic search of values of these parameters m = 2 , . . . τ = 2 , . . .
50, and for r = 0 . , . , . . . .
5, on a perfectly periodic test signal. This searchcan be considered as a scan for the optimum parameter values through all points on a three-dimensionalparameter cube with m , τ and r as the co-ordinates of that cube.All disordered voice speech signals will lie somewhere in between the extremes of perfect periodicity andcomplete randomness. Thus it will be useful to create, from the recurrence time probability density, astraightforward sliding scale so that normal and disordered voice signals can be ranked alongside eachother. This depends upon a simple characterisation of the recurrence probability density P ( T ). One simplemeasure of any probability density is the entropy which measures the average uncertainty in the value ofthe discrete-valued density p ( i ), i = 1 , . . . M [25]: H = − M X i =1 p ( i ) ln p ( i ) , (9)with units of nats (by convention, 0 ln 0 is taken to be zero). This measure can then rank disordered voicesignals by representing the uncertainty in the period of the disordered voice signal in the following way. Forperfectly periodic signals the recurrence probability density entropy (RPDE) is: H per = − T max X i =1 P ( i ) ln P ( i ) = 0 . (10)since P ( k ) = 1 and the rest are zero. Conversely, for the purely stochastic, uniform i.i.d. case (derived inthe appendix), the uniform density can be taken as a good approximation, so that the RPDE is: H iid = − T max X i =1 P ( i ) ln P ( i ) = ln T max , (11)in units of nats. The entropy scale H therefore ranges from H per , representing perfectly periodic examplesof Type I sounds, to H iid as the most extreme cases of noise-like Type III sounds. In practice, all soundswill lie somewhere in between these extremes.However, the entropy of a probability density is maximum for the uniform density, so that H iid is themaximum value that H can assume. Thus, in addition to ranking signals on a scale of aperiodicity, we canknow precisely the two extremes of that scale. For different sampling times ∆ t the value T max will change.Therefore, we can normalised the RPDE scale for subsequent usage: H norm = − P T max i =1 P ( i ) ln P ( i ) H iid , (12)a unit less quantity that assumes real values in the range [0 , H will also be an invariant of the system. We note thattraditional jitter and shimmer measurements do not share this invariance property, in this sense they donot give a reliable description for chaotic or Type III time series. Often, when stable phonation is initiatedin speech, the vocal folds will take some time to settle into a pattern of stable periodic or chaoticoscillation. The behaviour of speech signals during this “settling time” is similar to the transient behaviourtypically observed in nonlinear dynamical systems. In this study, we make use of voice signals which arestable phonations, and we discard any of these transient phases. Thus, to a reasonable approximation H can be considered as an invariant of the dynamics of the speech organs.Figure 3 shows the result of this recurrence analysis, applied to a synthetic, perfectly periodic signalcreated by taking a single cycle from a speech signal and repeating it end-to-end many times. It also showsthe analysis applied to a synthesised, uniform, i.i.d. random signal (on the signal range [ − , m , τ and r by gridded search. Even though exact results are impossible to obtain due to theapproximation inherent to the algorithm and only finite-length signals, the figure shows that a close matchis obtainable between the theoretical, predicted results and the simulated results. Detrended Fluctuation Analysis
In order to investigate the second aspect of disordered voices, that of increased breath noise, we require apractical approach to applying the scaling ideas introduced above. Detrended fluctuation analysis is onestraightforward technique for characterising the self-similarity of the graph of a signal from a stochasticprocess [54].It is designed to calculate the scaling exponent α in nonstationary time series (where the statistics such asmean, variance and autocorrelation properties might change with time), and has been shown to producerobust results when there are slowly moving trends in the data. These will naturally include low frequencyvibrations [54], which are similar in nature to the nonlinear vibrations of the vocal folds described by thefunction f in the model (1). Thus this technique can be used to characterise the properties of only thestochastic part ε ( t ) of the model (1). 15n this technique, the scaling exponent α is a measure of the ratio of the logarithm of the fluctuations orvertical height of the graph to the logarithm of the horizontal width of a chosen time window over whichthat vertical height is measured. The scaling exponent is calculated as the slope of a log-log graph of arange of different horizontal time window sizes against, the vertical height of the signal in those timewindows. Mathematically, for self-similar signals with positive scaling exponent α the self-similarityproperty of the graph of the signal s n should hold on all time scales, but we are limited by the finiteamplitude range of physical measurements to a certain maximum scale. Thus the signal is first integratedin order to create a new stochastic process that exhibits self-similarity over a large range of time scales(then, for example, a purely independent, stochastic process will result in a self-similar random walk).First, the time series s n is integrated: y n = n X j =1 s j , (13)for n = 1 , . . . N , where N is the number of samples in the signal s n . Then, y n is divided into windows oflength L samples. A least-squares straight line local trend is calculated by analytically minimising thesquared error E over the slope and intercept parameters a and b :arg min a,b E = L X n =1 ( y n − an − b ) . (14)Next, the root-mean-square deviation from the trend, the fluctuation, is calculated over every window atevery time scale: F ( L ) = " L L X n =1 ( y n − an − b ) . (15)This process is repeated over the whole signal at a range of different window sizes L , and a log-log graph of L against F ( L ) is constructed. A straight line on this graph indicates self-similarity expressed as F ( L ) ∝ L α . The scaling exponent α is calculated as the slope of a straight line fit to the log- log graph of L against F ( L ) using least-squares as above. For a more in-depth presentation and discussion ofself-similarity in signals in general, and further information about DFA, please see Kantz, Hu [25, 54].We are assuming that the signal, as the measured output of the new model, represents a combination ofdeterministic and stochastic dynamics. The deterministic part of the dynamics, dictated by the function f in equation (1) will result in slower changes in the signal s n taking place over a relatively long time scale.Similarly, the stochastic fluctuations in the signal indicated changes taking place over a much shorter timescale. Since the goal of DFA is to analyse the faster changing, stochastic properties of the signal, only alimited range of window sizes is investigated, over which the stochastic component of the signal exhibits16elf-similarity as indicated by a straight-line on the log-log graph of window size against fluctuation. As anexample, Type III would include some speech signals that are actually chaotic, where the chaos is due toslow, nonlinear dynamics in the vocal fold tissue and airflow, the characteristic time of this nonlinearoscillation will be much longer than the window sizes over which the scaling exponent is estimated. Thus,the nature of the chaotic oscillation will not affect the scaling exponent, which will respond only to anyrandom fluctuations occurring on a much shorter time scale.The resulting scaling exponent can assume any number on the real line. However, it would be moreconvenient to represent this scaling exponent on a finite sliding scale from zero to one, as we have done forthe RPDE measure. Thus we need a mapping function g : R → [0 , logistic function g ( x ) = (1 + exp( − x )) − [55],so that the normalised scaling exponent is: α norm = 11 + exp ( − α ) . (16)Therefore, each sound will lie somewhere between the extremes of zero and one on this scale, according tothe self-similarity properties of the stochastic part of the dynamics. As will be shown later, speech soundsfor which α norm is closer to one are characteristic of general voice disorder. Application to Examples of Normal and Disordered Voices
In this section we will apply these two new measures to some examples of normal and disordered voices, asa limited demonstration of characterising the extent of aperiodicity and breathiness in these signals.Figure 4 shows the RPDE value H norm calculated on the same two speech signals from the Kay database asshown in figure 1. Note that the second, disordered example is of Type III and shows significantly irregularvibration, which is detected by a large increase in H norm .Similarly, figure 5 shows two more speech examples, one normal and one disordered from the same databaseand the corresponding values of the scaling exponent α and α norm . In these cases, the disordered exampleis extremely “breathy”, and this turbulent noise is detected by an increase in the scaling exponent. Quadratic Discriminant Analysis
In order to test the effectiveness of these two measures in practice, one approach is to set up a classificationtask to separate normal control subjects from disordered examples using these measures alone. Here wechoose one of the simplest approaches, quadratic discriminant analysis (QDA) , which allows separation by17odelling the data conditional upon each class, here the normal (class C ) and disordered (class C ) cases,using joint Gaussian probability density functions [55]. For a J × K data matrix v = v jk of observationsconsisting of the measures j = 1 , k , these likelihooddensities are parameterised by the mean and covariance matrices of the data set: µ = E [ v ] , C = E h ( v − µ ) ( v − µ ) T i , (17)where E is the expectation operator, and µ is the mean vector formed from the means of each row of v .The class likelihoods are: f C ( w | C i ) = (2 π ) − J/ | C i | − / exp (cid:20) −
12 ( w − µ i ) T C − i ( w − µ i ) (cid:21) , (18)for classes i = 1 , w . It can be shown that, given these Gaussian classmodels, the maximum likelihood regions of the observation space R J are separated by a decision boundary which is a (hyper-)conic section calculated from the difference of log-likelihoods for each class, which is theunique set of points where each class is equally likely [55]. The maximum likelihood classification problemis then solved using the decision rule that l ( w ) ≥ w to class C , and l ( w ) < C , where: l ( w ) = − w T A w + A w + A , (19) A = C − − C − , A = µ T C − − µ T C − , (20) A = −
12 ln | C | + 12 ln | C | − µ T C − µ + 12 µ T C − µ . (21)In order to avoid the problem of overfitting, where the particular chosen separation model shows goodperformance on the training data but fails to generalise well to new, unseen data, the classifier resultsrequire validation.In this paper, we make use of bootstrap resampling for validation [55]. In the bootstrap approach, theclassifier is trained on K cases selected at random with replacement from the original data set of K cases.This trial resampling processes is repeated many times and the mean classification parameters E [ A ] , E [ A ] , E [ A ] are selected as the parameters that would achieve the best performance on entirelynovel data sets.Bootstrap training of the classifier involves calculating H k norm and α k norm (the observations) for each speechsample k in the database, (where the superscript k denotes the measure for the k -th subject). Then, K random selections of these values with replacement H ′ k norm and α ′ k norm form the entries of the vector18 k = H ′ k norm and v k = α ′ k norm Then the mean vectors for each class µ and µ and covariance matrices C , C for the whole selection are calculated. Next, for each subject, the decision function is evaluated: l ( w k ) = l ([ H k norm , α k norm ] T ) . (22)Subsequently, applying the decision rule assigns the subject k into either normal or disordered classes.Then the performance of the classifier can be evaluated in terms of percentage of true positives (when adisordered subject is correctly assigned to the disordered class C ) and true negatives (when a normalsubject is correctly assigned to the normal class C ). The overall performance is the percentage of correctlyclassified subjects, in both classes. This bootstrap trial process of creating random selections of themeasures, calculating the class mean vectors and covariance matrices, and then evaluating the decisionfunction on all the measures to obtain the classification performance is repeated many times. Assumingthat the performance percentages are normally distributed, then the 95% confidence interval of theclassification performance percentages can be calculated. The best classification boundary can then betaken as the mean boundary overall all the trials.Efficient implementations of the algorithms described in this paper written in C with Matlab MEXinterface accompany this paper: close returns [see Additional files 1 and 2] and detrended fluctuationanalysis [see Additional files 3, 4 and 5]. Algorithms for Classical Techniques
For the purposes of comparison, we calculate the classical measures of jitter, shimmer and HNR(Noise-to-Harmonics Ratio) [1]. There are many available algorithms for calculating this quantity, in thisstudy we make use of the algorithms supplied in the software package Praat [56]. These measures are basedon an autocorrelation method for determining the pitch period (see Boersma [57] for a detailed descriptionof the method).We also use the methods described in Michaelis [4]. This first requires calculating the measures EPQ(Energy Perturbation Quotient), PPQ (Pitch Perturbation Quotient), GNE (Glottal to Noise ExcitationRatio) and the mean correlation coefficient between successive cycles, measures which require theestimation of the pitch period using the waveform matching algorithm (see Titze [58] for a detaileddescription of this algorithm). The EPQ, PPQ, GNE and correlation coefficients are calculated oversuccessive overlapping “frames” of the speech signal. Each frame starts at a multiple of 0.26 seconds, andis 0.5 seconds long. For each frame, the EPQ, PPQ, GNE and correlation coefficients are combined into a19air of component measures, called Irregularity and Noise. We use the average of the Irregularity andNoise components over all these frames [4].
Classification Test Data
This study makes use of the Kay Elemetrics disordered voice database (KayPENTAX Model 4337, NewJersey, USA), which contains 707 examples of disordered and normal voices from a wide variety of organic,neurological and traumatic voice disorders. This represents examples of all three types of disordered voicespeech signals (Types I, II and III). There are 53 control samples from normal subjects. Each speechsample in the database was recorded under controlled, quiet acoustic conditions, and is on average aroundtwo seconds long, 16 bit uncompressed PCM audio. Some of the speech samples were recorded at 50kHzand then downsampled with anti-aliasing to 25kHz. Used in this study are sustained vowel phonations,since this controls for any significant nonstationarity due changes in the position of the articulators such asthe tongue and lips in running speech, which would have an adverse effect upon the analysis methods. Forcalculating the Irregularity and Noise components, the signals are resampled with anti-aliasing to 44.1kHz.
Results
Figure 6 shows hoarseness diagrams after Michaelis [4] constructed using the speech data and RPDE andDFA measures, the derived irregularity and noise components of Michaelis, along with the same diagramsusing two other combinations of the three classical perturbation measures for direct comparison. The threeclassical measures are jitter, shimmer and NHR (Noise-to-Harmonics Ratio) [1]. The normalised RPDE,DFA scaling exponents and derived irregularity and noise components are calculated for each of the K = 707 speech signals. Where the traditional perturbation algorithms did not fail to produce a result, thetraditional perturbation values were calculated for a smaller subset of the subjects, see [1] for details ofthese traditional algorithms. Also shown in figure 6 is the result of the classification task applied to thedataset; the best classification boundary is calculated using bootstrap resampling over 1000 trials. Table 1summarises all the classification performance results for the classification tasks on the hoarseness diagramsof figure 6. The RPDE parameters were the same as for figure 3, and the DFA parameters were the sameas for figure 5. 20 iscussion As shown in table 1, of all the combinations of the new and traditional measures, and derived irregularityand noise components, the highest overall correct classification performance of 91 . ± .
0% is achieved bythe RPDE/DFA pair. The combination of jitter and shimmer leads to the next highest performance. Theseresults confirm that, compared under the same, simple classifier approach, the new nonlinear measures aremore accurate on average than traditional measures or the derived irregularity and noise components. Wewill now discuss particular aspects of these results in comparison with traditional perturbation measures.
Feature Dimensionality
The curse of dimensionality afflicts all challenging data analysis problems [55]. In pattern analysis taskssuch as automated normal/disordered separation, increasing the size of the feature vector (in this case, thenumber of measures J in the classifier vector v ) does not necessarily increase the performance of theclassifier in general. This is because the volume of the feature space (the space spanned by the possiblevalues of the measures) grows exponentially with the number of features. Therefore, the limited number ofexamples available to train the classifier occupy an increasingly small volume in the feature space,providing a poor representation of the mapping from features to classes that the classifier must learn [55].Therefore the new measures help to mitigate this problem of dimensionality, since only these two newmeasures are required to obtain good separation performance. By comparison, we need to calculate fourdifferent measures in order obtain the irregularity and noise components [4]. Feature Redundancy – Information Content
It is also important to use as few features as possible because in practice, increasing the number of featurescauses excessive data to be generated that may well contain redundant (overlapping) information. Theactual, useful information contained in these vectors has a much smaller dimensionality. For clinicalpurposes, it is important that only this useful data is presented. This effect of redundant information forthe traditional measures can be clearly seen in figure 6, where combinations of pairs of (the logarithms of)these measures are seen to cluster around a line or curve in the feature space, indicating high positivecorrelation between these measures. Traditional measures occupy an effectively one-dimensional object inthis two-dimensional space. The irregularity and noise components occupy more of the area of the featurespace than traditional measures, and the new measures are spread evenly over the same space.21 rbitrary Parameters – Reproducibility
Minimising the number of arbitrary parameters used to calculate these measures is necessary to avoidselecting an excessively specialised set of parameters that leads, for example, to good normal/disorderedseparation on a particular data set but does not generalise well to new data.Many parameters are required for the algorithms used in calculating traditional perturbationmeasures [4, 5, 7]. For example, the waveform matching algorithm [1] requires the definition of roughmarkers, upper and lower pitch period limits, low-pass filter cutoff frequencies, bandwidth and orderselection parameters, and the number of pitch periods for averaging should these pitch period limits beexceeded [58]. Similarly, in just one of the noise measures (glottal-to-noise excitation ratio) used inMichaelis [4], we can determine explicitly at least four parameters relating to linear prediction order,bandpass filter number, order, cutoff selection, and time lag range parameters. There are two additionalparameters for the length and starting sample of the part of the signal selected for analysis.Our new measures require only five arbitrary parameters that must be chosen in advance: the length of thespeech signal N , the maximum recurrence time T max , and the lower value, upper value and increment ofthe DFA interval lengths L . We have also shown, using analytical results, that we can calibrate out thedependence upon the state space close recurrence radius r , the time-delay reconstruction dimension d andthe reconstruction delay τ . Interpretation of Results
We have found, in agreement with Titze [8] and Carding [2], that perturbation measures cannot beobtained for all the speech sounds produced by subjects (see table 1). This limits the clinical usefulness ofthese traditional measures. By contrast, the new measures presented in this chapter do not suffer from thislimitation and are capable of measuring, by design, all types of speech signals.Taking into account the classification performance achievable using a simple classifier, the number of thesemeasures that need to be combined to achieve effective normal/disordered separation, the number ofarbitrary parameters used to calculate the measures, and the independence of these measures, traditionalapproaches and derived irregularity and noise components are seen to be considerably more complex thanthe new measures developed in this paper. The results of the classification comparison with traditionalmeasures and the irregularity and noise components suggest that, in order to reach the classificationperformance of the new measures, we will either need much more complex classifiers, or need to combinemany more classical features together [5–7]. It is therefore not clear that traditional approaches capture22he essential biomechanical differences between normal and disordered voices in the most parsimoniousway, and an excessively complicated relationship exists therefore between the values of these measures andextent of the voice disorder. As a final comment, we note that the classical perturbation measures were, forthe majority of signals, able to produce a result regardless of the type of the signal. This is consistent withthe findings of other studies [4], where for Type II/III and random noise signals, the correct interpretationof these measures breaks down. Therefore, although it is no longer possible in these cases to assign acoherent meaning to the results produced by these measures, this does not per se mean that there is notsome as yet unknown connection between disorder and the these measures. For this reason, we do notdiscard the results of these measures for Type II/III and random cases.
Limitations of the New Measures
There are certain limitations to the new measures in clinical practice. These measures rely upon sustainedvowel phonation, and sometimes subjects experience difficulty in producing such sounds, which limits theapplicability. Also, at the beginning of a sustained vowel phonation, the voice of many subjects mayrequire some time to settle into a more stable vibration. As such, discarding the beginning of thephonation is usually a prerequisite (but this does not adversely affect the applicability of these methods).Nonetheless, the extent of breathiness in speech is not usually affected in this way. In practice we requirethat the subject maintains a constant distance from the microphone when producing speech sounds; thiscan be achieved, for example, with the use of head-mounted microphones. We note that these limitationsalso apply to existing measures.
Possible Improvements and Extensions
There are several improvements that could be made to these new measures. Firstly, every arbitraryparameter introduces extra variability that affects the reliability of the results. Much as it has beenpossible to calibrate out the dependence upon the RPDE parameters using analytical results, a theoreticalstudy of the DFA interval lengths based upon typical sustained phonation recurrence periods could revealvalues that would be found for all possible speech signals. These would be related to the sampling time ∆ t .The particular choice of normalisation function g for the scaling exponent might affect the classificationperformance, and better knowledge of the possible range of α values using theoretical studies of the DFAalgorithm would be useful for this. It should also be possible to increase the recurrence time precision ofthe RPDE analysis by interpolating the state space orbits around the times of close recurrence n , n . It23hould then be possible to achieve the same high-resolution as waveform matching techniques which wouldmake RPDE competitive for the detailed analysis of Type I periodic sounds. Conclusions
In this study, in order to directly characterise the two main biophysical factors of disordered voices:increased nonlinear, complex aperiodicity and non-Gaussian, aeroacoustic breath noise, we have introducedrecurrence and scaling analysis methods. We introduced a new, combined nonlinear/stochastic signalmodel of speech production that is capable, in principle, of producing the wide variation in behaviour ofnormal and disordered voice examples. To exploit the output of this model in practice, and hence all typesof normal and disordered voices, we explored the use of two nonlinear measures: the recurrence perioddensity entropy and detrended fluctuation analysis.Our results show that, when the assumptions of the model hold under experimental conditions (in that thespeech examples are sustained vowels recorded under quiet acoustic conditions), we can directlycharacterise the two main factors of aperiodicity and breath noise in disordered voices and thus construct a“hoarseness” diagram showing the extent of normality/disorder from a speech signal. The results also showthat on average, over all bootstrap resampling trials, these two measures alone are capable ofdistinguishing normal subjects from subjects with all types of voice disorder, with better classificationperformance than existing measures.Furthermore, taking into account the number of arbitrary parameters in algorithms for calculating existingperturbation measures, and the number of these existing measures that need to be combined to performnormal/disordered separation, we have shown that existing approaches are considerably more complex.We conclude that the nonlinearity and non-Gaussianity of the biophysics of speech production can beexploited in the design of signal analysis methods and screening systems that are better able characterisethe wide variety of biophysical changes arising from voice disease and disorder. This is because ultimatelythe biophysics of speech production generate the widely varying phenomenology.
Appendix – Mathematical Proofs
Periodic Recurrence Probability Density
We consider the purely deterministic case, i.e. when the model of equation (1) has no forcing term ε ( t ).Thus the measured time series is purely deterministic and points in the time series follow each other in anexactly prescribed sequence. When the measured, trajectory s n is a purely periodic orbit of finite period k { r n } , n ∈ Z in the reconstructed state space with r n = r n + k ,and r n = r n + j for 0 < j < k .Picking any point s in the reconstructed state-space, there are two cases to consider. In the first case, if s = r n for some n , then s is not the same as any other points in the periodic orbit except for r n + k , so thatthe trajectory returns with certainty for the first time to this point after k time steps. This certainty, withthe requirement the that probability of first recurrence is normalised for T = 1 , . . . implies that: P s ( T = r ) = (cid:26) r = k . (23)In the second case when s = r n for any n , the trajectory never intersects the point so that there are alsonever any first returns to this point. All the points in the reconstructed space form a disjoint partition ofthe whole space. Thus the probability of recurrence to the whole space is the sum of the probability ofrecurrence to each point in the space separately, appropriately weighted to satisfy the requirement that theprobability of first recurrence to the whole space is normalised However, only the k distinct points of theperiodic orbit contribute to the total probability of first recurrence to the whole space. Therefore, theprobability of first recurrence is: P ( T ) = 1 k k − X i =0 P r i ( T = r ) = (cid:26) r = k . (24) Uniform i.i.d. Stochastic Recurrence Probability Density
Consider the purely stochastic case when the nonlinear term f is in equation (1) is zero and the stochasticforcing term is an i.i.d. random vector. Then the measured trajectory s n is also a stochastic, i.i.d. randomvector. Since all the time series are normalised to the range [ − ,
1] then each member of the measurementtakes on a value from this range. Then the trajectories s n occupy the reconstructed state-space which is theregion [ − , m , and each co-ordinate s n is i.i.d. uniform. We form a equal-sized partition of this space into N m (hyper)-cubes, denoting each cubical region R . The length of the side of each cube R is ∆ s = 2 /N .Then the probability of finding the trajectory in this cube is P R = ∆ s m / m . Since the co-ordinates s n areuniform i.i.d., then the probability of first recurrence of time T to this region R is geometric [51]: P R ( T ) = P R [1 − P R ] T − = ∆ s m m (cid:20) − ∆ s m m (cid:21) T − . (25)This is properly normalised for T = 1 , . . . . However, we require the probability of first recurrence to allpossible cubes. The cubes are a disjoint partition of the total reconstruction space [ − , m . Thus the25robability of recurrence to the whole space is the sum of the probability of recurrence to each cubeseparately, appropriately weighted to satisfy the requirement that the probability of recurrence to thewhole space is normalised. Since the probability of first recurrence to each cube R , P R ( T ) is the same, theprobability of recurrence to all cubes is: P ( T ) = N m X i =1 ∆ s m m P R ( T ) = N m ∆ s m m P R ( T ) (26)= 2 m ∆ s m ∆ s m m P R [1 − P R ] T − = ∆ s m m (cid:20) − ∆ s m m (cid:21) T − . (27)For small cube side lengths ∆ s and close returns algorithm radius r , the first recurrence probabilitydetermined by the close returns algorithm is then: P ( T ) = ∆ s m m (cid:20) − ∆ s m m (cid:21) T − ≈ r m m (cid:20) − r m m (cid:21) T − . (28)Similarly, for small close returns radius r and/or for large embedding dimensions m , 1 − r m / m ≈ P ( T ) ≈ r m m . (29)Note that for fixed m and r this expression is constant. Since the close returns algorithm can only measurerecurrence periods over a limited range 1 ≤ T ≤ T max , and we normalise the recurrence histogram R ( T )over this range of T , then the probability of first recurrence is the uniform density: P ( T ) ≈ T max , (30)which is proportional to the expression r m / m above. Thus, up to a scale factor, for a uniform i.i.d.stochastic signal, the recurrence probability density is uniform. Authors’ Contributions
MAL lead the conceptual design of the study, developed the mathematical methods, wrote the softwareand data analysis tools, and prepared and analysed the data. PEM participated in the conceptual design ofthe study and the mathematical methods. SJR participated in the discriminant analysis of the data.DAEC participated in the data preparation. IMM contributed to the development of the mathematicalmethods. All authors read and approved the manuscript.
Acknowledgements
Max Little acknowledges the financial support of the Engineering and Physical Sciences Research Council,UK, and wishes to thank Prof. Gesine Reinert (Department of Statistics, Oxford University, UK), Mr26artin Burton (Radcliffe Infirmary, Oxford, UK) and Prof. Adrian Fourcin (University College London,London, UK) for valuable discussions and comments on early drafts of this paper.
References
1. Baken RJ, Orlikoff RF:
Clinical Measurement of Speech and Voice . San Diego: Singular Thomson Learning, 2ndedition 2000.2. Carding PN, Steen IN, Webb A, Mackenzie K, Deary IJ, Wilson JA:
The reliability and sensitivity tochange of acoustic measures of voice quality . Clinical Otolaryngology (5):538–544.3. Dejonckere PH, Bradley P, Clemente P, Cornut G, Crevier-Buchman L, Friedrich G, Van De Heyning P,Remacle M, Woisard V: A basic protocol for functional assessment of voice pathology, especially forinvestigating the efficacy of (phonosurgical) treatments and evaluating new assessmenttechniques. Guideline elaborated by the Committee on Phoniatrics of the EuropeanLaryngological Society (ELS) . Eur Arch Otorhinolaryngol (2):77–82.4. Michaelis D, Frohlich M, Strube HW:
Selection and combination of acoustic features for thedescription of pathologic voices . Journal of the Acoustical Society of America (3):1628–1639.5. Boyanov B, Hadjitodorov S:
Acoustic analysis of pathological voices . IEEE Eng Med Biol Mag (4):74–82.6. Godino-Llorente JI, Gomez-Vilda P: Automatic detection of voice impairments by means ofshort-term cepstral parameters and neural network based detectors . Ieee Transactions on BiomedicalEngineering (2):380–384.7. Alonso J, de Leon J, Alonso I, Ferrer M: Automatic detection of pathologies in the voice by HOSbased parameters . EURASIP Journal on Applied Signal Processing :275–284.8. Titze IR: Workshop on acoustic voice analysis: Summary statement
Jitter in sustained vowels and isolated sentences produced by dysphonic speakers . Speech Communication :61–79.10. Hirano M, Hibi S, Yoshida T, Hirade Y, Kasuya H, Kikuchi Y: Acoustic analysis of pathological voice.Some results of clinical application . Acta Otolaryngol (5-6):432–8.11. Quatieri TF:
Discrete-Time Speech Signal Processing: Principles and Practice . Upper Saddle River, NJ:Prentice Hall 2002.12. Box GEP:
Science and Statistics . Journal of the American Statistical Association (356):791–799.13. Herzel H, Berry D, Titze IR, Saleh M: Analysis of vocal disorders with methods from nonlineardynamics . Journal of Speech and Hearing Research (5):1008–1019.14. Akaike H: A new look at the statistical model identification . IEEE Trans. Automat. Contrl.
AC-19 (6):716–723.15. Maragos P, Potamianos A:
Fractal dimensions of speech sounds: computation and application toautomatic speech recognition . J Acoust Soc Am (3):1925–32.16. Banbrook M, McLaughlin S, Mann I:
Speech characterization and synthesis by nonlinear methods . IEEE Transactions on Speech and Audio Processing :1–17.17. Zhang Y, Jiang JJ, Biazzo L, Jorgensen M: Perturbation and nonlinear dynamic analyses of voicesfrom patients with unilateral laryngeal paralysis . Journal of Voice (4):519–528.18. Zhang Y, McGilligan C, Zhou L, Vig M, Jiang JJ: Nonlinear dynamic analysis of voices before and aftersurgical excision of vocal polyps . Journal of the Acoustical Society of America (5):2270–2277.19. Giovanni A, Ouaknine M, Triglia JL:
Determination of largest Lyapunov exponents of vocal signal:Application to unilateral laryngeal paralysis . Journal of Voice (3):341–354.20. Zhang Y, Jiang JJ, Wallace SM, Zhou L: Comparison of nonlinear dynamic methods and perturbationmethods for voice analysis . Journal of the Acoustical Society of America (4):2551–2560.
1. Hertrich I, Lutzenberger W, Spieker S, Ackermann H:
Fractal dimension of sustained vowel productionsin neurological dysphonias: An acoustic and electroglottographic analysis . Journal of the AcousticalSociety of America :652–654.22. Hansen JHL, Gavidia-Ceballos L, Kaiser JF:
A nonlinear operator-based speech feature analysismethod with application to vocal fold pathology assessment . Ieee Transactions on BiomedicalEngineering (3):300–313.23. Zhang Y, Jiang JJ: Nonlinear dynamic analysis in signal typing of pathological human voices . Electronics Letters (13):1021–1023.24. Jiang JJ, Zhang Y, McGilligan C: Chaos in voice, from modeling to measurement . Journal of Voice :2–17.25. Kantz H, Schreiber T: Nonlinear Time Series Analysis . Cambridge; New York: Cambridge University Press,new edition 1999.26. Behrman A, Baken RJ:
Correlation dimension of electroglottographic data from healthy andpathologic subjects . Journal of the Acoustical Society of America (4):2371–2379.27. Little M, McSharry P, Moroz I, Roberts S:
Nonlinear, biophysically-informed speech pathologydetection . In
Proc ICASSP 2006 , New York: IEEE Publishers 2006.28. McSharry PE, Smith LA, Tarassenko L:
Prediction of epileptic seizures: are nonlinear methodsrelevant?
Nat Med (3):241–2.29. Herzel H, Berry D, Titze I, Steinecke I: Nonlinear dynamics of the voice - signal analysis andbiomechanical modeling . Chaos :30–34.30. Jiang JJ, Zhang Y: Chaotic vibration induced by turbulent noise in a two-mass model of vocalfolds . Journal of the Acoustical Society of America (5):2127–2133.31. Krane MH:
Aeroacoustic production of low-frequency unvoiced speech sounds . Journal of theAcoustical Society of America :410–427.32. LaMar MD, Qi YY, Xin J:
Modeling vocal fold motion with a hydrodynamic semicontinuum model . Journal of the Acoustical Society of America :455–464.33. Ishizaka K, Flanagan JL:
Synthesis of Voiced Sounds From a Two-Mass Model of the Vocal Cords . ATT Bell System Technical Journal (6):1233–1268.34. Steinecke I, Herzel H: Bifurcations in an Asymmetric Vocal-Fold Model . Journal of the AcousticalSociety of America (3):1874–1884.35. Jiang JJ, Zhang Y, Stern J: Modeling of chaotic vibrations in symmetric vocal folds . Journal of theAcoustical Society of America (4):2120–2128.36. Dixit R:
On defining aspiration . In
Proceedings of the XIIIth International Conference of Linguistics , Tokyo,Japan 1988:606–610.37. Howe MS:
Theory of vortex sound . New York: Cambridge University Press 2003.38. Sinder D:
Synthesis of unvoiced speech sounds using an aeroacoustic source model . PhD thesis ,Rutgers University 1999.39. McLaughlin S, Maragos P:
Nonlinear methods for speech analysis and synthesis . In
Advances innonlinear signal and image processing , Volume 6 . Edited by Marshall S, Sicuranza G, Hindawi PublishingCorporation 2007:103–.40. Ockendon JR:
Applied partial differential equations . Oxford; New York: Oxford University Press 2003.41. Acheson DJ:
Elementary fluid dynamics . Oxford; New York: Oxford University Press 1990.42. Kinsler LE, Frey AR:
Fundamentals of acoustics . New York: Wiley, 2d edition 1962.43. Falconer KJ:
Fractal geometry: mathematical foundations and applications . Chichester; New York: Wiley 1990.44. Zhao W, Zhang C, Frankel SH, Mongeau L:
Computational aeroacoustics of phonation, part I:Computational methods and sound generation mechanisms . Journal of the Acoustical Society ofAmerica (5 Pt 1):2134–46.
5. Grimmett G, Stirzaker D:
Probability and random processes . Oxford; New York: Oxford University Press, 3rdedition 2001.46. Berry DA, Herzel H, Titze IR, Story BH:
Bifurcations in excised larynx experiments . J Voice (2):129–38.47. Little M, McSharry P, Moroz I, Roberts S: Testing the assumptions of linear prediction analysis innormal vowels . Journal of the Acoustical Society of America :549–558.48. Tokuda I, Miyano T, Aihara K:
Surrogate analysis for detecting nonlinear dynamics in normalvowels . Journal of the Acoustical Society of America (6):3207–17.49. Tokuda I, Tokunaga R, Aihara K:
A simple geometrical structure underlying speech signals of theJapanese vowel a . International Journal of Bifurcation and Chaos :149–160.50. Jackson P, Shadle C: Pitch-scaled estimation of simultaneous voiced and turbulence-noisecomponents in speech . IEEE Transactions on Speech and Audio Processing (7):713–726.51. Altmann EG, Kantz H: Recurrence time analysis, long-term correlations, and extreme events . Physical Review E (5):–.52. Ragwitz M, Kantz H: Markov models from data by simple nonlinear time series predictors in delayembedding spaces . Physical Review E (5):056201.53. Lathrop DP, Kostelich EJ: Characterization of an experimental strange attractor by periodic-orbits . Physical Review A (7):4028–4031.54. Hu K, Ivanov PC, Chen Z, Carpena P, Stanley HE: Effect of trends on detrended fluctuation analysis . Physical Review E :–.55. Bishop CM:
Neural Networks for Pattern Recognition . Oxford, New York: Clarendon Press; Oxford UniversityPress 1995.56. Boersma P, Weenink D:
Praat, a system for doing phonetics by computer . Glot International (9/10):341–345.57. Boersma P: Accurate short-term analysis of the fundamental frequency and theharmonics-to-noise ratio of a sampled sound . In
Proceedings of the Institute of Phonetic Sciences , Volume 17 , University of Amsterdam 1993.58. Titze IR, Liang HX:
Comparison of F0 extraction methods for high-precision voice perturbationmeasurements . Journal of Speech and Hearing Research (6):1120–1133. FiguresTables
Table 1 - Summary of disordered voice classification results
Summary of disordered voice classification task performance results, for several different combinations ofthe new measures, the derived irregularity (Irreg) and noise (Noise) components of Michaelis [4], andtraditional perturbation measures, jitter (Jitt), shimmer (Shim) and noise-to-harmonics ratio (NHR). TheRPDE parameters were the same as for figure 4, and the DFA parameters were the same as for figure 5.Combination Subjects True Positive True Negative
Overall
RPDE/DFA 707 95.4 ± ± ± Jitt/Shim 685 86.9 ± ± ± Shim/NHR 684 91.4 ± ± ± Irreg/Noise 707 78.4 ± ± ± Jitt/NHR 684 93.2 ± ± ±
00 400 600 800 1000 1200 1400-1-0.500.51 n (samples) s n (a) 200 400 600 800 1000 1200 1400-1-0.500.51 n (samples) s n (b) Figure 1: Discrete-time signals from (a) one normal (JMC1NAL) and (b) one disordered (JXS01AN) speechsignal from the Kay Elemetrics database. For clarity only a small section is shown (1500 samples).
Additional Files
Additional file 1 – close ret.c, 6K
Close returns algorithm implemented in C with Matlab MEX interface. Standard ASCII text file format.
Additional file 2 – close ret.dll, 53K
Close returns algorithm compiled as a DLL for Windows. Standard Windows DLL format.
Additional file 3 – fastdfa core.c, 9K
Efficient implementation of the detrended fluctuation analysis (DFA) algorithm written in C with MatlabMEX interface. Core C code. Standard ASCII text file format.
Additional file 4 – fastdfa core.dll, 53K
DFA algorithm core compiled as a DLL for Windows. Standard Windows DLL format.
Additional file 5 – fastdfa.m, 2K
Matlab function wrapper for above DFA algorithm implementation. Standard ASCII Matlab script fileformat. Type ’help fastdfa.m’ at the Matlab command prompt for usage instructions.30 n s n + (a) s n + -0.500.5-0.5 0 0.5-0.500.5 s n s n + (b) s n + Figure 2: Time-delay embedded discrete-time signals from (a) one normal (JMC1NAL) and (b) one disor-dered (JXS01AN) speech signal from the Kay Elemetrics database. For clarity only a small section is shown(1500 samples). The embedding dimension is m = 3 and the time delay is τ = 7 samples.31
00 400 600 800 1000-101 s n (a) n (samples)
200 400 600 800 1000-101 s n (c) n (samples)
200 400 600 800 100000.51 P ( T ) H norm = 0.06(b) T (samples)
200 400 600 800 100000.050.1 P ( T ) H norm = 0.91(d) T (samples)
Figure 3: Demonstration of results of time-delayed state-space recurrence analysis applied to a perfectlyperiodic signal (a) created by taking a single cycle (period k = 134 samples) from a speech signal andrepeating it end-to-end many times. The signal was normalised to the range [ − , P ( T )are zero except for P (133) = 0 . P (134) = 0 . P ( T ) is properly normalised. This analysisis also applied to (c) a synthesised, uniform i.i.d. random signal on the range [ − , P ( T ) is fairly uniform. For clarity only a small section of the time series (1000 samples) and therecurrence time (1000 samples) is shown. Here, T max = 1000. The length of both signals was 18088 samples.The optimal values of the recurrence analysis parameters were found at r = 0 . m = 4 and τ = 35.32
00 400 600 800 100000.51 P ( T ) H norm = 0.14(a) 200 400 600 800 100000.050.1 P ( T ) H norm = 0.89(b) T (samples)T (samples)
Figure 4: Results of RPDE analysis carried out on the two example speech signals from the Kay databaseas shown in figure 1. (a) Normal voice (JMC1NAL), (b) disordered voice (JXS01AN). The values of therecurrence analysis parameters were the same as those in the analysis of figure 3. The normalised RPDEvalue H norm is larger for the disordered voice. 33
00 400 600 800 1000-101 s n (a) n (samples)
200 400 600 800 1000-101 s n (c) n (samples) norm = 0.54 log L (b) l og F ( L ) (d) 1.699 1.7959 1.8751 1.942 20 = 1.75 norm = 0.85 l og F ( L ) log L Figure 5: Results of scaling analysis carried out on two more example speech signals from the Kay database.(a) Normal voice (GPC1NAL) signal, (c) disordered voice (RWR14AN). Discrete-time signals s n shownover a limited range of n for clarity. (b) Logarithm of scaling window sizes L against the logarithm offluctuation size F ( L ) for normal voice in (a). (d) Logarithm of scaling window sizes L against the logarithmof fluctuation size F ( L ) for disordered voice in (b). The values of L ranged from L = 50 to L = 100 in stepsof five. In (b) and (d), the dotted line is the straight-line fit to the logarithms of the values of L and F ( L )(black dots). The values of α and the normalised version α norm show an increase for the disordered voice.34 .2 0.4 0.6 0.80.550.60.650.70.750.8 H norm no r m -0.5 0 0.5 1-1-0.50 log (Jitt) l og ( NHR ) -0.5 0 0.5 100.51 log (Jitt) l og ( S h i m ) (Shim) l og ( NHR ) (a) (b)(c) (d)(a) (b)(c) (d)