P.-A Taillard
email: taillard@hispeed.ch

F Silva
email: silva@lma.cnrs-mrs.fr

Ph Guillemain
email: guillemain@lma.cnrs-mrs.fr

J Kergomard
email: kergomard@lma.cnrs-mrs.fr

Modal analysis of the input impedance of wind instruments. Application to the sound synthesis of a clarinet

Keywords: Wind Instruments, Modal Analysis, Sound Synthesis

This paper investigates the modal analysis of wind instruments as seen from the input of their air column. Beside the treatment of analytical models, a particular emphasis is given to the analysis of measured input impedances. This requires special care because the measurements cover only a limited frequency band and are affected by some unknown errors. This paper describes how the Prony analysis and the Least Squares Complex Exponential (LSCE) classical techniques can be used in this context and how the main pitfalls can be avoided in their application. A physically acceptable method of reconstruction of the low frequency band is proposed. A technique using fictitious points in the high frequency range is described in order to ensure the passivity of the resonator in the whole frequency band. The principles of a real-time synthesis of clarinet sounds based on the modal representation of the resonator is given as an application, with a method to efficiently handle the modal representation during the transition between fingerings. A musically relevant example finally illustrates the possibilities of the modal analysis applied to wind instruments.

Introduction

Analysis-synthesis of musical sounds produced by self-sustained oscillations in wind instruments is a difficult task, not yet fully attained for some synthesis models. This goal has been achieved for decades for some linear models, despite the difficulty to control such simple models. However, when the synthesis model relies on the physics of the functioning of the instrument, it may include strong nonlinearities that make the estimation of the parameters from the analysis of natural sounds challenging. The synthesis model incorporates playing parameters that generally remain unknown, unless they can be measured simultaneously with the sound. Many subsystems are present (excitation, nonlinear coupling by the input flow, resonator with a complex geometry, nonlinear effects at loud level) and the target is also difficult to define (radiated sound, mouthpiece pressure). The approach proposed in this paper aims at modeling precisely one of these subsystems, the resonator, from experimental data. The resonator is assumed to be linear (its features are independent of the sound level) and is coupled to a classical simple production model. Moreover this allows, by comparison with naturally produced sounds, to collect information on the shortcomings of both the production and the radiation models.

Following the pioneering work of Helmholtz [1], sound production in wind instruments is explained as the mutual coupling of a resonator and an exciter, such that it is often enough to know the behavior of the acoustical resonator as seen from the coupling point only. For this reason, the bore is often characterized, experimentally or theoretically, by its input impedance. For example, the latter has been extensively used for the prediction of the oscillation frequency of the full instrument based on the resonance frequencies [2], or transformed either into the impulse response [3] or into the reflection function [START_REF] Schumacher | Ab initio calculations of the oscillations of a clarinet[END_REF] to be used in numerical schemes for sound synthesis. Another usage relates to the prediction of the steady state periodic regime by means of the harmonic balance method [5]. To the best of our knowledge, Ref. [START_REF] Agulló | The reflection function r(t): a matrix approach versus F F T -1[END_REF] is the first to explore the modal analysis of woodwinds with the parametrization of the reflection function in damped exponentials, with an application to the simulation of some regimes of the tenora [7]. According to McIntyre et al [START_REF] Mcintyre | On the oscillations of musical instruments[END_REF], using the reflection function is allegedly more efficient than the impulse response for time-domain simulation. However, this is only true for schemes involving convolutions. Ref. [START_REF] Ph | On the use of the modal expansion for the sound synthesis in real-time: pitfalls and solutions[END_REF] shows that the modal representation of the input impedance replaces the convolution by lightweight IIR filters and enables synthesis schemes that are at least as efficient. The modal series can then be used directly for time-domain simulations of self-sustained oscillations (as in Ref. [START_REF] Silva | MoReeSC: a framework for the simulation and analysis of sound production in reed and brass instruments[END_REF]) or the estimation of oscillation thresholds [START_REF] Velut | Numerical Simulation of the Production of Pedal Notes in Brass Instruments[END_REF].

Experimental modal testing is a broad topic in acoustics and vibration. It usually benefits from the observation of the quantities of interest (e.g., displacement, velocity, acceleration, acoustic pressure) at several locations, but also from the use of multiple actuators (multiple-input-multiple-output, MIMO) to efficiently excite and identify the various modal shapes and frequencies [START_REF] Maia | Theoretical and experimental modal analysis[END_REF]. In fact, the modal poles are intrinsic characteristics of the tested system and are therefore common to all observations. As a consequence, the redundancy of the information improves the robustness of the identification techniques. Modal analysis using single-input-single-output (SISO) configurations entails particular difficulties, some of which relate to the loss of redundancy. This is the case of wind instruments characterized solely by their input impedance. Similarly to Ref. [START_REF] Gazengel | Time domain simulation of single-reed wind instrument. From the measured input impedance to the synthesis signal[END_REF] which showed the way to acceptable estimations of reflection function, the present paper aims at exposing the possible pitfalls of the modal analysis of wind instruments and how to avoid them.

In the case of an analytical model for the acoustic resonator, the modal parameters can be theoretically derived from the poles of the input impedance, but it generally accounts for an infinite number of modes due to the transcendental equations arising from the modeling of the wave propagation in the bore. The truncation of the modal series leads to a degraded reconstruction of the impedance, as explained in Ref. [START_REF] Ph | On the use of the modal expansion for the sound synthesis in real-time: pitfalls and solutions[END_REF] and again in Sec. 3. This also applies for the modal analysis of systems modeled by partial differential equations and discretized by the finite elements methods, such as the approach used by the Modalys software [START_REF] Ellis | Modalys demonstration[END_REF]. The structural dynamics community proposes several methods to compensate for the effect of neglected modes on the low-frequency range (known as static condensation techniques, see, e.g., Ref. [START_REF] Qu | Model Order Reduction Techniques[END_REF]) which is not typically relevant in the field of acoustics. Dynamic condensation methods try to restore the inertial effect of the deleted modes but result in a nonlinear problem in order to determine the modal frequencies and shapes.

Modal analysis is usually performed for wind instruments as an optimization process that minimizes the error between the measured impedance and the impedance reconstructed from the finite modal series. The optimized variables may be the poles and the modal coefficients (see Refs. [START_REF] Doc | A minimal model of a single-reed instrument producing quasi-periodic sounds[END_REF][START_REF] Velut | Numerical Simulation of the Production of Pedal Notes in Brass Instruments[END_REF]) or the poles only (Ref [START_REF] Maestre | Design of recursive digital filters in parallel form by linearly constrained pole optimization[END_REF], where the modal coefficients are explicitly determined within the evaluation of the cost function). Both these methods rely on an iterative procedure that requires a decent initialization to guarantee feasibility and a fast convergence.

The main goal of this paper is to propose an efficient method for the parametrization of theoretical or measured input impedances of acoustic waveguides in order to make them suitable for, e. g., the real-time physics-based synthesis of the sound of wind instruments or for the analysis of their functioning using methods from the dynamical systems community [START_REF] Velut | Numerical Simulation of the Production of Pedal Notes in Brass Instruments[END_REF]. The parametrization must preserve the passivity of the resonator over the full frequency range, while the above mentioned applications additionally require the models to introduce a low number of parameters (here a reasonable number of modes to be taken into account).

This paper is organized as follows: Sec. 2 reviews the basis of acoustic waveguides and of the modal analysis. Sec. 3 exposes the problems related to the truncation of the modal series of a transfer function. Sec. 4 treats the case of fitting a given transfer function to a modal series. Sec. 5 describes an application of the fitting method to the modeled or measured acoustic resonator. Sec. 6 shows how the results of the modal analysis can be used for sound synthesis and, in particular, how the transition between fingering can be handled. It also includes a musical example based on measurements of the input impedance of a clarinet, followed by a conclusion in Sec. 7.

Modal analysis in acoustic waveguides

Input impedance and reflection coefficient

Within linear acoustic theory, when the acoustic wavelength is large compared to the cross section dimensions of the bore, it is usual to consider one-dimensional models where acoustics is described in terms of the flow rate u(t, x) and the pressure p(t, x), or on their frequency-domain counterparts U (ω, x) and P (ω, x), respectively. These only depend on the angular frequency ω and on the axial abscissa x in the waveguide.

At the input of the duct (x = 0), the flow rate U (ω) = U (ω, 0) and the pressure P (ω) are related by the dimensionless input impedance Z in (ω) or the input admittance Y in (ω), both being frequency-dependent

Z in (ω) = 1 Y in (ω) = P (ω) Z c U (ω) . (1)
The characteristic impedance Z c = ρc/S in depends on the density of air ρ, the speed of sound c and the cross section area at the input of the pipe S in . The reflection coefficient R(ω) is defined as follows

R(ω) = Z in (ω) -1 Z in (ω) + 1 ⇔ Z in (ω) = 1 + R(ω) 1 -R(ω) . (2)
The inverse Fourier transforms of the input impedance and the reflection coefficient are called the impulse response h(t) and the reflection function, respectively. Traditional wind instruments are passive device as there is no energy production within the bore. Acoustical energy can only be dissipated. This implies that the energy flux at the input is positive, which is ensured by the constraint Re[Z(ω)] ≥ 0, or, equivalently, |R(ω)| ≤ 1, over the full frequency range. In addition, some other properties are commonly found in waveguides. At very low frequencies, the duct is mainly resistive: the flow behaves as a slowly-varying laminar one and the flow rate responds in phase with the input pressure in a way similar to that of the Poiseuille flow. The input impedance may even be assumed to vanish at zero frequency, as is the case in this paper. Conversely, at very high frequencies, dissipation becomes strong enough to damp any resonance, so that Z in approaches 1 (i.e., R approaches 0). This applies approximately for woodwind instruments and for brass instruments (when one removes the mouthpiece that would otherwise lead to an instantaneous reflection at the input of the bore.)

Example: the cylindrical pipe

Consider a cylindrical pipe with radius r and length . Following Ref. [START_REF] Chaigne | Acoustics of musical instruments[END_REF], the visco-thermal boundary layer losses are taken into account by means of the first-order approximation of the propagation constant Γ

Γ(ω) = jω c + 3.10 -5 r jω π , (3)
where R and the ω are expressed in MKS units, j = √ -1 is the imaginary unit. The value of the characteristic impedance is assumed to be the lossless value Z c = ρc/(πr 2 ) as explained in [START_REF] Chaigne | Acoustics of musical instruments[END_REF]. In Eq. ( 3), the term √ jω could lead to difficulties related to its branch cut [START_REF] Th | Representations with poles and cuts for the time-domain simulation of fractional systems and irrational transfer functions[END_REF]. However for visco-thermal losses, the use of the complete formula by Zwikker and Kosten leads to functions having only poles (see [START_REF] Kergomard | General equivalent circuits for acoustic horns[END_REF][START_REF] Bilbao | Passive models of viscothermal wave propagation in acoustic tubes[END_REF]). Therefore Eq. ( 3) is an excellent approximation except at extremely low frequencies, so we assume that the calculation of the poles at higher frequencies is satisfactory.

At the open end, the radiation is modeled by the dimensionless radiation impedance Z R given by Silva et al [START_REF] Silva | Approximation of the acoustic radiation impedance of a cylindrical pipe[END_REF] (see also [START_REF] Kergomard | Matching of fundamental modes at a junction of a cylinder and a truncated cone; Application to the calculation of some radiation impedances[END_REF]) as a fraction of two polynomials N R (ω)/D R (ω). Then, accounting for the transfer matrix between the input and the output of the pipe, the dimensionless input impedance can be expressed as the ratio of two functions Z(ω) = N (ω)/D(ω), N and D having zeros only (i.e., no poles)

N (ω) D(ω) = cosh (Γ(ω) ) sinh (Γ(ω) ) sinh (Γ(ω) ) cosh (Γ(ω) ) N R (ω) D R (ω) . (4)
This can easily be extended to other geometries like pipes with tone holes or branched tubes. Fig. 1 shows an example of the input impedance representation for a clarinet-like bore. It is evident that the waveguide is passive (| arg (Z(ω)) | ≤ π/2 ⇒ Re[Z(ω)] ≥ 0) and that the modulus of the reflection coefficient monotonically decreases with ω. A strong feature of wind instrument bores is the existence of weakly damped resonances visible on the modulus of the input impedance, that enable the musician to produce stable and clear tones. For most woodwind bores, the resonances are approximately harmonic, which influences the tuning and timbre of the self-sustained oscillations (see, e.g., [START_REF] Doc | A minimal model of a single-reed instrument producing quasi-periodic sounds[END_REF]).

Modal analysis

As mentioned in the introduction, a parametrization of the input impedance or of the reflection coefficient is required for numerical studies such as real-time simulation, or bifurcation analysis. A classical method is based upon the identification of the resonances in terms of angular frequencies ω m , quality factors Q m and magnitude Z m . We prefer the more general parametrization in terms of poles s m and modal coefficients C m to reconstruct the impedance as a series of modal contributions

Z in (ω) = m Z m 1 + jQ m ω ωm -ωm ω = m C m jω -s m . ( 5
)
For a bore with a high number of resonances M , the parametrization relies on 2M complex quantities, C m and s m , which may become computationally prohibitive. A model reduction can be achieved by truncating the series to account only for the resonances that are significant in the frequency range of interest.

Resonance frequencies above some cutoff frequency F cut are discarded, and the frequency range above F cut is referred to as the stop band, similar to the stop band of lowpass analog and digital filters. This definition has nothing to do with the distinction introduced by Benade [START_REF] Benade | Fundamentals of musical acoustics[END_REF] about the effects of a tone hole lattice on the input impedance of wind instruments. In addition, we denote F min the lower bound for which the model or measurement is valid, and the frequency range between F min and F cut is referred to as the 4)). = 1 m, r = 1 cm. pass band. Finally, from the zero frequency to the lower bound F min , the model or the measurement has to be extrapolated according to some physically acceptable rules. This frequency range is thus named the reconstructed band.

Effect of the truncation of a modal series

The purpose of this section is to highlight the influence of the truncation of the modal series. This truncation is required as only a finite number of modes can be dealt with numerically. To remove any bias from measurements, the analysis is performed on an input impedance calculated by classical analytical models. It is thus possible to compute a theoretically exact modal expansion using complex analysis. This is of limited utility as the models are not valid across the entire frequency range, and cumbersome in that it requires the model to be extended to the Laplace domain, i.e., as a function of the Laplace variable s. However it is useful to demonstrate the consequences of the truncation.

Considering the input impedance Z in (s) = N (s)/D(s) of the cylindrical pipe (see Sec. 2), there are an infinite number of poles s m , i.e., of roots of D(s). The application of the Cauchy's integral theorem shows that the poles s m have negative real parts, so that the corresponding time-domain functions exp (s m t) decrease with time. The application of the residue calculus [START_REF] Markushevich | Theory of functions of a complex variable[END_REF] to the input impedance yields

Z in (s) = m∈N C m s -s m with C m = N (s m ) D (s m ) , (6)
where the C m are the residues of Z in at the poles s m , assuming the poles to be simple (i.e., of order one).

The denominator D (s m ) of the residue is the derivative of D with respect to the variable s at the pole s m .

A finite number of poles are estimated numerically using the Newton-Raphson method. For each of the first M poles with a positive imaginary part, the iterative algorithm is initialized using one of the resonance frequency of the input impedance Z in (s = jω) evaluated on the frequency axis. The resulting poles are shown in Fig. 2, and evidence that eigenfrequencies Im[s m ] are almost odd multiples of the first one, and that damping increases with frequency. The residues C m can be analytically derived using the expressions of the radiation impedance and of the propagation constant, and the modal expansion of the input impedance is computed using Eq. ( 6) considering the M first poles (indexed from zero)

Z modal (s) = M -1 m=0 C m s -s m + C * m s -s * m = 2 M -1 m=0 s Re[C m ] -Re[C m s * m ] s 2 -2s Re[s m ] +|s m | 2 . ( 7
)
The resulting impedance and the related reflection coefficient evaluated on the frequency axis (i.e., for s = jω) are displayed in Fig. 3. Several drawbacks can be mentioned. First the modal impedance restores the resonant behaviour in the vicinity of the first resonances and then slowly decays to 0 with a -π/2 argument. This induces a reflection coefficient that asymptotically tends to unity above the last resonance frequency taken into account. The resulting resonator thus poorly dissipates power at high frequencies. Furthermore, the poles/residues expansion is exact when considering the infinite summation. Its truncation still preserves the resonances of the input impedance but the anti-resonances are perturbed due to the missing inertial contribution of the ignored higher frequency poles. This may have strong consequences (for example on the magnitude of the even harmonics of the mouthpiece pressure) and results from the choice of the elementary functions in Eq. ( 5) that essentially focuses on the resonances, the reconstruction out of the resonances being a by product that is only correct when accounting for the infinite set of poles. This is also visible on the modulus of the reflection coefficient that is correct at the resonance frequencies but oscillates in-between with a deviation that comes close to 1 near the truncation limit.

As a conclusion, the theoretical modal expansion is a powerful mathematical tool but it shows practical limitations due to the truncation of the series. In regard to the constraint of the finite number of poles that can be represented numerically, the poles/residues decomposition has to be adapted, for example, by allowing the poles s m to deviate from the analytical ones. This reverts to a fitting procedure as the one described in the next sections.

Other options are possible, for example, adding a corrective term such as the one introduced by Guillemain and Silva [START_REF] Ph | On the use of the modal expansion for the sound synthesis in real-time: pitfalls and solutions[END_REF] in the specific case of a clarinet-like bore. It is designed as a high-pass filter that is able 4), thin black line, same as Fig. 1) and from the modal expansion truncated to 20 modes (see Sec. 3, Eq. ( 7), thick red line). = 1 m, r = 1 cm.

|R| Reflection Coefficient R 0 1 2 3 4 5 f (kHz) -π/2 0 π/2 arg (Z) (rad) 0 1 2 3 4 5 f (kHz) 0 2 /c 4 /c τ g Analytical Modal
to restore the contribution of the ignored higher order poles in the low or mid frequency range and, as a consequence, the position of the anti-resonances. It also tends to the characteristic impedance above the frequency of truncation and thus ensures that the reflection coefficient decays to zero in the high frequency range. The parameters are determined based only on the analytical expression. The generalization of this method is outside the scope of the current paper. In the next section, a fitting procedure is proposed.

Fitting a given transfer function to a finite modal expansion

Prony analysis

Prony analysis, from the name of the French mathematician who developed the method in the late 18th century, is the transposition of the Fourier analysis for damped oscillations. It relies on the assumption that a given causal signal h(t) expands on a series of damped sinusoids

h(t) = Heaviside(t) M -1 m=0 C m e smt , (8)
or, considering the Laplace transform H(s) of the signal

H(s) = M -1 m=0 C m s -s m , (9)
where the modal frequencies s m and the modal coefficients C m are such that the frequency response H(s = jω) is hermitian symmetric, or, equivalently, that h(t) is a real function. This applies when H(s) is an input impedance and h(t) the related impulse response.

A measured signal, sampled with a time step T and assumed to vanish for negative time is expressed as

∀n ≥ 0, h[n] = M -1 m=0 C m y n m ≡ M -1 m=0 C m e nsmT , (10)
for y m = e smT or, equivalently, s m = (log |y m | + j arg (y m )) /T . This expansion implies that the signal h[n] is autoregressive, i.e., there exists a set of coefficients

(β 0 , β 1 . . . β M -1 , β M = 1) such as ∀n ≥ M, h[n] + M k=1 β M -k h[n -k] = 0. (11)
In fact it is possible to build an M -order polynomial with roots y m . Its coefficients are denoted as β k (with

leading coefficient β M = 1), so that ∀n ≥ M M k=0 β M -k h[n-k] = M k=0 β k h[n+k-M ] = M k=0 β k M -1 m=0 C m y n+k-M m = M -1 m=0 C m y n-M m M k=0 β k y k m =0 = 0. ( 12
)
The evaluation of Eq. ( 11) for n ∈ [M, 2M -1] defines the matrix equation

     h[M ] h[M + 1] . . . h[2M -1]      h = -      h[0] h[1] . . . h[M -1] h[1] h[2] . . . h[M ] . . . . . . . . . h[M -1] h[M ] . . . h[2M -2]      A •      β 0 β 1 . . . β M -1      β . (13)
The basic idea of the Prony analysis is that the poles s m derive from the roots y m of a polynomial whose coefficients are the β k which in turn can be obtained by numerically solving Eq. ( 13).

Least squares complex exponentials (LSCE)

The Prony method can deal with noisy data by considering more than M evaluations of the autoregressive equation ( 11), e.g. for n ∈ [M, M + N -1] with N > M . The matrix A is then rectangular and the matrix equation Aβ = -h is overdetermined. It is solved using the least squares Moore-Penrose pseudoinverse of A. This extension of the Prony analysis is the so-called Least squares complex exponentials method (LSCE, see, e.g., Ref. [START_REF] Maia | Theoretical and experimental modal analysis[END_REF]).

Evaluation of modal coefficients

The coefficients C m can be estimated from Eq. ( 10) for 0

≤ n < M        h[0] h[1] h[2] . . . h[M -1]        h =        1 1 . . . 1 y 0 y 1 . . . y M -1 y 2 0 y 2 1 . . . y 2 M -1 . . . . . . . . . y M -1 0 y M -1 1 . . . y M -1 M -1        B •      C 0 C 1 . . . C M -1      C ( 14
)
where B is a Vandermonde matrix and is invertible for simple poles. This equation can be problematic if h is obtained using a discrete Fourier transform (see last paragraph in Sec. 4.4). An evaluation in the spectral domain is usually preferable. The coefficients C m can be obtained from Eq. ( 9) by using the frequency response function (FRF) H n = H(s = jn∆ω) measured on a regularly sampled frequency grid

B • C = H (15)
where

H = H -N +1 . . . H N -1
T and B is a (2N -1) × M matrix with generic term B nm = 1/ (j(n -N )∆ω -s m+1 ). This includes the evaluation of the FRF for both positive and negative frequencies, spanning from low frequencies to the cutoff frequency, so that ∆ω = 2πF cut /(N -1). In order to make the estimation robust with respect to noise, the number N of computed or measured frequencies is generally much larger than M . The matrix B is thus rectangular and the system has to be solved again in the least squares sense, using the pseudo inverse of B .

Taking the z-transform of Eq. ( 10)

H(z) = n≥0 h[n]z -n = M -1 m=0 C m 1 -y m z -1 , (16)
the same procedure as Eq. ( 15) can be applied considering the frequency-domain characterization of the discrete-time transfer function, i.e., H(z) evaluated on the unit circle (z = e j(n-N )∆ωT for n

∈ [1, 2N - 1]
). The generic term of the matrix B is

B nm = 1 1 -y m e j(n-N )∆ωT . (17)

Numerical considerations

The methods proposed in the previous sections need to be treated carefully from a numerical perspective. First, the calculations involve exponentiation and high orders polynomials (see, e.g., Eqs. ( 12) and ( 14)) which is known to be problematic. Furthermore, when some modes are highly damped or when the C m span on many orders of magnitude, rounding errors may degrade the numerical accuracy of the computations. The use of high precision arithmetic is necessary essentially in the evaluation of the poles. Once this evaluation is done, the remaining operations do not need to be performed in high precision.

A second point concerns the fact that unstable poles (poles s m with positive real part) may emerge as roots of the polynomial defined by the coefficients β m . They may result from an overestimate of the number of modes M : the Prony analysis then uses the additional degrees of freedom to overfit noisy data. They have no physical meaning and can lead to growing exponentials in the modal synthesis. Once the poles have been estimated using Eq. ( 12), values of y m such that |y m | > 1 are discarded.

There is also a visible violation of the Nyquist-Shannon sampling theorem as the modal expansion in Eq. ( 9) has an unlimited bandwidth. Although the Prony's method is somewhat resilient to moderate violations of the sampling theorem, the estimation is improved if the transfer function H(ω) is guaranteed to converge towards zero instead of a finite non zero value. This may require a shift of the transfer function H → H -H(ω → ∞) which corresponds to an additional pole s m → -∞ such that y m = 0, i.e., a Dirac impulse in h(t).

Another important issue concerns Eqs. ( 15) and ( 17): the frequencies considered in the estimation of the modal coefficients only span the pass band. Without any additional constraint, the least squares fit may lead to unwanted behaviors in the stop band, for example violating the passivity property at high frequencies. This can be problematic when using the modal series for synthesis with a sampling frequency F s is much higher than the cutoff frequency F cut . It is thus necessary to enforce fictitious target points between F cut and F s /2 in order to guarantee the passivity, without degrading the fit in the pass band. This is done by extending the matrix B and the vector H by N fake values Ĥn for angular frequencies ωn > 2πF cut and their N corresponding hermitian symmetric counterparts Ĥn+N = Ĥ * n and ωn+N = -ω n .

        H Ĥ1 . . . Ĥ2N         H =         B 1 -y 0 e j ω1 T -1 . . . 1 -y M -1 e j ω1 T -1 . . . . . . 1 -y 0 e j ω2N T -1 . . . 1 -y M -1 e j ω2N T -1         B •C. ( 18
)
The choice of the fictitious points (ω n , Ĥn ) for the case of the modal analysis of wind instruments is discussed in Sec. 5. Finally, in musical acoustics, the measurements are usually performed in the spectral domain. As extensively explained in Ref. [START_REF] Gazengel | Time domain simulation of single-reed wind instrument. From the measured input impedance to the synthesis signal[END_REF], when the impedance Z(ω) is not purely resistive at the maximum measured frequency f cut (i.e., the impedance does not coincide with a resonance or an antiresonance at that frequency), strong ripple appears on the time-domain signals r(t) and h(t) obtained by the discrete inverse Fourier transform. This is known to lead to simulations with no physical sense. In the present problem, the ripple can also interfere with the evaluation of the modal coefficients C m by means of the Eq. ( 14). Conversely, the evaluation in the spectral domain as in Eqs. ( 15) and ( 18) does not suffer from this phenomenon.

Modal analysis of an acoustic resonator

The method described in the previous section is now applied to the input impedance of an acoustic resonator, first in the case where the impedance is known over the full frequency range (Sec. 5.1), and then when the measurements needs to be preprocessed (Sec. 5.2). The passivity over the stop band is constrained using fictitious points as described in Sec. 5.3, before an example of application is given in Sec. 5.4.

Analysis of an input impedance known over the pass band

When the input impedance is sampled from DC to the desired sampling frequency, we apply the LSCE method to the shifted variant of the impedance H = Z in -1 in agreement with the requirement of a transfer function decreasing to 0 when the frequency increases. The procedure is as follows:

It is also possible to apply the procedure to the reflection coefficient R without the shift (steps 2 to 7 and 9). Even if the reflection coefficient does not have visible resonances, the method can parametrize it as a modal series (as in Eq. ( 9)), generally with highly damped poles. As this is a smoother function of the frequency, the matrix B should be better conditioned. Furthermore, a better balance is maintained between impedance and admittance peaks than when applying LSCE to the shifted input impedance (which attaches more importance to the impedance peaks). However, there is a drawback as the results must be converted back to modal parameters for the impedance. This is done by searching the roots of the high order polynomial R -1 (see Eq. ( 2)) which requires high precision arithmetics again.

Analysis of an input impedance with missing data in the low frequency range

The case of measured input impedances requires special care because the measurements may be noisy and/or available and reliable only on a partial range of the frequency domain. The proposed procedure is as follows:

Extrapolation steps. It is known that the input impedance measurements are have limited validity at very low frequencies and the LSCE analysis appears to be very sensitive to the physical plausibility of the reconstruction in the low range. Therefore an extrapolation is required below the lower bound F min of the acceptable data. According to our experience, the reflection coefficient is a smoother function than the input impedance, and is the most suitable for its modulus and its argument to be approximated by low order polynomials.

1. Compute the reflection coefficient R meas from the measured input impedance, using Eq. (2). 2. Select a frequency range [F min , F max ] where R meas has an acceptable signal-to-noise ratio (typically from 120 to 230 Hz, depending on the signal-to-noise ratio of the measurement). 3. Fit an even polynomial q even (ω) to the modulus |R meas | on the selected frequency range. The order is typically set to 10, with the lowest coefficient set to 1 to enforce an input impedance vanishing at zero frequency. 4. Fit an odd polynomial q odd (ω), typically of order 11, to the unwrapped argument of (-R meas ).

Merging measured and extrapolated data. The two previous steps lead to the function R extrap (ω) = -q even (ω)e jq odd (ω) [START_REF] Chaigne | Acoustics of musical instruments[END_REF] that is hermitian symmetric and is used in the [0, F min ] frequency range. In order to overcome noisy measurements in the pass band and a possible too high frequency resolution (i.e., ∆ω too small, leading to huge dimensions of the matrix B ), a smoothing downsampling is desirable in the pass band.

5. The reflection coefficient is reconstructed in the low frequency range as follows:

R(ω) =            R extrap (ω) for 0 ≤ ω < 2πF min ; crossfade (R extrap , R meas ) for 2πF min ≤ ω < 2πF max ; R meas (ω) for 2πF max ≤ ω ≤ 2πF cut ; R(-ω) * for ω < 0; ( 20
)
where the crossfade is performed using the sigmoid function:

crossfade (R extrap , R meas ) = R extrap + R meas -R extrap 1 + exp [-2a (ω -π(F min + F max ))] . ( 21
)
a is chosen so that the sigmoid equals 1% and 99% at 2πF min and 2πF max , respectively.

Downsampling data.

6. Define a downsampling factor d s so that the frequency resolution becomes typically ∆ω/2π = 0.2 Hz after downsampling. 7. The smoothing downsampling filtering is applied on the samples of R(ω) in the pass band: consider 4d s + 1 samples, fit a low-order polynomial (typically order 3) on these samples and retain the evaluation of the polynomial at the center frequency of the frame. Advance d s samples, and repeat the procedure. It can be described as a moving average process applying on frequency-domain frames of length 4d s + 1 with a 3d s + 1 overlap.

Modal analysis on the reconstructed data.

8. Compute Z in using Eq. ( 2). 9. Apply Sec. 5.1 on Z in .

Modal analysis can again be performed on the reflection coefficient R with the drawback mentioned above. In our experience, the analysis of measured reflection coefficient requires an additional pole y 0 = 0 as for the input impedance (step 6 in Sec. 5.1) even without any shift of the data.

Choice of the fictitious points

In most cases, this issue is not complicated to solve for clarinet-like instruments. According to our experience, the fictitious points in the stop band should be spaced about 20∆ω to 200∆ω apart from another.

A crude initialization is Ĥ(0) n = 0, for which the least squares solution of Eq. ( 18) usually is a convenient solution C (0) . The stop band constraint is weakened in a second step by choosing new fictitious values Ĥ( 1)

n from the last 2N values of B • C (0)
. Again solve Eq. ( 18) for vector C (1) , this time with source term H (1) . A known difficult case is that of brass instruments measured with their mouthpieces. Because of the small cavity of the mouthpiece, the damping of the acoustic waves is low and R exhibits a modulus near to unity, even at high frequency. Thus more care is required for the initialization in order to preserve the passivity in the stop band. Notice that B -1 has to be computed only once, which means that many "guesses" can quickly be tested in a trial and error process. Providing a general methodology for such a search is beyond the scope of this paper.

Example of application

The input impedance of the fingering F 3 of a professional clarinet was measured with a CTTM impedance head [START_REF] Roux | A new impedance sensor for industrial applications[END_REF]. The measurements are processed according to Sec. 5.2. Typical values provided in previous sections were adopted for this analysis.

Fig. 4 shows the reconstruction of the reflection coefficient in the low frequency range. Fig. 5 depicts the result of the modal analysis performed on the reconstructed data, according to Sec. 5.1. The modulus of the reconstructed impedance deviates from the measurement with a standard error of 0.19 dB, the maximal error being obtained for the first antiresonance (1.5 dB local error). This is probably related to a lower signal-to-noise ratio of the measurements at the deepest antiresonance as also visible in Fig. 4. Concerning the argument, the standard deviation is 0.02 rad with a maximal error of 0.1 rad at the first and third antiresonance. The corresponding reflection coefficient R is illustrated in Fig. 6, including fictitious points, extrapolated up to F s /2. Figure 4: Reconstruction of the reflection coefficient in the low frequency range. Modulus of the reflection coefficient (top) and group delay (bottom) of the measured data (dots), the extrapolated data (thick dashed line), the reconstructed downsampled data (thin line) and of the result of the modal analysis (dash-dotted line). The bounds of the frequency range used for the extrapolation step are also shown (vertical dotted lines).

Application to sound synthesis

Principles of modal synthesis

A direct application of modal analysis of the input impedance (or reflection coefficient) of wind instruments lies in the design of algorithms for real-time sound synthesis. In fact, the functioning of wind instruments is classically viewed as the coupling of an acoustic resonator and a nonlinear excitation system, and this physics-based modeling paradigm can be used for sound synthesis (see Ref. [START_REF] Välimäki | Discrete-time modelling of musical instruments[END_REF] for a review). The numerical scheme then relies on two discrete-time models, one for the exciter and one for the acoustic resonator.

Accounting for the bore, the transposition of the continuous time model to the discrete time is performed in order to preserve the impulse response. The continuous time impulse response h(t) is obtained by the inverse Fourier transform of the modal series in Eq. ( 5) and the discrete time impulse response h[n] equals h(t) at the positive sampled times

h(t ≥ 0) = m C m e smt ⇒ h[n ≥ 0] = h(nT s ) = m C m e nsmTs = m C m x n m (22)
where x m = e smTs and T s = 1/F s are related to the sampling frequency F s . The latter is in general higher than the cutoff frequency F cut . Typical values for F s and F cut are 4 kHz and 44.1 kHz, respectively. The z-domain impedance is then n of the fictitious points determined by the least squares fit are depicted with red dots. These points were set to 0 at the beginning of the procedure.

H(z) = m C m 1 -x m z -1 .
Because h(t) is real, the coefficients s m must be either real (monopole) or complex conjugates (bipole), in order to respect the Hermitian symmetry. The same applies to the discrete-time response h[n] and the coefficients x m . Nevertheless, since a monopole and a discrete Dirac are degenerate bipoles, we can use the generic 2nd order IIR digital filter [START_REF] Smith | Introduction to Digital Filters with Audio Applications[END_REF]

H(z) = m Im[sm]≥0 H m (z) with H m (z) = b m,0 + b m,1 z -1 1 + a m,1 z -1 + a m,2 z -2 (23)
where the values of the coefficients b m,0 , b m,1 , a m,1 and a m,2 are computed from the poles s m and the modal coefficients C m according to Table 1. This filter design combines the contributions of complex conjugate poles s m and s * m into a single bipole. It is also important to note that the time steps used in measurements (T in Eq. ( 10)) and in the synthesis (T s in Eq. ( 22)) usually differ. This implies that the roots y m = exp (s m T ) of the former have to be converted into x m = exp (s m T s ). There may exist roots y m on the negative real axis (Re[y m ] < 0 and Im[y m ] = 0, i.e., a contribution oscillating at the cutoff frequency), and they have to be split into a pair of complex conjugate values x m and x * m , increasing the number M of poles. As a consequence, when applying the procedure described in Sec. 5, the modal coefficients C m must be evaluated (step 7. in Sec. 5.1 after converting y m into x m .

Type Root b m,0 b m,1 a m,1 a m,2 Dirac x 0 = 0 C 0 0 0 0 Monopole x m > 0 C m 0 -x m 0 Bipole (x m , x * m ) 2 Re[C m ] -2 Re[C m x * m ] -2 Re[x m ] |x m | 2
The filter H(z) enables the evaluation of the pressure p[n] at the input of the waveguide when it is excited by the flow rate u[n]. The pressure is the sum of the partial pressures p m [n], each one accounting for a real pole s m or a pair of complex conjugate poles (s m , s * m ) and obeying the difference equation

p m [n] = b m,0 u[n] + b m,1 u[n -1] -a m,1 p m [n -1] -a m,2 p m [n -2]. (24)
It is possible to arrange the relations so that the (total) pressure at the input of the bore at the current time depends on the current flow rate value and on some previous values

p[n] = m p m [n] = V 1 u[n] + V 2 with      V 1 = m b m,0 , V 2 = m b m,1 u[n -1] -a m,1 p m [n -1] -a m,2 p m [n -2]. (25)
This representation of the pressure at the input of the resonator is more convenient to couple with the classical representation of the exciter. The latter is generally described by a nonlinear time-domain relationship

u[n] = F (p[n]
) that may account for unsteady effects. In the general case, the coupling would require an iterative solution to the following problem

p[n] = V 1 u[n] + V 2 and u[n] = F (p[n]) (26)
at each time step n. However, some discrete-time models of the exciter, such as the one given in Appendix A, have an explicit solution. The current pressure p[n] and flow rate u[n] can be efficiently obtained and real-time implementation is achievable.

Transitions between fingerings

In real instruments, when the instrumentalist changes the pitch from some note A to another note B, he opens and/or closes one or many holes at the same time. According to Ref. [START_REF] Terroir | A simple dynamic tonehole model for real-time synthesis of clarinet-like instruments[END_REF], a simple crossfade model between the input impedances at the beginning and at the end of the transition is shown to lead to perceptually satisfactory simulations of the transitions between fingerings. This model is only required to be able to simulate the response of the bore for the initial and the final states during the N r time steps of the transition. The partial pressures p m,A and p m,B are computed simultaneously for the two configurations. At time step n ∈ [N 0 , N 0 + N r ], the pressure p[n] is the result of the linear crossfade

p[n] = m p m,A [n] + n -N 0 N r (p m,B [n] -p m,A [n]) . ( 27
)
as if the two pipes were simultaneously excited with the same flow rate. This can still be written as p

[n] = V 1 u[n] + V 2 .
Consequently the synthesis scheme remains unmodified. At the beginning of the transition, the pipe B is assumed to be silent (i.e. p m [N 0 -1] = 0 and p m [N 0 -2] = 0 for all m). Coupling through the flow transfers some energy from the pipe A to the pipe B.

We tested the refinement proposed by Guillemain and Terroir [START_REF] Terroir | A simple dynamic tonehole model for real-time synthesis of clarinet-like instruments[END_REF] on the basis of the input impedances of a clarinet measured with a tone hole being progressively closed. Their measurements show that the first resonance frequency varies gradually, like a glissando, while the amplitude of the peak decreases to a minimal value before increasing to the final value, approximately following the shape of a parabola. This behavior can be simulated in real-time for the most important peaks, such as the one that sustains the oscillation (so-called master mode, indexed as m) and, optionally, the peaks that are strongly excited in the playing situation. The modal coefficient C m follows a parabola-like evolution between the initial value C m,A and the final value C m,B while the pole s m varies linearly from s m,A to s m,B . The transition for less important peaks is still simulated using the simple crossfading algorithm at a lower computational cost. This algorithm is acceptable for simulating relatively small musical intervals (e.g., opening a single tone hole, as in Ref. [START_REF] Terroir | A simple dynamic tonehole model for real-time synthesis of clarinet-like instruments[END_REF]). During big jumps, the glissando-like transition sometimes sounds strange [START_REF] Smith | Pitch bending and glissandi on the clarinet: Roles of the vocal tract and partial tone hole closure[END_REF]. Such transitions are known by clarinetists to be difficult to achieve with a good legato. In some cases, different playing techniques are used to overcome this difficulty on the real instrument (such as: opening some holes earlier than others; half closed holes; glissando-like opening and closing of some holes working acoustically as "speaker keys"; subtle adaptations of the embouchure and the air pressure during the jump; special tuning of some resonances of the vocal tract; modification of the glottis using the crico-thyroid muscle; and others, according to the practice of the first author as a professional clarinetist). However, the MIDI wind controllers used to interface with the synthesizer do not allow such subtle controls.

A pragmatic solution consists of injecting some energy in the mode m of note B in order to initiate a bifurcation. At the beginning of the transition, all modes of the pipe B are silent, except mode m whose history is set to p

m [N 0 -1] = c jump p[N 0 -1] and p m [N 0 -2] = c jump p[N 0 -2].
The constant c jump controls the amount of injected energy. With c jump = 0, the simulated instrument sounds like the playing of an inexperienced beginner unable to control the jumps (but the legato between neighboring notes sounds good), whereas with c jump = 1 the jumps are much more secure (but the legato sounds strange because every note starts with a small accent). The good balance is about 0.4 < c jump < 0.7, depending on the piece being played.

The second concern deals with the computational cost, a critical factor for real-time implementations when the number of modes M is large. The simulation of low pitched wind instruments is particularly sensitive to this problem. In Eq. ( 27), both pipes are played in parallel. Therefore, the computational cost doubles during the transition. This issue is resolved using the algorithm depicted in Fig. 7. Prior to the simulation, the modes of each fingering are classified as primary modes (from 0 up to about M/2) and secondary modes (from about M/2 + 1 up to M ). The contribution of the latter to the input impedance is marginal, while the primary modes have either a big amplitude or are necessary to ensure the passivity of the filter. The transition begins with a decrescendo al niente of the secondary modes of fingering A, followed by the transition Eq. ( 27) applied to the primary modes of both fingerings and ends with a crescendo of the secondary modes of fingering B. The decrescendo and the crescendo are realized by multiplying the partial pressure of the secondary modes by a coefficient progressively varying respectively from 1 to 0 and from 0 to 1. Typically, the duration of each phase of the transition is 2 ms, 20 ms and 2 ms, respectively. With this algorithm, the computing cost remains approximately constant at every time step.

An example: sound synthesis of a clarinet

This section shows how the modal analysis performed on measured input impedances of a real instrument can be used to achieve the real-time synthesis of a musically relevant illustration, an excerpt of the Brahms Sonata Op. 120 no. 2 (see Fig. 8a).

The input impedance of a professional clarinet has been measured for each fingering used in the excerpt using the device described in Ref. [START_REF] Roux | A new impedance sensor for industrial applications[END_REF], and the data was processed according to Sec. 5.2. The motion of the reed tip (see Ref. [START_REF] Taillard | Comparison of two clarinets separated by 200 years of evolution (original title: Comparaison de deux clarinettes séparées par deux cent ans dévolution : tentative dhybridation amusante et instructive entre facture instrumentale, modèles physiques et synthèse sonore). 13ème Congrès Français d'Acoustique[END_REF]) and the pressure in the mouth and in the mouthpiece were then measured while an experienced musician played the excerpt on the same instrument. These signals are used to provide the rough estimation of the control parameters γ and ζ, as defined in Appendix. A. This is known to be a difficult task as many quantities involved in the computation of the control parameters are still unavailable (see, e.g., Ref. [START_REF] Perichon | Control parameters inversion using genetic algorithms applied to numerical impedance synthesis for woodwinds[END_REF] for clarinet and Refs. [34][35] for brass instruments). However, their time evolutions are estimated using the short-term average (10 ms frames) of the measured signals and then arbitrarily scaled to the dimensionless control parameters γ and ζ shown in Fig. 8b. This simplistic procedure should provide acceptable estimates for steady-state oscillations, but decreases the bandwidth of the signals which results in slower attack transients. The timings of the transition between fingerings had also to be retrieved from the recorded sound, which may introduce some mismatch between the measured and the synthesized sounds.

For real instruments, the transfer function between the mouthpiece pressure and the external pressure is different for each fingering. This would introduce another series of modal coefficients, thereby increasing the memory storage requirements and computing time. Appendix B gives the details of the design of an empirical transfer function. It is applied to the mouthpiece pressure measured on the real instrument and on the synthesized mouthpiece pressure signal, in order to produce comparable external sound estimates. The real-time implementation in Max/MSP of the sound synthesis scheme (accounting for the resonator, the exciter and the radiation) requires no more than 25% of CPU time on a modern personal computer. Many aspects of the simulations are satisfactory (see Fig. 8, and sound files available as supplementary material 1 ). The proposed algorithm performs the excerpt noticeably well, without any unwanted squeaks, wrong notes, or absent notes, despite the difficult bindings and jumps in the excerpt. The spectral content of the synthesized external sound (Fig. 8d) is also quite similar to the real one (Fig. 8c).

Fig. 8e provides some insights on the oscillation frequency and the deviation in cents. The musician performance has a mean deviation close to zero (0.5 ± 10 cents), with smooth variations of the playing frequency within almost all notes. The modal-based sound synthesis exhibits a higher overall deviation (15 ± 11 cents) without modulation within the notes. This has implications for the expressiveness of the musical performance [START_REF] Barthet | The effect of timbre in clarinet interpretation[END_REF] and relates to the artificial nature of the sound that can be perceived by experienced listeners. This is also noticeable in Fig. 8f, where the spectral centroid of the synthesized signal exhibits the same note-to-note trends as the reference, but lacks the variations within the notes.

The synthesis inherits problems relating to measurement inaccuracies such as the difference in temperature between the measurement and playing situations, or the use of an adapter that replaces the mouthpiece for the the input impedance measurements. These elements, in addition to the approximate control parameters, could explain the differences between the oscillation frequencies of the original and the synthesized sounds. The simulation of the nonlinear behavior of the tone holes, a more realistic model of the reed, the inclusion of a vocal tract and the treatment of the sound radiation will likely make the synthesized model more closely resemble the real one.

Conclusions

This paper investigates many aspects related to the modal representation of the input impedance of wind instruments. The LSCE method is a powerful tool in this context, but a careful treatment of its application to measured or analytical input impedance is necessary. The possible issues include: physically acceptable reconstruction of the impedance in the low range, implementation of fictitious points in the spectral domain ensuring the passivity of the digital filter at every frequency, use of high precision arithmetics depending of the required number of simulated modes, for instance. The described method allows to derive valuable digital filters for real-time synthesis. Though the usefulness of the method is only demonstrated on a clarinet model, the authors have investigated the relevance of the modal analysis on other wind instruments, notably on saxophone, french horn, trumpet and trombone.

4

 4

Figure 1 :

 1 Figure 1: Modulus and argument of the input impedance Z (top left and bottom left, respectively), modulus of the reflection coefficient R (top right) and group delay τg = -d(arg (R))/dω (bottom right) from the analytical model (Eq. (4)). = 1 m, r = 1 cm.

Figure 2 :

 2 Figure 2: Location of the poles sm of the input impedance of a cylindrical pipe (r = 1 cm and = 1 m) with radiation and boundary layer losses in the complex plane.

6

 6

Figure 3 :

 3 Figure 3: Modulus and argument of the input impedance Z (top left and bottom left, respectively), modulus of the reflection coefficient R (top right) and group delay τg = -d(arg (R))/dω (bottom right) from the analytical model (Eq. (4), thin black line, same as Fig.1) and from the modal expansion truncated to 20 modes (see Sec. 3, Eq. (7), thick red line). = 1 m, r = 1 cm.

Figure 5 :

 5 Figure 5: Input impedance Zin of the fingering F 3 of a professional clarinet (top: modulus, bottom: argument). The raw measurements are shown by the red points (∆w = 2π × 0.047 rad s -1) and the results of the modal analysis as used for sound synthesis (following the procedure in Sec. 5.2) are denoted by the thin dashed-dotted line.

|R|Figure 6 :

 6 Figure 6: Modulus of the reflection coefficient of the fingering F 3 from the modal analysis, including the extrapolation in the stop band necessary for the sound synthesis, up to Fs/2. The final values Ĥ(1)n of the fictitious points determined by the least squares fit are depicted with red dots. These points were set to 0 at the beginning of the procedure.

Figure 7 :

 7 Figure 7: Diagram illustrating the processing of the primary and secondary modes of the initial (A) and final (B) configurations during a transition between two fingerings.

Figure 8 :

 8 Figure 8: a: Excerpt of the first edition (1895) of Brahms Sonata Op 120 no. 2 that is played on a professional clarinet made in 1981. b: Control parameters γ and ζ estimated from measurements of the mouth pressure and of the reed motion (with an optical sensor) while playing the excerpt. c-d: Spectrograms of the external pressure obtained by applying the radiation filter (see Appendix B) to the mouthpiece pressure (c) measured on the real instrument or (d) produced by the synthesis model described in Sec. 6 using the estimated signals γ and ζ. e-f: Deviation from the tempered scale (e) and spectral centroid (f) calculated on the reconstructed external pressures for the real instrument (black thin markers) and for the synthesis (red thick markers). The parameters used in the synthesis are ωr = 2π × 1340 rad s -1 , qr = 0.3, and cjump = 0.7.

Table 1 :

 1 Coefficients of the digital filters with respect to xm and Cm.

Note to compositor: please add here a link to supplementary material

Acknowledgments

The authors thank Joël Bensoam, Louis Delebecque, Christophe Vergez for useful discussions, Anneke Scott and Erik Petersen for proofreading, and the He-Arc Ingénierie (Neuchâtel, Switzerland) for the facilities granted. This work is supported by the French ANR research project CAGIMA (ANR-11-BS09-0022).

Appendix A. Excitation model used in the sound synthesis scheme

The reed displacement is modeled as a single degree of freedom oscillator with natural angular frequency ω r and damping q r , driven by the dimensionless difference ∆p(t) = γ(t) -p(t) between the mouth pressure γ(t) and the mouthpiece pressure p(t)

when the reed channel is open, and with an additional unilateral contact force ensuring

when the channel is closed (y ≤ -1). The flow rate which excites the acoustic resonator takes into account the Bernoulli flow u b (t), proportional to a dimensionless embouchure parameter ζ(t) (related to the maximum flow rate that can enter the mouthpiece, see, e.g., Ref. [START_REF] Chaigne | Acoustics of musical instruments[END_REF]), as well as the flow rate u r (t) induced by the reed motion (with a factor λ)

and

A discrete implementation is proposed by Coyle et al [START_REF] Coyle | Predicting playing frequencies for clarinets: a comparison between numerical simulations and simplified analytical formulas[END_REF]. It uses a finite difference method for the reed displacement and the induced flow rate that are explicit:

) and ν = q r /(2T s ω r ). In Eq. (A.4), the current displacement y[n] is independent of p[n]. This implies that the current Bernoulli flow u b [n] and pressure p[n] are related by

where

) is the smoothed channel opening accounting for the regularization function

where ε is typically set to 0.04. This function is similar to the one provided in Ref. [START_REF] Kergomard | Idealized digital models for conical reed instruments, with focus on the internal pressure waveform[END_REF] (Sec. V.E) and operates when the reed beats against the lay. It follows that the Bernoulli flow u b [n] solves a simple 2nd order polynomial such that

and, finally, Eq. (A.5) can be used to compute p[n].

Appendix B. Filter used to estimate the external pressure

Sound radiation is certainly important for the musical quality of an instrument. Until now, this aspect was largely neglected by the musical acoustics community due to the complexity of the radiated field and the fact that it exerts only a minor influence on the functioning of the instrument. The literature offers expressions for the transfer function between the velocity at the open end of a waveguide and the radiated pressure in free space only for academic configurations: the plane and spherical wave approximations (see, e.g., Refs. [START_REF] Chaigne | Acoustics of musical instruments[END_REF] and [START_REF] Th | Radiation of a pulsating portion of a sphere: application to horn radiation[END_REF], respectively), the radiation of a semi-infinite circular pipe that is either unflanged [START_REF] Levine | On the radiation of sound from an unflanged circular pipe[END_REF], infinitely flanged [START_REF] Norris | Acoustic radiation from a circular pipe with an infinite flange[END_REF], or partially flanged [START_REF] Dalmont | Radiation impedance of tubes with different flanges: numerical and experimental investigations[END_REF]. The case of waveguides with several orifices is even more complex as the orifices radiate in a common space and external interaction can then not be ignored. Ref. [START_REF] Chaigne | Acoustics of musical instruments[END_REF] (Sec. 14.5) explains the role of the tone hole lattice and Ref. [START_REF] Lefebvre | External tonehole interactions in woodwind instruments[END_REF] provides an illustration of real instruments with finite elements computations.

In the context of real-time synthesis, we seek a very simplified model of radiation characterized by low computational cost; a single model common to all the fingerings; and one that creates perceptually realistic sounds. Simple models of radiation are available in the literature, but fail to satisfy the last requirement.

We made an attempt to derive such a model from experimental data: for all the fingerings of a professional clarinet, we measured the internal pressure in the mouthpiece p(t) and the external pressure p ext (t) on a microphone situated about 1.5 m in front of the player, in the practice room of a musician (i.e., not in an anechoic room [START_REF] Benade | From instrument to ear in a room: direct or via recording[END_REF]).

The signals p(t) and p ext (t) were averaged over all fingerings. The global transfer function was computed by dividing the discrete Fourier transforms of both averaged signals. The phase was eliminated in the computation of the impulse response by taking the discrete inverse Fourier transform of the modulus of the transfer function. The human ear is known to be practically insensitive to the phase [START_REF] Plomp | Effect of phase on the timbre of complex tones[END_REF]