# Playability of Self-Sustained Musical Instrument Models: Statistical Approaches

Martin Pégeot, Tom Colinot, Jean-Baptiste Doc, Vincent Fréour, Christophe Vergez

► To cite this version:

Martin Pégeot, Tom Colinot, Jean-Baptiste Doc, Vincent Fréour, Christophe Vergez. Playability of Self-Sustained Musical Instrument Models: Statistical Approaches. Acta Acustica, 2024, 8, pp.74. 10.1051/aacus/2024075. hal-04777809v3

HAL Id: hal-04777809

https://hal.science/hal-04777809v3

Submitted on 30 Dec 2024

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
SCIENTIFIC ARTICLE

OPEN ACCESS

# Playability of self-sustained musical instrument models: statistical approaches

Martin Pégeot<fnref n="1" />,*, Tom Colinot<fnref n="2" /> ID, Jean-Baptiste Doc<fnref n="3" /> ID, Vincent Fréour<fnref n="4" />, and Christophe Vergez<fnref n="1" /> ID

Received 1 March 2024, Accepted 25 October 2024

**Abstract** – Self-sustained musical instruments, such as wind or bowed string instruments, are complex nonlinear systems. They admit a wide variety of regimes, which sometimes coexist for certain values of the control parameters. This phenomenon is known as multistability. With fixed parameters, the selection of a regime and the shape of the transient depend not only on the values of the control parameters, but also on the initial conditions. In this article, we focus on the statistical influence of initial conditions on regime selection and transient duration. An existing sample-based method called basin stability is presented to calculate the probability of occurrence of each regime. A second sample-based method is proposed for the calculation of the probability density function of transient durations. Additionally, a study taking into account specific control scenarios is presented to highlight the influence of the distribution of initial conditions considered for the statistical methods. These methods are presented on a Van der Pol oscillator seen as a prototypical musical instrument model. They are then applied to a physical model of trumpet, to demonstrate their potential for a high dimensional self-oscillating musical instrument. Finally, their interest regarding questions of playability is discussed.

**Keywords:** Self-sustained musical instruments, Multistability, Basin stability, Transient duration, Playability

## 1 Introduction

Musical instruments are complex dynamical systems. Some of them are self-sustained oscillators, meaning that a continuous energy supply can make them oscillate. Wind instruments and bowed string instruments belong to this category. These instruments admit a wide variety of regimes, which sometimes coexist for certain values of control parameters. This is known as multistability. Multistability is a common phenomenon in self-sustained oscillators. It has been observed experimentally, theoretically and numerically on a wide variety of instruments, e.g., single and double reed instruments [1–3], flutes [4], brass instruments [5], vocal folds [6] and bowed string instruments [7]. There are multiple works focusing on mapping the operating regimes of these instruments [8]. We find in [9, 10] among the first curves representing the evolution of the amplitude of a solution as a function of a control parameter. The stability of these solutions is also studied few years later [11]. A seminal work of Dalmont et al. [12] represents bifurcation diagrams of self-sustained musical

instruments. These diagrams show the evolution of known solution properties as a function of one or more control parameters. The stability of these solutions is given, as well as some key features, such as amplitude and fundamental frequency for periodic solutions. Thanks to these diagrams, one can predict the behavior of an instrument whose control parameters vary in a quasi-static manner – i.e., slowly in relation to the system’s rate of evolution. In the case of multistability, bifurcation diagrams predict hysteretic behaviors, such as those of reed and brass instruments [13] or flute-like instruments [14]. Some unexpected regimes are also predicted with this method, such as the ghost note in brass instruments [15] or the wolf tone of the cello [16]. Nevertheless, it is much more difficult to predict the regime of a multistable instrument when the control parameters vary rapidly. In particular, the blowing pressure dynamics are shown to have an influence over the regime selection in flute-like instruments [17]. Similar results are shown in [18] for the saxophone.

Similarly, once the control parameters are fixed, the time required to reach the steady state depends on the previous evolution of the control parameters. This phenomenon was observed by [19, 20] on the clarinet. Generally

<footnote n="1">¹Aix Marseille Univ, CNRS, Centrale Med, LMA UMR7031, 4 impasse Nikola Tesla, CS 40006, 13453 Marseille Cedex 13, France</footnote>

<footnote n="2">²Buffet Crampon, 5 rue Maurice Berteaux, 78711 Mantes-La-Ville, France</footnote>

<footnote n="3">³Laboratoire de Mécanique des Structures et Systèmes Couplés, 2 rue Conte, 75003 Paris, France</footnote>

<footnote n="4">⁴Yamaha Corporation, Research and Development Division, 10-1 Nakazawa-cho, Naka-ku, Hamamatsu, Shizuoka 430-8650, Japan</footnote>

<footnote>*Corresponding author: pegeot@lma.cnrs-mrs.fr</footnote>
speaking, the duration of attack transients [7, 21] or transi-
tions between notes [22] are important topics, since the
quality of a transient perceived by a musician or a listener
depends greatly on its characteristics and duration [23]. The
duration of the transients that gives rise to a sound is thus
studied on multiple self-sustained instruments, from the
trombone [24] to the cristal Baschet [25], using linear stabil-
ity analysis. Like the prediction of the playing regime, the
duration of the transient is estimated for control parameters
varying in a quasi-static manner. Linear stability analysis
cannot predict dynamic phenomena such as those observed
by [19, 20]. Finally, despite their high efficiency in the quasi-
static regime, current methods for analyzing self-oscillating
instruments are limited when considering rapidly varying
control parameters. By taking these rapid variations into
account, the dimensionality of the problem becomes infinite
since the parameters can follow any time series. Conse-
quently, it motivates researchers to limit their studies to
specific transient scenarios. For example, several authors
consider pressure ramps of variable slope and study the
influence of the pressure rise rate on the selected regime
[26, 27] or the transient duration [20]. In this article, we
study rapid variations in control parameters through
stochastic initial conditions, in cases of multistability
between an equilibrium and a periodic regime. We are
interested in situations where the musician moves quickly
from one quasi-static configuration to another, as is the case
during an attack [21, 28], a break or a change of note
[22, 29, 30]. Quasi-statistical analysis is not valid during this
transition and the new quasi-static configuration is studied
under unknown initial conditions. With this in mind, we
study the statistical influence of initial conditions on the
selected regime and on the duration of the transient. The
probability of obtaining each regime is estimated with an
already existing method called basin stability [31], and an
original approach is proposed to determine the probability
of obtaining each transient duration. Finally, the distribu-
tions of initial conditions considered for these statistical
methods are discussed. In particular, initial conditions
generated by specific transient control scenarios are studied.

The remainder of the document is organized as follows. The system chosen to introduce the methods is given in Section 2. Section 3.1 focuses on the basin stability metric and its use for musical instruments. The transient duration is investigated in Section 3.2 and the influence of the distribution of initial conditions is explored in Section 4. To highlight the interest of the methods on more sophisticated and high dimensional systems, an application to a physical model of trumpet is presented in Section 5. The different methods and results are discussed in Section 6. Perspectives and conclusions are drawn in Section 7.

2 Minimal multistable system

Most physical models of self-sustained musical instru-
ments have no exact analytical solutions for transient or
steady-state regimes. Therefore, numerical methods are
often used to study them, such as finite differences [32] or

harmonic balance method [33]. A number of authors have
proposed less detailed models with exact analytical solu-
tions and only few state variables. These simplified models
help to investigate the mechanisms behind self-oscillations
[34]. In order to compare our statistical approach with exact
analytical results, we have chosen to illustrate our methods
on such a system. In this article, we focus on a Van der Pol
oscillator with fifth-order nonlinearity, as described in [35].
This oscillator is characterized by a region of multistability,
where both an equilibrium and a periodic solution are
stable. This feature can also be observed on several musical
instruments, such as saxophones [27] or brass instruments
[5]. Its phase space is of dimension two and can therefore
be displayed simply, enabling more detailed dynamic
analysis. A technical publication [36] presents the musical
interest of this system and provides a demonstrator solving
it in real time (https://zenodo.org/records/8413627). Note
that even though the global behavior of musical instru-
ments is varied and rich, the system studied here can be
obtained by change of variables and polynomial expansions
from any system presenting this type of multistable behav-
ior. It is only used to illustrate the interest of a statistical
approach to predict the behavior of a multistable system
under transient excitations.

The motion of the considered oscillator is given by the
following equation:

$$
\ddot{x} + [-\mu + \sigma(x^2 + \dot{x}^2) + \nu(x^2 + \dot{x}^2)^2]\dot{x} + x = 0. \quad (1)
$$

We impose $\sigma = -1.5$ and $\nu = 0.1$ in order to obtain a multistable saxophone-like or brass-like behavior. The parameter $\mu$ is the control parameter of this simple “musical instrument” and its effect is analogous to the blowing pressure for wind instruments [36]. A harmonic balance method with only one harmonic gives the amplitude $X$ of solutions of the form:

$$
x = X \cos(t + \varphi). \tag{2}
$$

There are up to three solutions, depending on the value of $\mu$:

$$
\left\{
\begin{array}{ll}
X_{Eq} = 0, & \\
\\
X_{2\pi+} = \sqrt{\frac{-\sigma + \sqrt{\sigma^2 + 4\mu v}}{2v}} & \text{if } \mu \ge \mu_{SN}, \\
\\
X_{2\pi-} = \sqrt{\frac{-\sigma - \sqrt{\sigma^2 + 4\mu v}}{2v}} & \text{if } \mu_{SN} \ge \mu \ge \mu_H,
\end{array}
\right.
\qquad (3)
$$

with $\mu_H = 0$ and $\mu_{SN} = -\frac{\sigma^2}{4v}$. These solutions consist of an
equilibrium ($X_{Eq}$) and two $2\pi$-periodic solutions ($X_{2\pi+}$ and $X_{2\pi-}$), which are represented on a bifurcation dia-
gram in Figure 1. This oscillator shows an inverse Hopf
bifurcation at $\mu = \mu_H$ and a saddle-node bifurcation at
$\mu = \mu_{SN}$. Between these two bifurcations, the system is
multistable: both the equilibrium $X_{Eq}$ and the largest peri-
odic solution $X_{2\pi+}$ are stable. Saxophones [27] and brass
instruments [5] are likely to present this bifurcation
sequence and therefore this multistability behavior.

A traditional bifurcation diagram only depicts the
steady state solutions. The influence of initial conditions
on the transient and the steady-state behavior is concealed.
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Bifurcation diagram of a fifth order Van der Pol oscillator given in equation (1) with respect to the control parameter $\mu$, with $\sigma = -1.5$ and $\nu = 0.1$. Stable solutions are represented with continuous lines, unstable solutions in broken lines and bifurcations are indicated with star markers. The region of multistability is shaded in grey.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>**Figure 2.** Illustration of the sensitivity to the initial conditions, for the Van der Pol oscillator, regarding the steady state regime and the transient duration. On the left side, the phase plane of the system for $\mu = -3$, with trajectories coming from three close initial states (indicated with star markers). On the right side, the evolution of the position $x$ with respect to time (dimensionless).</figcaption></figure>

This influence can be seen on a phase plane for a given value of the control parameters. In Figure 2, we show the phase plane of the fifth order Van der Pol oscillator for $\mu = -3$, which is a situation of multistability. For the remainder of the paper, we use this value each time we need to choose a specific $\mu$ for illustration. The equation of motion (1) is numerically integrated with the function ode45 of Matlab [37], for three initial states very close to one another. In that case, the slight change in initial conditions leads the system either to the equilibrium (in blue) or to the stable limit cycle (in red). The two oscillating trajectories have different transient durations: one is almost instantaneous while the other lasts approximately one period.

<figure><img src="image_4.png" /><figcaption>**Figure 3.** Basins of attraction of the fifth order Van der Pol oscillator ($\mu = -3$). Initial conditions taken in the red shaded area give rise to the periodic solution whereas initial conditions taken in the blue shaded area give rise to the equilibrium solution.</figcaption></figure>

In short, Figures 1 and 2 illustrate that the steady state and the transient duration of this minimal self-sustained musical instrument are sensitive to the initial conditions and to the control parameter value around a bifurcation. This sensitivity is merely a mark of non linearity and is shared by most non linear dynamical systems. In the Sections 3.1 and 3.2, we present two statistical approaches to predict the selected regime and the transient duration for such a system.

# 3 Methods

## 3.1 Basin stability

Some self-sustained musical instruments are multi-stable. In such a case, once their control parameters are settled, their steady state depends exclusively on their initial state. The subset of initial conditions leading to a specific attractor is called its basin of attraction [38, 39]. In Figure 3, we depict the basins of attraction of the fifth order Van der Pol oscillator, for $\mu = -3$. The basin of the equilibrium is shaded in blue and the basin of the periodic solution is shaded in red. In that case, the two basins are separated by the unstable limit cycle of amplitude $X_{2\pi-}$, plotted in broken line and given in equations (1) and (3). In the general case, however, the geometry of a basin of attraction can be more complex and its boundaries difficult to find. We discuss this problem in Section 6. Moreover, in the scope of this article, we consider that the initial conditions resulting from a quick change in the control parameters are unknown. They are randomly taken within a subset of the initial state space. Since the initial state is unknown, the precise geometry of the basins of attraction does not suffice to predict its playing regime.
To take into account the uncertainty of the initial state, we propose to use the notion of *basin stability*, a probabilistic metric introduced by [31]. The principle is to evaluate the size of the basins of attraction with respect to the total size of the considered space of initial states. The result gives the probability to reach each attractor, considering a random initial state. Formally, considering a subset of interest of the space of initial conditions $Q$, an attractor $A$, its basin of attraction $B$ and a probability density function $\rho$, the basin stability $S_B(A)$ is defined as follows:

$$S_B(A) = \int_Q \mathbf{I}_B(\mathbf{X})\rho(\mathbf{X})d\mathbf{X}, \quad (4)$$

where

$$\mathbf{I}_B(\mathbf{X}) = \begin{cases} 1, & \text{if } \mathbf{X} \in B, \\ 0, & \text{otherwise.} \end{cases} \quad (5)$$

The probability density function $\rho$ can follow a specific distribution if the random choice of initial conditions is not uniform. This aspect is discussed in Sections 4 and 6. Due to the nature of $\rho$ as a probability density function over $Q$, the basin stability metric is always comprised between 0 and 1. $S_B(A) = 0$ indicates that the basin of attraction of the solution is statistically never reached. $S_B(A) = 1$ means that the basin of attraction of $A$ occupies all of the region of interest $Q$, and that $A$ is reached by every trajectory. All the illustrations of Section 3 are realized with a uniform distribution. The volumic integral in the expression of $S_B$, in equation (4), is rarely calculable exactly and we rather compute the corresponding discrete sum. In other words, we take $N$ random samples in $Q$ and we evaluate the proportion of samples belonging to the studied basin of attraction. The basin stability estimation, which is noted $\hat{S}_B$, then writes:

$$\hat{S}_B(A) = M/N, \quad (6)$$

where $M$ is the number of samples belonging to $B$, the basin of attraction of $A$. As pointed out in [31, 40, 41], this computation corresponds to $N$ independent trials with probability of success $S_B$. The resulting standard error due to sub sampling writes:

$$\text{err}(\hat{S}_B) = \sqrt{\frac{S_B(1 - S_B)}{N}}. \quad (7)$$

In practice, $S_B$ is unknown so it is replaced by its estimation $\hat{S}_B$ in equation (7). This error only accounts for sub-sampling. It is interesting to notice that the error does not depend on the dimension of the state space. As a result, the basin stability is a metric suitable for low dimension problems as well as for high dimension ones. In general, for more complex systems, the basin boundaries are unknown and the classification method would consist in time integrating the system and observing toward which attractor it would converge. Such strategy is illustrated in Figure 4. However, for this minimal multistable system, the basin boundary is analytically known: the boundary is the unstable periodic solution whose $L_2$-norm $X_{2\pi-}$ is given

<figure><img src="image_5.png" /><figcaption>**Figure 4.** Illustration of the classification method in the general case. The plotted trajectories are obtained with time integrations starting from the stars as initial conditions. Red trajectories are inside the basin of attraction of the stable periodic solution $X_{2\pi+}$ whereas the blue ones belong to the basin of the equilibrium $\hat{X}_{Eq}$.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 5.** Basin stability of the two stable solutions computed from 100 samples. In red: basin stability of the periodic solution. In blue: basin stability of the equilibrium. Error bars correspond to the standard error of the size estimation of the Monte-Carlo method due to sub-sampling equation (7). In black: basin stability computed with the analytic expression of the basin boundary (circle of known diameter $X_{2\pi-}$ given in (3)). The region of multistability is shaded in grey.</figcaption></figure>

in equation (3). As a result, it is much faster to classify the samples with this analytical expression rather than with time integrations.

Figure 5 gives the basin stability of the Van der Pol oscillator. To compute it, we arbitrarily choose a subset of the phase space $Q$ which includes all the attractors and which is independent of $\mu$. For $\mu \in [-7, 1]$ as in Figure 1, the subset $Q = ([-4, 4] \times [-4, 4])$ meets these conditions. To uniformly sample $Q$, we use a Latin Hypercube Sampling of 100 samples. The basin stability of the equilibrium is plotted in blue and the one corresponding to the periodic solution in red. The error bars correspond to the
standard error given in equation (7). In addition, within the multistability region, i.e., for $\mu_{SN} \ge \mu \ge \mu_H$, the exact basin stabilities write:

$$ \left\{ \begin{array}{l} S_B(X_{Eq}) = \frac{\pi X_{2\pi-}^2}{\text{area}(Q)}, \\ S_B(X_{2\pi+}) = 1 - S_B(X_{Eq}). \end{array} \right. \qquad (8) $$

These exact solutions are represented with black lines. They mainly remain inside the error bars estimated with $\hat{S}_B$ (i.e., $S_B$ replaced by $\hat{S}_B$ in Eq. (7)). This validates the error estimation of the sample-based method given in equation (7). Outside the multistability area, the only stable solution has a basin stability equal to 1 whereas the other one has a basin stability equal to 0. Regarding this metric, non existing solutions and unstable solutions are identical and have a zero basin stability. Inside the multistability area, both solutions have a non zero basin stability. The stability of the periodic solution increases with $\mu$ whereas the stability of the equilibrium decreases. One can observe that $S_B$ is continuous at the Hopf bifurcation but discontinuous at the saddle-node bifurcation. Indeed, in the first case, the basin of attraction of the equilibrium shrinks progressively until it becomes a dot in the phase space for $\mu = \mu_H$. In the second case, the periodic solution suddenly disappears at $\mu_{SN}$, even though its basin of attraction had previously occupied a large part of $Q$. The value of $S_B$ corresponds to the probability to obtain the related solution, given a random initial state inside $Q$. If a musician were to impose uniformly random initial conditions in $Q = ([-4, 4] \times [-4, 4])$, Figure 5 indicates that he/she would have at least a 60% chance of reaching the periodic regime in the steady state. However, the assumption of a musician peeking a random initial condition uniformly seems oversimplified. This topic is discussed in Sections 4 and 6.

Overall, the basin stability indicates which values of the parameters (here, the only parameter is $\mu$) fosters one solution or another. It gives the probability to produce a specific regime, considering a subset $Q$ of the phase space and a probability density function $\rho$ of $Q$. The result depends greatly on the choice of $Q$ and $\rho$. The standard error of this metric depends on the number of samples and on the relative size of the studied basin of attraction. It does not depend on the phase space dimension however and thus is very well adapted to high but finite dimensional systems.

## 3.2 Transient duration

In music, transient phenomena play a crucial role in sound and instrument discrimination [42–44]. Musicians may also have expectations concerning the characteristics of these transients, and in particular concerning their duration. For example, Guettler and Askenfelt [23] highlight the importance of the type and duration of violin attacks on their quality as perceived by professional string players. Several authors [21, 22, 25, 45] assume that musicians and instrument makers look for short transitions, whether for the onset or for sequences of linked notes. In this paper, we thus consider the reaction time of an

<figure><img src="image_7.png" /><figcaption>Figure 6. (a) Phase trajectories and (b) time evolution of the associated norm of the system’s state for 300 random initial conditions. Blue trajectories lead to the equilibrium and red trajectories lead to the periodic solution. The neighborhoods of the two stable solutions are shaded in grey. They are defined at a distance $\epsilon = 0.1$ in $L_2$-norm from the stable solutions.</figcaption></figure>

instrument – i.e., its tendency to produce short transients – as a playability criterion. In self-sustained instruments, the transient depends on the parameters of the system but also on the initial conditions, as can be seen on Figure 6. In this figure, the evolution of the fifth order Van der Pol oscillator has been computed for a large number of random initial states, all other things being equal. In Figure 6a, the trajectories are represented in the phase plane whereas in Figure 6b, the same trajectories are represented with respect to time. The transient part of a trajectory is defined between its initial state and the moment it reaches the neighborhood of an attractor. The transient duration is denoted $\tau$ hereafter. The neighborhoods of the two attractors are shaded in grey and are defined at a distance $\epsilon$ in $L_2$-norm from the attractor. Although the size of these neighborhoods has an impact on the transient duration results, for now it is arbitrarily set to $\epsilon = 0.1$. Figure 6b highlights the diversity of transient durations that can be obtained only
<figure><img src="image_8.png" /><figcaption>**Figure 7.** (Top) Cumulative probability function of the transient duration distribution. In ordinate, the proportion of trajectories with a transient duration below the value given in abscissa. (Bottom) Corresponding violin plots. The box plots indicate the extreme values, the median, the upper and lower quartiles whereas the curves around represent the probability density function of the transient duration distribution. A latin hypercube of 1000 samples is used and the transient durations are defined with $\epsilon = 0.1$.</figcaption></figure>

by picking different initial conditions. Some tendencies seem to appear. For instance, for this distribution of initial conditions, trajectories leading to the equilibrium are usually longer than those leading to the periodic solution.

These tendencies stand out when the distributions of the transient durations are represented graphically. Information on a distribution can be displayed as in Figure 7. The upper panel gives the cumulative distribution function, i.e., the proportion of trajectories with a transient duration shorter than the value $\tau_i$ represented in abscissa. The two lower panels depict the corresponding violin plots [46, 47], where the box plots give the extreme values (tips of the whiskers), the median and the two quartiles, and the shaded curves represent the probability density function of the distributions. To compute the latter, we used the kernel density estimator of Matlab (ksdensity). This statistical information is calculated for each attractor separately. As in Section 3.1, we use a Latin Hypercube Sampling to uniformly explore the space of initial states. With this sampling, we do not control the exact number of samples in each basin of attraction, but it can be estimated with the basin stability $S_B$ previously computed. As a result, for lowest values of $S_B$, the statistical transient analysis relies on only few samples and may not be very representative. It should be remembered that such misrepresentations occur only on highly improbable regimes. If a specific

<figure><img src="image_9.png" /><figcaption>**Figure 8.** Violin plots of the transient duration with respect to the control parameter $\mu$. In blue: equilibrium. In red: periodic solution. Number of samples: $N = 1000$. Panel (b) is a magnified version of the transients leading to the periodic solution. The region of multistability is shaded in grey.</figcaption></figure>

case calls for in-depth study of these regimes, the region of interest $Q$ can be adapted, leaving the method otherwise unchanged. The cumulative distribution function lends itself well to interpretation in musical context, where one could define a longest acceptable transient duration. The probability of the transient being shorter than this upper limit can be read directly on the cumulative distribution function, under the hypothesis of random initial conditions following a uniform law.

We use the violin plot representation in Figure 8 to depict the evolution of the transient duration distribution with respect to the control parameter $\mu$. Plotting the probability density function over the boxplot is particularly interesting when there are multiple maxima in the distribution since multiple maxima are not visible on a box plot. For instance, some examples of bi-modal distributions can be observed around the saddle-node bifurcation (for $\mu \approx \mu_{SN}$).

Figure 8 shows that most transients leading to the periodic solution are very short compared with its period ($T = 2\pi$). Indeed, the violin plots are centered around low values of $\tau$ (a magnified view for $\mu = -3$ is visible in Fig. 8b). However, the maximum transient durations are markedly larger than the median duration, up to 800 times longer for $\mu = -0.26$. This indicates that, for most values of $\mu$, a few initial conditions have extremely long transients.
<figure><img src="image_10.png" /><figcaption>**Figure 9.** Control scenarios used to compute the transient-informed initial conditions. The multistability region is shaded in gray. $\mu_0$ is the initial value of $\mu$, $\mu_{max}$ is the overshoot, $\mu_{end}$ is the final value of $\mu$, at which the basin stability is computed, $\mu_H$ and $\mu_{SN}$ are the value of $\mu$ at the Hopf and at the saddle-node bifurcations and $\tau_{att}$ is the attack duration.</figcaption></figure>

In fact, these particular initial states are very close to the unstable limit cycle. In that case, the trajectory follows the unstable limit cycle for some time while slowly diverging from it, which lengthens the transient duration. These rare and very long durations have particularly large values in the vicinity of the Hopf bifurcation. For the equilibrium, transient durations are longer and distributed more uniformly for all the values of $\mu$. Outside the vicinity of the two bifurcations, the transient durations leading to both solutions statistically decrease as $\mu$ increases. This is coherent with the fact that $-\mu$ is the linear part of the damping. As $\mu$ increases, the system becomes less and less damped until $\mu$ eventually becomes positive and energy is added to the system. Near the saddle-node bifurcation, this trend does not apply. The equilibrium becomes the only stable solution but its transient duration globally increases around $\mu_{SN}$. This is a usual phenomenon around a saddle-node bifurcation. It is interpreted by Strogatz in [38] as the remaining *ghost* of a nearby attractor that does not exist anymore. Around the Hopf bifurcation, the periodic solution becomes the only stable solution but its maximum transient duration significantly increases. However, apart from these very few long transients, Figure 8b shows that the overall trend does not change significantly for the periodic solution around the Hopf bifurcation. Considering initial conditions chosen from a uniform distribution, these long transients are highly improbable.

## 4 Influence of the distribution of initial conditions

In order to provide a more meaningful analysis in a context of musical performance, we consider another distribution for the initial conditions, informed by temporal evolutions of the control parameter. More precisely, we now consider initial conditions resulting from a transient

variation of the control parameter $\mu$, just before it reaches a constant target value. Our goal is therefore to consider what we call hereafter a “transient-informed” distribution and to analyze how it affects the basin stability results. To obtain a transient-informed distribution, we compute the time evolution of the system under several transient variations of $\mu$. These transient-informed distributions could also be used to compute the transient duration statistics, but we only focus on the basin stability here, since the results are more remarkable.

### 4.1 Control scenarios

Systems such as the fifth order Van der Pol oscillator presenting an inverse Hopf bifurcation followed by a saddle-node bifurcation produce a sound of non-zero amplitude when the control parameter exceeds the Hopf bifurcation (this is the case for the trumpet and the saxophone for certain sets of parameters). Nevertheless, it is possible to reduce the sound amplitude afterward by decreasing the control parameter under the Hopf bifurcation value, into the multistability region. To produce a low-volume sound at the onset, one can attempt to realize this gesture quickly by reducing the control parameter, and move the system into the multistable region before the limit cycle is attained. However, this strategy may fail if the system remains too close to the equilibrium during the transient. It would become trapped in the equilibrium’s basin of attraction when the parameter stops varying, which would result in no stable sound being produced. We consider piecewise linear control scenarios, as depicted in Figure 9. These scenarios are defined by four parameters: their initial value $\mu_0$, their maximal value $\mu_{max}$ that corresponds to what can be called the overshoot value, their final value $\mu_{end}$ at which the basin stability is computed and their duration $\tau_{att}$ that corresponds to the “attack duration”. For simplicity, we choose the same slope (in absolute value) for the increasing and the decreasing phases (before and after the overshoot).

The order of magnitude of $\mu_0$ is based on the bifurcation diagrams of the trumpet models presented in [5]. The objective is to give as much importance to the multistability as in these trumpet models. Therefore $\mu_0$ is set to have a similar ratio between the size of the multistability region [$\mu_{SN}$, $\mu_H$] and the size of the stability region of the equilibrium [$\mu_0$, $\mu_H$] as in the models of [5]. $\mu_0$ is then set to $\mu_0 = -25$. At the end of the attack, $\mu = \mu_{end}$ and the basin stability is computed for this value of $\mu$. The overshoot $\mu_{max}$ and the duration of the attack $\tau_{att}$ are the two remaining parameters of the control scenarios. In the following, we study the influence of a given overshoot $\mu_{max}$ on the regime selection, while the attack duration $\tau_{att}$ remains free to take different values. We arbitrarily choose $\mu_{max} \in [1, 10]$ and $\tau_{att} \in [1, 20]$ (dimensionless values). In [21], Ernoult and Fabre measured on three recorder players pressure rise times going approximately from 5 to 150 ms, for all fingerings between C5 (1046.5 Hz) and C6 (2093 Hz). Therefore, by choosing $\tau_{att} \in [1, 20]$, we keep a ratio between the longest and the shortest attacks with the same order of magnitude as in [21].
<figure><img src="image_11.png" /><figcaption>**Figure 10.** Example of a transient-informed distribution generated with 10 control scenarios of different durations. (•) Pre-initial Latin Hypercube of size $N_0 = 10$; (●) states of the system when $\mu = \mu_{max}$ (●) transient-informed distribution at the end of the attack. The black dotted line is the fast eigendirection of the equilibrium at $\mu = \mu_{max}$ and the red one is its slow eigendirection at $\mu = \mu_{end}$.</figcaption></figure>

## 4.2 Distribution generation

To generate a transient-informed distribution of initial conditions for a specific $\mu_{end}$, we apply 10 control scenarios with different durations ($\tau_{att} \in [1, 20]$) and with a unique value of $\mu_{max}$, to $N_0 = 10$ pre-initial conditions chosen close to the equilibrium (Latin Hypercube Sampling taken inside $x \in [-0.1, 0.1]$ and $\dot{x} \in [-0.1, 0.1]$). As a result, we obtain a transient-informed distribution of $N = 100$ initial conditions depicted in Figure 10. Firstly, it can be observed that this distribution is far from being uniform; most transient-informed initial conditions are very close to either the equilibrium or the periodic solution. It is expected for such a system to leave an unstable equilibrium through its fast eigendirection (i.e., associated with the eigenvalue of highest modulus) and to converge towards a stable equilibrium through its slow eigendirection (i.e., associated with the eigenvalue of smallest modulus). Such behavior is explained by Strogatz in example 5.2.3, p.133 of [38]. Precisely, around the equilibrium in Figure 10, the states of the system at $\mu = \mu_{end}$ are gathered along a specific direction whereas the states of the system at $\mu = \mu_{max}$ are regrouping along another direction. In addition, we plot on the same figure the fast eigendirection of the equilibrium at $\mu = \mu_{max}$ and its slow eigendirection at $\mu = \mu_{end}$. At the overshoot, the system clearly leaves the equilibrium through its fast eigendirection. At the end of the attack, some of the trajectories converge toward the equilibrium. If the simulation was performed over a longer duration, they would eventually gather along the slow eigendirection of the equilibrium.

<figure><img src="image_12.png" /><figcaption>**Figure 11.** Basin stability of the periodic solution computed with different distributions. The red curves are computed with transient-informed distributions of initial conditions, each curve corresponding to a specific overshoot $\mu_{max}$. The black curve is computed with a uniform distribution (it is the same as in Fig. 5 without the error bars).</figcaption></figure>

## 4.3 Basin stability with transient-informed distribution

The basin stability is then computed using these transient-informed distributions. The results are shown in Figure 11 for three different overshootvalues. For the sake of readability, only the basin stability of the periodic solution is represented (we recall that the basin stability of the equilibrium is the complementary to one).

The result computed with the uniform distribution is represented in black. Firstly, this Figure 11 highlights that the basin stability depends greatly on the considered distribution, even though some features are preserved: $S_B(\mu)$ is discontinuous at $\mu = \mu_{SN}$, continuous at $\mu = \mu_H$ and monotonously increasing for the cases where $\mu_{max} = 3$ and 4. In the case where $\mu_{max} = 2$, the overshoot scenario does not suffice to leave the basin of the equilibrium and no oscillations are observed in the multistable region. Secondly, for the situation considered, this Figure 11 shows that the large-overshoot attacks are more likely to make the system end up on the periodic solution. Moreover, the basin stability increases faster when $\mu_{end}$ gets closer to $\mu_H$. In terms of playability, a slowly varying $S_B$ can be interpreted as a region where the difficulty to produce a sound does not depend much on the final value of the control parameter. For the situation under study, this means that for a given value of $\mu_{max}$, all the attacks with $\mu_{end}$ in this region of slow varying $S_B$ have almost the same probability to converge towards the periodic solution. As a result, if we ignore the difficulty for a musician to maintain $\mu > \mu_{SN}$, there is no additional difficulty to reducing the sound amplitude to its minimal value.

# 5 Application to a trumpet model

In order to evaluate the interest of the methods introduced in previous sections to a more advanced and high
dimensional model of musical instrument, the analysis of
basin stability and transient durations is applied in this
section to a physical model of B, trumpet.

**5.1 Presentation of the model**

The model used is identical to the one considered in [5]. It is detailed in Appendix A. The lips are represented by a one-degree-of-freedom oscillator. The air column inside the body of the instrument, also known as the resonator, is modeled by a modal truncation of its input impedance. These two elements are coupled by a nonlinear function that models the air jet through the lip channel. In this article, we consider an 11-mode truncation of the resonator and fix the resonance frequency of the lips around the second impedance peak ($f_L = 200$ Hz and $f_2 = 232.7$ Hz). This peak corresponds to the lowest note of the instrument, in "normal" playing conditions (the first peak being associated with the "pedal note"). The variables of this model are the position and velocity of the lip as well as the real and imaginary parts of the 11 modal pressures. This results in a system of dimensionality 24. The parameter values used in this article are given in Appendix A.

**5.2 Bifurcation diagram**

Unlike the fifth-order Van der Pol oscillator, this trum-
pet model has no analytical solutions. To determine its solu-
tions, numerical continuation can be used. Here, we use the
Manlab software [33]. The bifurcation diagram of the model
is shown in Figure 12. The amplitude of the mouthpiece
pressure *p* and the frequency of the sound produced *f*<sub>play</sub>
are plotted with respect to the blowing pressure *p*<sub>0</sub>. In this
zone, the system exhibits an inverse Hopf bifurcation fol-
lowed by a saddle-node bifurcation, like the Van der Pol
oscillator presented in previous sections. The periodic
regime is associated with the second mode of the resonator,
represented on the bottom panel of Figure 12 by a horizon-
tal dotted line.

5.3 Basin stability and transient durations computed
with bSTAB

To calculate the system’s basin stability and transient
durations, we use the bSTAB toolbox [48]. For our applica-
tion, we have made a number of modifications to this tool-
box, which are detailed in Appendix B. These modifications
rely on convergence criteria in order to identify when the
system has reached steady state. The numerical integration
lasts until these criteria are met, rather than for an arbi-
trary time which would be identical for all trajectories.
These modifications allow us both to reduce the duration
of the time integrations required to compute basin stability,
and to calculate the duration of transients (which is not
provided in the original toolbox). These criteria must be
defined carefully, in order to avoid samples missclassifica-
tion, leading to errors in the estimation of basin stability
and transient durations.

<figure><img src="image_13.png" /><figcaption>**Figure 12.** Bifurcation diagram of the 11-mode trumpet model, calculated with Manlab. Top pane: amplitude of the oscillations $||p||_{L_2}$ versus blowing pressure $p_0$. Bottom pane: fundamental frequency of the oscillations $f_{play}$ versus blowing pressure $p_0$. Stable solutions are shown as solid lines, unstable solutions as dotted lines. Bifurcations are represented by stars.</figcaption></figure>

**5.3.1 Initial conditions from a uniform distribution**

As for the Van der Pol oscillator, we first consider a
uniform random distribution of initial conditions Q_u
selected inside a hyper-rectangle. This hyper-rectangle is
chosen large enough to contain all the solutions, whatever
the value of p₀. Its bounds are set as the extrema of the
stable periodic solution at p₀ = 0.8 kPa, rounded up to
the nearest integer value.

The basin stability calculated with Q_u for N = 300
samples is given in Figure 13a, the transient duration distri-
butions are given in Figure 13b and a zoomed-in view is
given in Figure 13c. The basin stability shares some features
with the Van der Pol oscillator. Indeed, it is discontinuous
at the saddle-node bifurcation and remains continuous at
the Hopf bifurcation. The observations and comments
made about this feature on the simple Van der Pol oscilla-
tor translate directly to this more complicated system. In
addition, S_B is monotonous. However, unlike the Van der
Pol oscillator, for dynamical systems of dimensionality
three or more, the basins' boundaries are different than
the unstable periodic solutions. These objects even have
different dimensionalities (the basins' boundaries are hyper-
surfaces and the limit cycles are hyperlines). Consequently,
the basin stability of the equilibrium cannot be inferred
from a size estimate of the unstable periodic solution and
the sampling approach is then required.

Concerning transient durations, the Hopf bifurcation
induces a few long transients in its vicinity, leading both
to the equilibrium and to the periodic solution. This is prob-
ably related to the increasing timescale of the slowest
eigendirection whose stability is reversed at the Hopf
<figure><img src="image_14.png" /></figure>

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

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

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

<figure><img src="image_18.png" /><figcaption>**Figure 13.** Basin stability (a) and transient duration distributions (b) of the trumpet model calculated with a uniform $Q_u$ distribution of $N = 300$ initial conditions. The dotted lines in (b) indicate the limit of the scale used for the zoomed view in (c). Error bars are given by the equation (7).</figcaption></figure>

bifurcation. Eventually, the timescale of this eigendirection goes to the infinity at the Hopf bifurcation. The ghost effect around the saddle-node bifurcation, where the periodic solutions do not exist yet, is also visible but less than for the Van der Pol oscillator. This may be due to the fact that in a phase space of dimensionality two, for a value of the control parameter slightly below the value of the saddle-node, the trajectories initiated far away from the equilibrium necessarily have to enter the region of influence of this ghost limit cycle. On the contrary, in a phase space of higher dimensionality, these trajectories can pass away from this region. Moreover, the median and the interquartile range of both solutions increase around the Hopf bifurcation (even though it is hardly noticeable for the periodic solution with that scale).

Finally, this trumpet model sometimes exhibits transients of more than several seconds, which is extremely long compared to usual musical timescales. Beyond a certain duration, these asymptotic regimes can be considered unplayable or at least difficult to play. However, the regime that is heard during one of these extremely long transients can be treated as a playable regime by the musician, even if it is unstable. In this respect, one could see a percussion instrument as a dynamic system with only one asymptotic solution, the equilibrium, but whose transients are long enough to be used to make music. We could then consider a “transient” category for samples with a transient duration exceeding a certain “musical duration”, rather than classifying them according to their asymptotic solution.

### 5.3.2 Initial conditions from a distribution based on archetypal time variations of $p_0$

In the following, we study the influence of attack transients on basin stability. We apply the same method as in Section 4 to generate distributions $Q_{ib}$ based on transient evolutions of $p_0$. Again, we consider pressure rises with overshoot (see Fig. 9), in the manner of a strong tonguing attack or a sforzando, as can be found in [49]. The configuration and the parameters chosen here differ from [49]. Consequently, we consider slightly different control transients than those measured in that article, so that $p_{0,max} \in [1, 10]$ kPa and $\tau_{att} \in [10, 500]$ ms. Ten attack transients are applied to each sample of a uniform distribution of $N_0 = 100$ pre-initial conditions. This leads to a final distribution of $N = 1000$ initial conditions. Each variable of the pre-initial distribution is selected in ranges 0.3 times smaller than in Section 5.3.1. These ranges are chosen arbitrarily and have a significant influence on the results.

In Figure 14a, we present the influence of the attack duration $\tau_{att}$ on the behavior of the system. Each curve has a fixed value of $\tau_{att}$ and 10 different values of $p_{0,max} \in [1, 10]$ kPa. It is striking to see that the basin stability does not evolve monotonously with $\tau_{att}$. Indeed, one would expect that the longer the system stays in the region of monostability of the periodic solution, the higher the chances it has to end up in its basin of attraction.
<figure><img src="image_19.png" /><figcaption>**Figure 14.** Basin stability of the trumpet model calculated with distributions $Q_{tb}$ of $N = 1000$ initial conditions based on attack transients (10 scenarios applied to each of the $N_0 = 100$ pre-initial conditions). (a) Each curve corresponds to a specific attack duration $\tau_{att}$. (b) Each curve corresponds to a specific overshoot value $p_{0,max}$. Only the basin stability of the periodic regime is shown.</figcaption></figure>

However, when $\tau_{att}$ increases from 10 to 50 ms, $S_B$ decreases. $S_B$ does not evolve much for $\tau_{att} \in [50, 100]$ ms, and it then increases, as expected, for $\tau_{att} > 100$ ms. These observations can be analyzed by considering recent results concerning dynamical bifurcations [50]: how close the system has come to the equilibrium when the Hopf bifurcation is crossed is crucial to predict the dynamics above the bifurcation. This closeness is the result of the time spent by the system below the bifurcation but also of the number of significant digits used in the simulation (numerical noise is therefore inevitable). This makes direct interpretation not straightforward. In Figure 14b, we study instead the influence of the overshoot amplitude $p_{0,max}$ for various attack durations $\tau_{att} \in [10, 500]$ ms. In that case, $S_B$ increases with the overshoot value, until it reaches its maximal value for $p_{0,max} \in [5, 7]$ kPa. Hence, it seems that there are optimal values of transient control parameters to maximize $S_B$. Further study would be needed to understand the dynamical reasons leading to this non trivial behavior.

# 6 Discussion

## 6.1 Computational cost

The methods presented in this paper allow to compute the basin stability and the transient duration distributions of a system. These metrics describe the global, statistical behavior of a system and could be very useful to study musical instruments, notably their playability. However, in order to take into account all the potential behaviors of a system, their computation can be costly. As a result, the computation cost limits the possible applications of these methods, for parametric studies for example.

Initially, [31] proposed to use Monte-Carlo techniques to numerically estimate the basin stability. Some limitations of these techniques have been drawn [40] when applied on strange basins of attraction (e.g., fractal and riddled/intermingled). However, such basins remain to be exhibited on self-sustained musical instruments. Some authors also proposed enrichment to the metric, to take into account uncertainties and variations of the parameters [51, 52]. Additionally, improvements on classical Monte Carlo approaches, notably using machine learning tools, show promise to accelerate basin stability computation. This could be interesting for more complex musical instrument models. For example, [41] proposed a method based on support vector machines. The idea is to find the basins' boundaries with a limited number of samples. Work still needs to be done to determine which method should be used depending on the situation and on the classification cost of a single sample. Indeed, if the classification of each sample is fast, it may be slower to find basins' boundaries rather than to apply a simple Monte-Carlo size estimation. For instance, it is the case for the Van der Pol oscillator presented in this article.

If the transient duration statistics and the basin stability are both studied, the two metrics can be computed simultaneously, using the same time integrations (as in Sect. 5). The cost needed to compute the basin stability with Monte-Carlo methods only depends on the size of the basin as indicated by the absolute standard error given by equation (7). For systems with higher dimensions, that cost could increase if the classification cost does. Moreover, the relative standard error, which writes $err_{rel}(\hat{S}_B) = \sqrt{\frac{1}{N} \sum_{i=1}^{N} (S_B - \bar{S})^2}$, increases as $S_B$ decreases. Higher dimension systems may have more stable solutions which lead to smaller basin stability values. More samples would thus be needed to keep the relative standard error low. This standard error formula stands for independent trials with only two outcomes – the sample is either inside or outside the studied basin of attraction. The calculation of the transient duration distributions does not enter into this category since the outcome is a continuous random variable: a transient duration. With a sparse sampling, some transient behaviors that appear for specific initial conditions might be missed. Consequently, it could be interesting to estimate the convergence of the probability density functions in order to choose an adequate number of samples.

The statistical methods presented here can be applied applied at low computational cost on the fifth order
Van der Pol oscillator, because the two attractors and the boundary of the basins of attraction are analytically known and characterized by their $L_2$-norm. On more complex systems, the attractors and basins' boundaries usually do not have analytic expressions. The main difficulty encountered in carrying out the statistical analysis of the trumpet model presented in Section 5 was to define adequate classification and stopping criteria. Indeed, non-linear systems sometimes evolve with very varied dynamics, making it difficult to distinguish between a steady state and a transient regime with slow dynamics. This kind of behavior leads to significant differences in transient duration between different initial conditions, especially at a bifurcation, where one eigendirection becomes infinitely slow. These discrepancies in transient durations are particularly visible in Figure 13b, around the Hopf bifurcation. In this regard, the Manlab solutions were very helpful to define efficient classification and convergence criteria (see Appendix B). For systems with a wide range of transient durations, we would strongly recommend to set a convergence criteria rather than to use a fixed integration time, as this is the case in the version of bSTAB currently online. This would improve performances and limit missclassification.

## 6.2 Initial conditions and musical gestures

In Section 3, the presented statistical approach is based on a uniform distribution of initial conditions. This is a general distribution and it allows different instrument configurations to be compared. However, the resulting probabilities probably fail to translate the musician's playing experience. As a first attempt to tackle this issue, we study in Section 4 initial conditions resulting from an archetypal control gesture. This approach highlights the importance of the considered distribution of initial conditions and it can be interpreted in musical terms more directly than random initial conditions. We would like to draw the reader's attention to the fact that the results obtained with this transient-based approach not only depend on the transient control parameters, but also on the choice of the pre-initial conditions. Overall, applicative studies on playability would greatly benefit from a preliminary description of the control scenarios and, where relevant, initial conditions distributions. In that perspective, it could be of great interest to measure transient blowing scenarios provided by humans (musicians and non-musicians) in an experimental study. Moreover, as illustrated by Figure 10, dynamical systems are more likely to cross specific regions of their phase space during transients. Choosing distributions of initial conditions based on the slow and fast eigendirections (depending on the solutions that exist during the transient) could also be an interesting idea.

In this paper, we focus on initial conditions induced by the fast evolution of one control parameter. However, musicians may also impact the initial conditions by acting directly on a state variable. One example among many others is the use of the tongue in reed instruments to impose initial reed positions and velocities. Measurements of tongue-induced reed positions in a clarinet are presented in

[53] for example. To improve the statistical approaches, these initial conditions should also be taken into account.

## 6.3 Stability, transient duration and playability

Other studies were interested in predicting the transient duration of a self-sustained musical instrument [24, 25]. A common approach is to linearize the system around an unstable equilibrium and to consider the exponential growth of a perturbation. This exponential growth depends on the positive real part of the eigenvalues of the matrix of derivatives (i.e., the jacobian matrix). This method is better adapted for studying initial conditions close to an equilibrium solution. The Floquet theory might be used in the same way for initial conditions next to a periodic solution, although it has yet to be applied in that manner to a musical instrument model (for more details on the Floquet theory, refer to [54, 55]). In this article, we consider initial conditions not necessarily close to any solution. Hence, these linear analyses might not be valid for every sample. On more complex systems with several competing multi-stable solutions, its local nature prevents it from giving a complete analysis of the transients. However, the systematic and fast nature of linear analysis makes it a good candidate to bolster certain sections of a complete statistical study.

Finally, the statistical study of the trumpet model presented in Section 5 highlights several interests of the methods. The long transients observed in Figure 13 raise the question of the playability of a regime as a function of its transient duration. In a standard musical context, we argue that a regime that takes several seconds to establish cannot really be considered playable. Therefore, certain stable regimes could actually be unplayable due to transient behaviors, which shows the relevance of nuancing the notion of stability using transients when studying playability. The transient-based analysis of Section 5.3.2 also indicates that the transient control parameters of a musical gesture may have nontrivial optimal values. Future works could analyze these types of transient control scenarios into more details, relying on the statistical approaches as playability guidelines.

## 7 Conclusion and perspectives

In this article, we present a sample-based approach that can be used to enrich a bifurcation diagram. In this study we use it to investigate two playability issues: the prediction of the steady state in multistable situations and the prediction of the transient duration. The method used for the first issue is called the basin stability [31], whereas the method proposed for describing the responsiveness of a dynamical system is new, to the authors' knowledge. These two methods rely on distributions of initial conditions that represent the transient action of a musician. The results strongly depend on this choice of distribution and we proposed a method based on time integration to generate transient-informed distributions.
The methods are then applied to a physical model of
trumpet with the aim to evaluate their interest in a more
practical case. The transient-based basin stability of this
model present non trivial tendencies regarding some tran-
sient control parameters. Moreover, very long transients
are highlighted, which raises the question of the playability
of such asymptotic regime. Overall, these statistical methods
show interest for the analysis of musical instruments, but
they can also be extended to a large variety of systems.

Future works will focus on two aspects. First, these
methods will be applied on musical instrument models in
configurations showing multiple oscillating regimes and rich
transients. Then, measurements on musicians will be con-
ducted to evaluate the initial conditions that they induce,
depending on the desired musical effect.

Conflicts of interest

The authors declare that they no conflicts of interest in
relation to this article.

Data availability statement

No new data were created or analyzed in this study.

References

1. T. Idogawa, M. Shimizu, M. Iwaki: Acoustical behaviors of an oboe and a soprano saxophone artificially blown. Some problems on the theory of dynamical systems in applied science, in: Proceedings of the Institute for Mathematical Analysis, Institute for Mathematical Analysis, Kyoto University, 1992, pp. 71–93.

2. T. Idogawa, T. Kobata, K. Komuro, M. Iwaki: Nonlinear vibrations in the air column of a clarinet artificially blown. Journal of the Acoustical Society of America 93, 1 (1993) 540–551.

3. K.Y. Takahashi, H.A. Kodama, A. Nakajima, T.A. Tachi-
bana: Numerical study on multi-stable oscillations of wood-
wind single-reed instruments. Acta Acustica united with
Acustica 95, 6 (2009) 1123–1139.

4. J. Bocanegra, D. Borelli: Review of acoustic hysteresis in
flute-like instruments, in: Proceedings of 26th International
Congress on Sound and Vibration, Canada, Montreal, 7–11
July, 2019.

5. V. Fréour, L. Guillot, H. Masuda, E. Tominaga, Y. Tohgi, C. Vergez, B. Cochelin: Numerical continuation of a physical model of brass instruments: Application to trumpet comparisons. Journal of the Acoustical Society of America 148, 2 (2020) 748–758.

6. J.J. Jiang, Y. Zhang: Chaotic vibration induced by turbulent noise in a two-mass model of vocal folds. Journal of the Acoustical Society of America 112, 5 (2002) 2127–2133.

7. R.T. Schumacher, J. Woodhouse: The transient behaviour of models of bowed-string motion. Chaos: An Interdisciplinary Journal of Nonlinear Science 5, 3 (1995) 509–523.

8. T.A. Wilson, G.S. Beavers: Operating modes of the clarinet. Journal of the Acoustical Society of America 56, 2 (1974) 653–658.

9. N.H. Fletcher: Nonlinear interactions in organ flue pipes. Journal of the Acoustical Society of America 56, 2 (1974) 645–652.

10. B. Lawergren: On the motion of bowed violin strings. Acta Acustica united with Acustica 44, 3 (1980) 194–206.

11. C. Maganza, R. Caussé, F. Laloë: Bifurcations, period doublings and chaos in clarinetlike systems. Europhysics Letters 1, 6 (1986) 295.

12. J.P. Dalmont, J. Gilbert, J. Kergomard: Reed instruments, from small to large amplitude periodic oscillations and the Helmholtz motion analogy. Acta Acustica united with Acustica 86, 4 (2000) 671–684.

13. J. Gilbert, S. Maugeais, 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 (2020) 27.

14. S. Terrien, C. Vergez, B. Fabre: Flute-like musical instruments: a toy model investigated through numerical continuation. Journal of Sound and Vibration 332, 15 (2013) 3833–3848.

15. R. Mattéoli, J. Gilbert, S. Terrien, J.-P. Dalmont, C. Vergez, S. Maugeais, E. Brasseur: Diversity of ghost notes in tubas, euphoniums and saxhorns. Acta Acustica 6 (2022) 32.

16. E. Gourc, C. Vergez, P.-O. Mattei, S. Missoum: Nonlinear dynamics of the wolf tone production. Journal of Sound and Vibration 516 (2022) 116463.

17. S. Terrien, R. Blandin, C. Vergez, B. Fabre: Regime change thresholds in recorder-like instruments: influence of the mouth pressure dynamics. Acta Acustica united with Acustica 101, 2 (2015) 300–316.

18. S. Terrien, B. Bergeot, C. Vergez: Dynamic basins of attraction in a toy-model of reed musical instruments, in: Proceedings of Forum Acusticum, Turin, Italy, 11–15 September, 2023.

19. B. Bergeot, A. Almeida. C. Vergez, B. Gazengel: Measure-
ment of attack transients in a clarinet driven by a ramp-like
varying pressure, in: Acoustics 2012, April, Nantes, France,
2012. https://hal.science/hal-00810986.

20. B. Bergeot, A. Almeida, C. Vergez, B. Gazengel: Prediction of the dynamic oscillation threshold in a clarinet model with a linearly increasing blowing pressure. Nonlinear Dynamics 73 (2013) 521–534.

21. A. Ernoult, B. Fabre: Temporal characterization of experi-
mental recorder attack transients. Journal of the Acoustical
Society of America 141, 1 (2017) 383–394.

22. S. Logie, S. Bilbao, J. Chick, M. Campbell: The influence of transients on the perceived playability of brass instruments, in: Proceedings of 20th International Symposium on Music Acoustics, Sydney and Katoomba, 26–27 August, 2010.

23. K. Guettler, A. Askenfelt: Acceptance limits for the duration of pre-Helmholtz transients in bowed string attacks. The Journal of the Acoustical Society of America 101, 5 (1997) 2903–2913.

24. L. Velut, C. Vergez, J. Gilbert, 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, 1 (2017) 132–148.

25. A. Couineaux, F. Ablitzer, F. Gautier: Minimal physical model of the cristal Baschet. Acta Acustica 7 (2023) 49.

26. B. Bergeot, S. Terrien, C. Vergez: Predicting transient dynamics in a model of reed musical instrument with slowly time-varying control parameter. Chaos: An Interdisciplinary Journal of Nonlinear Science 34, 7 (2024) 073146.

27. T. Colinot, C. Vergez, P. Guillemain, J.-B. Doc: Multistability of saxophone oscillation regimes and its influence on sound production. Acta Acustica 5 (2021) 33.

28. M. Castellengo: Acoustical analysis of initial transients in flute like instruments. Acta Acustica united with Acustica 85, 3 (1999) 387–400.

29. A. Almeida, R. Chow, J. Smith, J. Wolfe: The kinetics and acoustics of fingering and note transitions on the flute.
Journal of the Acoustical Society of America 126, 3 (2009) 1521–1529.

30. A. Almeida, W. Li, E. Schubert, J. Smith, J. Wolfe: Recording and analysing physical control variables used in clarinet playing: A musical instrument performance capture and analysis toolbox (MIPCAT). Frontiers in Signal Processing 3 (2023) 1089366.

31. P.J. Menck, J. Heitzig, N. Marwan, J. Kurths: How basin stability complements the linear-stability paradigm. Nature Physics 9, 2 (2013) 89–92.

32. F. Silva, P. Guillemain, J. Kergomard, C. Vergez, V. Debut: Some simulations of the effect of varying excitation parameters on the transients of reed instruments. Proceedings of Meetings on Acoustics 19 (2013) 035058.

33. B. Cochelin, C. Vergez: A high order purely frequency-based harmonic balance formulation for continuation of periodic solutions. Journal of Sound and Vibration 324, 1–2 (2009) 243–262.

34. R. Bader: Musical instruments as synchronized systems, in: R. Bader (ed.), Springer handbook of systematic musicology, Springer, Berlin, Heidelberg, 2018, pp. 171–196

35. D. Dessi, F. Mastroddi, L. Morino: A fifth-order multiple-scale solution for Hopf bifurcations. Computers & Structures 82, 31–32 (2004) 2723–2731.

36. T. Colinot, C. Vergez: How to build a MATLAB demonstrator solving dynamical systems in real-time, with audio output and MIDI control. Acta Acustica 7 (2023) 58.

37. L.F. Shampine, M.W. Reichelt: The matlab ode suite. SIAM Journal on Scientific Computing 18, 1 (1997) 1–22.

38. S.H. Strogatz: Nonlinear dynamics and chaos: with applications to physics, biology, chemistry, and engineering, Addison-Wesley Publication, Reading, MA, 1994.

39. A.H. Nayfeh, B. Balachandran: Applied nonlinear dynamics: analytical, computational, and experimental methods, Wiley-VCH GmbH, Weinheim, Germany, 2008.

40. P. Schultz, P.J. Menck, J. Heitzig, J. Kurths: Potentials and limits to basin stability estimation. New Journal of Physics 19, 2 (2017) 023005.

41. Y. Che, C. Cheng, Z. Liu, Z.J. Zhang: Fast basin stability estimation for dynamic systems under large perturbations with sequential support vector machine. Physica D: Nonlinear Phenomena 405 (2020) 132381.

42. K. Siedenburg: Specifying the perceptual relevance of onset transients for musical instrument identification. Journal of the Acoustical Society of America 145, 2 (2019) 1078–1087.

43. K. Siedenburg, M.R. Schädler, D. Hülsmeier: Modeling the onset advantage in musical instrument recognition. Journal of the Acoustical Society of America 146, 6 (2019) EL523–EL529.

44. K. Siedenburg, S. McAdams: Four distinctions for the auditory “wastebasket” of timbre. Frontiers in Psychology 8 (2017) 1747.

45. P.M. Galluzzo: On the playability of stringed instruments, PhD thesis, University of Cambridge, 2004.

46. J.L. Hintze, R.D. Nelson: Violin plots: a box plot-density trace synergism. American Statistician 52, 2 (1998) 181–184.

47. S. Weöglarczyk: Kernel density estimation and its application. ITM Web of Conferences 23 (2018) 00037.

48. M. Stender, N. Hoffmann: bSTAB: an open-source software for computing the basin stability of multi-stable dynamical systems. Nonlinear Dynamics 107, 2 (2022) 1451–1468.

49. T. Bianco, V. Freour, I. Cossette, F. Bevilacqua, R. Caussé: Measures of facial muscle activation, intra-oral pressure and mouthpiece force in trumpet playing. Journal of New Music Research 41, 1 (2012) 49–65.

50. B. Bergeot, C. Vergez: Analytical prediction of delayed Hopf bifurcations in a simplified stochastic model of reed musical instruments. Nonlinear Dynamics 107, 4 (2022) 3291–3312.

51. P. Brzeski, M. Lazarek, T. Kapitaniak, J. Kurths, P. Perlikowski: Basin stability approach for quantifying responses of multistable systems with parameters mismatch. Meccanica 51 (2016) 2713–2726.

52. C. Mitra, J. Kurths, R.V. Donner: An integrative quantifier of multistability in complex systems based on ecological resilience. Scientific Reports 5, 1 (2015) 16196.

53. M. Pàmies-Vilà, A. Hofmann, V. Chatzioioannou: Analysis of tonguing and blowing actions during clarinet performance. Frontiers in Psychology 9 (2018) 617.

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

55. L. Guillot, A. Lazarus, O. Thomas, C. Vergez, B. Cochelin: A purely frequency based Floquet-Hill formulation for the efficient stability computation of periodic solutions of ordinary differential systems. Journal of Computational Physics 416 (2020) 109477.

Cite this article as: Pégeot M. Colinot T. Doc J.-B. Fréour V. Vergez C. 2024. Playability of self-sustained musical instrument models: statistical approaches. Acta Acustica, xx, xx. https://doi.org/10.1051/aacus/2024075.

# Appendix A

## Trumpet model

We consider the brass instrument model described in [5]. The convention used to represent the lip position is given in Figure A1. The dimensioned and unregularized equations are as follows:

$$ 
\left\{ 
\begin{array}{l}
\ddot{x} + \frac{\omega_L}{Q_L}\dot{x} + \omega_L^2(x - x_0) = \frac{p_0 - p}{\mu_L}, \\
\dot{p}_n - s_n p_n = Z_c C_n u \quad \forall n \in [1, N_m], \\
p = 2 \sum_{n=1}^{N_m} \mathfrak{R}(p_n), \\
u = W \sqrt{\frac{2p_0 - p}{p}} \cdot \text{sign}(p_0 - p) \cdot \max(x, 0),
\end{array}
\right. 
\quad (A1) 
$$

with

$$ 
\left\{ 
\begin{array}{l}
\omega_L = \sqrt{\frac{k}{m}}, \\
Q_L = \frac{\sqrt{km}}{c}, \\
\mu_L = \frac{m}{S}.
\end{array}
\right. 
\quad (A2) 
$$

$S$ is the surface area of the lip to which the pressures $p_0$ and $p$ are applied, $\mu_L$ is the surface mass of the lip, $Q_L$ its quality factor and $\omega_L$ its angular eigenfrequency. The resonator is represented by its input impedance, which is treated as a sum of $N_m$ modes with poles $s_n$ and residues $C_n$. The lip allows air to pass through a rectangular surface of height $x$ and width $W$. The system given in (A1) is scaled as follows:
<figure><img src="image_20.png" /><figcaption>**Figure A1.** Schematic diagram of the one-degree-of-freedom lip model. Assume $x=0$ when the lip is in the closed position and note $x_0$ the position of the lip at rest (the situation shown here).</figcaption></figure>

$$ \left\{ \begin{array}{l} \ddot{\tilde{x}} + \frac{\omega_L}{Q_L} \dot{\tilde{x}} + \omega_L^2 (\tilde{x} - 1) = \omega_L^2 (\gamma - \tilde{p}), \\ \dot{\tilde{p}}_n - s_n \tilde{p}_n = C_n \tilde{u} \quad \forall n \in [1, N_m], \\ \tilde{p} = 2 \sum_{n=1}^{N_m} \Re(\tilde{p}_n), \\ \tilde{u} = \zeta \sqrt{|\gamma - \tilde{p}|} \cdot \operatorname{sign}(\gamma - \tilde{p}) \cdot \max(\tilde{x}, 0). \end{array} \right. \quad (A3) $$

Variables affected by this scaling are marked with a ~ symbol and are defined as follows:

$$ \left\{ \begin{array}{l} \tilde{x} = \frac{x}{x_0}, \\ p_M = \mu_L \omega_L^2 x_0, \\ \tilde{p} = \frac{p}{p_M}, \\ \gamma = \frac{p_0}{p_M}, \\ \tilde{p}_n = \frac{p_n}{p_M}, \\ u_M = \frac{p_M}{Z_c}, \\ \tilde{u} = \frac{u}{u_M}, \\ \zeta = Z_c W \tilde{x}_0 \sqrt{\frac{2}{pp_M}}. \end{array} \right. \quad (A4) $$

Irregular functions appearing in the flow term are regularized:

$$ \left\{ \begin{array}{l} |\bullet|_r = \sqrt{|\bullet|^2 + \epsilon}, \\ \operatorname{sign}(\bullet)_r = \frac{\bullet}{|\bullet|_r}, \\ \max(\bullet, 0)_r = \frac{\bullet + |\bullet|_r}{2}. \end{array} \right. \quad (A5) $$

The regularization term $\epsilon$ is arbitrarily set to $\epsilon = 10^{-6}$.

Time is scaled by the first modal angular frequency $\omega_1 = \mathfrak{I}(s_1)$ and dimensionless time derivatives are denoted'. The time dimensionless poles, residues and lip angular frequency respectively write $\tilde{C}_n, \tilde{s}_n$ and $\tilde{\omega}_L$. Finally, we separate the real $R_n$ and imaginary $I_n$ parts of the modal pressures $p_n$, and the system we solve is as follows:

<table><caption>**Table A1.** Resonator modal parameters. Complex residues $C_n$ and poles $s_n$ and resonance frequencies $f_n$.</caption><thead><tr><th>n</th><th>C<sub>n</sub> (s<sup>-1</sup>)</th><th>s<sub>n</sub> (rad · s<sup>-1</sup>)</th><th>f<sub>n</sub> (Hz)</th></tr></thead><tbody><tr><td>1</td><td>744.6</td><td>-13.98 + i522.5</td><td>83.15</td></tr><tr><td>2</td><td>954.5</td><td>-22.42 + i1462</td><td>232.7</td></tr><tr><td>3</td><td>1335</td><td>-28.64 + i2187</td><td>348.1</td></tr><tr><td>4</td><td>2582</td><td>-37.64 + i2907</td><td>462.6</td></tr><tr><td>5</td><td>3140</td><td>-45.82 + i3658</td><td>582.1</td></tr><tr><td>6</td><td>4191</td><td>-49.82 + i4339</td><td>690.6</td></tr><tr><td>7</td><td>4013</td><td>-58.42 + i5029</td><td>800.4</td></tr><tr><td>8</td><td>2602</td><td>-66.77 + i5705</td><td>908.1</td></tr><tr><td>9</td><td>1278</td><td>-72.24 + i6459</td><td>1028</td></tr><tr><td>10</td><td>909.7</td><td>-94.40 + i7211</td><td>1148</td></tr><tr><td>11</td><td>620.7</td><td>-128.6 + i7931</td><td>1262</td></tr></tbody></table>

$$ \left\{ \begin{array}{l} \tilde{x}'' + \frac{\hat{\omega}_L}{Q_L} \tilde{x}' + \hat{\omega}_L^2 (\tilde{x} - 1) = \hat{\omega}_L^2 (\gamma - \tilde{p}), \\ \tilde{R}'_n = \mathfrak{R}(\hat{s}_n) \tilde{R}_n - \mathfrak{I}(\hat{s}_n) \tilde{I}_n + \mathfrak{I}(\hat{C}_n) \tilde{u} \\ \phantom{\tilde{R}'_n =} \quad \forall n \in [1, N_m], \\ \tilde{I}'_n = \mathfrak{I}(\hat{s}_n) \tilde{I}_n + \mathfrak{I}(\hat{s}_n) \tilde{R}_n + \mathfrak{I}(\hat{C}_n) \tilde{u} \\ \phantom{\tilde{I}'_n =} \quad \forall n \in [1, N_m], \\ \tilde{p} = 2 \sum_{n=1}^{N_m} \mathfrak{R}(\tilde{p}_n), \\ \tilde{u} = \zeta \sqrt{|\gamma - \tilde{p}|_r} \cdot \operatorname{sign}(\gamma - \tilde{p})_r \cdot \max(\tilde{x}, 0)_r. \end{array} \right. \quad (A6) $$

The parameter values chosen for this article are as follows:

$$ \left\{ \begin{array}{l} \omega_L = 2\pi \times 200 \text{ rad} \cdot \text{s}^{-1}, \\ \mu_L = 2 \text{ kg} \cdot \text{m}^{-2}, \\ Q_L = 3, \\ x_0 = 1 \cdot 10^{-4} \text{ m}, \\ W = 8 \cdot 10^{-3} \text{ m}, \\ Z_c = 1.83 \text{ kg} \cdot \text{s}^{-1} \cdot \text{m}^{-4}, \\ \epsilon = 10^{-6}. \end{array} \right. \quad (A7) $$

## Appendix B

### Modifications to bSTAB

The bSTAB toolbox [48] has been designed to calculate the basin stability of any dynamic system as automatically as possible. It generates initial conditions, performs a time integration for each of them, then classifies the obtained regime by comparing it to reference signals. In our case, we compare the peak-to-peak amplitude of $R_2$, the real part of the second modal pressure, to that of the MANLAB solutions. Since classification must be performed on the steady part of the regime, the time integration continues until the following convergence criteria are met.

During the last duration $t^*$:

1. The amplitude of $R_2$ is within $\epsilon$ of that of a reference solution;

2. The difference between the amplitude of $R_2$ and that of the reference decreases.

The transient duration is defined as the duration after which the envelope of $R_2$ stays within $\epsilon$ of the amplitude of the reference solution. In this article, we choose $\epsilon = 0.1$.

The computation of the transient duration is illustrated in Figure B1a. The corresponding trajectory is also represented in Figure B1b, with its projection along $x$, $R_2$ and $R_4$.
<figure><img src="image_21.png" /><figcaption>**Figure B1.** (a) Computation of the transient duration (represented by the dotted line) and (b) projection of the corresponding trajectory along *x*, *R*₂ and *R*₄. The red dot is the initial state, the blue cross is the equilibrium, the red line is the stable periodic solution and the black broken line is the unstable periodic solution.</figcaption></figure>