Amplitude-dependent modal coefficients accounting for localized nonlinear losses in a time-domain integration of woodwind model

This article develops the design of a sound synthesis model of a woodwind instrument by modal decomposition of the input impedance, taking into account viscothermal losses as well as localized nonlinear losses at the end of the resonator. This formalism has already been applied by Diab et al. (2022) to the study of forced systems. It is now implemented for self-oscillating systems. The employed method extends the denition of the input impedance to the nonlinear domain by adding a dependance on the RMS acoustic velocity at a geometric discontinuity. The poles and residuals resulting from the modal decomposition are interpolated as a function of this velocity. Thus, the pressure-ow relation dened by the resonator is completed by new equations which account for the dependence with the velocity at the end of the tube. To assess the ability of the model to reproduce a real phenomenon, comparisons with the experimental results of Atig et al. (2004) and Dalmont et al. (2007) were carried out. Simulations show that the model reproduces these experimental results qualitatively and quantitatively.


Introduction
Sound synthesis by modal decomposition is based on an input impedance measurement, which captures the passive acoustical response of a real instrument. This method thus has the advantage of capturing the acoustical subtleties that differentiate two clarinets by means of a single measurement. Other methods, such as delay lines [1], waveguides [2,3] or spatial discretization [4,5] would require more effort on the geometrical description of the resonator to render such acoustical subtleties. Furthermore, sound synthesis by modal decomposition requires little RAM compared to such other methods, which can be beneficial for embedded devices. Indeed, in modal decomposition synthesis, the temporal integration scheme requires to keep only a few previous iterations (precisely two for [1]) to compute a new one. For the delay line synthesis, it is necessary to keep in memory all the temporal iterations during a round trip, i.e. during 2LF s /c 0 iterations for an academic cylindrical resonator without lateral holes. Although this is not a problem for most embedded processors, fine modeling of viscothermal losses in the waveguide formalism also requires additional computations involving fractional derivatives [6]. Finally, the ODE formalism is particularly well suited to the bifurcation analysis, which is not the case of the PDE [6] or of the DDE formalisms [7].
The input impedance is a measure of the linear frequency response of the resonator, for a low-amplitude excitation. However, in the case of woodwind instrument playing, the measured acoustic pressure and velocity in the resonator can be very high. For instance, [8] measure a pressure in the mouthpiece of 4 kPa, i.e. 163 dB for a monochromatic wave. According to [9] of Chap. 8.4.5, from an acoustic speed of roughly 1 m/s, jet separation phenomena appear. At the level of a geometrical discontinuity, such as the open end of a pipe [10] or a lateral orifice [11], vortices are observed. A part of the kinetic energy of the jet is absorbed by these vortices and is dissipated as heat by friction. These losses are accounted for through a nonlinear relationship derived from Bernoulli's law, as demonstrated by the implementation in waveguide modeling or delay lines of Refs. [7,[12][13][14].
Following [11], the works from Atig et al. [12,15] and Dalmont et al. [13,16] model localized nonlinear losses through a resistive impedance Z t . This impedance is in series with the radiation impedance Z R , and depends on the amplitude of the acoustic velocity v 0 at the location of the discontinuity: where Z c = q 0 c 0 /S is the characteristic impedance of the medium for plane waves, and c d 2 [0; 3] is a parameter depending on the geometry of the termination. The larger the radius of curvature at the output, the smaller this coefficient is. The coefficient c d is difficult to predict theoretically, and easier to determine experimentally [12]. As shown by Atig et al. [12], nonlinear losses have a significant influence on the playing range of a clarinet, hence models of sound production in woodwinds should include such effects. Models of woodwind instruments taking into account localized nonlinear losses at the end of the resonator have never been developped in the framework of modal synthesis. This is the main contribution of this work.
First, a method proposed by Diab et al. [17] to account for nonlinear losses by modal decomposition is presented. It is then applied to a cylindrical tube with nonlinear losses located at the open end. From this resonator, a selfoscillating clarinet-like system is defined and simulated by time integration. These simulations are finally compared to the experimental results published by Atig et al. [12] and Dalmont and Frappé [13].

Method of modal decomposition accounting for localized nonlinear losses
The definition of a "nonlinear impedance" (i.e. adjusting the impedance so that the link between acoustic flow and pressure stays valid in nonlinear conditions) opens the way to consider localized nonlinear losses by modal decomposition. In Diab et al. [17], the evolution of the surface impedance of perforated plates is computed by temporal simulation for a broadband excitation, with respect to the RMS amplitude of the acoustic velocity at the level of the hole, noted v RMS . From the decomposition of the impedance in the linear domain into a sum of N modes, two methods are studied for increasing values of v RMS : the interpolation of the impedance, on one hand, and the regression of the poles and complex residues (s n , C n ) on the other hand. This second method will be considered hereafter.
Nonlinear effects are taken into account by allowing the poles and residues to vary with respect to v RMS . In the Laplace domain, the modal decomposition of the input impedance is written: where s is the Laplace variable and * denotes the complex conjugation operation. The relationship between poles and residues and v RMS is assumed to be a rational fraction [17]. In the present work, the approximation is limited to a polynomial of degree N p : where the index (0) refers to the linear part of the modal coefficient. In Diab et al. [17], the regression coefficients of the poles and residues are then adjusted on a nonlinear surface impedance model Z NL (v RMS ), given by Laly et al. [18]. It is worth noting that in the present paper, the poles and residues are fitted on the dimensioned quantity of the RMS velocity. Therefore, the unit of coefficients C ðkÞ n and s ðkÞ n are expressed in rad Á s À1 Á ½m Á s À1 Àk . The last step consists in computing v RMS , given by: where v 0 (t) = u 0 /S is the acoustic velocity at the discontinuity, of cross section S. This definition is rather impractical in the framework of a time-domain simulation. In order to obtain a state formalism, it is replaced by: which is equivalent. By considering the right-hand side of equation (6) as an input for the ODE solver, tv 2 RMS can be computed, hence v RMS .

Application to a cylindrical tube with nonlinear losses at the open end
The impedance regression technique presented by Diab et al. [17] is applied to the case of a closed-open cylinder of radius R and length L, which is a similar case study to Ref. [12]. Viscothermal losses are taken into account in the propagation: where g = 3 Â 10 À5 s 1/2 . The boundary condition at x = L is given by the radiation impedance Z R , which includes in series the nonlinear impedance Z t defined by equation (1): and jk = s/c 0 , Dl % 0.6R ( [9], Chap. 12.6.1.3), K NL = 4c d /(3p). In the following, v RMS (L, t) will be denoted as v RMS for reading comfort. The input impedance is defined by: introducing the dimensionless notation z = Z /Z c . A linear input impedance z ðlinÞ in is also defined, such that which is equivalent to z in for v RMS = 0. Figure 1 illustrates the evolution of z in (computed through equation (9) and equation (10), in which the Laplace variable is substituted by j2pf ) as a function of the RMS velocity at the open end. It can be observed that taking into account nonlinear losses has a consequence on the amplitude, but not on the frequency of the resonance peaks. Figure 1b shows that this impact is particularly accentuated on the first resonance peak: its amplitude is 43% lower for v RMS = 24 m/s than in the case without nonlinear losses. The difference is reduced to 30% for the second peak, and to 16% for the sixth peak. Near the antiresonances, the negative relative discrepancies are high, since the impedance values are all close to zero. The absolute deviation remains very low. From the lack of frequency variation of the resonance peaks with v RMS , we could anticipate that nonlinear losses at the end of the pipe would rather impact dynamics than intonation.

Modal decomposition of the nonlinear input impedance
In order to take into account nonlinear losses in temporal simulation, the modal coefficients of the input impedance must be determined. In the same way that the nonlinear input impedance depends on the acoustic velocity at the open end, the coefficients C n and s n derived from the modal decomposition of z in also depend on v RMS . To formulate the relationship between the modal coefficients and v RMS , a minor adaptation of the method presented in Ref. [7] of Chap. 5.5.3 is used. The definition of z R used in [9] is substituted by the expression given by equation (8), taking into account nonlinear losses at the end of the tube.

Expression of the poles s n
Equation (9) can also be written as: The poles s n ðv RMS Þ are the solutions of cosh[CL + h] = 0, i.e.: The solutions of equation (12) give s n , for different input values of v RMS . They are plotted on Figure 2a. The value of c d has been set to 13/9, which corresponds to a radius of curvature of 0.3 mm, according to [29]. From top to bottom, from left to right: modulus of the input impedance; relative gap between nonlinear and linear definition of z in ; detail view on the first peak; detail view on the fourth peak.
The consideration of nonlinear losses at the end of a cylindrical pipe has almost no influence on Iðs n Þ, which is the resonance angular frequency of peak n. However, a larger v RMS increases jRðs n Þj. This shift remains almost the same regardless of the index of the pole. However, in terms of relative deviation, the ratio jRðs n À s ðlinÞ n Þ=Iðs n Þj is decreasing as n increases. This is signaled by the amplitude of peaks of higher index appearing less affected by nonlinear losses, as shown in Figure 1b. 3.1.2 Expression of the residues C n It remains to calculate residue C n (v RMS ). In the vicinity of a pole s n , the denominator D of z in can be written through a first-order series expansion of cosh: introducing notation 0 = o/os. According to Chaigne and Kergomard ( [7], Chap. 5.5.3), after application of the residues theorem, the modal decomposition of the input impedance for the nth peak can be written as The evolution of coefficients C n with respect to v RMS is plotted on Figure 2b. Real and imaginary parts of C n slightly decrease when v RMS increases. Following Diab et al. [17], it remains to fit the modal coefficients with respect to v RMS , as in the example of equation (3). Figure 3a shows the mean relative fitting error for polynomials of different degrees, for the same data as in Figure 2. For linear regression, the mean relative error stays around 10 À7 for each mode n, both for C n and s n . Further simulations reveal that choosing an excessive values of c d tends to increase this error. For a tube with sharp edges at the open end, Atig et al. [12] estimated a maximum value of c d = 2.8. For c d = 5, for instance, the mean relative error of linear regression is around 10 À5 . Although this error has increased, it remains very low.

Polynomial fitting of the modal coefficients
Moreover, Figure 3b shows that choosing an polynomial regression degree larger than 1 has no consequence on the error over z in between modal decomposition with fitted coefficients and the definition given by equation (9). This error is therefore inherent to the modal decomposition approximation (Eq. (2)) of the input impedance (Eq. (9)).
According to these previous results, in the rest of this article, the modal coefficients will both be fitted by a polynomial of order 1.

Application to the time-integration simulation of a clarinet-like system
The cylindrical tube considered in Section 3 is now treated as the resonator of a clarinet-like instrument. The physical model governing the self-sustained oscillations of the instrument is presented in Section 4.1. It is then simulated by time integration, imposing a linear increase of the blowing pressure (crescendo) followed by a linear decrease (diminuendo). Details and results from the simulations are given in Section 4.2.  [29]. As in Figure 1, c d = 13/9. (a) Real and imaginary parts of poles s n (v RMS ), calculated by solving equation (12). Detailed view on s 1 . (b) Real and imaginary parts of residues C n (v RMS) , calculated with equation (14). Detailed views on C 1 and C 4 .

Physical model for the self-sustained oscillations
The self-sustained system consists of three main parts: reed dynamics, reed channel and resonator. This last subsystem is enriched compared to the classical modal decomposition formalism [19][20][21] by the addition of nonlinear losses at the end of the pipe.

Reed dynamics
The dynamic behavior of the reed is modeled for its dimensionless displacement x(t) by the following equation: where x r = 2p Â 2200 rad/s is the reed resonance angular frequency [22] and q r = 0.4 is the reed damping [23]. c(t) = P m (t)/p M is the dimensionless blowing parameter and p(t) = p dim (t)/p M is the dimensionless pressure at the input of the resonator. The beating pressure p M is the value of the blowing pressure P m for which the reed closes the reed channel, in quasistatic regime. F c ðx; _ xÞ is the contact force function. In the following, a model of "ghost reed" [24] is considered, i.e. F c = 0.

Reed channel
The characteristic of the flow u(t) through the reed channel is described by the following equation, according to [25]: where k = 5.5 Â 10 À3 /c 0 is the reed flow parameter [26], f is the embouchure parameter, and c(t) is the dimensionless blowing pressure. The operator [] + = (+||)/2 refers to the positive part function. For the simulations, the absolute value functions are regularized to smooth the irregularities of the quasi-static flow characteristic, without altering its behavior, according to [27]. Thus, j j7 ! ffiffiffiffiffiffiffiffiffiffiffiffi ffi 2 þ g p , with g = 0.001.

Resonator
The resonator is described in the frequency domain by its input impedance z in , under its modal decomposition form given by equation (2). The modal poles and residues s n and C n are chosen to be linearly dependent on v RMS , i.e.: The modal decomposition of z in allows to write the relation between the dimensionless pressure p and flow u in the temporal domain: _ p n ðtÞ À s n ðv RMS Þp n ðtÞ ¼ C n ðv RMS ÞuðtÞ; where pðtÞ ¼ 2 The values of s n and C n must be updated at each time step by computing v RMS (L, t). The method for calculating the mean velocity at nonlinearity is detailed in the following section.

Computation of v RMS at each time step
To compute v RMS (L, t), the pressure at the termination L must be calculated. First, a linear problem is considered. The value of s n is therefore chosen for v RMS = 0 m/s. The modal components of the pressure are related to the total pressure at L through the following equation: and n 2 [0, L] is the distance along the longitudinal axis of the resonator. For reading comfort, the expression p(t) will be reserved to denote p(0, t).
The dimensionless acoustic velocity at the open end is then calculated from its definition which is related to the pressure field through the dimensionless Euler equation: The dimensionless acoustic velocity v(L, t) is thus obtained by numerical integration of the following equation: The RMS velocity is finally obtained by double time integration of equation (20), according to equation (6). Since the equations of the system characterizing the self-sustained oscillations are dimensionless, it is necessary to resize v RMS , because the modal coefficients C ðkÞ n and s ðkÞ n presented in equation (3) have been regressed from the dimensioned velocity. To do so, v ðadimÞ RMS is multiplied by p M /(q 0 c 0 ). In the following simulations, a value of p M = 8.5 kPa is set, according to experimental results from [12].
Regarding the calculation of v(L, t), it is not numerically guaranteed that the mean value of _ vðL; tÞ remains zero.
Consequently, the integral vðL; tÞ ¼ Z t 0 _ vðL; tÞdt may diverge. To avoid these divergence problems when integrating _ vðL; tÞ and v 2 (L, t), a short-memory term s is added within each integral. Thus, v(L, t) is now defined by the following convolution product: which is written, in the Laplace (L) domain, as: The same method is applied for the computation of v RMS :

Complete system of equations
The complete self-sustained system is governed by the equations defined in Sections 4. 1.1, 4.1.2, 4.1.3 and 4.1.4. This article shows the resolution of this problem using an ode solver. The Cauchy problem _ Y ¼ F ðY; tÞ writes, according to equation (28), as: xðtÞ þ x 2 r pðtÞ À cðtÞ À xðtÞ ½ where u(t) is computed by using equation (16), C k (v RMS ) are computed with equation (14), and s k (v RMS ) are computed with equation (12). The initial conditions Y 0 will be set to 0 in the following simulations.

Simulation including nonlinear losses at the open end
The model including nonlinear losses at the open end as represented by the system of equation (28) is now simulated by time integration, using the solver ode45 from Matlab. The input data of the problem are based on experimental results from [12].

Simulation parameters
The parameters related to the resonator, the reed and the reed channel are presented in Table 1. Six values of c d have been taken from [12], which are between c d = 0 and c d = 2.8. The value of f = 0.28 was calculated using data corresponding to a "loose embouchure" configuration. For the calculation of v RMS , the short-memory term was set to the period of the first impedance peak, i.e. s ¼ 2p=Iðs 1 Þ. The evolution of c is first linearly ascending, from c = 0 to c = 3 in 8 s, then linearly descending, from c = 3 to c = 0 in 8 s. The total simulation time is therefore 16 s. Simulations were performed for N = 4 modes. A convergence study showed that the dynamic behavior of the system did not change for a higher number of modes.

Results
The bifurcation diagram of the input pressure p (dim) = p Á p M over the blowing pressure P m = c Á p M is represented on Figure 4, with the parameters detailed in Section 4.2.1.
During a crescendo (Fig. 4a), the oscillation threshold is located near 5.7 kPa (c th % 0.67). This threshold remains the same, whether nonlinear losses are taken into account (c d 6 ¼ 0) or not (c d = 0). However, nonlinear losses have a significant influence on the extinction threshold c extup , i.e. the blowing pressure from which the reed stops oscillating and is pressed completely against the mouthpiece. This influence is shown in Table 2. When c d is increased (i.e. nonlinear losses increase), c extup diminishes. Similarly to the experiments conducted by Atig et al. [12], nonlinear losses have an important influence on the dynamic playing range of the musician. Figure 4 shows that when nonlinear losses are low, the range of stable oscillation amplitude that can be obtained is larger.
In the diminuendo phase, the inverse threshold is the same (c inv % 0.89) for each geometry at the open end of the pipe. Around this threshold, the amplitude of the input pressure decreases slightly as the losses increase. This behavior is also observed in the experimental curves of Ref. [12] of Figures 12a, 12c and 12e.
This comparaison with the experimental results of Atig et al. is limited to a qualitative study. Parameters such as the reed resonance angular frequency x r and reed damping q r have a strong influence on the dynamics of the system, as demonstrated by Silva et al. [28] for the oscillation threshold. The complete recalibration of the model on experimental results is beyond the scope of this article. Dalmont and Frappé (2007) An attempt to quantitatively validate the model of nonlinear losses is carried out. To do so, the values of the saturation threshold c sat measured by Dalmont and Frappé [13] are employed (see Fig. 5). This threshold corresponds to the blowing pressure c for which p is maximum during a crescendo. The authors measured this threshold for different reed openings, which are transcribed here in different values of f. The dimensioning pressure values p M used here correspond to the beating pressure estimated by the authors for a diminuendo. Furthermore, the measurements are performed on a cylindrical tube of dimensions L = 50 cm and R = 8 mm. The geometry at the opening has been estimated by the authors at c d = 2.8.

Comparison with experimental results from
The experimental data are compared to simulations using the same parameters L and R as in the experiment. Simulations are performed for two cases: the first one does not take into account nonlinear losses (c d = 0), the other one takes them into account (c d = 2.8). The simulations  are performed for ascending ramps of blowing pressure varying from c ¼ 0 to c ¼ 3 in 10 s. The other parameters are the same as those given in Table 1.
The model presented in this paper is also compared to the Raman model including nonlinear losses proposed by Dalmont and Frappé [13]. In this model, viscothermal losses are simplified by a frequency-independent coefficient a which replaces RðCÞ. The simplicity of this model allows the authors to obtain an analytical expression of the saturation threshold c sat as a function of the parameter c d .
In the present work, the value of a was first adjusted at RðCðjx 1 ÞÞ, where x 1 is the angular frequency of the impedance peak supporting the oscillation. However, this value produced excessively high estimates of c sat compared to experimental results, especially for high values of f. To better match the experimental data, a ¼ 2:7RðCðjx 1 ÞÞ % 0:147 was chosen. It should be noted that the analytical expression of c sat exhibits a linear dependence in 1=a, for low values of a. Thus, c sat is extremely sensitive to the value of a in Raman's model. Figure 5 illustrates the evolution of the saturation threshold as a function of f in the experimental case, through modal decomposition simulation, and through the analytical expression based on the Raman model. It appears first that the saturation threshold is overestimated when nonlinear losses are not taken into account ðþÞ. In comparison, the two models taking into account nonlinear losses produce results much closer to the experiment. This shows the importance of taking this phenomenon into account in wind instrument simulations.
Furthermore, for the analytical solution given by the Raman model ðÃÞ, the evolution of c sat follows a straight line whose slope depends on 1=a. Although the results are close to the experimental ones, the high sensitivity of c sat to a as well as the linear evolution of c sat reflect the excessive simplicity of the Raman model to describe the dynamics of a wind instrument. Finally, the model presented in this paper ðÂÞ has a good agreement with the experimental data, except for very low values of the mouthpiece parameter (f < 0:12). For these low values of f, the model overestimates c sat compared to the experimental results. This overestimation may be related to the transient behavior at extinction caused by the temporal evolution of the control parameter c. A characterization of the system by continuation could give results independent of the evolution rate of cðtÞ.
In conclusion, the correspondence between the results from the model presented in this article and the experimental data from [13] highlights its ability to describe the dynamic behavior of a simplified wind instrument at saturation.

Conclusion
This article develops the design of a sound synthesis model of a reed instrument by modal decomposition of the input impedance, taking into account viscothermal losses as well as nonlinear losses at the end of the resonator. The input impedance now depends on the RMS acoustic velocity at a geometric discontinuity (here, the open termination). Poles and residues resulting from the modal decomposition are fitted as a function of this velocity. In a physical model of wind instrument, the pressure-flow relation defined by the resonator is then completed by new equations which account for this dependence with the velocity at the end of the pipe.
To evaluate the ability of the model to reproduce a real phenomenon, comparisons with the experimental results of [12] and [13] have been made. In the first case, simulations show a qualitatively similar behavior regarding the evolution of the extinction threshold depending on the geometry at the open end. In the second case, the model gives a good correspondence with the experimental results, in particular compared to a model without nonlinear losses.
This formalism could be promising in the sound synthesis of wind instruments in real time, which employs the modal decomposition formalism in order to keep a minimum of memory. The inclusion of nonlinear losses in a complete clarinet model could more accurately translate dynamic phenomena essential to the experience of the musician, as suggested by Dalmont and Frappé [13]: "any realistic model of the clarinet should include nonlinear losses in the side holes".

Conflict of interest
The authors declare no conflict of interest.