# Parameter identification of a physical model of brass instruments by constrained continuation

Vincent Fréour, Louis Guillot, Hideyuki Masuda, Christophe Vergez, B.

Cochelin

► To cite this version:

Vincent Fréour, Louis Guillot, Hideyuki Masuda, Christophe Vergez, B. Cochelin. Parameter identification of a physical model of brass instruments by constrained continuation. Acta Acustica, 2022, 6, pp.9. 10.1051/aacus/2022004. hal-03837108

HAL Id: hal-03837108

https://hal.science/hal-03837108

Submitted on 16 Nov 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Parameter identification of a physical model of brass instruments by constrained continuation

Vincent Fréour¹, Louis Guillot², Hideyuki Masuda¹, Christophe Vergez², Bruno Cochelin²

¹YAMAHA Corporation, Research and Development Division,
10-1 Nakazawa-cho, Naka-ku, Hamamatsu, Shizuoka 430-8650, Japan.

²Aix Marseille Univ., CNRS, Centrale Marseille, LMA UMR7031, Marseille, France

## Summary

Numerical continuation using the Asymptotic Numerical Method (ANM), together with the Harmonic Balance Method (HBM), makes it possible to follow the periodic solutions of non-linear dynamical systems such as physical models of wind instruments. This has been recently applied to practical problems such as the categorization of musical instruments from the calculated bifurcation diagrams [15]. Nevertheless, one problem often encountered concerns the uncertainty on some parameters of the model (reed parameters in particular), the values of which are set almost arbitrarily because they are too difficult to measure experimentally. In this work we propose a novel approach where constraints, defined from experimental measurements, are added to the system. This operation allows uncertain parameters of the model to be relaxed and the continuation of the periodic solution with constraints to be performed. It is thus possible to quantify the variations of the relaxed parameters along the solution branch. The application of this technique to a physical model of a trumpet is presented in this paper, with constraints derived from experimental measurements on a trumpet player.

## 1 Introduction

Physical modeling is a valuable strategy in order to study numerically the response of a musical instruments [1]. In the case of brass instruments, various studies have focused on the modeling and simulation of the brass player and his/her instrument [2–10]. Some techniques such as numerical continuation are particularly relevant in order to compute periodic solution branches of wind instrument models [11]. Such a technique has been applied to woodwind instruments [14] and brass instruments [12, 13], and recently to the objective comparison of trumpets, through the extraction of descriptors calculated from the solution branches obtained by numerical continuation [15].

Nevertheless, the question of the values of the reed parameters (lip parameters for brass instruments) remains particularly sensitive. In brass instruments, the lip-reed system consists in a sophisticated biological system

1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

26
made of several layers of biological tissues, whose mechanical properties can be controlled by the player by activating various orofacial muscles. The question of the equivalent parameters of a biomechanical lip-reed model is made complicated by the difficulty in estimating these parameters experimentally, as well as by the uncertainty on the variations of these parameters in the course of a given musical task.

Some experimental studies on artificial player systems [6, 16], or human performers [2, 17–20], have shed some light on the dynamics of the lip-reed valve and proposed some values for the mechanical parameters of the lips, which are derived from frequency response measurements. Other strategies, based on numerical optimization techniques applied to the internal pressure, have been proposed in woodwind instruments in order to derive the mechanical parameters of a reed model from experimental measurements [21, 22]. In physical models of brass instruments, the mechanical parameters of the lip-reed model are often assumed constant for a given note. Velut et al. [23] proposed a review on the values used in different studies focused primarily on the trombone.

In this context, the question of matching the output of a physical model with experimental measurements becomes an important topic of interest. Furthermore, the strategy of numerical continuation applied to a physical model of brass instrument offers a global view of the response of the system (evolution of the state of the system with respect to one parameter of the model), likely to be used as a basis for fitting the model with experimental observations. The idea we propose in this study then comes from a dynamical system approach, where the modeled system can be augmented with additional equations that impose some constraints on the outputs of the model. In this article, these constraints are established from experimental observations on a musician. Adding these constraints then requires relaxation of some parameters of the model. The parameters which values and evolution are relatively uncertain along a periodic branch will be prioritized: the mechanical parameters of the lip model in the case of brass instruments.

This approach has two great benefits: 1- it enables identification of the evolution of the lip parameter values along a solution branch, based on some constraints established from experimental observations. This then provides information about the “gesture” required by the player’s model to fulfill the additional constraints. 2- it provides new information about the instrument, by allowing the variation of lip parameters required to achieve a given task (defined by the constraints) to be quantified. This is then suitable for instrument comparison, which is one of the motivations behind this work.

In this paper, application of constrained continuation to a physical model of the trumpet is presented, where the constraints are established from measurements on a trumpet player: frequency and/or amplitude of the sound produced during a crescendo. This paper is organized as follows: the physical model and basic continuation results are recalled in Section 2, measurements on a trumpet player are described in Section 3, the constrained continuation approach is then detailed in Section 4. Results and conclusions are presented in Sections 5 and 6.
# 2 Numerical continuation of physical model of brass instruments

## 2.1 Model and dimensionless equations

We consider a one-dimensional lip model, coupled to the resonator impedance described by a series of complex modes similar to what is proposed in [15]. The coupling between the mechanical oscillator and the acoustic resonator is achieved by a stationary Bernoulli flow equation, considering turbulent mixing in the mouthpiece with no pressure recovery. The mechanical and acoustic equations are given in system 1, where $y$ is the vertical lip position ($y_0$ is the lip position at rest), $\omega_l$, $Q_l$, $\mu_l$ and $b$ the lip mechanical parameters (resonance angular frequency, quality factor, mass per surface area and lip opening width respectively), $s_k$ and $C_k$ with $k \in [1, N]$ the modal parameters (poles and residues respectively) of the $N$ resonances of the acoustic impedance of the instrument, $Z_c$ the characteristic impedance, $u$ the volume flow, $p$ the downstream pressure at the input of the instrument (in the mouthpiece), and $p_0$ the upstream (mouth) static pressure.

$$ \begin{cases} \ddot{y}(t) + \frac{\omega_l}{Q_l} \dot{y}(t) + \omega_l^2 (y(t) - y_0) = \frac{1}{\mu_l} (p_0 - p(t)) \\ \dot{p}_k(t) = Z_c C_k u(t) + s_k p_k(t), \forall k \in [1, N], \end{cases} \quad (1) $$

with $p(t) = 2 \sum_{k=1}^{N} \Re(p_k(t))$ and $u = \sqrt{\frac{2|p_0-p|}{\rho}} b \cdot \text{sign}(p_0 - p) \cdot \theta(y)$, where $\theta(y) = \frac{|y|+y}{2}$, $b$ is the lip width and $\rho$ is the air density.

The case of a negative opening of the lips is managed by introducing the function $\theta(y)$ which enforces $u = 0$ if $y < 0$. The modal parameters of the $N$ modes of the impedance are extracted from the measured input impedance, corrected to 27°C [24], using the high resolution method ESPRIT [25]. We remind that nonlinear propagation phenomena that may originate from large pressure levels (at high mouth pressure values) are not taken into account in this model.

Figure 2 represents a reconstruction of the input impedance of a Bb trumpet in open fingering (no valve pressed) from the superposition of 11 complex modes, against the measured input impedance. Despite some discrepancies at some anti-resonances, this representation shows overall a very satisfactory match in both magnitude and phase between the two curves. The corresponding poles and residues values are given in Table 1.

The choice is made to work with the Asymptotic Numerical Method (ANM) implemented in the software MANLAB [26]. Recently, this method has been associated with the Harmonic Balance Method (HBM) for the search of periodic solutions of oscillating systems [28]. One requirement of MANLAB relies on the recast of nonlinearities of the model into, at most, quadratic nonlinearities. This system of equations can be made dimensionless by introducing the following variables:
<figure><img src="image_2.png" /><figcaption>Figure 1: (Color online) Measured input impedance (dashed black) and reconstructed input impedance from 11 complex modes (red) for a B♭ trumpet with open fingering.</figcaption></figure>

<table><caption>Table 1: Values of poles and residues extracted from the input impedance of a B♭ trumpet with open fingering using the ESPRIT method.</caption><thead><tr><th>peak index</th><th>s<sub>k</sub></th><th>C<sub>k</sub></th></tr></thead><tbody><tr><td>1</td><td>-1.3979.10<sup>1</sup> + 5.2247.10<sup>2</sup>i</td><td>1.3601.10<sup>9</sup> + 2.5960.10<sup>7</sup>i</td></tr><tr><td>2</td><td>-2.2419.10<sup>1</sup> + 1.4621.10<sup>3</sup>i</td><td>1.7437.10<sup>9</sup> + 3.6566.10<sup>7</sup>i</td></tr><tr><td>3</td><td>-2.8642.10<sup>1</sup> + 2.1870.10<sup>3</sup>i</td><td>2.4381.10<sup>9</sup> + 3.4217.10<sup>7</sup>i</td></tr><tr><td>4</td><td>-3.7641.10<sup>1</sup> + 2.9066.10<sup>3</sup>i</td><td>4.7175.10<sup>9</sup> - 6.5252.10<sup>6</sup>i</td></tr><tr><td>5</td><td>-4.5031.10<sup>1</sup> + 3.6577.10<sup>3</sup>i</td><td>5.7356.10<sup>9</sup> + 8.7953.10<sup>5</sup>i</td></tr><tr><td>6</td><td>-4.9819.10<sup>1</sup> + 4.3390.10<sup>3</sup>i</td><td>7.6553.10<sup>9</sup> - 1.1749.10<sup>7</sup>i</td></tr><tr><td>7</td><td>-5.8425.10<sup>1</sup> + 5.0290.10<sup>3</sup>i</td><td>7.3308.10<sup>9</sup> + 1.3414.10<sup>8</sup>i</td></tr><tr><td>8</td><td>-6.6774.10<sup>1</sup> + 5.7054.10<sup>3</sup>i</td><td>4.7527.10<sup>9</sup> + 2.2058.10<sup>8</sup>i</td></tr><tr><td>9</td><td>-7.2239.10<sup>1</sup> + 6.4595.10<sup>3</sup>i</td><td>2.3351.10<sup>9</sup> + 2.8082.10<sup>8</sup>i</td></tr><tr><td>10</td><td>-9.4396.10<sup>1</sup> + 7.2109.10<sup>3</sup>i</td><td>1.6618.10<sup>9</sup> + 2.8424.10<sup>8</sup>i</td></tr><tr><td>11</td><td>-1.2862.10<sup>2</sup> + 7.9310.10<sup>3</sup>i</td><td>1.133810<sup>9</sup> + 4.9805.10<sup>8</sup>i</td></tr></tbody></table>
$$
\left\{
\begin{aligned}
x &= \frac{y}{y_0} \\
P_M &= \mu_l \omega_l^2 y_0 \\
\gamma &= \frac{p_0}{p_M} \\
\tilde{p} &= \frac{p}{p_M} \\
\tilde{R}_k &= \frac{R_k}{p_M} \\
\tilde{I}_k &= \frac{I_k}{p_M} \\
\tilde{u} &= u \frac{Z_c}{p_M} \\
\tilde{v} &= \frac{v}{\sqrt{2P_M/\rho}} \\
\zeta &= Z_c by_0 \sqrt{\frac{2}{\rho P_M}} \\
\tilde{t} &= t \Im(s_1) \\
\tilde{\omega}_l &= \frac{\omega_l}{\Im(s_1)} \\
\tilde{C}_k &= \frac{C_k}{\Im(s_1)} \\
\tilde{s}_k &= \frac{s_k}{\Im(s_1)}
\end{aligned}
\right.
\qquad (2)
$$

with $R_k$ and $I_k$ the real and imaginary parts of the pressure components $p_k$ with $k \in [1, N]$. The complete quadratic dimensionless model can then be written as follows.

88

89

$$
\left\{
\begin{aligned}
\tilde{R}_k &= \Re(\tilde{C}_k)\tilde{u} + \Re(\tilde{s}_k)\tilde{R}_k - \Im(\tilde{s}_k)\tilde{I}_k, \quad \forall k \in [1, N] \\
\tilde{I}_k &= \Im(\tilde{C}_k)\tilde{u} + \Im(\tilde{s}_k)\tilde{R}_k + \Re(\tilde{s}_k)\tilde{I}_k, \quad \forall k \in [1, N] \\
\tilde{x} &= \tilde{\omega}_l \tilde{z} \\
\tilde{z} &= \tilde{\omega}_l \left(1 - x - \frac{1}{Q_l} z + \gamma - \tilde{p}\right),
\end{aligned}
\right.
\qquad (3)
$$

with the auxiliary equations:

90

$$
\left\{
\begin{aligned}
0 &= 2 \sum_{k=1}^{N} \tilde{R}_k - \tilde{p} \\
0 &= x^2 + \epsilon_x - s^2 \\
0 &= \gamma - \tilde{p} - \tilde{v}w \\
0 &= \tilde{v}^2 + \epsilon_v - w^2 \\
0 &= \zeta \frac{(s+x)}{2} \tilde{v} - \tilde{u},
\end{aligned}
\right.
\qquad (4)
$$

with $0 < \epsilon_v \ll 1$ and $0 < \epsilon_x \ll 1$, the regularization constants such as $\epsilon_v = \epsilon_x = 10^{-3}$. The reader is invited to consult [15] for more details.

91

92
## 2.2 Continuation by ANM

The dynamical system described by Eqs. 3 and 4 can be analyzed by numerical continuation, especially using the Asymptotic Numerical Method (ANM) [27] implemented in the software MANLAB [26]. This method is based on the expansion of the solutions under the form of truncated Taylor series, providing analytical approximate formulations of the branch of solution. Associated with the Harmonic Balance Method (HBM), the ANM allows for search of periodic solutions of the dynamical system, the unknowns being the Fourier coefficients of each variable [29, 30] and the oscillation frequency. For more details, the reader is invited to refer to the specific literature on the subject [31, 32].

Figure 2 shows the results of continuation the system described by Eqs. 3 and 4: it represents the bifurcation diagram (peak to peak amplitude of *p* with respect to *p*₀) of the branch of periodic solution corresponding to a Bb4. This result is obtained with the lip parameters given in Table 2.

<table><caption>Table 2: Lip parameters used for the numerical continuation of the model.</caption><tr><td>$Q_l$</td><td>3</td></tr><tr><td>$\mu_l$</td><td>2 kg.m<sup>-2</sup></td></tr><tr><td>$y_0$</td><td>0.1 mm</td></tr><tr><td>$b$</td><td>8 mm</td></tr></table>

The lip natural frequency is set by Linear Stability Analysis (LSA) [15, 23] to $f_l = \omega_l/2\pi = 382.18$ Hz, so that the playing frequency is closed to a Bb4 ($f_0 \approx 470$ Hz).

The bifurcation diagram obtained in Fig. 2 sheds light on some behavior that brass players are familiar with, and that have been described in previous studies [12, 15]. These features include: 1- an inverse bifurcation and Hopf point around $p_0 = 2.2$ kPa; 2- an unstable section oriented towards the left until it reaches a fold; 3- at the fold, a change of stability towards a stable branch where the amplitude of the internal pressure *p* increases with the mouth pressure $p_0$. Note that for the lips parameters chosen (Table 2), the order of magnitude of $p_0$ is coherent with experimental observations in the trumpets [33].

## 3 Measurements on trumpet player

Whereas the numerical results obtained by continuation seem coherent with the behavior of a brass instrument, it is then relevant to compare more accurately these numerical outputs with some experimental data obtained on brass players. For this purpose, experimental bifurcation diagrams were collected by measurements conducted on a trumpet player using the trumpet whose impedance was measured in order to implement the model described in the previous section. For these measurements, the tuning slide was set to its default setting, as in the impedance measurements, in order to use exactly the same resonator in both conditions.
<figure><img src="image_3.png" /><figcaption>Figure 2: (Color online) Bifurcation diagram (peak-to-peak amplitude of *p* with respect to *p*₀) of the periodic branch of solution for a Bb4 (470Hz). The dotted line indicates the unstable part of the solution branch, while the solid line indicates the stable part of the branch. ○: Hopf bifurcation (*p*<sub>H</sub> point), ∇: fold, △: |*p*| at the Hopf bifurcation, □: |*p*| when *p*₀ = 5 kPa.</figcaption></figure>

**3.1 Experimental setup and protocol**

A National Instrument acquisition card and two ENDEVCO pressure sensors (model 8510B) are used to measure the upstream (mouth) pressure $p_{0e}$ and downstream (mouthpiece) pressure $p_e$ respectively (the sampling rate is set to $f_s = 51,2$ kHz). The mouthpiece sensor is set up so that the extremity is mounted flush in the mouthpiece cup. The mouth sensor is connected to a small tube that the musician is asked to insert in the mouth. An illustration of the experimental setup is shown in Fig. 3.

In order to collect experimental data comparable as much as possible with a bifurcation diagram obtained numerically, the player is instructed to produce a slow crescendo-decrescendo from *p*pp*pp* to *ff* (see Fig. 4). The objective is to be as close as possible to the hypothesis of a quasi-static variation of the mouth pressure. The player is then instructed to avoid using the tongue at the note onset (no tonguing at the attack) and to make the crescendo-decrescendo as regular and as slow as possible (recommended tempo: 70 bpm). No particular instructions were given to the players with regards to the stability of the fundamental frequency.

The collected data are treated as follows: a moving window is applied to the data (window size 20 ms with no overlap). The peak-to-peak amplitude of the mouthpiece pressure $|p_e|$ is extracted for each window, as well as the mean of the low-pass filtered mouth pressure $\bar{p}_{0e}$. The duration of the crescendo and decrescendo phase are then extracted from the peak-to-peak envelope (Fig. 5). In order to impose a relatively low variability between
<figure><img src="image_4.png" /><figcaption>Figure 3: (Color online) Trumpet player during experiments: two miniature pressure sensors are used to measure the upstream (mouth) and downstream (mouthpiece) pressures during performance.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 4: Musical task asked to the player. Recommended tempo: 70 bpm. Note that the musical staff shows written pitch C4, a tone higher than sounding pitch B♭4.</figcaption></figure>

the tasks used for analysis, only the maneuvers where the duration of the crescendo and decrescendo are within 3.5 and 4.5 seconds are considered. Out of eight repetitions of the task, three maneuvers are then discarded. The time-evolution of $|p_e|$ and $\bar{p}_{0e}$ for the five selected repetitions of the crescendo-decrescendo task is shown in Fig. 6

The experimental bifurcation diagrams ($|p_e|$ with respect to $\bar{p}_{0e}$, and $f_0$ with respect to $\bar{p}_{0e}$) associated with the data represented in Fig. 6 are then represented in Fig. 7.

The results of Fig. 7 show some common features with the numerical diagram of Fig. 2: 1- a minimum blowing pressure in the same pressure range as calculated by continuation, 2- an inverse bifurcation that appears under the form of an hysteresis in the experimental data, 3- a relatively linear increase of $|p_e|$ with $\bar{p}_{0e}$ as observed in numerical results, 4- relatively small variations of $f_0$ along the maneuver. These results then support our hypothesis regarding the ability of the 1-mode outward-striking lip model to generate behaviors that seem relatively coherent with human behaviors.

Figure 8 shows the five bifurcation diagrams of Fig. 7 overlapped on the same plot. This representation highlights the relatively high repeatability of the player across the five iterations. In the case of $|p_e|$, the average of the linear regression coefficients of each bifurcation diagram is calculated for $\bar{p}_{0e} > 1.5$ kPa (plotted as a red line in Fig. 8, top). In the case of $f_0$, the average of $f_0$ values for $\bar{p}_{0e} > 1.5$ kPa is calculated (plotted as a red line in Fig. 8, bottom).
<figure><img src="image_6.png" /><figcaption>Figure 5: (Color online) Crescendo-decrescendo maneuvers performed by the trumpet player. Mouthpiece pressure (gray), mouth pressure (blue) and peak-to-peak mouthpiece pressure envelope (red). The crosses indicate the onset, maximum amplitude and sound extinction locations. The green area corresponds to the crescendo phase, the magenta area to the decrescendo phase.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 6: (Color online) Five crescendo-decrescendo maneuvers performed by the trumpet player. Top: mouth pressure measured (black) and low-pass filtered (red). Bottom: mouthpiece pressure (black) and peak-to-peak envelope (red).</figcaption></figure>

## 3.2 Comparison with numerical solutions

Although some important similarities can be found between the bifurcation diagrams obtained numerically (Fig. 2) and experimentally (Fig. 8), it seems very difficult, with constant lip parameters, to generate a numerical bifurcation diagram with a stable branch of similar slope as the measured diagrams. Figure 9 shows results of 100 bifurcation diagrams calculated from random variations of the lip parameters. More precisely, the parameters

152

153

154

155

156
<figure><img src="image_8.png" /><figcaption>Figure 7: Experimental bifurcation diagrams obtained from one player playing a Bb4 with slow crescendo-decrescendo. Top: $|p_e|$ with respect to $\bar{p}_{0e}$. Bottom: $f_0$ with respect to $\bar{p}_{0e}$.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Figure 8: (Color online) Experimental bifurcation diagrams obtained from one player playing a Bb4 with slow crescendo-decrescendo (same data as in Fig. 7). Top: $|p_e|$ with respect to $\bar{p}_{0e}$. Bottom: $f_0$ with respect to $\bar{p}_{0e}$. The top red line is obtained by averaging the linear regression coefficients of each diagram for $\bar{p}_{0e} > 1.5$ kPa. The bottom red line is obtained by averaging $f_0$ values for $\bar{p}_{0e} > 1.5$ kPa.</figcaption></figure>

are varied by ±50% around the values given in Table 2 for $Q_l$ and $b$, and by ±90% for $\mu_l$ and $y_0$. For each calculation, the lip natural frequency $f_l$ is obtained from LSA. In this figure, the amplitude of $p$ is estimated by the $L^2$ norm, which choices is motivated by the mathematical form of this norm, as it will be explained in the next section.
Linear regressions applied to the stable parts of the calculated diagram (||p||L2 w.r.t. p0) allows extraction of a slope value for each diagram. On the top right plot of Fig. 9, these slope values and the mouth pressure at the Hopf point, pH, are represented for each diagram, along with the slope and intercept with the horizontal axis of the red line obtained from experimental data and represented on the top left plot of Fig. 9. This intercept with the horizontal axis corresponds indeed to the Hopf point of the target solution. Despite large variations of the lip parameters across calculations, the slopes of the calculated stable branches are lower than the slope of the experimental diagrams. Furthermore, on the bottom right plot, f0 values, as well as the variations of f0 along the solution branch (Δf0 = max(f0) - min(f0)), are represented for all diagrams, along with the values obtained from the experimental data. It can be noticed from fig. 8 that Δf0 is associated with an overall decrease of f0 with p0. The calculated fundamental frequencies are significantly above the values observed experimentally, which is usually expected with an outward striking lip model [34].

These observations suggest that it might be relatively difficult to find a combination of lip parameters allowing experimental diagrams to be replicated numerically. In the following, an approach is proposed to study if a numerical bifurcation diagram closer to the experimental one can be reached if some parameter values are allowed to vary with respect to the bifurcation parameter p0.

More precisely, by enabling some lip parameters to vary along the solution branch, and by imposing some constraints to the solution, the periodic solutions of this new extended system are derived, as well as the evolution of the relaxed parameters along the solution branch.

# 4 Constrained continuation

In order to generate numerical diagrams with similar features as the one obtained with the trumpet player, some constraints should then be defined from the bifurcation diagrams obtained experimentally.

## 4.1 Constraints

Two constraints are defined (red lines of Fig. 9), which represent the result of the musician's action. The first constraint concerns the relationship between the L2 norm of $\tilde{p}$ (||$\tilde{p}$||L2) and the dimensionless mouth pressure $\gamma$, along the stable part of the branch of the bifurcation diagram:

$$ ||\tilde{p}||_{L^2} = S\gamma + I, \qquad (5) $$

where $\gamma = p_0/P_M$ is the dimensionless mouth pressure with $P_M = \mu_l\omega_l^2y_0$, S and I are constant values, and $||\tilde{p}||_{L^2} = 2||\sum_{k=1}^N \Re(\tilde{p}_k)||_{L^2} = 2\sqrt{\sum_{k=1}^N \Re(\tilde{p}_k)^2}$. The value of $\omega_l$ is calculated by LSA as in Section 2 such as $f_l = \omega_l/2\pi = 382.18$ Hz. The choice of the L2 norm is motivated by the quadratic form of Eq. 5 w.r.t. the unknowns when elevated at the power of 2, quadratic nonlinear equations being a prerequisite of the ANM.
<figure><img src="image_10.png" /><figcaption>Figure 9: (Color online) Left: 100 bifurcation diagrams obtained from random variations of the lip parameters around the values given in table 2, the plain red lines are the targets obtained from human performance. The dashed line on the bottom plot corresponds to the modified target. Top right: mouth pressure at Hopf bifurcation and slope calculated from 100 linear regressions performed on the stable sections of the 100 bifurcation diagrams, the red diamond corresponds to the red line in the left plot. Bottom right: $\Delta f_0$ and $f_0$ obtained from the 100 bifurcation diagrams, the plain red diamond correspond to the red line in the left plot. The empty red diamond corresponds to the modified target.</figcaption></figure>

The second constraint simply writes as follows:

$$ \tilde{\omega} = F, \qquad (6) $$

where $F$ is a constant value.

## 4.2 Relaxed parameters and new system to solve

Equations 5 and 6 are added to the original system (systems 3 and 4). Adding two equations to the system has to be balanced by the introduction of two new unknowns. $Q_L$ and $\zeta = Z_c by_0 \sqrt{\frac{2}{\rho P_M}}$ are chosen as the relaxed parameters. The parameter $\zeta$ can be interpreted as an “embouchure” parameter as it depends only on lip and mouthpiece parameters, and in particular on the lip width $b$ that appears exclusively in the expression of $\zeta$.

This requires the system of equations to be recast in order to preserve the quadratic property of the model.
The main equations are unchanged compared to Eq. 3:

197

$$
\left\{
\begin{aligned}
\tilde{R}_k &= \Re(\tilde{C}_k)\tilde{u} + \Re(\tilde{s}_k)\tilde{R}_k - \Im(\tilde{s}_k)\tilde{I}_k, \forall k \in [1, N] \\
\tilde{I}_k &= \Im(\tilde{C}_k)\tilde{u} + \Im(\tilde{s}_k)\tilde{R}_k + \Re(\tilde{s}_k)\tilde{I}_k, \forall k \in [1, N] \\
\tilde{x} &= \tilde{\omega}_l \tilde{z} \\
\tilde{z} &= \tilde{\omega}_l (1 - x - q_l z + \gamma - \tilde{p}),
\end{aligned}
\right.
\qquad (7)
$$

with $q_L = 1/Q_L$. Two equations are introduced to account for the constraints:

198

$$
\left\{
\begin{aligned}
0 &= (S\gamma + I)^2 - \sum_{k=1}^{N} \tilde{R}_k^2 \\
0 &= 2\pi F - \tilde{\omega}.
\end{aligned}
\right.
\qquad (8)
$$

One additional auxiliary variable $m$, as well as the corresponding equation are introduced in order to preserve the quadratic property of the model. Indeed, since $\zeta$ is now an unknown, the last equation of system 4 is no longer quadratic. The system of auxiliary equations then becomes:

199

200

201

$$
\left\{
\begin{aligned}
0 &= 2 \sum_{k=1}^{N} \tilde{R}_k - \tilde{p} \\
0 &= x^2 + \epsilon_x - s^2 \\
0 &= \gamma - \tilde{p} - \tilde{v}w \\
0 &= \tilde{v}^2 + \epsilon_v - w^2 \\
0 &= \frac{(s+x)}{2} \tilde{v} - m \\
0 &= \zeta m - \tilde{u}.
\end{aligned}
\right.
\qquad (9)
$$

The system formed by Eq. 7, 8 and 9 is dimensionless and quadratic. It is then compatible with the application of the Asymptotic Numerical Method (ANM) using MANLAB.

202

203

The following vector of unknowns is then considered:

204

$$
X = [R_1, I_1, R_2, I_2, \dots, R_N, I_N, x, z, \tilde{p}, s, \tilde{v}, w, m, \tilde{u}, \zeta, q_L, \tilde{\omega}]. \qquad (10)
$$

# 5 Results of constrained continuation

205

## 5.1 Values of the constraints

206

The two constraints are defined based on the experimental observations of Section 3:

207

1. The values of $S$ and $I$ are calculated from the average of the linear regressions applied to the musician's data: $S = 1.5461$ and $I = -0.9724$. In Fig. 9, the red line crosses the beam of blue diagrams for $\|p\|_{L^2}$

208

209
(top left). This suggests that it should be possible to find a solution with relaxed parameters that will
follow the red constraint.

2. In Fig. 9, although the red line crosses the beam' of blue diagrams for ||p||_L² (top left), it is not the case for f₀ (bottom left), where the plain red line that comes from experimental results does not share any point with the numerical results. This mismatch is due to the outward striking lip model which, by construction, oscillates at frequencies slightly above human lips. This observation implies that no solution with relaxed parameter should be able to reach measured f₀ values. The value of F is then set to 480 Hz, slightly above the averaged value measured on the player, as shown in Fig. 9. This frequency remains in the vicinity of B₅: 42.8 cents above the theoretical B₅ at 468.28 Hz (with 442 Hz A₄ tuning reference).

## 5.2 Initial conditions

<figure><img src="image_11.png" /><figcaption>Figure 10: (Color online) Initial bifurcation diagram (blue) and target diagram (dashed red).</figcaption></figure>

In order to initialize the calculation with constraints, an initial point should be defined, preferably not too far
from the target solution branch. To this end, among the 100 bifurcation diagrams generated in Section 3, one
diagram, and the associated lip parameters, is selected, so that it crosses the target diagram within the p₀ range
of interest (between 1 kPa and 4 kPa). This diagram is represented in Fig. 10 and the associated lip parameters
given in Table 3. A point near the intersection with the target (around p₀ = 2.5 kPa) is extracted and is used
as the initial point for the constrained continuation. Priority is given to the intersection with the ||p||_L² target
since, after several trials, it appears to be more critical than f₀ for the initialization of the calculation and for
<table><caption>Table 3: Lip parameters used to initialize the constrained continuation. The lip natural frequency is then set by LSA to $f_l = 414$ Hz.</caption><tr><td>$Q_l$</td><td>5.25</td></tr><tr><td>$\mu_l$</td><td>2.46 kg.m<sup>-2</sup></td></tr><tr><td>$y_0$</td><td>0.078 mm</td></tr><tr><td>$b$</td><td>12.78 mm</td></tr></table>

the convergence of towards the constrained solution.

## 5.3 Results

Figure 11 shows the results of continuation with the constraints defined above. It can be seen that the constrained continuation can be performed successfully over a large range of $p_0$ (up to about 5 kPa). The variations of $\zeta$ can be observed along the solution branch. Within the $p_0$ range covered by the human player (up to about 3 kPa), these variations are such that $0.04 < \zeta < 0.1$ and $2.9 < Q_l < 5$. For $Q_l$, this interval falls within expected values for a brass player lip Q factors [23]. For $\zeta$, this is equivalent to a variation by about 60% of the maximum value. Among the lip parameters, the lip width $b$ is the only parameter that appears only in the expression of $\zeta$, and not in the definition of $P_M$ (System 2). Furthermore, $b$ is proportional to $\zeta$, assuming other parameters constant. The variations of $\zeta$ can then be interpreted as variations of $b$ by 60% of its maximum value, which we believe to be quite realistic: it is physically acceptable to consider that a player may vary the lip width by such amount within a crescendo-decrescendo maneuver, although this should be eventually confirmed by experimental measurements if possible.

Figures 12 (left plot) shows the waveform of $p$, obtained from experimental measurements with the trumpet player, for three values of $p_0$: 1 kPa, 2 kPa and 3 kPa. In addition, the waveform of $p$ obtained from the bifurcation diagram of Fig. 11 is represented for the same $p_0$ values in Fig. 12 (right plot). Despite some differences in peak-to-peak amplitudes ($\|p\|_{L^2}$ is not the peak-to-peak amplitude of $p$), both figures show a clear, and relatively similar, evolution of the waveform amplitude and shape with increase in $p_0$. This observation confirms that relevance of the model and of the constrained continuation in reproducing the behavior monitored on a human player.

To illustrate this statement, two-second duration sounds corresponding to the waveforms of Fig. 12 are generated (after normalization and application of a linear envelope at attack and release), from the experimental measurements for $p_0 = 1$ kPa, $p_0 = 2$ kPa, and $p_0 = 3$ kPa, as well as from the bifurcation diagram for $p_0 = 1$ kPa, $p_0 = 2$ kPa, and $p_0 = 3$ kPa. Despite some differences in intonation between the model and the experiment (as expected with the outward striking lip model), these sound excerpts confirm a similar evolution of the timbre with $p_0$ between the model and the musician's recording.

In order to further investigate the nature of the calculated solution, the constrained bifurcation diagram, along with a number of bifurcation diagrams where the lip parameters are set constant and taken at various points of
<figure><img src="image_12.png" /><figcaption>Figure 11: (Color online) Bifurcation diagram with constraints and evolution of the relaxed parameters $\zeta$ and $Q_t$ with respect to $p_0$. The target diagram is represented by the red dashed line.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Figure 12: Waveform of $p$ obtained from experimental measurements on a trumpet player (left) and from constrained continuation (right), at $p_0 = 1$ kPa (solid line), $p_0 = 2$ kPa (dashed-dotted line) and $p_0 = 3$ kPa (dashed line).</figcaption></figure>

the constraint solution, are represented in Fig. 13. The intersection points (circle markers) highlight the path followed by the constraint solution across constant-parameters solution branches. Note that to be considered
as an intersection point with a constant-parameter diagram, both $f_0$ and $\|p\|_{L^2}$ diagrams should cross a same constant-parameter diagram. In this figure, it appears that all the calculated diagrams are characterized by inverse bifurcations, and that below a certain $p_0$ value (around 1400 Pa), the constrained solution crosses unstable branches only. In other words, the constrained solution is formed by an ensemble of points that belong to bifurcation diagrams which emerge from inverse Hopf bifurcations. Some of these points are then stable and other are unstable, which makes the constrained solution stable above $p_0 = 1400$ Pa, and unstable below this value. Note also that the constrained solution thus does not emerge from any Hopf bifurcation, especially as it is produced through the joint modification of three parameters ($\zeta, Q_{lip}$ and $p_0$), where a Hopf bifurcation is defined trough the variation of a single parameter ($p_0$ in our case). That being said, the question whether this stability result depends on the initial conditions (lip parameters at the initial calculation point) remains to be clarified: it is possible that another set of lip parameter values at the initial calculation point will produce a different stability result.

<figure><img src="image_14.png" /><figcaption>Figure 13: (Color online) Bifurcation diagram with constraints (blue), overlapped with constant-parameter bifurcation diagrams where $\zeta$ and $Q_l$ are set to values obtained from constrained continuation. Circle markers indicate intersection points between the diagram with constraints and the constant-parameters diagrams.</figcaption></figure>
## 6 Conclusions

In this article, a method of continuation where constraints are included in the system of equations, and parameters of the model are relaxed, is proposed in order to force the model to follow a solution branch close to a human performance. Contrary to an optimization strategy, this method provides the solution of the constrained problem without any need of defining a cost function nor exploring the parameter space. Within some limitations regarding the fundamental frequency, that we may assume to be linked to the properties of the outward striking lip model, this method was successfully applied to constraints established from measurements on a trumpet player. The obtained trajectories of the relaxed parameter reveal significant but realistic variations along the solution branch, which confirms the great potential of the one-degree-of-freedom lip model in producing results close to human performances, as well as the importance to make lip model parameters vary to produce these outputs. Nevertheless, the obtained results also show some limitations of the present model in producing oscillations in the same range as human players, particularly a stable branch down to $||p||_{L^2} = 0$. Modifying the lip model parameters that define the initial conditions, may result in different behaviors (e.g. direct bifurcations), then allowing to compute stable constrained solutions down to lower pressure values. Moreover, applying the constraint from a given point of the initial bifurcation diagram (from the fold point for instance) may contribute to generate a constrained solution closer to experimental bifurcation diagrams. Closer investigations at the influence of the initial conditions on the constrained continuation results, as well as at the possibility to apply the constraints from specific landmarks of an initial solution, should be the objects of future investigations.

In the future, we also plan to apply this method to more sophisticated lip models, although uncertainties on some parameter values of the model (second mode, contact parameters) should be addressed carefully. More advanced constraints should also be defined, in order to attempt a better match with experimental observations.

Regarding potential applications, this method opens some interesting perspectives for the comparison of musical instruments, by providing new indicators related to the control of sound production. For instance, the variations of the lip parameters along the solution branch may provide some useful indications on corrections needed at the level of the embouchure in order to produce a constant-pitch crescendo. It may then bring some new basis for objective comparisons of brass instruments using physical modeling and numerical continuation.

## 7 Data availability statement

The sound files associated with this article are available in https://medihal.archives-ouvertes.fr/hal-03545981

## 8 Acknowledgments

The authors would like to thank the Yamaha employees who participated in the trumpet playing measurements.
References

299

[1] M. Campbell and J. Gilbert and A Myers. *The Science of Brass Instruments*. Springer International Publishing, 2021.

300

301

[2] S. J. Elliot and J. M. Bowsher. Regeneration in brass instruments. *Journal of Sound and Vibration*, 83(2):181–217, 1982.

302

303

[3] N. H. Fletcher. Autonomous vibration of simple pressure-controlled valve in gas flows. *J. Acoust. Soc. Am.*, 93(4):2172–2180, 1993.

304

305

[4] S. Adachi and M. A. Sato. Trumpet sound simulation using a two-dimensional lip vibration model. *J. Acoust. Soc. Am.*, 99(2):1200–1209, 1996.

306

307

[5] C. Vergez and X. Rodet. Comparison of real trumpet playing, latex model of lips and computer model. In *Proc. ICMC97*, 1997.

308

309

[6] J. S. Cullen, J. Gilbert, and D. M. Campbell. Brass instruments: Linear stability analysis and experiments with an artificial mouth. *Acta Acustica*, 86:704–724, 2000.

310

311

[7] T. Kaburagi, N. Yamada, T. Fukui, and E. Minamiya. A methodological and preliminary study on the acoustic effect of a trumpet player's vocal tract. *J. Acoust. Soc. Am.*, 130(1):536–545, 2011.

312

313

[8] V. Fréour, N. Lopes, T. Hélie, R. Caussé, and G.P. Scavone. In-vitro and numerical investigations of the influence of a vocal-tract resonance on lip auto-oscillations in trombone performance. *Acta Acustica united with Acustica*, (101):256–269, 2015.

314

315

316

[9] N. Lopes and T. Hélie. Energy balanced model of a jet interacting with a brass player's lip. *Acta Acustica united with Acustica*, 102:141–154, 2016.

317

318

[10] H. Benjamin, B. Lombard, C. Vergez, and E. Cottanceau. Time-domain numerical modeling of brass instruments including nonlinear wave propagation, viscothermal losses, and lips vibration. *Acta Acustica united with Acustica*, 103:117–131, 2017.

319

320

321

[11] S. Karkar. *Méthodes numériques pour les systèmes dynamiques non linéaires - Application aux instruments de musique auto-oscillants*. PhD thesis, Université de Provence Aix-Marseille I, Marseille, France, 2012.

322

323

[12] L. Velut. *Contrôle par le musicien des régimes d'oscillation des instruments de la famille des cuivres*. PhD thesis, Université de Provence Aix-Marseille I, Marseille, France, 2016.

324

325

[13] J. Gilbert, S. Maugeais, and C. Vergez. Minimal blowing pressure allowing periodic oscillations in a simplified reed musical instrument model: Bouasse-Benade prescription assessed through numerical continuation. *Acta Acustica*, 4(6):27, 2020.

326

327

328
[14] T. Colinot, C. Vergez, P. Guillemain, and J.B. Doc. Multistability of saxophone oscillation regimes and its influence on sound production. *Acta Acustica*, 5:33, 2021.

[15] V. Fréour, L. Guillot, H. Masuda, S. Usa, E. Tominaga, Y. Tohgi, C. Vergez, and B. Cochelin. Numerical continuation of a physical model of brass instruments: Application to trumpet comparisons. *J. Acoust. Soc. Am.*, 148(2):748–758, 2020.

[16] J. Gilbert, S. Ponthus, and J. F. Petiot. Artificial buzzing lips and brass instruments: Experimental results. *J. Acoust. Soc. Am.*, 104(3):1627–1632, 1998.

[17] M.J. Newton, M. Campbell, and J. Gilbert. Mechanical response measurements of real and artificial brass player lips. *J. Acoust. Soc. Am.*, 123(1):EL 14–EL20, 2008.

[18] V. Fréour, and G.P. Scavone. Acoustical interaction between vibrating lips, downstream air column, and upstream airways in trombone performance. *J. Acoust. Soc. Am.*, 134(5):3887–3898, 2013.

[19] T. Hézard, V. Fréour, R. Caussé, T. Hélie, and G.P. Scavone. Synchronous Multimodal Measurements on Lips and Glottis: Comparison Between Two Human-Valve Oscillating Systems. *Acta Acustica united with Acustica*, 100:1172–1185, 2014.

[20] H. Boutin, N. Fletcher, J. Smith, and J. Wolfe. Relationships between pressure, flow, lip motion, and upstream and downstream impedances for the trombone. *J. Acoust. Soc. Am.*, 137(3):1195–1209, 2015.

[21] V. Chatziioannou, and M. van Walstijn. Estimation of Clarinet Reed Parameters by Inverse Modelling. *Acta Acustica united with Acustica*, 98:629–639, 2012.

[22] A. Munoz Arancon, B. Gazengel, J-P. Dalmont, and Conan E. Estimation of saxophone reed parameters during playing. *J. Acoust. Soc. Am.*, 139(5):2754–2765, 2016.

[23] L. Velut, C. Vergez, J. Gilbert, and M. Djahanbani. How well can linear stability analysis predict the behaviour of an outward-striking valve brass instrument model? *Acta Acustica united with Acustica*, 103:132–148, 2017.

[24] J. Gilbert, L.M. Leboso Ruiz, and S. Gougeon. Influence de la température sur la justesse d'un instrument à vent. In *Proc. 8ème Congrès Français d'Acoustique*, pages 599–602, Tours, France, 2006.

[25] R. Roy and T. Kailath. Esprit: Estimation of signal parameters via rotational invariance techniques. *IEEE Trans. Acoust. Speech, Signal Process.*, 37(7):984–995, 1999.

[26] Manlab. http://manlab.lma.cnrs-mrs.fr/spip/. Accessed: 2021-07-26.

[27] B. Cochelin. A path following technique via an asymptotic-numerical method. *Computers and Structures*, 53(5):1181–1192, 1994.
[28] B. Cochelin and C. Vergez. A high order purely frequency-based harmonic balance formulation for continuation of periodic solutions. *Journal of Sound and Vibration*, 324:243–262, 2009.

[29] M.S. Nakhla and J. Vlach. A piecewise harmonic balance technique for determination of periodic response of nonlinear systems. *IEEE Transactions on Circuit Theory*, 23(2):85–91, 1976.

[30] J. Gilbert, J. Kergomard, and E. Ngoya. Calculation of the steady-state oscillations of a clarinet using the harmonic balance technique. *J. Acoust. Soc. Am.*, 86(1):39–41, 1989.

[31] A. Lazarus and O. Thomas. A harmonic-based method for computing the stability of periodic solutions of dynamical systems. *Comptes Rendus Mécanique*, 338(9):510–517, 2010.

[32] L. Guillot, B. Cochelin, and C. Vergez. A generic and efficient taylor series-based continuation method using a quadratic recast of smooth nonlinear systems. *International Journal of Numerical Methods in Engineering*, 119(4):261–280, 2019.

[33] A. Bouhuys. Sound-power production in wind instruments. *J. Acoust. Soc. Am.*, 37(3):453–456, 1965.

[34] M. Campbell. Brass instruments as we know them today. *Acta Acustica united with Acustica*, 90:600–610, 2004.