Pierre Margerit

Arthur Lebée

Jean-François Caron

Xavier Boutillon

High Resolution Wavenumber Analysis (HRWA) for the Mechanical Characterisation of Viscoelastic Beams

Keywords: Structure identification, Non-destructive testing, Signal Processing, Wavenumber extraction

The High-Resolution Wavenumber Analysis (HRWA) is presented. It identifies complex wavenumbers and amplitudes of waves composing the harmonic response of a beam. Based on the frequency dependence of these wavenumbers, experimental dispersion branches corresponding to various beam motions (e.g bending, torsion) can be retrieved. The HRWA method is compared to the Mc Daniel and the Inverse Wave Correlation (IWC ) methods. It overcomes some drawbacks of these methods: the wavenumber resolution is enhanced. Also, the wavenumber search problem is expressed as a linear problem, making the method computationally efficient. A number of wavenumbers can be identified automatically, thanks to a statistical criterion. First, the noise sensitivity of each method is investigated in the basis of synthesised measurements. For this criterion, the HRWA and Mc Daniel method performances are close and much better than IWC. Moreover, the HRWA is five to twenty times faster to compute than other methods, depending on the mesh size. Second, an experimental case is presented where bending and torsion waves are identified, yielding an apparent viscoelastic Young and shear moduli on a wide-frequency range.

Introduction

With the apparition of full-field contactless measurement devices (e.g. scanning doppler laser vibrometers or high-speed cameras), displacement or velocity measurements can be performed on fine meshes, with a large number of points. Thanks to this major improvement, non-destructive structural characterisation methods have been developed [START_REF] Avril | Overview of identification methods of mechanical parameters based on fullfield measurements[END_REF][START_REF] Pezerat | Identification of vibration sources[END_REF][START_REF] Leclère | Vibration source identification using corrected finite difference schemes[END_REF][START_REF] Wassereau | Characterization of materials and flaw detection using force analysis technique[END_REF]. Some of them focus on the identification of structural waves in the medium frequency domain [START_REF] Berthaut | Contribution à l'identification large bande des structures anisotropes[END_REF][START_REF] Berthaut | K-space identification of apparent structural behaviour[END_REF][START_REF] Mcdaniel | Estimation of structural wave parameters from spatially sparse response measurementsle[END_REF][START_REF] Palan | Characterization of an experimental wavenumber fitting method for loss factor estimation using a viscoelastically damped structure[END_REF][START_REF] Rak | Identification of structural loss factor from spatially distributed measurements on beams with viscoelastic layer[END_REF]. These methods aim at filling the gap between low-frequency methods, (Oberst [START_REF] Oberst | Über die dämpfung der biegeschwingungen dünner bleche durch fest haftende beläge. {On the damping of bending vibrations of thin sheets by adherent coverings[END_REF][START_REF] Oberst | Über die dämpfung der biegeschwingungen dünner bleche durch fest haftende beläge II. {On the damping of bending vibrations of thin sheets by adherent coverings II[END_REF], modal analysis [START_REF] Berthelot | Damping analysis of unidirectional glass and Kevlar fibre composites[END_REF][START_REF] Berthelot | Damping analysis of laminated beams and plates using the Ritz method[END_REF][START_REF] Berthelot | Damping analysis of composite materials and structures[END_REF][START_REF] Ege | High-resolution modal analysis[END_REF][START_REF] Rébillat | Measurement of relevant elastic and damping material properties in sandwich thick plates[END_REF] or Dynamical Mechanical Analysis (DMA) [START_REF] Jaouen | Elastic and damping characterizations of acoustical porous materials: Available experimental methods and applications to a melamine foam[END_REF]) and ultrasonic testing [START_REF] Hosten | Inhomogeneous wave generation and propagation in lossy anisotropic solids. Application to the characterization of viscoelastic composite materials[END_REF][START_REF] Audoin | Measurement of stiffness coefficients of anisotropic materials from pointlike generation and detection of acoustic waves[END_REF][START_REF] Kersemans | The quasi-harmonic ultrasonic polar scan for material characterization: Experiment and numerical modeling[END_REF].

The present work focuses on applications to beam, considered as waveguides. Waves travelling in these structures are representative of the beam section motions: for example bending, twist or longitudinal motion. In these unidimensionnal structures, the wavenumbers are characteristic of the local structural behaviour: far from singularities, boundary conditions and sources have an influence on the wave amplitudes only. These wavenumbers are given by reduced beam models (e.g. Euler or Timoshenko models) or more elaborated schemes like Wave Finite elements [START_REF] Duhamel | Finite element analysis of the vibrations of waveguides and periodic structures[END_REF]. The experimentally extracted wavenumbers can be used to identify apparent viscoelastic parameters of reduced models via an inverse problem [START_REF] Rak | Identification of structural loss factor from spatially distributed measurements on beams with viscoelastic layer[END_REF].

The existing Mc Daniel [START_REF] Mcdaniel | Estimation of structural wave parameters from spatially sparse response measurementsle[END_REF] and IWC [START_REF] Berthaut | Contribution à l'identification large bande des structures anisotropes[END_REF][START_REF] Berthaut | K-space identification of apparent structural behaviour[END_REF] (Inhomogeneous Wave Correlation) methods are two candidate for the extraction of waves in the harmonic response of beams. However, these methods have some drawbacks (i) in low frequency, the Fourier based IWC method suffers from resolution limitations; (ii) in practice, the beam response is composed of multiple wave types (e.g. twist, compression, shear) and both existing methods fail to separate the different wave contributions, as they postulate the number of waves present in the signal; (iii) the formulation of these methods leads to a computationally expensive non-linear problem that has to be solved, with a complex wavenumber as parameter.

The aim of the High Resolution Wavenumber Analysis (HRWA) presented in this paper is to overcome the limitations of the existing wavenumber extraction techniques. It makes use of the subspace-based identification algorithm ESPRIT [START_REF] Roy | ESPRIT-Estimation of Signal Paramters via Rotational Invariance Techniques[END_REF] (Estimation of Signal Parameters via Rotational Invariance Techniques). Subspace-based methods are widely-used in linear system identification, using for example the state-variable framework [START_REF] Ho | Effective construction of linear state-variable models from input/output functions 1 )[END_REF][START_REF] Van Overschee | Subspace identification for linear systems: Theory-Implementation-Applications[END_REF]]. Another example is the ERA [START_REF] Juang | An Eigensystem Realization Algorithm for Modal Parameter Identification and Model Reduction[END_REF] (Eigenvalue Realisation Algorithm), which is devoted to the identification of the modal parameters of a measured system. Thanks to the use of the ESPRIT algorithm, some limitations of the Mc Daniel and IWC methods are overcome: (i) the algorithm resolution is high, as it uses a recurrence property of the signal to identify the wave parameters; (ii) by using the subspace decomposition, the number of waves contained in the signal can be estimated automatically with the ESTER criterion [START_REF] Badeau | A new perturbation analysis for signal enumeration in rotational invariance techniques[END_REF]; (iii) the complex wavenumbers are the solution of an optimisation-free problem thus the computational cost is lightened.

With the HRWA, a discrete number of complex wavenumbers is identified in the harmonic response of a beam with a high resolution. Based on the dependence of the extracted wavenumbers on frequency, experimental dispersion branches are retrieved. From these branches are identified the beam viscoelastic properties. Thanks to the high resolution aspects of the HRWA, the low-frequency limit of wavenumber extraction is lowered. Also, with the automated identification of multiple wavenumbers, strain mechanisms can be well separated, extending the upper frequency limit.

The paper is organised as follows. Firstly, theoretical wavenumbers are derived from the Euler and the Timoshenko beam models. The ability to identify beam properties from the wavenumbers is illustrated. Then, a common framework of wavenumber extraction methods is given, and the existing IWC method and McDaniel method reformulated in this framework. Secondly, the HRWA is developed in details. A summary of the algorithm is given, that allows to identify viscoelastic properties of the beam. Thirdly, a numerical study based on synthesised harmonic responses of an Euler beam in bending motion only is developed, where the Mc Daniel method is taken as reference. The Mc Daniel method, IWC method and the proposed HRWA are compared in terms of sensitivity to noise ratio and computation time. Finally, an experimental result is given on the simultaneous identification of frequency-dependent apparent viscoelastic Young and shear modulus, for both beam models.

Natural wavenumbers in a beam

All along the paper, the beam section is considered as homogeneous and made of a linear viscoelastic isotropic material with density ρ, Young modulus E and shear modulus G. However, the present method is applicable to higher-order beam models and more complex material configurations. The axis and the geometry of the beam are presented in Fig. 1 (x being the beam's direction). Not considering the bending motion along y, the beam neutral axis is assumed to remain in the (O, x, z) plane.

The Euler model

For the sake of simplicity, the Euler beam model is used for the numerical case considered in this paper. The linearised displacement field u is given by:

u(x, y, z, t) = V (x, t) e x + W (x, t) e z -z W (x, t) e x + Θ(x, t) (y e z -z e z ) (1)
where V and W are respectively the longitudinal and transverse displacements, Θ is the rotation angle of the beam section with respect to the x axis, and • denotes the partial derivative with respect to x.

Along the paper, the harmonic dependence on time of real physical quantities is accounted for by making them complex according to the e iωt convention (where i denotes the imaginary unit). The local harmonic response of the beam at the angular frequency ω (except on sources or at boundaries locations) obeys the following uncoupled linear homogeneous equations:

E V = -ω 2 ρ V (2a) I b E W = -ω 2 M W (2b) I t G Θ = -ω 2 J Θ (2c) with A = b h, I b = b h 3 /12, J = ρ (b h 3 + h b 3 )/12, I t = κ J/ρ , M = ρ A
and κ being a correction factor for twist inertia momentum, related to the shape of the beam's section. By taking exponential solutions as function of x, the displacement components are expressed as a sum of the solutions of the previous equations:

V (x) = v 1 e ikvx +v 2 e -ikvx (3a) W (x) = w 1 e ik b x +w 2 e -ik b x +w 3 e k b x +w 4 e -k b x (3b) Θ(x) = θ 1 e iktx +θ 2 e -iktx (3c)
where the wave amplitudes v r , w r and θ r are given by boundary conditions. According to the e iωt convention specified above, the e -ikx waves are forward-propagating waves. The three following dispersion laws are obtained:

k 2 v = ω 2 ρ E , k 4 b = ω 2 M EI b and k 2 t = ω 2 J GI t (4)
It is of common practice to account for the viscoelastic behaviour of the material by taking complex moduli E = E(1 + i η E ) and G = G(1 + i η G ) with positive (and small) imaginary parts, according to the e iωt convention. The wavenumbers become complex too:

k = k(1 -i γ) (5)
generically, with k = (k). Specifically, the three groups of wavenumbers given by Eq. ( 4) become:

k v = ±k v (1 -i γ v ) (6a) k b = ±k b (1 -i γ b ) and k b = ±ik b (1 -i γ b ) (6b) k t = ±k t (1 -i γ t ) (6c)
where these expressions are physical (waves attenuating along their direction of propagation) only if the generic coefficient γ defined above is positive. The second group of bending waves

(k b = ±ik b (1 -i • γ b )
) correspond to fast decaying and slowly oscillating waves. Finally, Young and shear moduli can be identified from bending and twist wavenumbers composing the harmonic response of the beam. Considering the Euler model, the Young and shear moduli are given:

E eul (ω) = M I b ω 2 k 4 b and G(ω) = J I t ω 2 k 2 t ( 7
)
The coefficients that are obtained are apparent Young and shear moduli, in the sense that they refer to an underlying theory of the Euler beam behaviour and not directly to elasticity. Consequently, they may vary in frequency.

The Timoshenko model

One can define one equivalent slenderness ratio µ as the bending wavelength divided by the beam thickness:

µ(ω) = λ b (ω) h = 2π k b (ω)h (8)
When the frequency increases, µ decreases. As a consequence, out-of-plane shear and rotary inertia effects are no longer negligible and the beam behaves as a thick beam. The Timoshenko model aims at taking into account these effects. The governing equation for bending motion in the harmonic regime (Eq. (2b)) can be found in [START_REF] Bhaskar | Elastic waves in Timoshenko beams: the 'lost and found' of an eigenmode[END_REF]:

I b E W -ω 2 M W + ω 2 ρI b 1 + E ξG W + ω 4 ρ 2 I b ξG = 0 (9)
with ξ = π 2 /12 the shear correction factor. As in Eq. ( 4), the bending wave dispersion law is derived:

k 4 b I b E -ω 2 M -ω 2 k 2 b ρI b 1 + E ξG + ω 4 ρ 2 I b ξG = 0 (10)
As in the preceding case, this fourth order equation in k gives four solutions, of which two are almost real and two are almost imaginary. Finally, an apparent Young modulus E tim that takes both out-of-plane shear effects and rotary inertia into account can be identified from the bending wavenumber k b and the shear moduli G(ω) (Eq. ( 7)):

E tim (ω) = ω 2 M + ω 2 k 2 b ρI b -ω 4 ρ 2 I b ξG(ω) k 4 b I b -ω 2 k 2 b ρ I b ξG(ω) = E eul (ω)   1 1 -ω 2 k 2 b ρ ξG(ω) + k 2 b I b A   (11)
The next section presents the signal model of a measurement performed on the beam, which parameters are identified by means of the proposed HRWA and of two existing methods suited for wavenumber identification: the Mc Daniel method [START_REF] Mcdaniel | Estimation of structural wave parameters from spatially sparse response measurementsle[END_REF] and the Inhomogeneous Wave Correlation [START_REF] Berthaut | Contribution à l'identification large bande des structures anisotropes[END_REF].

Signal model

The harmonic response of a beam is measured on a domain away from sources and boundaries, along a regularly ∆-spaced mesh x of N points, aligned with the beam neutral axis:

x u = x 1 + (u -1)∆ , u ∈ [[ 1, N ]] (12)
A typical setup is schematised in Fig. 2. For a beam structure vibrating in the linear regime, the measurement of a component of u(x, y, z) along the mesh x is expected to provide the spatially discrete signal s formed by a sum p = u i (x, y 0 , z 0 ) of R damped exponentials and noise n : where a r are the amplitudes corresponding to the poles

s u = u i (x u , y 0 , z 0 ) + n u = p u + n u = R r=1 a r z u r + n u (13) x 0 L x 1 x 2 • • • x N ∆ F x F
z r = e ikr∆ , k r = k r (1 -i • γ r ).
The p part of the signal can be expressed as:

p = V N (z R ) • a R ( 14
)
where

a R = [a 1 , . . . , a R ] , z R = [z 1 , . . . , z R ] and V N (z R ) is the Vandermonde matrix: V N (z R ) =       1 • • • 1 z 1 • • • z R . . . . . . . . . z N -1 1 • • • z N -1 R       (15)
The main problem is to find the poles z R . Once the poles are known, the complex amplitudes a R can be estimated in the measured noisy signal s in the least-square sense:

a R = (V * N (z R ) • V N (z R )) -1 (V * N (z R ) • s) (16)
where • * denotes complex conjugate or hermitian transpose.

Two existing wavenumber identification methods

As well as the HRWA method presented in this paper, the Mc Daniel and IWC methods make use of a Vandermonde matrix, but with different R and N values.

One pole : Inverse Wave Correlation

The Inverse Wave Correlation method [START_REF] Berthaut | Contribution à l'identification large bande des structures anisotropes[END_REF][START_REF] Berthaut | K-space identification of apparent structural behaviour[END_REF] aims at finding the wavenumber k wich maximises the normalised scalar product C IW C (k) between the Vandermonde vector v N (k) = V N (e ik∆ ) and the noisy signal s:

C IWC (k) = s • v N (k) s v N (k) (17)
One can notice that the IWC correlation coefficient C IW C (k) can be seen as a normalised z-transform on the complex wavenumber k. Therefore, it suffers the same resolution limitations as any method based on a Fourier transform.

Bending waves's four poles: the Mc Daniel method

Mc Daniel suggested [START_REF] Mcdaniel | Estimation of structural wave parameters from spatially sparse response measurementsle[END_REF][START_REF] Palan | Characterization of an experimental wavenumber fitting method for loss factor estimation using a viscoelastically damped structure[END_REF] to identify simultaneously the four waves derived from the bending operator (Eq. (2b), Eq. (3b) and Eq. ( 4b)). The four poles are imposed to be z 4 (k) = [e ik∆ , e -ik∆ , e k∆ , e -k∆ ], where k only has to be determined. A signal reconstruction s(k) is made:

s(k) = V N (z 4 (k)) • a 4 (18)
with the amplitudes a 4 estimated by Eq. ( 16). If only bending waves are present in the measured signal s, its normalised correlation C McDan (k) with the reconstructed signal s(k):

C McDan (k) = s • s(k) s s(k) (19)
is maximum when k matches the natural bending wavenumber.

Both methods are based on the maximisation of a correlation coefficient (C IWC or C McDan ). This is a non-linear problem on two unknowns: the real and imaginary parts of k. For the comparisons given in the present paper, the procedure proposed by Rak [START_REF] Rak | Identification of structural loss factor from spatially distributed measurements on beams with viscoelastic layer[END_REF] was adopted, using a Nelder-Simplex research algorithm.

High Resolution Wavenumber Analysis (HRWA) method: implementation for the characterisation of beams

The HRWA method makes use of the ESPRIT algorithm [START_REF] Roy | ESPRIT-Estimation of Signal Paramters via Rotational Invariance Techniques[END_REF] (Estimation of Signal Parameters via Rotational Invariance Techniques) to extract wavenumbers from the harmonic response of a beam. It is inspired from the pioneer works of Prony [START_REF] Prony | Essai experimental et analytique. {experimental and analytical assay[END_REF], who stated that signals composed of complex poles follows a recurrence relation. This led, at the end of the twentieth century, to the development of a family of identification methods such as Matrix Pencil [START_REF] Hua | Matrix Pencil Method for Estimating Parameters of Exponentially Damped/Undamped Sinusoids in Noise[END_REF], Pisarenko method [START_REF] Pisarenko | The Retrieval of Harmonies from a Covariance Function[END_REF], MUSIC algorithm [START_REF] Schmidt | A Signal Subspace Approach to Multiple Emitter Location and Spectral Estimation[END_REF] and ESPRIT algortihm. The main goal of these methods consists in identifying, in a measured signal, complex poles with their associated amplitudes. These methods are widely used in array processing and radar applications [START_REF] Shahbazpanahi | Distributed source localization using ESPRIT algorithm[END_REF][START_REF] Rouquette | Estimation of frequencies and damping factors by two-dimensional ESPRIT type methods[END_REF], musical applications [START_REF] Emiya | A parametric method for pitch estimation of piano tones[END_REF][START_REF] Ege | High-resolution modal analysis[END_REF][START_REF] Badeau | Eds Parametric Modeling And Tracking Of Audio Signals[END_REF] and high-resolution spectral analysis [START_REF] Badeau | High-resolution spectral analysis of mixtures of complex exponentials modulated by polynomials[END_REF]. Some applications have been developed for the impact localisation on plates [START_REF] Qiu | Impact imaging of aircraft composite structure based on a modelindependent spatial-wavenumber filter[END_REF].

A Prony method has been already used in the past for the retrieval of dispersion relations in cylindrical shells [START_REF] Plona | Axisymmetric wave propagation in fluid-loaded cylindrical shells. II: Theory versus experiment[END_REF] [START_REF] Vollmann | High-resolution analysis of the complex wave spectrum in a cylindrical shell containing a viscoelastic medium .2. Experimental results versus theory[END_REF]. However, it was limited by the sensibility to noise, the unknown number of travelling waves and the number of measurement points required to retrieve evanescent waves [START_REF] Grosh | Complex wave-number decomposition of structural vibrations[END_REF]. The ESPRIT algorithm, used here, applies subspace decomposition to improve the resistance to noise. In addition, the ESTER criterion [START_REF] Badeau | A new perturbation analysis for signal enumeration in rotational invariance techniques[END_REF] (ESTimation or ERror) is used to estimate the signal order (number of poles contained in the signal). Finally, the procedure is here applied to full-field measurements, which allows to make measurement over fine meshes and overcome the limitations due to a limited number of points.

The ESPRIT algorithm 5.1.1. Signal and noise subspaces

The first step of the algorithm consists in the decomposition of the signal between signal and noise subspaces. First, a Hankel matrix H is formed with the measured signal s of length N :

H =       s 1 s 2 • • • s N -K s 2 s 3 • • • s N -K+1 . . . . . . . . . s K • • • s N       (20)
where the integer parameter K corresponds to the sum of the dimensions of signal and noise subspaces. The HRWA method makes use of the autocovariance matrix R ss :

R ss = HH * = W * DW [START_REF] Duhamel | Finite element analysis of the vibrations of waveguides and periodic structures[END_REF] which eigenvectors W spans the same subspace as the Hankel matrix singular vectors, and has the advantage to be an asymptotically non-biased estimator of the signal autocovariance, in presence of white gaussian noise. The eigenvectors matrix W spans the entire subspace of the noisy signal s. This total subspace can be decomposed into signal subspace W p and noise subspace W n , with W p built with the eigenvectors corresponding to the R dominant eigenvalues (with the number of poles R supposed known here):

W = [W p ] K×R [W n ] K×(K-R) K×K ( 22
)
It can be seen, with this last equation, that the parameter K has to be larger than R. Also, the dimension of the noise subspace is given by this parameter. In fact, the (K-R) additional poles contained in the noise subspace permits to sample an hypothetical correlated part of the noise [START_REF] Kumaresan | Estimating the parameters of exponentially damped sinusoids and pole-zero modeling in noise, Acoustics, Speech and Signal Processing[END_REF], to separate it from the signal pole extraction. Numerous studies investigates this parameter influence on the errors of pole estimation [START_REF] Laroche | The use of the matrix pencil method for the spectrum analysis of musical signals[END_REF][START_REF] Badeau | Méthodes à haute résolution pour l'estimation et le suivi de sinusoïdes modulées. application aux signaux de musique {High-resolution methods for the estimation and the tracking of modulated sinusoids[END_REF], and a common strategy consists in choosing this parameter so that K ≈ N/2.

Rotational Invariance

Thanks to the regularity of the uniform measurement mesh x of Eq. ( 12), a rotational invariance property is expressed via two Vandermonde matrices, related to two subparts of the signal:

V ↓ (N -1) (z R ) = V ↑ (N -1) (z R ) • Z ( 23
)
with Z = diag(z R ) and the two matrices V ↑ (N -1) and V ↓ (N -1) respectively corresponding to the (N -1) first and last samples of the signal. The ESPRIT algorithm strategy consists in estimating the signal poles z R via the signal subspace matrix W p . As this last matrix and the Vandermonde matrix spans close subspaces, they are related by a transfer matrix T : V N (z R ) = W p T . The rotational invariance property is then expressed as a function of the signal subspace matrix:

W ↓ p = W ↑ p F , with F = T Z T -1 and: W ↑ p = I (K-1) 0 (K-1)×1 W p W ↓ p = 0 (K-1)×1 I (K-1) W p ( 24
)
Consequently, the matrix F is estimated in the Least-Square sense:

F = (W ↑ p ) -1 W ↓ p ( 25
)
Finally, the poles z R are extracted from the diagonalisation of F .

Estimation of the number of poles

In the previous steps, the number of poles R, or signal order, has been supposed known. However, the number of dominant waves in the signal is unknown in most of the practical applications. As a wrong estimated signal order would drastically decrease the poles extraction accuracy, a signal order estimation method is needed. A number of signal order selection criterion have been proposed [START_REF] Stoica | Model-Order Selection: A review of information criterion rules[END_REF], as it is a common issue encountered in subspacebased identification methods. The HRWA method makes use of the ESTER criterion [START_REF] Badeau | A new perturbation analysis for signal enumeration in rotational invariance techniques[END_REF] (ESTimation or ERror), which is particularly adapted to the ESPRIT algorithm.

The estimation of matrix F (Eq. ( 25)) is sensitive to the assumed signal order R used to extract the signal subspace (Eq. ( 22)). The ESTER criterion consists in finding the signal order R that minimises the F estimation residual, so that :

R = min r ∈ [[ r min ,rmax ]] W ↑ p (r)F (r) -W ↓ p (r) 2 (26)
With this criterion, the signal order R can be estimated automatically in a given range of signal orders r ∈ [[ r min , r max ]].

HRWA Implementation

For applying the HRWA to the characterisation of beams, any type of excitation may be used (impulse, random noise or sweep), with any type of actuator (e.g. shaker at x F , exerting a point force F on the beam). The method has no requirement with regard to boundary conditions, provided that they are not applied too close to the measured zone. The displacement, velocity or acceleration S(x, t) is measured along a regularly spaced mesh x of Eq. ( 12), as in Fig. 2. A collection of harmonic responses is obtained by computing the Fourier Transform S(x, ω) of the measured signal over the time dimension for a given number of individual frequencies ω i . The HRWA procedure consists in applying the following steps to each obtained harmonic response s(ω i ) = S(x, ω i ):

1. A Hankel matrix H is built, based on the array s(ω i ) (see Eq. ( 20)).

2. The Covariance Matrix R ss is computed and diagonalised yielding the matrix of eigenvectors W (see Eq. ( 21)). 3. The ESTER method is applied to estimate the signal order, in other words the number of detectable waves in the (noisy) signal. Then, for each r ∈ [[ r min , r max ]] : (a) the r eigenvectors corresponding to the r dominant eigenvalues are extracted to form the approximated signal subspace matrix W p (Eq. ( 22)) (b) the truncated signal subspace matrices W ↑ p and W ↓ p are built (Eq. ( 24)) (c) the Least-Squares estimation of F is computed (Eq. ( 25)) (d) the ESTER criterion as function of r is evaluated (Eq. ( 26)) Finally, the signal order R is chosen to minimise the ESTER criterion. 4. The eigenvalues of the matrix F (see Eq. ( 25)) are computed, yielding the wavenumbers k r of all the waves which can be detected with the HRWA at the angular frequency ω i :

k r (ω i ) = ln(i • z r (ω i )) ∆ (27)
5. If a signal reconstruction is needed, the Vandermonde matrix V N (z R ) is computed (Eq. ( 15)), and the complex amplitudes a R are estimated (Eq. ( 16)).

Comparison of methods for a simulated cantilever bending-only beam

This section is devoted to a comparison between the Mc Daniel method, the IWC method and the HRWA, regarding the sensitivity to noise and the computation time. For this purpose, a simple case is investigated numerically. Using synthesised data for the comparison makes possible the exploration of various parameters (e.g. mesh size, material characteristics, etc.). A thin aluminium cantilever Euler beam of length L is harmonically excited at the pulsation ω with a force F = 1N at (x F , y F , z F ) = (L, 0, h/2) (Fig. 3), in order to be submitted to bending only. In this situation, the signal model of the Mc Daniel method is exact: the harmonic response of the beam is composed of four bending waves, with wavenumbers k b solutions of the fourth root of Eq. (4b). It is thus the reference method for this case study.

The characteristics of the beam are summarised in Tab. 1. The beam is excited between 10 Hz and 5 kHz. Eleven resonance peaks appear in this bandwidth (see Fig. 4, top frame). The first modal frequency is around 14 Hz. The response of the beam is derived analytically since the four complex amplitudes w j in Eq. ( 3) are determined by the boundary conditions. The response yields the virtual measurements on the mesh x of Eq. ( 12) along a sub-part of the beam (x u ∈ [0.1, 0.9] × L). The beam response is computed at 100 equally-spaced frequencies in the desired bandwidth, and for 100 regularly spaced points over the sub-part of the beam.

x 0 L x 1 x 2 • • • x N ∆ F

Noise Sensitivity

To study the sensitivity to noise of each wavenumber identification method, a white gaussian noise is added to the signal s. The amplitude of this noise is determined by a given Signal to Noise Ratio (SNR). For each of the 100 SNR values ranging from 10 0 to 10 4 and each of the 100 computed beam harmonic responses, 100 realisations of the noisy signal are computed. Altogether, the estimation of the bending wavenumber is done one million times for each method. Based on these 100 realisations of the noise the mean and variance of the wavenumber estimation error are computed, as a function of frequency and SNR.

The Mc Daniel method and IWC method identify only one wavenumber k. By contrast, the HRWA identifies a varying number of wavenumbers. Consequently a valid wavenumber selection strategy has to be chosen to perform the comparison. All along this numerical investigation, the maximum signal order r max of the HRWA procedure is given by the equilibrium (Eq. ( 3)). The algorithm is then free to identify from one to eight waves of different wavenumbers (R ∈ [ [START_REF] Avril | Overview of identification methods of mechanical parameters based on fullfield measurements[END_REF][START_REF] Palan | Characterization of an experimental wavenumber fitting method for loss factor estimation using a viscoelastically damped structure[END_REF]]). The choice is automated with the ESTER criterion. The chosen strategy, which was observed to work the best among several options, consists in isolating the P propagative wavenumbers k p so that 0 < γ p < 1. The mean of these wavenumbers is then computed as follows:

k = 1 P   P p=1 Re(k p ) 2 -i P p=1 Im(k p ) 2   ( 28
)
the square power allowing to take into account both forward and backward propagating waves. This strategy eliminates the identified evanescent (fast decreasing) waves, the identification of which is sensitive to noise. In Fig. 4 and Fig. 5, average relative errors on the identification of the real part k of the wavenumber and the spatial decay γ are represented for the IWC method, the Mc Daniel method and the HRWA. In Fig. 4, three contour lines are plotted, corresponding on three levels of performance (relative error of 10%, 1% and 0.1%). The wavenumber identification turns out to be sensitive to the number of spatial periods contained in the measurement mesh. As the wavelength depends on frequency, it has been chosen to represent the relative error on the complex bending wavenumber estimation as a function of frequency and of Signal to Noise Ratio.

Regarding the real part of the bending wavenumber, both HRWA and Mc Daniel method give an estimation with less than 1% of relative error, above 500 Hz and for small signal to noise ratios (SNR > 3). The Mc Daniel method provides a relative error of 0.1% in the estimation, up to a 10 2 SNR, whereas the HRWA delivers such an estimation above 3kHz. The HRWA and the Mc Daniel methods display comparable abilities to provide an estimation of the spatial decay with less than 10% of relative error. For high signal to noise ratios (SNR> 10 3 ), the Mc Daniel method can give an estimation of the spatial decay with less than 1% of relative error, above 2.5 kHz, whereas this frequency bound is higher for the HRWA (4kHz).

In contrast, the estimate of the IWC method depends on the frequency. This behaviour is related to its resolution limitations. Except for few high frequencies, no estimation of the real part of the wavenumber with IWC can be made with less than 0.1% or relative error. Moreover, the estimation of the spatial decay γ is not possible to the IWC method in this SNR and frequency ranges. These results agree with those of Rak [START_REF] Rak | Identification of structural loss factor from spatially distributed measurements on beams with viscoelastic layer[END_REF].

The HRWA method exhibits a frequency-dependent behaviour for the fine estimation of real part of wavenumber (0.1% of relative error) and spatial decay (1% of relative error): both real and imaginary parts of the wavenumber are better identified at modal frequencies. In order to explain this particular behaviour, results of the study for two given frequencies are plotted in Fig. 5. It shows the results of complex bending wavenumber identification as a function of SNR, for the three methods. First, Fig. 5a shows the results at 1000 Hz, between two modal frequencies (see Fig. 4, top frame, marker a ). For this particular frequency, the HRWA converges shortly to a first plateau at 0.5% of relative error in k identification, for SNR higher than 10 1 . Then, one can observe a jump at SNR of 2000, HRWA finally reaching a very fine estimation of the real part k of the wavenumber, comparable to the Mc Daniel method. The level and the extension of the plateau are due to the variable estimation of the number of waves in the signal by the ESTER criterion, an effect which appears to be stronger between modal frequencies. The spatial decay estimation is affected by this plateau effect too, the relative error converging to a value of 3%, below a SNR of 2000. By comparison, the error of the estimation of k with the IWC method reaches a plateau too, which is higher (0.9% of relative error) and related to its resolution limitation. Furthermore, the Mc Daniel method shows a uniform convergence rate, with a relative error fifty times smaller than the noise ratio on k (∆k/k ≈ (50 SNR) -1 ) and twenty times larger on the spatial decay (∆γ/γ ≈ 20 (SNR) -1 ). In Fig. 5b, the same data has been plotted, but for a higher frequency of 4400Hz which coincides with a modal frequency (see Fig. 4, top frame, marker b ), where the three methods provides a good estimation of the real part of the wavenumber. Here, the plateau reached by the HRWA error is much lower (≈ 6 × 10 -5 of relative error) than in the previous case. The HRWA, before this plateau, shows the same convergence rate as Mc Daniel method, both for k and γ estimation. It has to be noted that IWC method, in this case, estimates the spatial decay with approximately 15% of relative error.

The plateau effect of HRWA shown previously is due to the possibility given to the ESPRIT algorithm, to select the number of waves in the signal (thanks to the ESTER criterion), and by making no assumption on a relation between each wavenumber. This is the main difference with Mc Daniel method's signal model, which confers to this last method its stability against the noise level, when only bending waves are present. However, this plateau appears below 1% of relative error on the real part of the wavenumber and 10% of relative error on the spatial decay, which is are good estimations, altogether.

Comparison of Computation Times

In addition to noise sensitivity, the computation time of each method has been investigated. As these methods extract wavenumbers from one single frequency component, the computation time is independent from frequency. It is only dependent on the size N of the measurement mesh. From 10 to 1000 measurement points have been exported from the synthesised forced response's shape. The computation time has been averaged on 20 frequency points and 10 iterations. The results of the study are plotted in Fig. 6.

This figure illustrates the advantages of HRWA in terms of computational cost: it is, depending on the mesh size, five to twenty times faster than the two other methods. This makes wavenumber analysis suitable for real-time applications, such as structural health monitoring. The computational cost of Mc Daniel and IWC methods is mainly caused by the amount of iterations needed in the two dimensional non-convex minimisation problem in k.

Number of Measurement Points

Experimental Results

An experimental application is performed on a simple rigid PVC (polyvinyl chloride) beam. The beam's section is square of side 15mm, and its length of 1m. The velocity is measured with a Doppler laser vibrometer along a regular mesh of 500 points over the top surface of the beam. Consequently, the mesh is shifted from the beam neutral axis (y = 0 in Eq. ( 1)). Hence both bending and twist motion contributes in the measured velocity. A shaker is fixed at one end of the beam ((x F , y F , z F ) = (0, -b/2, -h/2), see Fig. 1), whereas the other end is supported by a silicon block.

A first measurement is performed with an impulse excitation, achieved by sending an electrical pulse to the shaker. A second measurement is made with a random noise excitation filtered between 5 kHz and 22 kHz. In both cases, the frequency response is estimated from the average of 20 realisations, respectively below 6 kHz and between 5 and 22 kHz. For both experiments, a sampling frequency of 50 kHz was used.

The HRWA is applied to the computed frequency transfer estimator H1 between the measured velocity and the electrical excitation signal (no force measurement). As the measurement mesh is shifted from the neutral axis, it contains both the bending and twist motion's contributions. Fig. 7a shows the HRWA results, for both excitation type. It shows that the HRWA is able to identify both bending (blue) and twist (red) wavenumbers, at each frequency. It can be observed that both the real and imaginary parts of the identified wavenumbers are independent of the excitation type: the wavenumbers extracted in the overlapping frequency band (between 5Hz and 6kHz) are identical. By combining different excitations, a wavenumber extraction can be achieved in a very wide frequency domain. The Loss Factor 2 (%) )) and the Timoshenko model (E tim (green) from Eq. ( 11))).

2 E eul 2 G 2 E tim (b)
minimum identified wavelength λ = 2π/k is approximately 3.7 cm, which means that at least 18 measurement points are contained in each wavelength. Following the Nyquist criterion, a coarser measurement mesh could have been taken (≈ 8 times less points). However, taking a finer mesh improves the identification accuracy, as the parameter identification is made in the least-square sense (Eq. ( 25)).

Observing the wavenumber dependence as a function of frequency, two dispersion branches can be built, related to the bending and to the twist motions. Within the frame of the Euler beam theory, one can connect wavenumbers to the Young and shear moduli according to Eq. ( 7). The results of the identification of the apparent Young modulus E eul (in blue, from bending branch) and shear modulus G (in red, from twist branch) are shown in Fig. 7b (with κ = 0.845 for a square section and the PVC density ρ = 1380 kg/m 3 ). The low-frequency estimation is noisier because of the plateau effects of the HRWA related to the modal behaviour of the beam (see preceding section). As the frequency increases, the beam's response tends to be smoother and the number of spatial periods increases, leading to a better wavenumber estimation.

The dependency on frequency of the identified complex apparent Young modulus E eul is mainly due to the transverse shear effects generated by the bending motion, that increases with the frequency. To explain these effects, one can use the equivalent slenderness ratio defined in Eq. ( 8). At low frequencies, this ratio is large (µ(100 Hz) 42) and the beam behaves as a thin beam. As the frequency increases, µ decreases (µ(5 kHz) 6.5, µ(20 kHz) 2.7) and the beam behaves as a thick beam. In order to take these effects into account, the Timoshenko model is used to better estimate the Young modulus. By injecting the identified shear modulus G in the equation [START_REF] Oberst | Über die dämpfung der biegeschwingungen dünner bleche durch fest haftende beläge II. {On the damping of bending vibrations of thin sheets by adherent coverings II[END_REF] leads to the estimation of an other Young modulus, E tim (in green). The result is shown on Fig. 7b in yellow. It can be observed that the identified modulus is almost constant ( 5.5 GPa) in the considered frequency range. The resulting Poisson ratio is estimated as ν = E/(2G) -1 0.4. In addition, the identified Young loss factor η E tim seems very close to the shear loss factor η G .

Conclusion and Perspectives

A common framework for wavenumber extraction methods has been presented and the High Resolution Wavenumber Analysis has been proposed. This new method aims at overcoming three following drawbacks of the IWC and the Mc Daniel methods: (a) IWC suffers from resolution limitations (b) the nonlinear wavenumber search problem is computationally heavy and (c) the number of identified waves is fixed and small.

The HRWA procedure includes a linear search of wavenumbers. Thus, the computation cost is reduced and the wavenumber identification is free from any initialisation. Moreover, it allows an automated determination of wavenumbers in a noisy signal without a priori assumption on their number. The detected waves are not limited to bending waves but can be associated with twist, compression, shear, etc.

Firstly, numerical investigations and comparisons of HRWA with the two other methods has been led in simple cases, including that of a thin prismatic Euler beam. When only bending waves are present, the Mc Daniel method can be taken as a reference method.

In this simple case, this method is the least sensitive to noise. The IWC method fails to identify a spatial decay, in the frequency band studied. The HRWA exhibits a convergence rate as function of the signal to noise ratio comparable to the Mc Daniel method (∆k/k ≈ 0.02 SNR -1 ), with certain fluctuations for a very precise determination of the real part of the wavenumber (0.1% of relative error). These fluctuations are mainly due to the fact that the method has no a priori assumptions with regard to the number of wavenumbers present in the signal, as opposed to the Mc Daniel method. The spatial decay has a reduced influence on the measured beam harmonic response, compared to the real part of the wavenumber. Therefore, its determination is more sensitive to noise. An estimation with a relative error of 10% can be achieved with both Mc Daniel and HRWA methods, for good signal to noise ratios (SNR> 100). Moreover, the comparison of computational times shows that HRWA can be twenty times faster than the two other methods.

Secondly, the HRWA has been applied to experimental data. It emphasises the ability of the HRWA to simultaneously identify apparent complex Young and shear modulus, as a function of frequency, independently from boundary conditions or excitation type. Identified with the Euler model, the Young modulus vary in frequency. To account for the out-ofplane shear and rotary inertia effects, the Timoshenko model is finally used and leads to the identification of a Young modulus which is almost constant in the frequency range of interest.

It has been observed that most of the identification errors of the HRWA results are related to a wrong signal order estimation. Indeed, the ESTER criterion performs well only if the signal model is correct. In order to improve the robustness of the HRWA, other signal order selection criterion could be investigated. For example, the use of a stabilisation diagram [START_REF] Bakir | Automation of the stabilization diagrams for subspace based system identification[END_REF][START_REF] Reynders | Fully automated (operational) modal analysis[END_REF] to estimate the order of the signal could improve the accuracy of the wavenumber estimation, at the price of a reasonable increase of the computational burden.

At this point, it appears clearly that no hypothesis has been made on the nature of the waves: the HRWA is not limited to the simple beam models presented here. For example, a different section or material configuration may induce coupled strain mechanisms and coupling the characteristic in equations ( 2) and (3). The detected wavenumbers would be different, but still be representative of the structural behaviour, reasonably far from the boundary conditions.

The present paper is focused on material characterisation but the HRWA can also be used for the identification of boundary conditions (not presented here). Moreover, multidimensional versions of ESPRIT [START_REF] Rouquette | Estimation of frequencies and damping factors by two-dimensional ESPRIT type methods[END_REF] and ESTER [START_REF] Liu | Subspace techniques for multidimensional model order selection in colored noise[END_REF] have been developed and can be applied to the experimental identification of wavevectors in the harmonic response of plates. This will be the object of a forthcoming paper.

Figure 1 :

 1 Figure 1: Beam geometry and coordinate axis

Figure 2 :

 2 Figure 2: Typical setup for wavenumber identification methods.

Figure 3 :

 3 Figure 3: The simulated cantilever beam.

  40 mm 70 GPa 0.5% 0.3 2500

Figure 4 :

 4 Figure 4: Synthesised Mean Square Velocity without noise (top frame). Contour plots of the relative errors in the determination of the real part of the wavenumber k (middle frame) and the spatial decay γ (bottom frame) for different SNR levels (ordinates) as a function of the frequency. Comparison of the HRWA (black), Mc Daniel (red) and IWC (blue) methods. Contour lines display the 10% (continuous lines), 1% (dashed lines) and 0.1% (dotted lines) limits for the relative identification errors. The results are obtained as the mean of the estimation over 100 virtual tests.

Figure 5 :

 5 Figure 5: Relative error in the identification of the complex wavenumber k = k(1 -i • γ) as a function of the SNR for two fixed frequencies. Mean ± Standard Deviation. Comparison between HRWA (black), Mc Daniel (red) and IWC (blue) methods. (a) at 1 kHz, between two modal frequencies. (b) at 4.5 kHz, on a modal frequency.

Figure 6 :

 6 Figure 6: Computation time of each wavenumber identification method as a function of the measurement mesh size, averaged on 200 wavenumbers extractions. Comparison between methods.

Figure 7 :

 7 Figure 7: Experimental application of the HRWA to a PVC beam of length L = 1 m and square section of side h = b = 15 mm, for frequencies between 0 and 22 kHz. (a) Identification of wavenumbers on the measured beam with the proposed HRWA, for two excitation types: impulse (dashed lines) and steady noise (solid line). Bending wavenumber k b (blue) and twist wavenumber k t (red). (b) Apparent complex moduli identification from the extracted wavenumbers: With the Euler model (E eul (blue) and G (red) from Eqs. (7)) and the Timoshenko model (E tim (green) from Eq. (11))).

Table 1 :

 1 Characteristics of the beam.

Acknowledgments

This work is part of the ANR founded (french National Research Agency) project MAESSTRO (ModélisAtion Et Synthèse Sonore pour Tables d'haRmonie de pianO, ANR-14-CE07-0014). The authors thank Kerem Ege for his precious help on high-resolution signal processing methods.


Francesco Foscarin

Florent Jacquemard

Philippe Rigaux

Masahiko Sakai

A Parse-based Framework for Coupled Rhythm Quantization and Score Structuring

We present a formal language-based framework for MIDI-toscore transcription, the problem of converting a sequence of symbolic musical events with arbitrary timestamps into a structured music score. The framework aims at solving in one pass the two subproblems of rhythm quantization and score production. It relies, throughout the process, on an apriori hierarchical model of scores given by generative grammars. We show that this coupled approach helps to make relevant and interrelated decisions, and we present an algorithm computing transcription solutions optimal with respect to both the fitness of the quantization to the input, and a measure of complexity of music notation.

Introduction

Music transcription is the act of converting a music performance into music notation (i.e., a score). Several aspects of this problem are studied in the literature, according to the variety of music representation format considered. One of the most studied is the conversion of an audio recording into an unquantized MIDI file (audio-to-MIDI, A2M) [START_REF] Benetos | Automatic music transcription: challenges and future directions[END_REF], i.e. the extraction, from an audio signal, of a symbolic representation with explicit event descriptors such as pitch, onset and offset, expressed in a real-time unit (seconds). We focus on the complementary problem of converting unquantized MIDI into a music score (MIDI-to-score, M2S). M2S transcription can itself be divided into two subproblems: (i) Rhythm Quantization (RQ) is the conversion of time values from real-time units to musical-time units (beats, or fraction of bars) [START_REF] Raphael | A hybrid graphical model for rhythmic parsing[END_REF][START_REF] Cogliati | Transcribing human piano performances into music notation[END_REF][START_REF] Nakamura | Rhythm transcription of polyphonic piano music based on merged-output HMM for multiple voices[END_REF]. RQ alone is generally achieved via the manipulation of linear data structures (e.g., sequences of messages in MIDI files). (ii) Music score production (MSP) involves the determination of higher level information: voices, bars, grouping of events with tuplets and beams, encoding of durations with ties and dots etc, see [START_REF] Gould | Behind Bars: The Definitive Guide to Music Notation[END_REF]. A salient feature of music notation is its hierarchical nature, already advocated by the models such as Rhythm Trees [START_REF] Agon | Representation and rendering of rhythm structures[END_REF]. An accurate MSP procedure thus requires the manipulation of hierarchical data structures. Traditionally, subproblems (i) and (ii) are considered independently, in sequence. This allows to delegate subproblem (ii) to the MIDI import module of a score editor. Such an approach to M2S transcription might be satisfying in simple cases. It strongly depends, however, on how the quantized MIDI input fits the specifics of the music notation language. To put it differently, the linear structure produced by step (i) might be hardly compatible with the hierarchical notational structures of music scores used by step (ii). Even if each step yields quite satisfying results regarding its specific goal, their combination migh therefore exhibit a discrepancy, and possibly yield poor transcription results. This is particulalry true for complex rhythms with e.g. deep nesting, mixed tuplets, rests, grace notes, etc.

In the present paper, we propose a framework for M2S where subproblems (i) and (ii) are tightly coupled: the structural information needed for score construction in (ii) is built during step (i), and takes into consideration an apriori music notation model. More precisely, (i) is solved by a parsing algorithm, and the parse tree defines a rhythmic structure of the output score, similar to Rhythm Tree representations [START_REF] Agon | Representation and rendering of rhythm structures[END_REF][START_REF] Jacquemard | A Structural Theory of Rhythm Notation based on Tree Representations and Term Rewriting[END_REF]. This framework offers several distinctive advantages that contribute to improve the result quality. First, it makes it possible to jointly consider all the decisions made along the transcription process, and to model it as a multicriteria optimization problem, for two criterias: fitness between output to input, and of rhythmic complexity, in the sense of [START_REF] Nauert | A theory of complexity to constrain the approximation of arbitrary sequences of timepoints[END_REF]. Second, we can leverage on expressive and powerful computational formalisms. On the one hand, we rely on weighted context-free grammars (WCFG), a standard formalism for modeling and ranking hierarchical constituent structure in computational linguistics. WCFG are used to describe an a priori music notation language. On the other hand, optimal tree representation can be obtained by efficient parsing algorithms [START_REF] Huang | Advanced Dynamic Programming in Semiring and Hypergraph Frameworks[END_REF] using Dynamic Programming. One of the main contribution of the paper is to show that this formal machinery can be adapted to solve M2S accurately and efficiently.

We expose the formal background of M2S in Section 2. The algorithmic part is developed in Section 3. Section 5 concludes the paper by describing experiments and further work.

Framework Definition and Objective

M2S takes as input a sequence of (unquantized) events and returns a music notation of this sequence. This section presents the formalisms used to represent these input and output and to model the framework.

Time Units and Tempo Curves. Timestamps can be expressed either in real-time unit (rtu), used for unquantized events, or musical-time unit (mtu) for music score events. In both cases, the temporal domain is Q+. In the rest of the paper, we assume a rtu of 1 second and a mtu of 1 bar. Given a time signature, every time value in mtu can be converted to a value in beats. A tempo curve is monotonically increasing function θ : Q+ → Q+, converting rtu values into mtu values see e.g. [START_REF] Honing | From time to time: The representation of timing and tempo[END_REF]. Let τ = τ0, . . . , τm be a sequence of rtu values with τ0 = 0 (typically the timestamps of input events). A tempo model M compatible with τ is a set of tempo curves piecewise linear, with slope changing at the timestamps of τ . More precisely, every θ ∈ M is defined by a sequence T0, . . . , Tm-1 such that for all 0 ≤ i ≤ m-1, the restriction of θ to [τi, τi+1[ is a line of slope Ti ∈ Q+ (expressed in bars per second in our case, instead of bpm). Typically M expresses restrictions on the changes of Ti, according to the real durations τi+1 -τi, in order to ensure a certain smoothness of the tempo curves. We skip unecessary details about the specification of M.

Timelines. A time interval is a right-open interval I = [τ, τ [⊂ Q+ .
I is called unbounded when τ = +∞ and bounded otherwise. The left bound τ is called the start of I and denoted start(I). We call partition of a time interval I a sequence of disjoint time intervals I1, . . . , I k , with k ≥ 1, such that k j=1 I k = I. We also write I = I1 + . . . + I k in this case. The k-split of a bounded time interval I (for k > 0) is a partition of I of size k such that the duration of each component is dur (I) k . In the case of a 2-split I1, I2 of I, we write left(I) for I1 and right(I) for I2. We assume given a notational alphabet E, which is a finite set of symbols to encode musical artifacts.

Example 1. A possible choice for E is the set of MIDI message symbols extended for explicit representation of rests. Pitch values are in [0, 128], 128 being for rests, velocity in [0, 127], and a flag distinguishes onsets/offsets. This flag, useful for polyphonic music, can be skipped in monophonic case. ♦

An event e is a pair η, τ made of a symbol η ∈ E, denoted symb(e), and a timestamp τ ∈ Q, denoted date(e). A timeline I is a pair I, ē , where I is a time interval denoted carrier (I) and ē = e1, . . . , em is a finite sequence of events denoted events(I), with increasing timestamps and such that date(ei) ∈ I for all 1 ≤ i ≤ m. A timeline with timestamps in rtu (resp. mtu) is called a real -timeline (resp. musical -timeline). Operations on time intervals, like e.g. +, are extended to timelines as expected.

Example 2 (Toy running example). Let I1 and I2 be timelines defined by: carrier (I1) = [0, 1[ and events(I1) = e1, e2, e3 , with respective timestamps 0.07, 0.72, 0.91; carrier (I2) = [1, 2[ and events(I2) = e4, e5, e6 with respective timestamps 1.05, 1.36, 1.71. ♦

Semirings. Domains of weight used to rank solutions to transcription are abstractly defined as semirings, that can be instantiated into several concrete domains (e.g. probabilities or costs). A semiring S = S, ⊕, 0, ⊗, 1 is a structure with a domain S = dom(S), two associative binary operators ⊕ and ⊗, and neutral elements 0 and 1; ⊕ is commutative, ⊗ distributes over ⊕: ∀x, y, z ∈ S, x ⊗ (y ⊕ z) = (x ⊗ y) ⊕ (x ⊗ z), and 0 is absorbing for ⊗: ∀x ∈ S, 0 ⊗ x = x ⊗ 0 = 0. Components of a semiring S may be subscripted by S when needed. We simply write x ∈ S to mean x ∈ S.

Intuitively, in the application presented below, ⊕ selects an optimal value amongst two values and ⊗ combines two values into a single value.

A semiring S is commutative if ⊗ is commutative. It is idempotent if for all x ∈ S, x ⊕ x = x. It is monotonic wrt a partial ordering ≤ iff for all x, y, z, x ≤ y implies x ⊕ z ≤ y ⊕ z, x ⊗ z ≤ y ⊗ z and z ⊗ x ≤ z ⊗ y.
Every idempotent semiring S induces a partial ordering ≤S called the natural ordering of S and defined by: for all x and y, x ≤S y iff x⊕y = x. It holds then that S is monotonic wrt ≤S . S is called total if it is idempotent and ≤S is total i.e. when for all x and y, either x ⊕ y = x or x ⊕ y = y.

In practice, we use two kinds of total semirings: Viterbi semiring defining probability models, whose domain [0, 1] ⊂ R+, ⊕ is max, 0 = 0, ⊗ is real product, and 1 = 1, and tropical semirings, defining cost models whose domain R+ ∪ {+∞}, ⊕ is min, 0 = +∞, ⊗ is sum, and 1 = 0.

Weighted Context-Free Grammars. A WCFG over a semiring S and an alphabet E is a tuple G = Q, init, P, weight, mus where: Q is a finite set of non-terminal symbols (nt), init ∈ Q is an initial nonterminal, P is a set of production rules in one of the following forms: (k-div) q → q1, . . . , q k with q, q1, . . . , q k ∈ Q, and rank k > 1, or (term) q → ē with q ∈ Q and ē ∈ E * (ē is called terminal symbol). weight assigns to each production rule in P a weight value in S, mus assigns to each (k-div) production rule in P a function associating to a musical-time interval O a partition O1, . . . , O k of O.

The components of a WCFG G may be subscripted by G when needed. We use the respective notations q --→ w q1, . . . , qn and q --→ w ē for (k-div) and (term) productions rules of weight w ∈ S. The (k-div) rules (for k ≥ 2) define the possible divisions of musical time intervals, e.g. the division of a quarter note into 2 eight notes or into a triplet. The weight associates a rhythmic complexity in S to each division. The recursive application of (k-div) rules represents nested divisions. Their complexity values will be composed using ⊗.

Example 3. The following (2-, 3-div) production rules, with weight values in a tropical semiring, define two possible divisions of a bounded time interval represented by the nt q0, into respectively a duplet and a triplet.

ρ1 : q0 0.06 --→ q1, q2 , ρ2 : q0 0.12 --→ q1, q2, q2 .
In those rules, q1 represents the first event in a division, and q2 the others. Further binary divisions of time sub-intervals are possible with:

ρ3 : q2 0.1 --→ q3, q3 , ρ4 : q3 0.11 --→ q4, q4 . ♦
The (term) production rules specify the musical symbols of E that can occur in a time interval. An empty sequence in the right-hand-side of such rule represents the continuation of an event started before and not yet released -notated with a tie or a dot in a score. In practice, in order to keep G reasonably small, we use as set of terminal symbols a finite abstraction of E * , like in the following example.

Example 4. In the case of monophonic input, simultaneous events are interpreted as grace notes: a singleton sequence η1 represents a single note, η1, η2 represents a grace note η1 followed by a note η2, η1, η2, η3 represents two grace notes η1, η2 followed by a note η3, etc. The set F = {0, 1, 2, 3} is a finite abstraction of E * where the symbols 0, 1, 2 represent resp. a continuation, one and two symbols of E and 3 represents a sequence of E * of length ≥ 3. In the following, we assign a weight value (in a tropical semiring) to (term) productions rules of a grammar G, depending on the number of grace notes.

ρ5 : q0

0.15 --→ 0, ρ6 : q0 0.01 --→ 1, ρ7 : q0 0.79 --→ 2, ρ8 : q0 1.02 --→ 3, ρ9 : q1 0.02 --→ 0, ρ10 : q1 0.01 --→ 1, ρ11 : q1 0.25 --→ 2, ρ12 : q1 0.64 --→ 3, ρ13 : q2 0.02 --→ 1, ρ14 : q3 0.04 --→ 0, ρ15 : q3 0.01 --→ 1, ρ16 : q4 0.01 --→ 1
The nt q0 represents a whole bar, with a single event in ρ6 (e.g. a whole note in a 4 4 bar), as a tied note in ρ5 or preceded by 1 or 2 grace notes in ρ7 and ρ8; q1 represents a first note in a division with a rule of Example 3 (preceded by 1 or 2 grace notes in ρ11 and ρ12); q2 represents the next notes in the same divisions and q3 and q4 further levels of divisions (grace notes are not allowed for q2, q3, q4, ties are not allowed for q4).

♦ Symbols in abstractions of E * may embed more information, like e.g. pitch-contour for sequences of grace notes, or velocity, or on/off flag.

Parse Trees and Serialization. Given a WCTG G over a semiring S, the set T (G) of parse trees is the smallest set of trees labelled by production rules of G such that:

for all (term) rule ρ in G, ρ ∈ T (G), with root ρ, for all (k-div) rule ρ = q w -→ q1, . . . , q k in G, and all t1, . . . , t k ∈ T (G) whose respective roots have heads q1, . . . , q k , ρ(t1, . . . , t k ) ∈ T (G) with root ρ. In the second case, we call head of a rule its left-hand-side nt. We write T (G, q) for the subset of parse trees of T (G) whose root is headed with nt q. The weight of a parse tree is obtained by recursively applying ⊗.

weight(t) := weight(ρ) when the label of t is of type (term), weight ρ(t1, . . . , t k ) := weight(ρ) ⊗ k i=1 weight(ti) . We associate to every parse tree t of T (G) and every (output) mtu interval O a musical-timeline denoted t O and defined by:

-

q w -→ a O = O, η1, start(O) , . . . , ηp, start(O) , when a = η1, . . . , ηp , -ρ(t1, . . . , t k ) O = t1 O 1 +. . .+ t k O k when O1, . . . , O k = mus(ρ)(O).
This mapping, called serialization, defines the timestamps for the symbols of E labeling the leaves of t. Moreover, t also yields a grouping structure for the resulting events. In other terms, t is a consistent representation of music events with respect to the notation defined by the grammar, and a music score can be constructed straightforwardly from it.

Example 5. Taking the rules of Examples 3, 4, assuming that mus(ρ1), mus(ρ2), mus(ρ3), mus(ρ4) return respectively 2-, 3-, 2-and 2-splits of bounded intervals, Figure 1 presents parse trees and their serialization in a 1 bar mtu-interval. Note that t1 has 3 leaves, but t1 [0,1[ contains only 2 events, because its second leaf is a continuation with 0 event. ♦ We extend δ to sequences of events ē = e1, . . . , em and ē = e 1 , . . . , e n by δ ē, ē = m i=1 δ(ei, e i ) if m = n and 0 otherwise, and to timelines by δ(I, O) = δ events(I), events(O) .

Example 7. For a tropical semiring, let δ0(τ1, τ2) = |τ1-τ2|. Its extension δ is a measure of the accumulation of time-shifts of events, in rtu. ♦

We also use a measure γ of tempo variations defined by a tempo curve, based on a given γ0 : Q × Q+ → S. Given a real-timeline I, and a tempo curve θ in a model M compatible with events(I), γ is defined as

γ(θ) = m-1 i=0 γ0(Ti+1 -Ti, τi+1 -τi)
where T0, . . . , Tm-1 is the sequence of slope values defining θ (see page 2), and τ0, . . . , τm are the timestamps in events(I).

Example 8. For a S tropical, we can define γ0 as the ratio between the variation of slopes γ0(dT, dτ ) = dT dτ when dτ = 0, and γ0(dT, 0) = 0. ♦ Altogether, we define the fitness of a quantized musical-timeline O (a score) to the real-timeline I (a performance), wrt a tempo curve θ, as

fit(I, O, θ) = δ θ(I), O ⊗ γ(θ).
In our settings, the output timeline O will be the serialization t O of a parse tree t of a WCFG G over S (for a given mtu time interval O).

Transcription Objective. Assuming an alphabet E, a commutative, idempotent and total semiring S and a fitness measure based on δ0 and γ0 as above, the M2S problem is defined as follows.

Input:a real-timeline I, non-empty (i. Output:a tempo curve θ ∈ M defined on carrier (I), a parse tree t ∈ T (G), such that weight(t) ⊗ fit(I, t O , θ) is minimal wrt ≤S . Therefore, the objective of M2S is to find a parse tree t representing a score that optimizes a combination of its rhythmic complexity (weight wrt G), and its fitness to the input I. The two criteria are antinomic, in the sense that improving the fitness generally increases the complexity of the tree and viceversa. Let us discuss the relevance of the above combination with ⊗ by reviewing two concrete domains used for for S. (i) If S is a Viterbi semiring, then the weight wrt G and the fitness are probability values and we want to maximize their product with ⊗. (ii) If S is a tropical semiring, then the weight and the fitness can be seen as two unrelated quality criteria, and one can resettle M2S as a multicriteria optimization problem [START_REF] Jahn | Vector Optimization. Theory, Applications, and Extensions[END_REF]. Let P = T (G)×M be the solution space for M2S, and let us consider the two objective functions c and d of P into the tropical semiring S defined, for p = (t, θ) ∈ P, by c(p) = weight(t) and d(p) = fit I, t O , θ (for the given I and O). By monotonicity of S, we can restrict our search to so-called Pareto-optimal points p ∈ P, i.e. such that there is no p ∈ P with c(p ) <S c(p) and d(p ) <S d(p). M2S is expressed as the minimization of the combination c(p) ⊗ d(p), where ⊗ is interpreted as a sum in R+. This is similar to a technique of scalarization by weighted sum, selecting a point p with minimal α.c(p) + d(p) (in R+) called scalar optimal. Intuitively, α can be seen as a user parameter setting how much one want to favour the rhythmic complexity against the fitness. In practice, one can apply the coefficient α to weight(t) by multiplying by α all the weight values in productions of G. This approach is correct and complete in the following sense: every scalar optimal point p ∈ P (for some α) is Pareto-optimal and for all Paretooptimal point p ∈ P there exists a coefficient α such that p is a scalar optimal for α (Theorem 11.17 and Corollary 11.19 of [START_REF] Jahn | Vector Optimization. Theory, Applications, and Extensions[END_REF], chapter 11.2.1).

Transcription Algorithm

We now present a transcription algorithm that works in two steps: it computes first a WCFG K by augmenting G with some information from the input I (Sections 3.2, 3.3 describe two examples of construction of such K, corresponding to different use cases) and then it solves M2S by computing an optimal parse tree for K, using a Dynamic Programming algorithm presented in next Section 3.1.

Viterbi 1-best algorithm

Let K be a WCFG with nt set K over a total semiring S. The following recursive function bestK associates to every nt k ∈ K a parse tree t ∈ T (K, k) with a weight optimal wrt ≥S . By abuse, we make no distinction below between a parse tree t ∈ T (K) and its weight value weight(t) ∈ S. Since S is total, it means that ⊕, applied to parse trees of T (K), selects the tree with minimal weight wrt ≤S .

bestK(k) = ρ 0 =k w 0 --→ K a ρ0 ⊕ ρ=k w -→ K k 1 ,...,kn
ρ bestK(k1), . . . , bestK(kn)

If K is acyclic, then the above definition is well founded and the following results holds (e.g. by induction on k following a topological sort of K).

Lemma 1. For all k ∈ K, bestK(k) = t ∈T (K,k) t = min ≤ S {t ∈ T (K, k)}.
Remember that the weight of a parse tree is the product of all the weights of the productions rules labeling its nodes. Therefore, the above formula can be understood as an alternation of sums with ⊕ (selection of one parse tree of optimal weight) and products with ⊗ (of weights of all subtrees). The function best can be computed by a straightforward adaptation of a Viterbi-like Dynamic Programming algorithm returning the best derivations for weighted acyclic hypergraphs ([9], Section 5.1).

With a tabulation over K, in order to avoid recalculation of solution for subproblems, this algorithm runs in time linear in the size of K.

In the case where K is not acyclic, one can use a generalization by Knuth of the Dijkstra shortest path algorithm (Algorithm 6 of [9]) 5 . We apply this algorithm to a WCFG K built on the top of G from an input timeline I. Two particular computations of K corresponding to different case studies are presented below.

Constant Tempo

In this first case study, we assume a constant tempo. This case study, illustrated in Example 9, corresponds to performances recorded with a metronome. The tempo model is M = {θ0}, where a single tempo curve θ0 represents the constant tempo value T which, for the sake of simplicity, we assume below to be the identity (T = 1).

In this case, the purpose of M2S transcription is essentially to correct local time-shifts of events. A parse tree t defines a partition of a given time interval O, by recursive application of the division rules labeling the nodes of t (see the definition of t O ). This partition can be seen as a "grid '' containing the time positions of the bounds of the sub-intervals. M2S then consists in the realignment of the events of I to the nearest bound in the grid. The cost of this alignment is then the distance, with δ, between I and the score represented by t. K = E * , K, init, P, weight, mus is built to represent all the time subintervals defined from carrier (I) by recursive divisions, according to the rules of G. For this purpose, every nt k ∈ K contains two components accessible with the following attributes:

k.nt : a nt of G, k.car : a real-time interval embedded in carrier (I). (div) productions rules of P are of the form ρ = k w -→ k1, . . . , kn , where 1. ρ = k.nt --→ G w k1.nt, . . . , kn.nt is a rule of G, mus(ρ) = mus(ρ ), 2. k1.car , . . . , kn.car is the application of mus(ρ ) to k.car , 3. events(I| k.car ) = ∅, meaning that k is inhabited6 .

The last condition drastically reduces the size of K in practice, wlog since it is useless to divide empty intervals. (term) rules of P deal with the realignement of the events of I to the bounds defined by a parse tree of K, and compute the input-output distance as in Section 2. For a nt k, we know the events of I inside C = k.car . Some of these events (those in the first half of C, called early) will be aligned to the left bound of C. The others (those in the second half of C, called late) will be aligned to the right bound of C, i.e. actually to the left bound of the interval defined by the next leaf in a parse tree. To deal with this situation correctly, we add two components to every nt k ∈ K, accessible with the following attributes:

k.post : the late events of I in k.car (i.e. those in its second half), k.pre : a buffer memorizing the post for the previous leaf.

And we add the following conditions for every (div) production rule: This can be showed by induction on t, using the above conditions (2-4) for inductions steps. It follows from Proposition 1 and the fact that I| k+0 = I that πG(bestK(k0)), θ0 is a solution of M2S for G and I.

Coupled Tempo Inference and Rhythm Quantization

In this second case, we sketch a construction of grammar K that permits a joined evaluation of a tempo curve θ and a parse tree t. Instead of realigning late events to the right bound of the current real-time interval, we move the right bound of the interval, such that the modified interval contains only early events. The new interval length induces the definition of a piecewise linear tempo curve θ like in Section 2.

In addition to the attributes k.nt and k.car (Section 3.2), every nt k of K has a third attribute k.mod for a real-time interval modified from k.car . The latter is defined in (term) production rules k --→ w ē, such that events ē in k.mod are all early in this interval. The propagation of mod values for (div) rules is similar to that of pre and post (Section 3.2). Moreover, the nts of K also store the slope values induced by the k.mod , and the (div) rules also verify that the corresponding tempo curve θ (defined as in Section 2) belongs to the given model M.

Experiments

The above framework for M2S transcription has been fully implemented 7 in C++. We ran manual experiments with a dataset composed of monophonic extracts from the classical repertoire, chosen from a progressive textbook for learning rhythmic reading [START_REF] Lamarque | D'un rythme à l'autre[END_REF]. The extracts have a length of about 15 bars on average, different time signatures (assumed known in the experiments), and contain rhythmic notations of various style and complexity. They were recorded with a MIDI keyboard by several players of different levels (non-pianist, semi-professional and professional pianists) and processed using a command line executable. The WCFGs used for experiments, ranged from generic grammars crafted manually, to specific grammars learned from the score [START_REF] Foscarin | Modeling and Learning Rhythm Structure[END_REF]. These experiments gave promising results, especially in cases that are traditionally complex to handle (Figure 2), e.g. complex rhythms with mixed tuplets or alternations between short and long notes, grace notes, etc. The transcription time for each extract was below 100ms. Several optimisations (out of the scope of this paper) on the internal representation of grammars (with attributes) and scores (by binary trees instead of sequences) permitted important efficiency gains. We are currently developing a general framework for automated evaluation, that stores MIDI inputs along references scores, trains WCFGs [START_REF] Foscarin | Modeling and Learning Rhythm Structure[END_REF] from style-consistent corpus, and identifies errors in transcription output.

Conclusion

We have presented a modular framework for M2S transcription. It relies on a parsing algorithm for a given weighted grammar G, and handles the problem of complex rhythms and grace notes. Grammar models contrast with the linear Markovian models used in [START_REF] Raphael | A hybrid graphical model for rhythmic parsing[END_REF][START_REF] Cogliati | Transcribing human piano performances into music notation[END_REF][START_REF] Nakamura | Rhythm transcription of polyphonic piano music based on merged-output HMM for multiple voices[END_REF]. Using hierarchical models is a trend sucessfully explored for rhythmic notation processing e.g. [START_REF] Nauert | A theory of complexity to constrain the approximation of arbitrary sequences of timepoints[END_REF][START_REF] Agon | Representation and rendering of rhythm structures[END_REF][START_REF] Jacquemard | A Structural Theory of Rhythm Notation based on Tree Representations and Term Rewriting[END_REF], meter detection [START_REF] Mcleod | Meter Detection in Symbolic Music using a Lexicalized PCFG[END_REF], melodic search [START_REF] Bernabeu | Melodic identification using probabilistic tree automata[END_REF], and music analysis [START_REF] Granroth-Wilding | Statistical parsing for harmonic analysis of jazz chord sequences[END_REF][START_REF] Rohrmeier | Towards a generative syntax of tonal harmony[END_REF][START_REF] Nakamura | Treestructured probabilistic model of monophonic written music based on the generative theory of tonal music[END_REF][START_REF] Marsden | No Longer 'Somewhat Arbitrary': Calculating Salience in GTTM-style Reduction[END_REF]. The approach is founded on the conviction that music structure complexity exceeds linear models. An application of the same framework to polyphonic inputs is also under study, applying our framework to note-on and note-off input events to couple the voice separation problem with rhythm quantization.

7 See https://qparse.gitlabpages.inria.fr for sources and complete examples.

Fig. 1 .Example 6 .

 16 Fig. 1. Some parse trees and their linearization (Example 5).

  e. with |events(I)| > 0), a WCFG G = E * , Q, init, P, weight, mus over S, a musical-time interval O = [0, N [ with N ∈ N+ ∪ {+∞}, a tempo model M compatible with events(I).

Example 9 .

 9 Let I = I1 + I2 (see Example 2), G be the WCFG from the Ex. 3,4,6, a musical-time interval O = [0, 1[ and a tempo model containing a single tempo curve θ0 mapping 1 seconds to 1 measure. Two possible solutions to M2S are the parse trees t5 = ρ0(t2, t3) (Ex. 6) and t6 = ρ0(t1, t4). Their serialization t5 O and t6 O have both six events with respective timestamps (0, and t5 and t6 respectively correspond to and . Using the distance defined in the Example 7 we obtain weight(t5) ⊗ fit(I, t5 O , θ0) = 0.51+0.265 = 0.775 and weight(t5)⊗fit(I, t6 O , θ0) = 0.63 + 0.32 = 0.95. That means that t5 is preferred over t6, and actually the algorithm of Section 3.2 will return t5 as optimal solution. The reason is that in G the weight for having a grace note (rule ρ11) is quite high compared to the other rules. If we lower the weight of ρ11 to 0.07, then weight(t6) ⊗ fit(I, t6 O , θ) = 0.77 and t6 becomes the optimal solution. This illustrates the notation preferences defined in G, e.g. for favouring grace-notes or precise rhythmic notation. ♦

4 .

 4 k1.pre = k.pre, ∀1 ≤ i ≤ n ki+1.pre = ki.post, kn.post = k.post. This ensures a correct bottom-up propagation of pre and post values. The other rules of P have the form (term) k --→ w ē, where ē = k.pre + events I| left(k.car ) (concatenation of the pre buffer and the vector of events of I early in k.car ), events I| right(k.car ) = k.post and w = w0 ⊗ |ē| i=1 δ0 start(k.car ), date(ei) , with q --→ w 0 ē in G. The weight of the above (term) rule combines with ⊗ the weight w0 of the corresponding rule in G (i.e. a complexity value) with the distance for the realignment of the points of ē from their positions in I to new positions defined by the bounds of k.car . Moreover, k0 = initK is defined by k0.nt = initG, k0.car = carrier (I), k0.pre = k0.post = 0.

Fig. 2 .

 2 Fig. 2. Two extracts of transcription experiments; the output of our framework (qparse) is compared with the output of commercial softwares 7 .

  Correctness and Completeness of Construction. Let us now sketch a proof that from bestK(k0), one can build a solution for M2S conform to the definition in Section 2. First, one can observe that every parse

tree t ∈ T (K) can be projected onto a parse tree πG(t) ∈ T (G). Indeed, by condition (1.) for (div) rules of K, it is sufficient to replace, in every label of t, every nt k by k.nt and, for (term) rules, to replace the rule's weight by the weight defined in G. Next, for k ∈ K, let use us define I| k by carrier (I| k ) = carrier (I| k.car ) and events(I| k ) = k.pre + events I| left(k.car ) (using the above notations). Proposition 1. For all k ∈ K and t ∈ T (K, k), it holds that weight(t) = weight(t ) ⊗ fit I| k , t k.car , θ0 , where t = πG(t).

Acyclic WCFG are however sufficient for our purpose.

The restriction I|C of I is such that by carrier (I|C ) = C and events(I|C ) is the sub-sequence of events of events(I) inside C.