Alexandre Guillet

Alain Arneodo

Pierre Argoul
email: pierre.argoul@email.

Françoise Argoul

Quantifying the Rationality of Rhythmic Signals

published or not. The documents may come   L'archive ouverte pluridisciplinaire

Introduction

Scientific approaches of natural systems have been revolutionized in the last part of the XXth century with the advent of miniaturized electronic and computer systems. Beyond their impressive beauty, it was offered to human beings to demonstrate that nature is constructed from multi-scale intertwined networks, (in time and in space) and that these networks are the field of highly complex nonlinear dynamics (non linear and/or non stationary rhythms) [START_REF]Encyclopedia of Ecology[END_REF][START_REF] Ivanov | [END_REF]3]. Even though apparently distinct biological rhythms (endogenous and exogenous) have been recognised as universal features of all organisms (neural signals, heart, hormone secretion, metabolism, tidal, circadian, lunar, seasonal, annual clocks, life cycle, ....) [START_REF] Goldbeter | Au Coeur des Rythmes du Vivant. La Vie Oscillatoire[END_REF], the variability of these rhythms and their spatio-temporal interplay is still considered as incidental or ignored. Despite the fact that we can concretely demonstrate that the frequencies of these rhythms pave more than 10 decades, still, time (and frequency) is considered as varying linearly in living systems. In particular the presence of strong nonlinearities can give us greater sensing resolution to less intense stimuli. These mechanisms are ubiquitous across animal species and across all sensory modalities. Interestingly, the mappings between an external stimuli and the internal perception (psychophysical) of scales and laws are rather logarithmic than linear. A simple and more commonly encountered example for the non-specialist is the perception and emission of acoustic vibrations (sounds) by living species, these processes occur in logarithmic scales in time and frequency domains [START_REF] Oppenheim | [END_REF]. It has also been demonstrated experimentally that the cochlear filters of the inner ear are not spaced at linear frequency intervals but that their spacing is approximately logarithmic [START_REF] Schnupp | Auditory Neuroscience: Making Sense of Sound[END_REF].

The emission of sound (speech, songs) by human cord tract (larynx, pharynx, mouth) is a complex nonlinear process that combines both muscles and tissues with different temporal and spatial scales, and the entire autonomic and central nervous systems. In this study, we analyse human voice signals (a single note maintained for a few seconds) that characterise the physiology of the vocal organ (larynxpharynx-mouth) in healthy and pathological situations. To compare different signals and their spectral composition, we define a log-frequency correlation measure on spectral densities that gives the distribution of frequency ratios (rational or irrational numbers) between two signals. Using the wavelet transform formalism we extend this measure to a time-frequency correlation measure, that offers the possibility to estimate the temporal variability of this log-frequency correlation. We introduce reference spectral expansions as sums of Dirac terms that resume the characteristic property of these voice signals (harmonics as integer multiples of a fundamental frequency). Finally, we define a new integral cross-correlation of the previously defined measure which quantifies the rationality of the rhythms of two compared signals. We call it sonance, by analogy with the term consonance (resp. dissonance) that counts the perceived affinity or agreement (resp. disagreement) between different sounds. We validate this method on numerical model and voice signals collected from different sources. The first section is this introduction. The second section describes the mathematical methodology for log-frequency correlations (or spectrum of frequency ratios) and its generalization to time-frequency expansions in terms of analytic wavelet transforms. The third section illustrates these concepts on numerical signals (sums of sine functions) and voice recordings from the Voice-Icar-Federico II database, introduces the sonance measure and illustrates it on the previously computed log-frequency correlation measures of voice signals. Finally, we leave the medical application of voice dysphonia diagnosis with the comparison of an untrained voice with a singer voice that have similar spectral envelopes.

Spectrum of frequency ratios. Formalism and time-frequency generalization 2.1 Correlation functions for signal comparison

Let us consider two signals x and y of finite total energy L 2 (R) : x, x < +∞ and y, y < +∞ where •, • is the ordinary inner product of L 2 (R) and x is the complex conjugate of x:

x, y = +∞ -∞

x(u)y(u)du.

(

) 1
The comparison of these two signals x and y is usually performed through a deterministic correlation function R[x, y](ξ) constructed from a time shift (translation) operator T ξ : R[x, y](ξ) = x, T ξ y = +∞ -∞

x(u)y(u + ξ)du.

(

) 2
This definition, given for energy signals or square-integrable functions, can be extended to power signals. Thus, for signals which can be described by sums of periodic functions (stochastic signals with finite power), the cross-correlation function reads:

C[x, y](ξ) = lim T →∞ 1 2T T -T x(u)y(u + ξ)du, (3)
When x = y, we get the auto-correlation function C[x, x](ξ), that characterises the similarity between observations of a same signal as a function of the lag ξ between them. The auto-correlation function is Hermitian: C[x, x](-ξ) = C[x, x](ξ). The absolute value of C[x, x](ξ) is maximum at the origin, where the auto-correlation function is real, positive and equal to the power of the signal x. When the signal x is real, this implies that the auto-correlation function is real and even.

Note that when u = t, t being the time variable, the function C[x, y](ξ) is the cross-correlation function commonly used for time signals but u could be replaced by any other type of variable, and in particular the frequency (or log-frequency) when comparing spectral signals, as will be discussed below.

Translation-based correlation functions are very important for physics. They turn functions of a relative quantity (such as time or space position whose value depends on a translation from an arbitrary origin) into a function of an absolute quantity (such as time or space interval). However, the value of absolute quantities that have a physical dimension still depends on its comparison with an arbitrary standard: the physical unit. Since a scaling is involved, the unit plays the role of an arbitrary origin for the logarithm of these quantities. That is the reason why dilation-based correlation functions can be of interest for physics, as long as they compare functions of an absolute quantity: a new variable made of the ratio of two absolute quantities with the same physical unit neither depends on an origin nor on the unit; it is a pure proportion.

To extend the concept of correlation functions to absolute physical quantities, we need first to revisit the definition of the inner product. We make use of the logarithm to change from the translation-invariant group (R, +) to the dilation-invariant one (R + , ×). The change of variable u = log v applied on Eq. (1) yields:

X, Y = ∞ -∞ X (u) Y (u)du = ∞ 0 X (log v)Y (log v)d log v. (4)
The change from the function X (u) and the measure du to the function X • log(v) and the measure d log v = dv/v means, for numerical computations, that we replace linearly sampled functions by geometrically sampled ones (of positive variable).

In the following, we choose to make explicit the composition with the logarithm in each function. The previous translation operator T ξ is naturally replaced by a dilation operator D q :

T log q [X](log v) = D q [X • log](v) = X (log(qv)). (5)
Combining Eqs ( 4) and ( 5), we obtain from Eq. ( 2) a similar correlation function adapted to geometrically sampled signals:

R[X, Y ](log q) = ∞ 0 X (log v)Y (log(qv))d log v, (6)
where q is positive. For functions X and Y , the finite energy condition for the validity of this integral takes the form X, X < +∞. It can also be reformulated for finite power signals in a similar way as in Eq. ( 3). The dilation correlation function in Eq. ( 6) inherits the following symmetry and linearity properties from Eq. ( 2):

R[Y, X](log q) = R[X, Y ](-log q), (7) R[X, Y + Z](log q) = R[X, Y ](log q) + R[X, Z](log q). (8)
Note that the logarithm does not allow to study functions of a negative absolute quantity (for instance negative delays or frequencies), nor negative ratios q < 0.

Spectrum of frequency ratios: a frequency ratio distribution

For the application of interest here, the unfolding of rhythms from real signals (their spectral "timbre"), we concentrate on "geometric" spectral densities that we define as real and positive functions S(log f ) ≥ 0 of the logarithm of the frequency. The log-frequency correlation function between two such densities

R[S 1 , S 2 ](log q) = ∞ 0 S 1 (log f )S 2 (log(q f ))d log f , (9)
captures all the spectral relations between frequency modes of S 1 (log f ) and S 2 (log f ). R[S 1 , S 2 ](log q) is positive and gives the distribution of frequency ratios q of S 1 (log f ) and S 2 (log f ), hence the notation R for ratio distribution. Similarly to standard correlation function of linearly sampled variables, the existence of this integral R[S 1 , S 2 ](log q) requires that both distributions S 1 (log f ) and S 2 (log f ) be square integrable with the geometric measure of f (linear measure for log f ).

Both the log-frequency distribution S(log f ) and the frequency ratio distribution R[S 1 , S 2 ](log q) can be normalised as probability density functions:

∞ 0 R[S 1 , S 2 ](log q)d log q = ∞ 0 S 1 (log f )d log f ∞ 0 S 2 (log f )d log f = 1 . (10)
Frequency ratio distributions can be written in analytic form from spectral densities defined as isolated or sum of Dirac δ functions. For example, the two spectral densities S j (log f ) = δ(log f f j ), j = 1, 2 have a single frequency ratio f 2 f 1 , and give a frequency ratio distribution R[S 1 , S 2 ](log q) = δ(log q f 1 f 2 ). If we define S(log f ) as a doublet of Dirac deltas S(log f ) = S 1 (log f ) + S 2 (log f ), from the linearity property Eq. ( 8), we can write the ratio distribution R[S, S](log q) = δ(log q f 1 f 2 ) + 2δ(log q) + δ(log q f 2 f 1 ). This simple analytic case is illustrated in Fig. 1, where we distinguish from R[S, S](log q) three peaks, corresponding to the frequency pairs: (4:4) and (8:8) for log q = 0, (4:8) for log q = log 2, and (8:4) for log q =log 2. Fig. 1 (a) Ideal distribution S(log f ) in log-frequencies of a doublet of Dirac deltas such that the highest frequency is twice the lowest. (b) Representation of the spectrum of self-relations R[S, S](log q) in logarithmic scale (base 2). The peak of ratio log 2 q = 0 represents the selfrelation of each frequency peak, whereas the ratios log 2 q = -1, 1 represent their cross-relations.

The log-frequency spectral distributions S(log f ) cannot be assimilated to linfrequency spectral densities defined the from Fourier transform of s: ŝ( f ) = +∞ -∞ s(t)e -2πi f t dt, because their computation from linear measures in time and frequency faces some difficulties. The main one is practical, power spectral densities estimated with Fast Fourier Transform (FFT) algorithms are sampled linearly, whereas the integral of Eq. ( 9) requires a geometric frequency sampling. Re-sampling strategies of the Fourier spectra have been proposed in the literature [START_REF] Haines | [END_REF], and could be used for stationary signals, however they require greater memory size and are computer time-consuming. Importantly, in the context of physiological signals which are often non-stationary, the extension of time-averaged spectral quantities to time-frequency distributions is mandatory. The wavelet transform answers to both issues, it provides not only a time-frequency representation of the spectral quantities, but also allows a geometric sampling in frequency. Using time-frequency decompositions we can straightforwardly extend our definition of log-frequency ratio distributions Eq. ( 9) to time-log-frequency ratio distributions for the analysis of non-stationary signals.

Wavelet transform formalism

Time-frequency analysing tools based on the wavelet transform have been introduced in the second half of the twentieth century and applied to many scientific domains for characterising and modelling non-stationary processes [8,9,[START_REF] Combes | Wavelets: Time-Frequency Methods and Phase Space[END_REF][START_REF] Delprat | [END_REF][START_REF] Carmona | Wavelets and Statistics[END_REF][START_REF] Carmona | [END_REF][START_REF] Carmona | Practical Time-frequency Analysis: Gabor and Wavelet Transforms with an Implementation in S[END_REF]. The wavelet transform of a finite energy signal s(t) ∈ L 2 (R) is defined as its inner product with the shifted copies of an analysing absolute integrable and finite energy wavelet ψ(t) ∈ L 1 (R) ∩ L 2 (R) [9,[START_REF] Carmona | Practical Time-frequency Analysis: Gabor and Wavelet Transforms with an Implementation in S[END_REF][START_REF] Torresani | Analyse Continue par Ondelettes. Collection Savoirs Actuels[END_REF][START_REF] Flandrin | Time-Frequency Time-Scale Analysis[END_REF]:

W (p) ψ [s](a, b) = ψ a, b , s = a -1 p +∞ -∞ s(t)ψ t -b a dt, (11)
b ∈ R and a ∈ R + are the shift and scaling parameters. ψ is the complex conjugate of the analysing wavelet ψ, p is a parameter which defines the normalisation of the wavelet. Two values of p are usually found in the literature: p = 1, corresponding to the L 1 (R) norm and p = 2, corresponding to the L 2 (R) norm, respectively. p = 1, often used for time-localized signals with different amplitudes, is appropriate when the magnitude of the modulus wavelet transform is wished to reflect the amplitude of the analysed signal s(t). p = 2 is appropriate when the modulussquared wavelet transform is wished to reflect the energy of the analysed signal s(t).

In the frequency domain, the expression of the wavelet transform reads:

W (p) ψ [s](a, b) = a 1-1 p +∞ -∞ ŝ( f ) ψ(af )e 2iπ f b d f , (12)
where ŝ, ψ denote the Fourier transforms of the signal and the wavelet. This time-scale representation is quite suited for non-stationary signals since it localizes the analysis around time b and operates a band-pass filtering scaled by the parameter a. Importantly, a can be sampled arbitrarily, in our case we will sample it geometrically. It is common practice to consider the scale a as proportional to an inverse frequency 1 f a :

a = f ψ f a , (13)
where f ψ is a characteristic frequency of the mother wavelet ψ. Three meaningful frequencies are classically used for f ψ [17]: the peak frequency f 0 ψ where the frequency domain mother wavelet magnitude ψ( f ) is maximum, the energy (norm 2) frequency f * ψ which is the mean of ψ( f ) 2 and the norm 1 frequency fψ , that can be interpreted as an instantaneous frequency for progressive wavelets. An asymmetry in the frequency domain of the mother wavelet leads to distinct values for the previous frequencies f ψ . For the computation of the log-frequency correlation functions, the expression W ψ [s]( f ψ / f a , b) for the wavelet transform given in Eq. ( 12) can be turned to a time-frequency analysis by using Eq. ( 13) for a given characteristic frequency f ψ :

W (p) ψ [s] f ψ f a , b = a 1-1 p +∞ -∞ ŝ( f ) ψ f ψ f a f e 2iπ f b d f . ( 14
)
For our applications to physiological signals, the Banach space L 1 (R, dt) norm corresponding to p = 1 will be preferred for the wavelet transform definition. The main reason is due to the fact that when rescaling time in the input signal as s t ρ , with ρ > 0, both the time and the scale of the wavelet transform are rescaled, but without changing its magnitude. Thus as the Fourier transform of s t ρ is: ρ s( ρ f ), Eq. ( 12) when p = 1 leads to W (1) ψ [s] ( a ρ , b ρ ). The (1) is dropped in the following. The peak frequency f 0 ψ will be then adopted for the characteristic frequency f ψ in Eqs [START_REF] Carmona | [END_REF], [START_REF] Carmona | Practical Time-frequency Analysis: Gabor and Wavelet Transforms with an Implementation in S[END_REF]. The admissibility condition for an analysing wavelet ψ ∈ L 1 (R) ∩ L 2 (R) establishes that the number

c ψ = +∞ 0 | ψ(u)| 2 du u ( 15
)
must be finite, nonzero and independent of f ∈ R + . If this admissibility condition is fulfilled, then every s ∈ L 2 (R) can be reconstructed from the convergent integral:

s(t) = 1 c ψ +∞ -∞ ∞ -∞ W ψ [s] (a, b) ψ t -b a da |a| db . (16)

Time and frequency window for the analysing wavelet

The time-frequency window can be computed from the expression of the analysing wavelet ψ, assuming that ψ and ψ verify tψ(t) ∈ L 2 and f ψ( f ) ∈ L 2 (R) [START_REF] Chui | An Introduction to Wavelets[END_REF]. If the center and the radius (with the norm 2) of the window function ψ are respectively t * ψ and ∆ ψ , ψ((tb)/a) is a window function with center b + at * ψ and radius equal to

a∆ ψ : [b + at * ψ -a∆ ψ , b + at * ψ + a∆ ψ ] . (17)
This windows narrows (respectively widens) for small (resp. large) values of a. In the frequency domain, the window of ψ is defined similarly, assuming that the center and width of ψ are f * ψ and ∆ ψ , ψ(a f ) is centered around f * ψ /a and has a radius ∆ ψ /a:

      f * ψ a - 1 a ∆ ψ, f * ψ a + 1 a ∆ ψ       . ( 18
)
In the following discussion, the center f * ψ of ψ is assumed to be positive. There are different ways of defining the wavelet resolution, called the quality factor of the wavelet. A first definition, given in [START_REF] Le | [END_REF], uses the bandwidth and the norm 2 frequency as follows:

Q * = f * ψ /a 2∆ ψ /a = f * ψ 2∆ ψ , (19)
which is independent of the scale parameter a. Alternatively, we could also use the full width at half maximum height of | ψ( f )| 2 instead of ∆ ψ . We thus define another quality factor, Q, such as

Q = f 0 ψ | f 2 -f 1 | (20)
where

| ψ( f 1 )| 2 = | ψ( f 2 )| 2 = | ψ( f 0 ψ )| 2 /2 and f 1 < f 0 ψ < f 2 .
This factor is usually computed to characterise the qualitative damping behavior of simple damped oscillators [START_REF] Rocard | Dynamique Générale des Vibrations[END_REF].

The choice of the quality factor is essential to obtain an adapted time-frequency resolution and consequently a "good" analysis of the processed signals. The authors in [START_REF] Erlicher | [END_REF] propose three bounds to obtain a range of acceptable values. When the signal is composed of several frequency components, the proximity of their characteristic frequencies provides a lower bound. The exponential decay rate of the amplitude imposes another upper bound. Eventually, the length of the signal determines yet another upper bound.

Choice of the analysing wavelet: the Grossmann wavelet

In the absence of a suitable unifying theory for wavelet behaviors, the choice of a particular wavelet for a particular problem may often appear arbitrary. For rhythmic signals, complex analytic analysing wavelets are preferred, leading to: ψ( f ) = 0, ∀ f ≤ 0. In that case, the measure appears naturally in these integrals (Eq. ( 16)), as in Eqs ( 4) and ( 6), because the analysing wavelet is scale invariant (under dilations).

In the following, we choose a single-parameter progressive wavelet, introduced for the decomposition of Hardy functions by Grossmann and Morlet [8]:

ψQ ( f ) =        ψ 0 e -1 2 Q log f f 0 2 ∀ f > 0 ; 0 ∀ f ≤ 0 , (21)
of peak frequency f 0 ψ = f 0 , for which the maximum value is ψ 0 . This wavelet is symmetric in log-frequencies about log f 0 . The other characteristic frequencies are

f * ψ = f 0 e 3 4Q 2 and fψ = f 0 e 3 2Q 2
. The Grossmann wavelet is also centered in time (t * ψ = 0), with a width (radius in norm 2):

∆ ψ = 1 + 2Q 2 4π f 0 . ( 22
)
Both previously defined quality factors depend on Q only:

Q * = 1 2 e 1 2Q 2 -1 -1 2 , (23)
Q = 2 sinh log 2 Q -1 . ( 24
)
When Q is large enough, the leading term in the expansions gives Q * Q/ √ 2 and Q Q/ log 2 respectively, followed by a term of order 1 Q . Consequently, we will refer to the parameter Q as the quality factor for this wavelet.

When choosing the value ψ 2 0 = Q √ π , the admissibility constant is one and | ψ( f )| 2 can be considered as a probability density function in log-frequencies:

c ψ = ∞ 0 | ψQ ( f )| 2 d log f = 1 . ( 25
)
In Fig. 2, we plot the Grossmann wavelet for two values of Q, respectively Q = 8 (top plots) and Q = 64 (bottom plots). For larger Q values, the number of oscillations of ψ Q (t) and its width increases whereas ψQ ( f ) narrows. We can observe in Fig. 2(c,f) that the wavelet in the log-frequency domain is symmetric around log f 0 = 0 whereas it is asymmetric around f 0 = 1 in linear frequencies (Fig. 2(b,e)). An important aspect of the oscillating progressive wavelets is how many oscillations are fitting inside their time-window [START_REF] Le | [END_REF]17]. This number of oscillations determines the acuteness of the local frequency detection of a given rhythm and is of order Q. If this number is too large, the wavelet averages over too much oscillations and cannot provide a correct estimation. Conversely, if the number of oscillations is insufficient (less that ∼ 3) the detection of a local rhythm will not be possible. The choice of this parameter is particularly important if the signal presents sharp transitions or close frequencies, as will be illustrated in the following figures.

The authors in [17] showed that the Grossmann wavelet can be seen as a scaling limit of a general family of progressive wavelets with two parameters, the Morse wavelet [22,23,24,25,26,27,28]. The Cauchy-Paul wavelet, intensively used in quantum mechanics and in the context of analytic functions [START_REF] Paul | Wavelets and Their Applications[END_REF], as well as the analytic version of the derivative of Gaussian wavelet or the Airy wavelet all belong to the Morse family.

Extension of frequency ratio distributions to time-frequency ratio distributions

From the Grossmann progressive wavelet transform defined in the previous section, we define a time-frequency distribution for non-stationary signals:

S (Q) (log( f a ), b) = W ψ Q [s] f 0 f a , b 2 . ( 26
)
Note that the integral of the wavelet transform definition in Eq. ( 12) is sampled linearly in f , but that the values of the frequencies f a (or scale a) can be chosen arbitrarily, for our purpose we will select them geometrically distributed. In the following, b = t and f a = f are considered as time and frequency parameters, which simplifies the notation of S (Q) (log f , t). This distribution is computed for strictly positive values of f and we can extend the definition of the cross-correlation function to time-frequency distributions:

R[S (Q) 1 , S (Q) 2 ](log q, t) = ∞ 0 S (Q) 1 (log f , t)S (Q) 2 (log(q f ), t)d log f (27) = ∞ 0 W ψ Q [s 1 ] f 0 f , t 2 W ψ Q [s 2 ] q f 0 f , t 2 d log f . (28)
The log-frequency autocorrelation function is defined as R[S (Q) , S (Q) ](log q, t). The temporal mean of S (Q) (log f , t): S (Q) t (log f ) that can be seen as a power spectral density based on the wavelet transform

W ψ Q [s].

Computation of the log-frequency correlation function

Using the convolution theorem, R[S (Q) 1 , S (Q) 2 ]
can be computed quite efficiently using the fast Fourier transform (FFT) several times (that discretizes the Fourier transform here denoted F ): on a first step with respect to the time variable the signal (noted F ), and on a second step with respect to the log-frequency variable (noted F log f ) and the computation step is an inverse FFT in log-frequency space (noted

F -1 log f ). R[S (Q) 1 , S (Q) 2 ](log q, t) = F -1 log f F log f |W ψ Q [s 1 ](., t)| 2 F log f |W ψ Q [s 2 ](., t)| 2 (log q) (29)
where

W ψ Q [s]( f , t) = F -1 ψQ f f F [s]( f ) (t) . (30)
This supposes that the frequency f (or scale a) parameter of the CWT is sampled geometrically. The slowest operations consist in matrix multiplications. The fact that the second step requires Fourier transforms of the distributions S on log-frequency scale implies that the computed range of log-frequency values is enlarged, and padded with zeros to avoid extra-ratios arising from the FFT computation by an artificial periodisation of the S distribution.

Computation of log-frequency distributions from numerical and real signals 3.1 Model signals constructed from sine functions

In Fig. 3, we construct an artificial non-stationary signal from the sum of two sine functions: s(t) = sin(φ 1 (t)) + sin(φ 2 (t)), with φ 2 (t) = 4πt linear in time, and (8) , S (8) ](log q, t ). (f) R[S (64) , S (64) ](log q, t ). R[S (Q) , S (Q) ](log q, t ) is defined in Eq. (27).

φ 1 (t) = 2πt H (-t) + 3πt H (t)
With this signal we estimate a lower bound of Q that is suitable for a frequency discrimination according to [START_REF] Erlicher | [END_REF]: for t < 0, Q 10 and for t > 0, Q 14. Moreover, the signal length gives the constraint Q 285. For t < 0, the signal possesses two frequencies, highlighted on the colour-coded image of S (Q) (log f , t) by two horizontal bands ( f 1 = 1 and f 2 = 2), the width of which depends on the quality factor Q, (Q = 8 near the lower acceptable Q bound in Fig. 3(b) and Q = 64 in Fig. 3(c)). For t > 0, we can again recognise the two bands f 1 and f 2 , and, as for t < 0, their narrowing for the larger Q values. The transition zone of this two bands, below and above t = 0, needs to be discussed. Fig. 4 highlights this transition with sections of S (Q) (log f , t) performed for remarkable values of f ; 1, 3/2 and 2. From the sections of Fig. 4(a), we estimate the width of this transition ∼ 4.8s for Q = 8, and ∼ 39s for Q = 64. Another interesting phenomenon emerges in the t > 0 regime, where the two frequency bands become closer. A low-frequency modulation of the wavelet transform squared modulus in the intermediate frequency range [ f 1 , f 2 ] with period 2s appears, corresponding to frequency f m = f 2f 1 (0.5Hz in this example). The matrix of the wavelet transform modulus is not simply the superimposition of the wavelet transform squared moduli of the sine alone,

|W ψ Q [s 1 + s 2 ](a, b)| 2 |W ψ Q [s 1 ](a, b)| 2 + |W ψ Q [s 2 ](a, b)| 2 , but extra terms such as 2 ψ(a f 1 ) ψ(a f 2 ) cos(2π( f 2 -f 1 )b)
are also involved and are not negligible when f 1 and f 2 become too close (which is the case in Fig. 3(b)). This effect disappears quite completely for larger Q values because the product ψ(a f 1 ) ψ(a f 2 ) vanishes. We conclude that the choice of Q is a compromise between two objectives, (i) discriminating close frequencies (in which case larger Q values will be preferred), (ii) affording a correct temporal resolution for the detection of steep frequency changes (in which case smaller Q values will be more efficient).

Fig. 4 (a,c,e) Sections of the log-frequency spectral distributions S (8,64) (log f i , t ) selected from Fig. 3(b,c) for the three frequencies f 1 = 1, f 2 = 3/2 and f 3 = 2 Hz, and the same values of the quality factor Q (8 and 64). (b, d, f) Sections of R[S (Q) , S (Q) ](log q i , t ) selected from Fig. 3(e,f) for three values of q: q 1 = 1, q 2 = 4/3, q 3 = 2, corresponding to local maxima of R[S (Q) , S (Q) ]. For each selected q value, R[S (Q) , S (Q) ] was scaled by its maximum in the time interval. Fig. 3(e,f) shows the corresponding colour-coded maps R[S (Q) , S (Q) ](log q, t) for the same signal and the same values of Q (8 and 64). We recognise for t < 0 three horizontal bands of constant q, corresponding respectively to frequency ratios q = 1/2, 1, 2. The intensity of the middle band (q = 1) is more contrasted (×2) because it corresponds to the sum of self-relations ( f 1 : f 1 ) and ( f 2 : f 2 ). The two symmetric weaker bands correspond to cross-frequency ratios ( f 1 : f 2 ) and ( f 2 : f 1 ). For t > 0, the three bands become closer, and similarly to the maps of S (Q) in Fig. 3(b), a slow temporal modulation of R[S (Q) , S (Q) ](log q, t) superimposes to the bands, due to coupling terms in the wavelet transform modulus. As expected, and similarly to what was observed on S (Q) (log f , t) maps, increasing Q from 8 to 64 produces a strong narrowing of the bands and a strong reduction of the low-frequency modulation. The sections at fixed q of these R[S (Q) , S (Q) ](log q, t) maps are shown in Fig. 4 to highlight similarly the transition zone around t = 0, its widening for larger Q values, and the slow temporal modulations observed for Q = 8. Due to the use of the Grossmann wavelet, sections at fixed t of both S (Q) (log f , t) and R[S (Q) , S (Q) ](log q, t) are Gaussian of widths ( √ 2Q) -1 and Q -1 respectively when the bands are not interfering (independent of t).

Another family of model signals (Fig. 5(e)), particularly interesting with respect to the applications to voice signals, is defined as the sum of sine functions n i=1 sin(2π f i t) with f i = i f 1 , i positive integer. We use here a simple model with discrete and constant frequency components, that does not pretend to account neither for the intrinsic randomness nor the nonstationarity of physiologic signals. To improve the matching of model equations with real signals, we suggest two recent works based on the implementation of stationarity-breaking operators on Gaussian stationary random signals [START_REF] Meynard | Spectral analysis for nonstationarity audio[END_REF]31]. In Fig. 5, we take n = 6 and perform the same time-frequency analysis with a Grossmann analysing wavelet with two values of Q, respectively Q = 8 (b,c) and Q = 128 (f,g). We note again that the larger Q, the finer and distinguishable the peaks of both S (Q) (log f , t) and R[S (Q) , S (Q) ](log q, t). The already noticed low-frequency modulations in the previous example again appear in this example (Fig. 5(b,f)) for Q = 8. Amazingly the frequency of this slow mode is precisely the fundamental frequency of this signal, and this modulation is the most intense for the highest harmonic ( f 6 = 6 f 1 ), this effect is due to the ordering of these 6 frequencies as integer multiples of f 1 , giving a constant frequency step between successive harmonics f i+1f i = f 1 . This f 1 = 1 Hz slow modulation mode appears when two frequencies of the list are too close (in log-scale) for being separated properly by the analysing wavelet. (8) , S (8) ](log q, t ). (g) R[S (128) , S (128) ](log q, t ). (h) R[S (Q) , S (Q) ](log q, t = 0). In (d) and (h) the plots for Q = 8 (resp. 128) are coloured in black (resp. red). R[S (Q) , S (Q) ](log q, t ) is defined in Eq. ( 27).

We observe that in the time-frequency distribution of s shown in Fig. 5(b) for Q = 8, the higher harmonics cannot be distinguished, and their separation requires increasing markedly the value of Q (for instance Q = 128 in Fig. 5(c)). Fig. 5(d) illustrates two sections of these distributions for t = 0 (black curve, Q = 8 and red curve, Q = 128). This phenomena is even more visible on the ratio distribution R[S (Q) , S (Q) ](log q, t) in Fig. 5(h). R[S (Q) , S (Q) ](log q, t) presents an odd number of peaks, it is symmetric around the central peak (q = 1). In that example, each of the sixth frequency components contributes to this central peak. 11 lateral peaks emerge for q > 1, and accumulate closer to the central peak. The positions of these peaks correspond to all the possible distinguishable frequency ratios of the signal, and the amplitude of these peaks is proportional to the number of combinations of frequencies that produces a given ratio. To distinguish all the peaks in Fig. 5(g,h) it was necessary to increase Q to 128. The total number of frequency ratios (for q > 1) is n-1 i=1 i = n(n -1)/2 if n > 2 , in this example it is equal to 15. When there is no redundancy in the frequency ratios, for instance if harmonic frequencies are prime multiples of the fundamental frequency, each frequency ratio occurs once in R[S (Q) , S (Q) ](log q, t).

Physiological signals: voice recordings

The voice signals reported in this manuscript were selected from the VOice-ICar-fEDerico II (VOICED) database [32] recorded by the "Institute of High Performance Computing and Networking of the National Research Council of Italy (ICAR-CNR)" and the Hospital University of Naples "Federico II" during 2016 and 2017. This database can be downloaded from the PhysioNet website [33]. It has been proposed lately as a new element in research on automatic voice disorder detection and classification. Together with medical phonetic examinations of a set of 208 individuals, among which 73 male and 135 female, voice signals, proportional to a local sound emission intensity, were acquired for about 4-5 s and sampled at 8000 Hz at 32 bit, vocal folds were examined by laryngoscopy and two medical questionnaires were collected at the ambulatories of Phoniatrics and Videolaryngoscopy of the "Federico II" Hospital of the University of Naples or at the medical room of the ICAR-CNR. The protocol description is reported in [32]. Dysphonia is a quite common voice disorder (1/3 of adults will suffer from it once in their lifetime), it may originate from a functional or organic alteration of the vocal apparatus and its mechanics and may not systematically be considered as pathologic [START_REF] Le Huche | La Voix, Tome 2: Pathologies Vocales d'Origine Fonctionnelle[END_REF][START_REF] Le Huche | La Voix, Tome 3: Pathologies Vocales d'Origine Organique[END_REF]. On the one side, laryngoscopy is an invasive technique that gives a direct view of the physical alterations of the vocal tract [START_REF] Berke | [END_REF]. On the other side, the analysis of the voice acoustic signal is not intrusive and, thanks to the improvement of signal analysis methods, it can nowadays be used to guide or assist the recognition of the origin of a suspected dysphonia. Voice classification methods from voice recordings by the recognition quantification of the voice timbre (or tone color) has rapidly attracted the interest of electronic and computer science engineers. Globally, one can classify these methods in three groups [START_REF] Jouvet | 2017 25th European Signal Processing Conference (EUSIPCO) (IEEE[END_REF], (i) the time-domain methods which use autocorrelation functions or their variants [START_REF] Ross | [END_REF]39] to search for repeatability between a temporal waveform and its time lagged version, (ii) frequency domain methods which locate characteristic frequencies and conclude to a spectral "coloraturas" for the voice, these methods meet rapidly their limitations if the signal is not stationary, (iii) time-frequency domain techniques [START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF][START_REF] Cohen | Time-frequency Analysis[END_REF][START_REF]Hilbert-Huang Transform and its Applications[END_REF], that we have chosen for this study.

The voice signal s(t), numbered #008, is that of a female of 51 years without deep vocal impairment at the time of the test, ranked in the group of reflux laryngitis (Fig. 6). This example was chosen because it has marked peaks which can be detected by thresholding the signal (this is quite rare, because it requires both a particular shape of the signal and a global stationarity of its amplitude). The Fourier spectra of this signal (reported in log-log and log-lin scales in Fig. 6(b) and (c) respectively) weight the power (in log-scale) of its spectral components; a fundamental mode with frequency f 1 ∼ 188 Hz and higher modes (harmonics), ranked as integer multiples of f 1 : i f 1 with i = 2, 3, 4, 5, ... with different powers. This simple frequency decomposition was observed in most of the signals provided in the VOICED database, this is a conspicuous characteristic of the human voice. These voice signals appear as the alternance of quite regular large and sharp peaks (which give the fundamental mode) and smaller oscillations which may be very irregular. In some cases these smaller oscillations may be difficult to discriminate from the noise produced by some friction of the vocal tracts. Even though this type of signal can be compared to the sum of sine functions introduced in Fig. 5(e), the higher number of harmonics of this signal and their different power means that it could be reproduced by a nonlinear dynamical system (ruled by nonlinear ordinary differential equations) where the different frequency components follow nonlinear rules [START_REF] Travieso | [END_REF]. Our purpose in this paper is not to discuss the physical and biological mechanisms or the modelling of voice signals, we have selected these examples as illustrations for our log-frequency correlation method because their spectral decomposition is very rich in harmonics (overtones) of the fundamental frequency.

The temporal change of the fundamental mode frequency f 1 (t i ) and the largest peak amplitude A(t i ) can be extracted from the #008 voice signal by thresholding its largest amplitude peaks (maxima: s P (t i ) and minima s p (t i )) as depicted in Fig. 6(e). Fig. 6(c) shows that f 1 (t) is modulated in time, suggesting an irregularity of the rhythm coming from some difficulty of the patient to maintain a constant value of f 1 . In this example, a similar temporal modulation is also visible on the largest peak amplitude A(t i ) (Fig. 6g). If these temporal variations were solely produced by instrumental noise, the first return scatter plots of f 1 and A at successive peaks would give a symmetric cloud of points around the diagonal. In Fig. 6(d) for the fundamental frequency modulation and in Fig. 6(h) for the amplitude modulation these first return scatter plots are anisotropic, meaning that the dispersion of these values extends beyond instrumental noise. This conclusion is also confirmed by the temporal evolution of f 1 (t i ) (Fig. 6(c)) and A(t i ) (Fig. 6(g)), we notice that, in the first second, the modulations of f 1 (t i ) have the largest amplitude and are quasiperiodic, this first regime can also be recognised from the modulations of A(t i ). This patient has a rather mild dysphonia (classified as produced by reflux laryngitis) which can be recognised by an important set of harmonics and a rather low vocal fold noise. (20 ms) showing the local maxima s P (t i ) (black dots) and minima s p (t i ) (black stars) which are used to compute both the local frequency: f 1 (t i ) = 1/∆T i and the amplitude A(t i ) of each larger amplitude peak:

A(t i ) = s P (t i ) -s p (t i ). (f) | ŝ | (
in base 10 log-scale) plotted versus f in Hz (linear scale). (g) Amplitude of the largest peaks A(t i ) versus time (see (e) for their detection). (h) First return peak amplitude (A(t i+1 ) vs A(t i )) scatter plot.

The second voice signal illustrated here is that of a female of 62 years (#169 in the VOICED database), with hyperkinetic dysphonia. In that case, the quasiperiodicity observed in signal #008 is so much disrupted that it is impossible to use the previous threshold method for extracting the largest signal peaks; the time-frequency analysis is required to check to which extent we can find a timbre for this voice, and how it changes with time. Fig. 7(b,c) reports the colour-coded images of the #008 (s 1 ) and the #169 (s 2 ) time-frequency distributions S (64) (log f , t). With a temporal averaging of these time-frequency distributions S (Q) (log f , t) t , we get a smooth estimate of power spectrum distributions for these two examples (Fig. 7(d)). From the logfrequency filtering by the Grossmann analysing wavelet, the shapes of the averaged peaks shown in Fig. 7(d) are in general different from those which could be obtained with a Welch estimator [START_REF] Hayes | Statistical Digital Signal Processing and Modeling[END_REF]. The fundamental band frequency of #169 is much broader than that of #008, and shifted to greater values f 1,1 ∼ 188.8Hz (voice #008) and f 1,2 ∼ 268Hz (voice #169), and we also note that it is quite impossible in #169 to discriminate more than one harmonics from the averaged frequency spectrum. The time frequency distributions in Fig. 7(b,c) highlight these differences. Whereas the fundamental mode band and its harmonics are weakly modulated in time for #008, that of #169 are very irregular, the third and fourth harmonics can be mixed up and indistinguishable, the harmonics above five are no longer visible. The vocal folds of #169 can no longer maintain their tight contact that is essentiel for a correct sound emission and the resulting effect, when hearing the voice, is that of a scratching noise which covers completely the expected tone. This person is quite unable to sing a melody.

Tuning voice pitches via the computation of correlation functions

Reference frequency distribution S 0, j

For each signal s j , an "ideal" frequency distribution S 0, j is first introduced in order to compare the real frequency distribution S (Q) j with a reference through the crosscorrelation defined in Eq. ( 27). Let us consider the vibrating string model used to represent the sounds emitted by stringed instruments. When the string is plucked at its ends, its natural frequencies are integer multiples of the fundamental frequency depending on the square root of the force of tension of the string [START_REF] Gérardin | Mechanical Vibrations. Theory and Application to Structural Dynamics[END_REF]. By analogy to this model, the reference frequency distribution S 0, j (log f ) is a Dirac comb model in log scale defined as a sum of integer multiples of the fundamental frequency f 1, j of the studied signal

S 0, j (log f ) = n c n δ log f n f 1, j , (31)
with weights c n ≥ 0 and possible cut-off (c n = 0, ∀n > N). We assume that the series of c n is bounded: n c n < +∞.

We compare then the spectrum of the comb reference model to itself. We build an ideal ratio distribution by computing the auto-correlation of S 0, j derived from Eq. ( 28):

R 00 (log q) = R[S 0, j , S 0, j ](log q) = n m c n c m δ log q n m , (32)
this ratio distribution depends neither on time, nor on the signal under study. When log q = 0 (q = 1), R 00 (0) = n c 2 n < +∞.

Cross-correlation

R[S 0, j , S (Q) j ]
The time distribution cross-correlation function between S 0, j and S (Q) j is deduced from Eq. ( 28):

R[S 0, j , S (Q) j ](log q, t) = n c n S (Q) j (log(qn f 1, j ), t). (33)
Different degrees of "frequency matching" can also be captured by high peaks in R[S 0 , S (Q) j ](log q, t), especially when q is a simple frequency ratio of harmonics of f 1, j and f 1,0 , for example 1:2 -octave, 2:3 -fifth, 3:4 -fourth would give a perfect consonance, 3:5 -major sixth and 4:5 -major third would give a medial consonance, 5:6 -minor third and 5:8 -minor sixth would give imperfect consonance. "Unmatched" frequency configurations would be obtained if a couple of frequency ratio of harmonics belong to the dissonance list: 8:9 -major second, 8:15 -major seventh, 9:16 -minor seventh, 15:16 -minor second, 32:45 (∼ 1/ √

2) -tritone [START_REF] Helmholtz | Theorie Physiologique de la Musique Fondée sur l'Etude des Sensations Auditives[END_REF]. In the following, we will denote R i j = R[S (Q) i , S (Q) j ], the time-frequency window is fixed (Q = 64).

Application to two voice signals from the VOICED data base

For each of the two voice signals #008 and #169, we construct a Dirac comb model as reference distribution. These distributions are such that their lowest frequency peak matches the signal fundamental frequency (for instance for the signal #008: f 1,1 = 188.8Hz and for signal #169: f 1,2 = 268Hz). The frequency of the highest harmonic of the comb model is limited by the sampling frequency F s : n f 1, i F s /2 (F s = 8000 Hz). We take c n = 1, ∀n ≤ 15:

S 0, j (log f ) = 15 n=1 δ log f n f 1, j , j = 1, 2. (34)
For numerical computations, the frequency f is discretized in

f k = f min • α k-1
, with k = 1, 2, 3, ...N and N the size of the frequency vector f , f min its minimum value and α the geometric factor determined from N, f min = 100 Hz (fixed by the voice database), and f max = F s /2. The correlation function R 0 j (log q, t) is computed by combining this comb distribution with that of the voice signal as in Eq. ( 29), using the analytic expression of the log( f )-Fourier transform (F log f ) of the comb distribution. We do not take its Fourier transform numerically because it is the source of numerical artefacts. For the comb model aligned to the fundamental frequency f 1, j of signal s j , it reads:

F log f S 0, j (u) = 15 n=1 exp(-i2πu log(n f 1, j / f min )). (35)
u is the conjugated variable (through Fourier transformation) of log( f ). In that space, we take the scalar product of F log f S 0, j with the conjugate of

F log f S (Q) j (computed numerically from W ψ Q [s j ](log f , t))
and compute its inverse Fourier transform to recover the correlation function R[S 0, j , S (Q) j ](log q).

R 0i = R[S 0, j , S (Q) j t ](log q) ( j = 1, 2, Q = 64)
with their reference comb frequency distributions (defined in the text). These correlations have been normalised to their maximum for ease of comparison. (b, c) Plot of R[S 0, j , S (Q) j t ] versus q for q > 1 (b) and 1/q for q < 1 (c).

The time-averaged frequency ratio distributions R[S 0, j , S (Q) j t ](log q) of each voice signals #008 and #169 with its "best-fitted" Dirac comb model are presented in Fig. 8. If the signals were regular and quasi-stationary, these ratio distributions should pinpoint ratios corresponding to the multiples of the fundamental frequencies. The plot of these correlation functions in linear q and 1/q scales in Figs 8(b) and 8(c) highlights a strong asymmetry, it is due to different amplitudes of the fundamental mode and its harmonics compared to the constant coefficients in the comb model. Again, as for frequency distributions, we note a strong difference of the ratio distributions for signals #008 (s 1 ) and #169 (s 2 ). Confronting R[S 0 , S (Q) ](log q, t) with S (Q) (log f , t) for voice signal #169 (Fig. 9) unveils important features which were not visible from the time averaged ratio distribution R[S 0 , S (Q) t ](log q) (Fig. 8). Even if the fundamental mode frequency and its harmonics vary a lot during these 3s record, their ratios do not change dramatically, as a characteristic property of the mechanics of the vocal folds. In the middle of this signal (1.4s< t < 1.55s) (see Fig. 9(a) for a zoom in this interval), four flat ratio bands can be noticed, suggesting that this person put sufficient effort to recover for a short period of time a "mild sensation" of timbre. How this intermittent loss and recovery of the voice timbre occurs, the time range of these alternating sequences could be used as diagnosis criteria or aftercare follow-up (invasive intervention is necessary if soft or hard nodules are detected on the vocal cords (stage III), or voice exercises for earlier stages). It has been recently shown for patients with neurodegenerative diseases, that not only the patient's ability to speak and formulate sentences is altered but also their voice purely acoustic features [START_REF] Al-Hameed | [END_REF]. t ](log q). (f) Colour-coded map of the time-ratio distribution R[S 0 , S (Q) ](log q, t ).

Cross

-correlation R[S (Q) i , S (Q) j

] of two voice signals

There are two interpretations for the log-frequency cross-correlation function R[S (Q) i , S (Q) j ](log q, t), leading to different possible applications. Either we see it as a distribution of the ratios q between the frequencies of S (Q) i (log f , t) and S (Q) j (log f , t), or we view it for each q as a measure of how well S (Q) i (log f , t)

and S (Q) j (log(q f ), t) match. For instance, in the example reported here, S (Q) 1 , S (Q) 2 are obtained from two different persons holding a pitch in their vocal range. The peaks in the correlation function R[S (Q) 1 , S (Q) 2 ](log q) indicate the importance of the corresponding frequency ratios between the voices, in accordance with its first interpretation as a ratio distribution. When these ratios are close to simple rational numbers m/n, they indicate the presence of a m : n-synchronisation, that corresponds to the consonance of the voices simply sung together. This outlines our strategy to assess how rational the spectral relations of the voices are. The other interpretation is as follows: assuming S (Q) 2 (log q f ) models the second voice transposed by q to a different pitch, the peaks in R[S (Q) 1 , S (Q) 2 ](log q) also indicate for which pitch transpositions the second voice would best match the first voice. This allows us to tune one voice with the other.

The possibility to match a real voice signal with reference model signals is very interesting because it can limit the maximum harmonics frequency for this cross-correlation, as an "intelligent" low pass filtering. This would not be possible by computing directly R[S (Q) i , S (Q) j ]. We compare in Figs 10 and 11 two different "normal" voices ("a" vowel) from the clinic research in the speech therapy laboratory UNADREO in Toulouse (France). The frequency distributions S (Q) i (log f , t) and S (Q) j (log f , t) (i j) plotted in Fig 10, have the same characteristic frequency peaks structure as the voice #008, but we notice that the frequency distribution of the voice s 1 has greater energy in the harmonics around 1000 Hz, which is characteristic of the emission of trained singer voices.

The cross-correlation ratio distribution R

i j = R[S (Q) i , S (Q) j ] is quite different from the auto-correlation ratio distributions R 00 , R (Q)
ii and R (Q) j j . A common reference Dirac comb is chosen for both s 1 and s 2 (i = 1 and j = 2) and is aligned to the fundamental frequency f 1,1 = 307.2Hz. The highest central peak indicates the ratio of the fundamental frequencies, it is centered for R 00 , R 11 and R 22 and shifted to q = f 1,2 / f 1,1 = 247/307.2 ∼ 2 -0.315 for R 12 . R 00 shows very sharp and narrow peaks which line up symmetrically on either sides of q = 1. The amplitude of these peaks recapitulates the weighting of the frequency ratios q for simple comb models and gives us which distribution would be obtained if all the frequency components of the signals had exactly the same power. We finish with the question of comparing the ratio distributions of the voices to the ideal reference.

Frequency distribution matching and sonance

We propose to find the best match between R 12 and R 00 by computing yet another cross-correlation:

R[R 00 , R 12 ](log x, t) = ∞ 0 R 00 (log q)R 12 (log xq, t)d log q (36) = ∞ 0 R 00 log q x R 12 (log q, t)d log q = n m R 12 log m n x, t . (37)
This new quantity can be obtained by two equivalent paths:

R[R 00 , R 12 ](log x, t) = R[R 01 , R 02 ](log x, t), (38)
each corresponding to two uses of the newly introduced parameter x that we call pitch transposition. Either R 12 (log xq, t) is seen as the distribution of ratios q between the first voice S (Q) 1 (log f , t) and the second voice of transposed pitch S (Q) 2 (log x f , t). Or the x in R 00 (log(q/x)) is seen as a varying ratio between the fundamental frequency of the ideal distribution S 0 .

Indeed, the best matching is expected when the pitch of the second voice is transposed to match the first one, thus when the voices are sung at unison, or equivalently when the fundamental frequency ratios are matched between the pairs of distributions: x = f 1,2 / f 1,1 . As a result, for two voices of fundamental frequency f 1,1 and f 1,2 /x, the quantity R[R 00 , R 12 ](log x, t) as a function of the pitch transposition x has the following interpretation: it measures how "ideal" (similar to the model R 00 ) the spectral relations are between the voices. Extrema of this curve appear directly related to the musical property of consonance or dissonance of certain fundamental frequency ratios. For this reason, we call this quantity the sonance between the two voices, and we rewrite its definition (36) using the reference ratio distribution as the density of a measure dG(log q) = R 00 (log q)d log(q): G[R 12 ](log x, t) = R[R 00 , R 12 ](log x, t) = ∞ 0 R 12 (log xq, t)dG(log q), (39) that we could denote equivalently G[S (Q) 1 , S (Q) 2 ](log x, t). This sonance measure is a geometric function of a pitch transposition quantity x, its maxima indicate the optimum relative pitch transpositions for which the two voices sung together would match best. This term sonance bears some analogy with the concepts of consonance and dissonance which were first suggested by Pythagoras (sixth century BC), hence our choice of the symbol G, but this similarity of terms must be nuanced. Dissonance and consonance are not mathematical quantities since they have been used to describe an empirical sensation of human beings (combination of cochlea physiology and cognitive training) when hearing a mixture of sounds (two or more) [49,50].

In Fig. 11, we compute the sonance of the two voice records shown in Fig. 10, from time-averaged frequency distributions (for convenience) referenced to the same comb model. The comparison of the sonance profiles in (c) to the ones of the generating auto-and cross-correlation functions R 01 and R 02 in (a), and R 12 , R 11 , R 22 and R 00 in (b) draws our attention to important features. Apart from the asymmetry (symmetry) of the cross-(auto-)sonances, expected for any correlation function, the sonance profiles have a strong positive baseline and are much less peaked than their simpler counterparts R i j . This can be linked to the combined effects of a dense forest of peaks in R 00 and their width in the voices ratio distributions. Nevertheless, this baseline can be subtracted, as shown in Fig. 11(d), for further comparison. The highest peaks of G[R i j ], pointing to the interval changes between the two voices that would lead to the best consonances, are very similar to the ones already present in R i j . For instance, the global maxima of G[R 12 ] to the left of x = 1 corresponds to the ratio between the fundamental frequencies (voice s 2 has a lower pitch than voice s 1 ). However, their prominence are different, and sonance profiles contain a more detailed landscape of smaller peaks and wells, tracking the gains and losses of consonance of the voices sung simultaneously at the corresponding relative pitch transposition. The almost inexistent peak of G[R 12 ](0) indicates that these two voices sung together without any tone adjustment would be quite dissonant because their fundamental and harmonic frequencies have few commensurability: the frequency ratios are not close to simple rational numbers m/n. In agreement with the common intuition, the central value of the auto-sonance is also higher for the singer voice than for the untrained one, G[R 11 ](0) > G[R 22 ](0).

The implementation of a frequency reassignment would be particularly beneficial to this time-dependent application, improving successively the discrimination of higher harmonics in the CWT, the resulting ratio distributions and the resolution of the sonance landscape. In particular, recent developments in this area based on the synchrosqueezing transform [51,52] could be important for further developing the sonance concept proposed herein.

As a last remark, the sonance profile of the voices is directly influenced by two choices: first, the quality factor of the wavelet Q, which determines the distinguishability of the frequencies and their ratios, and second, the choice of the number of harmonics and their amplitude in the reference comb model S 0 . We believe that, for a realistic sonance profile, Q should be related to the critical band of the ear [49] and is, together with the design of the reference ratio distribution R 00 , representative of the musical training of the ear.

Conclusion

We have introduced time-log-frequency ratio distributions based on analytic wavelets that we have applied to model and physiological signals (voice records). We found that the Grossmann wavelet is a natural shape for this task. A second correlation operation was defined to compare the matching of the voices ratio distribution with an ideally rational one, called sonance. This function of a pitch transposition estimates, in a sense, the "harmony" produced by the two voices sung together. This work has shown that a geometric correlation function, in log-frequency is best suited to uncover characteristic frequency ratios between different signals. The application to voice records has been selected not only for its simplicity to perform and reproduce, but also because it gives credit to the concept of frequency ratios in voiced sounds. This method is presently generalized to physiological signals recorded from different organs or tissues, such as the heart and the breath, extending the application of these ratio distributions.

Fig. 2

 2 Fig. 2 Grossmann analysing wavelet (log-normal in frequency). ψ (t ) is computed by inverse Fourier transform of ψQ ( f ) (Eq. (21)). (a) ψ Q (t ). (b) ψQ ( f ) in linear-frequency scale. (c) ψQ ( f ) in a (base 10) logarithmic frequency scale. (a,b,c) are computed for Q = 8. (d,e,f) Same as (a,b,c) for Q = 64. In (a) and (d) |ψ Q (t ) | (respectively ψ Q (t ) ) are plotted in thick (resp. light) black lines.

  with the Heaviside step function H, and we compare the wavelet transform analysis for the two quality factors Q = 8 and Q = 64.

Fig. 3

 3 Fig. 3 Analysis of a model signal defined as the sum of two sine functions s(t ) = sin(2π f 1 (t )t ) + sin(2π f 2 t ), with f 2 = 2 constant, and f 1 (t ) = H (-t ) + 3 2 H (t ) with the Heaviside step function. (a) Plot of the frequencies f 1 (t ) and f 2 (t ) in (base 2) logarithmic scale. (b) S (8) (log f , t ), computed for Q = 8. (c) S (64) (log f , t ), computed for Q = 64. (d) Temporal signal s(t ) in the time window [-40s, 40s]. (e) R[S(8) , S(8) ](log q, t ). (f) R[S (64) , S (64) ](log q, t ). R[S (Q) , S (Q) ](log q, t ) is defined in Eq. (27).

Fig. 5

 5 Fig. 5 Analysis of a model signal defined as the sum of 6 sine functions s(t ) = 6 1 sin(2π f i t ), with f i = i f 1 constant. (a) Fourier spectra of s: | ŝ |( f ) on which the two axes have been inverted. (b) S (8) (log f , t ), computed for Q = 8. (c) S (128) (log f , t ), computed for Q = 128. (d) S (Q) (log f , t = 0). (e) Temporal signal s(t ) in the time window [-6.5s, 6.5s]. (f) R[S(8) , S(8) ](log q, t ). (g) R[S(128) , S (128) ](log q, t ). (h) R[S (Q) , S (Q) ](log q, t = 0). In (d) and (h) the plots for Q = 8 (resp. 128) are coloured in black (resp. red). R[S (Q) , S (Q) ](log q, t ) is defined in Eq. (27).

Fig. 6

 6 Fig.[START_REF] Schnupp | Auditory Neuroscience: Making Sense of Sound[END_REF] Analysis of the #008 voice signal s(t ) from the voice-icar-federico-ii database[32]. (a) Zoom of the signal during 0.2 s. (b) | ŝ | plotted versus f in logarithmic scales (base 10 and base 2 respectively). (c) Local frequency f 1 (t ) computed from the detection of the extrema of larger prevalence from the signal (see (e)). (d) First return inverse of interpeak intervals f 1 (t i ) = 1/∆T i scatter plot (these large amplitude peaks are marked with black dots in (e)). (e) Zoom of the signal on the short interval (20 ms) showing the local maxima s P (t i ) (black dots) and minima s p (t i ) (black stars) which are used to compute both the local frequency: f 1 (t i ) = 1/∆T i and the amplitude A(t i ) of each larger amplitude peak: A(t i ) = s P (t i )s p (t i ). (f) | ŝ | (in base 10 log-scale) plotted versus f in Hz (linear scale). (g) Amplitude of the largest peaks A(t i ) versus time (see (e) for their detection). (h) First return peak amplitude (A(t i+1 ) vs A(t i )) scatter plot.

Fig. 7

 7 Fig. 7 Comparison of the time-frequency analysis of voice signals #008 (s 1 ) and #169 (s 2 ). (a) Zooms of s 1 and s 2 in a 0.1s window. (b,c) Associated time-frequency distributions (Eq. (26)) S (64) [s 1 ](log f , t ), and S (64) [s 2 ](log f , t ) computed with a Grossmann analysing wavelet and a quality factor Q = 64. The horizontal bands highlight the fundamental and harmonic frequencies. (d) Corresponding temporal averages of the frequency distributions reported in panels (b) black line and (c) red line. The ordinate of (d) (here the horizontal axes) is arbitrary and the frequency distributions are normalised.

Fig. 8

 8 Fig. 8 Ratio distributions of voice signals ((1) #008 and (2) #169) and their frequencymatched Dirac comb models (R 0 j ). (a) Correlations of the frequency distribution R 0i = R[S 0, j , S (Q)

Fig. 9

 9 Fig. 9 Comparing R[S 0 , S (Q) ](log q, t ) with S (Q) (log , t ) for voice signal #169. (a) Zoom of s(t ) in the [1.4s, 1.5s] interval. (b) Plot of a middle selection of 3s out of the 4.5s recorded voice signal. (c) Temporal average of the frequency distribution S (Q) t (log f ) computed with a Grossmann analysing wavelet with quality factor Q = 64. (d) Colour-coded map of the timefrequency distribution S (Q) (log f , t ). (e) Ratio distribution of the averaged frequency distribution: R[S 0 , S (Q)t ](log q). (f) Colour-coded map of the time-ratio distribution R[S 0 , S (Q) ](log q, t ).

2 Fig. 10 1 ( 2 (

 21012 Fig. 10 Comparison of the time-frequency analysis for two "normal" voice signals: s 1 is a sung vowel, s 2 is simply a maintained vowel. (a) Zooms of s 1 and s 2 in a 0.1s window. (b,c) Associated time-frequency distributions (Eq. (26)) S (Q) 1 (log f , t ), and S (Q) 2 (log f , t ) computed with a Grossmann analysing wavelet with quality factor Q = 64. The horizontal bands highlight the fundamental and harmonic frequencies. (d) Corresponding temporal averages of the frequency distributions reported in panels (b) black line and (c) red line. The ordinate of (d) (here the horizontal axes) is arbitrary and the frequency distributions are normalised.

Fig. 11

 11 Fig. 11 Comparing ratio distributions and sonance of the two voice signals of Fig. 10. (a) Plots of the normalised ratio distributions R 0 j = R[S 0,1 , S (Q) j ](log q) with j = 1, 2. For both signals we use the same Dirac comb model with the fundamental mode frequency of voice s 1 . (b) Plots of the normalised ratio distributions R 11 , R 22 , R 12 with R i j = R[ S (Q) i , S j ](log q) computed from the two voices, and R 00 = R[S 0,1 , S 0,1 ](log q) computed from the comb Dirac model. (c) Plot of the sonance (cross-and auto-) of the two voices with respect to the reference comb model ratio distribution for an arbitrary pitch transposition log(x): G[R i j ](log x) = R[R 00 , R i j ](log x) with i, j = 1, 2 green line, i, j = 1, 1 blue line and i, j = 2, 2 red line. (d) Sonance curves of Fig. 11(c) are corrected by subtracting their lower envelope.

Acknowledgements

This work has been supported by ANR under contract: ANR-18-CE45-0012-01, EURLight S&T funding. We are very indebted to L. Delmarre, E. Harte, and S. Polizzi for fruitful discussions.


Fabrice Ofridam

E Gagnière

Noureddine Lebaz
email: noureddine.lebaz@univ-lyon1.fr

Denis Mangin

Abdelhamid Elaïssari

Emilie Gagniere

Abdelhamid Elaissari

pH-sensitive polymethacrylate derivatives Eudragit E100 & L100 nanoparticles preparation in aqueous medium

come

Introduction & Context ❑

Polymers in general and polymethacrylate derivates are largely used in numerous applications including drug encapsulation and in controlled oral drug delivery [1]. ❑ pH-sensitive polymers are macromolecules that contain potential ionizable groups or hydrolyzable compounds and which react to pH variation by undergoing structural and property changes such as surface activity, chain structure, conformation, solubility and configuration [2] .

❑

The use of pH-sensitive polymers to prepare pH-sensitive capsules or particles has not been widely reported probably due to the few potential applications of pH sensitive environment.

❑

Various pH-sensitive carriers have been reported but the literature is free from any deep characterization of the used polymers. ❑ Charges Screening between SDS 4 and 8 mM (SDS/E100 charge ratio between 0,4 and 0,8) ❑Polymer precipitation as E100-SDS particles ❑Solubilization of complexes in micellar regimen ( > 8mM SDS) ❑ Charges Screening in micellar regimen between CTAB 4,5 and 6 mM (CTAB/L100 charge ratio between 0,3 and 0,4) ❑Polymer precipitation as L100-CTAB particles ❑Solubilization of complexes at high CTAB concentration ( > 8mM SDS) Anionic surfactant: E100/SDS Complexes at pH=2,5 [5] Cationic surfactant: L00/CTAB Complexes at pH=7 [5] Volume-based hydrodynamic particle size distribution of particle formed in presence of polyethylenimine (PEI) as a function of w ratio (mass ratio polyelectrolyte/Eudragit polymer) [4].

Conclusion & Perspectives

❑

  Solubilization and precipitation domains of oppositely charged polyelectrolytes Eudragit L100 & E100 were clearly established ❑ Eudragit L100 & E100 micrometric particles were prepared in buffer solutions ❑ Stabilization of submicrometric Eudragit particles in aqueous solution was able through formation of PECs (200 nm spherical Eudragit L00polyethylenimine particles) or PSCs (stable oppositely charged E100/SDS & L100/CTAB complexes dispersions). ❑ Eudragit particles preparation by acidification or basification in confined systems (using the two microemulsions process) and based on results obtained in aqueous solutions was not successful. ❑ It would be interesting to study the encapsulation efficiency of complexes prepared in aqueous solutions and their ability as potential nanocarriers. ❑ A better understanding of confined systems is needed at this point to continue the exploration of those systems for preparation of Eudragit L100 & E100 nanoparticles Objectives ❑ Investigation of the effect of different physicochemical parameters such as pH, salinity, polymer concentration, buffer concentration, and incubation time on the solubility and precipitation of two oppositely charged polymethyl methacrylate derivatives. ❑ Investigation of the effect of added secondary polymer or surfactant on the precipitation of both polymers ❑ Effect on size, size distribution, zeta potential and morphology of the particles. Solubilization of Eudragit E100 in HCl medium as a function of polymer amount and time [3] Solubilization of Eudragit L100 in NaOH medium as a function of polymer amount and time [3] Solubilization studies in buffer solutions confirms that Eudragit L100 was soluble above pH = 6.50 and Eudragit E100 up to pH = 5.20 By pH titration, it was found that Eudragit L100 contains around 4.8 nmol/mol of carboxylic groups and Eudragit E100 contains around 8.5 nmol/mol of ammonium groups. L100 and E100 Eudragit polymers act respectively as polyacid and poly base In addition to the initial pH of the medium, the solubility conditions of the two investigated polymers are closely related to the polymer content and its incubation time. Grey bars on the 3D graphics refer to solubilized samples. Polymers precipitation: Particle preparation and characterization ❑ Precipitated particles were micrometric in size with mean size between 30 and 50 μm, and large size distributions. ❑ TEM images show that all obtained dispersions are spherical polymer-based particles and highly polydisperse, which agrees with diameters deduced from diffraction light analysis. ❑ Stabilization of submicronic Eudragit particles was made through complexation either as Polyelectrolytes complexes (PECs) or Polyelectrolyte-Surfactant complexes (PSCs) ❑Narrow size distribution around 200 nm was obtained when Eudragit L100 is precipitated in presence of polyethylenimine (PEI) with particles zeta potential found to be +45mV. No aggregation and no sedimentation occurred. ❑Results suggest the formation of well defined spherical and monodisperse Eudragit particles stabilized by polyethylenimine.