# DYNAMICS-BASED CRITERIA FOR THE LIMITS OF VALIDITY OF CYLINDRICAL REED INSTRUMENT MODELS USING MODAL DECOMPOSITION

Tom Colinot, Etienne Gourc, Christophe Vergez

► To cite this version:

Tom Colinot, Etienne Gourc, Christophe Vergez. DYNAMICS-BASED CRITERIA FOR THE LIMITS OF VALIDITY OF CYLINDRICAL REED INSTRUMENT MODELS USING MODAL DECOMPOSITION. Forum Acusticum 2023, Sep 2023, Turin, France. hal-04799995

HAL Id: hal-04799995

https://hal.science/hal-04799995v1

Submitted on 23 Nov 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
# DYNAMICS-BASED CRITERIA FOR THE LIMITS OF VALIDITY OF CYLINDRICAL REED INSTRUMENT MODELS USING MODAL DECOMPOSITION

Tom Colinot<fnref n="1" />*

Étienne Gourc<fnref n="2" />

Christophe Vergez<fnref n="2" />

## ABSTRACT

Self-sustained instrument models such as clarinets or reed organ pipes include two main elements : a strongly non-linear exciter and a resonator, usually considered linear. In the past, the resonator has been modelled with either delay lines or modal decompositions, which are vastly different approaches. This work formulates objective criteria based on the self-sustained dynamics of the model, particularly the nature of its bifurcations, which allow to circumscribe the domain where the modal and delay lines approaches yield agreeing results. Notably, artifact fold bifurcations due to modal decomposition are exhibited in a specific region of the control parameter space. The dynamical characteristics of the models are studied through numerical continuation techniques (Harmonic Balance Method and Shooting) for periodic regimes and bifurcation analysis, as well as time-domain integration to include transient and aperiodic phenomena. The criteria introduced can also be applied when comparing the models to experimental data, or as tools to assess the performance of a physical model as a synthesizer.

**Keywords:** Modal decomposition - Numerical continuation - Bifurcation analysis - Time-domain integration

## 1. INTRODUCTION

The physical models of self-oscillating musical instruments, due to their strongly nonlinear nature, can produce a great variety of signals depending on the actions of the musician. Potentially, this diversity corresponds to that of a real instrument, which is among the features that make physical-based modelling a promising candidate to numerically emulate a real playing experience. Categorizing and describing the variety of playing regimes attainable with the physical model is therefore an essential task to bring it as close as possible to real instruments. However, comparing the physical model to a real instrument poses numerous challenges mostly unrelated to the dynamics themselves, such as a robust and trustworthy estimation of the real system's parameters. For this reason, and as a stepping stone towards the ultimate goal of reproducing a real instrument's dynamics, this paper compares the dynamics of two physical models differing only in the way that they model the resonator.

There are several approaches to the physical modeling of self-oscillating instruments, and notably their resonator (the vibrating string or air column which is put in motion by the action of the musician). This paper focuses on two types of methods: modal methods that represent the resonator as a combination of oscillators representing its modes [1,2], and methods which rather model wave propagation and reflection along the resonator using delays and filters [3]. A particular case of the second category uses the reflection function linking the backward and forward wave as the response of the resonator [4]. Each method has its advantages and drawbacks in terms of implementation and memory requirements, and with the correct pa-

<footnote n="1">¹ Buffet Crampon, 5 Rue Maurice Berteaux, 78711 Mantes-La-Ville, France</footnote>

<footnote n="2">² Aix Marseille Univ, CNRS, Centrale Marseille, LMA UMR7031, 4 impasse Nikola Tesla, 13453 Marseille Cedex 13, France</footnote>

<footnote>*Corresponding author: tom.colinot@buffetcrampon.com.
Copyright: ©2023 Colinot et al. This is an open-access article distributed under the terms of the Creative Commons Attribution 3.0 Unported License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.</footnote>
rameters can produce signals that reproduce experimental results (most commonly, static regimes [5, 6] or specific playing techniques [7, 8]), rendering them relevant as musical instrument substitutes in certain contexts.

To represent the dynamics of a system along one dimension of the parameter space, a bifurcation diagram is used. It is obtained by continuation using the shooting method (see section 2.2), which is applied to a discrete version of the equations of the model. Therefore, the system on which the continuation is performed can be used with no modifications in a real-time signal synthesizer.

## 2. PHYSICAL MODEL AND CONTINUATION METHOD

The model is that of a clarinet-like resonator (lossy cylinder) attached to an exciter system composed by a ghost reed [9] and a regularized version of the nonlinear characteristic [10]

$$u = \zeta \frac{h + |h|}{2} \text{sign}(\gamma - p) \sqrt{|\gamma - p|}, \quad (1)$$

where dimensionless parameters $\gamma$ and $\zeta$ represent the action of the musician through their blowing pressure and action on the reed respectively (see [11] for more detail). The dimensionless Kirchoff variables are the pressure $p$ in the mouthpiece and the flow $u$ entering the mouthpiece. They relate to the dimensionless forward and backward waves in the resonator as

$$p^{+} = \frac{p+u}{2} \quad (2)$$

$$p^{-} = \frac{p-u}{2}. \quad (3)$$

### 2.1 Time-domain integration

The discretization scheme is based on [12]. In the reflection function formalism, the resonator is taken into account using its reflection function $r$, leading to the backward waves $p^-$ being expressed through the discrete convolution

$$p^{-}[n] = \sum_{i=1}^{D} r[i] p^{+}[n-i]. \quad (4)$$

Note that with the simplified resonator used here, this expression is strictly causal ($r[0] = 0$).

In the modal formalism, the resonator is represented by $N_m$ oscillators whose modal pressures $p_k$ evolve as a

function of the flow rate $u$ and their past values:

$$p_k[n] = b_{0,k}u[n] + b_{1,k}u[n-1] \quad (5)$$

$$-a_{1,k}p_k[n-1] - a_{2,k}p_k[n-2] \quad (6)$$

Their evolution is decomposed into a strictly causal part $V_k$ and an instantaneous part depending on $u[n]$, such that

$$V_k = b_{1,k}u[n-1] - a_{1,k}p_k[n-1] - a_{2,k}p_k[n-2] \quad (7)$$

$$p_k[n] = b_{0,k}u[n] + V_k \quad (8)$$

Both formalisms lead to the apparition of the strictly causal part of the pressure $p$, which is denoted $V$ and is expressed in each case as

$$V = \sum_{k=1}^{N_m} V_k \quad \text{or} \quad V = 2p^{-}[n]. \quad (9)$$

The reed is considered to be a simple spring whose displacement is unaffected by the mouthpiece lay (according to the so called 'ghost reed' hypothesis). The reed position is

$$x = p - \gamma, \quad (10)$$

The height of the reed channel is given by the regularized positive part of the reed displacement

$$h = \frac{x+1+\sqrt{(x+1)^2+\eta}}{2}, \quad (11)$$

where $\eta$ is a regularization parameter set here at $\eta = 10^{-2}$. This parameter makes the reed closure phenomenon smoother which has two separate effects: it mitigates the grazing phenomena around $\gamma = 0.5$ which simplifies the model dynamics, and it limits extremely high harmonics generation from the contact. The fact that the contact between reed and mouthpiece is not rigid is coherent with the physics of the instrument, where the reed deforms when it comes in contact with the curved lay of the mouthpiece. The reed channel height $h$ is used afterwards to determine the flow $u$.

The flow entering the mouthpiece is determined following the method presented in [12], meaning that the nonlinear characteristic [10] is recast into

$$u[n] = 0.5S_{\Delta}W(-V_{1}W^{2} + \sqrt{V_{1}^{2}W^{2} + 4A_{\Delta}}), \quad (12)$$
where $W = \zeta h$ and in the non regularized version,

$$A_{\Delta} = |\gamma - V| \\ S_{\Delta} = \text{sign}(\gamma - V)$$

which are replaced here by

$$A_{\Delta} = \sqrt{(\gamma - V)^2 + \eta} \qquad (13)$$

$$S_{\Delta} = \frac{\gamma - V}{A_{\Delta}}. \qquad (14)$$

The pressure is then computed as

$$p[n] = Bu[n] + V, \qquad (15)$$

where in the modal formalism $B = \sum_{k=1}^{N_m} b_{0,k}$ and in the reflection function formalism $B = 1$.

## 2.2 Continuation method

In order to do a systematic exploration of the regimes available in the $\gamma$ parameter range, a continuation method is used [13]. A shooting method based on the method proposed in [14] allows to retain the time integration formalism exactly as it is used in a real-time formalism. This ensures that the regimes found in continuation are those accessible to the end user. Shooting is limited to periodic regimes, and requires an initialization based on a time-domain signal which is provided using time-domain integration. The retained variables are the initial conditions (in the phase space) and the period of the signal for each point. The phase space is of dimension $2N_m$ for the modal formalism, which accounts for the modal pressures and their derivatives. It is generally much larger for the reflection function formalism: although a continuous-time formalism would entail an infinite-dimensional phase space, the discretized version studied here is entirely determined by $D+1$ dimensions (all forward wave $p^+$ samples in the resonator and the current pressure $p$).

# 3. COMPARISON OF BIFURCATION DIAGRAMS

## 3.1 Studied system

For the compared resonators to be as close as possible, we study a single cylinder with a simplified loss model and perfect radiation condition [11]. The reflection coefficient is

$$R(\omega) = \exp(-\Gamma L) \qquad (16)$$

where $\Gamma = j\omega/c(1 + \alpha_1\sqrt{-2j}/r_v - j\alpha_2/r_v^2)$ and $r_v = R\sqrt{\omega/(c\ell_v)}$, $\ell_v = \mu/(\rho c)$, with the viscosity of air $\mu = 1.8071.10^{-5}$ kg.m$^{-1}$.s$^{-1}$ and the density $\rho = 1.2047$ kg.m$^{-3}$, and the length of the cylinder is $L = 0.3348$ m, which corresponds to 86 samples at 44.1 kHz sample rate. By inverse Fourier transform, one gets

$$r(t) = H(t - 2L/c) \frac{D_r}{\sqrt{\pi}} \exp\left(\frac{D_r^2}{t - 2L/c}\right) / (t - 2L/c)^{3/2}, \qquad (17)$$

where $H$ denotes the Heaviside step function, $D_r = \alpha_1 \frac{L}{R} \sqrt{\frac{2\ell_v}{c}}$. The discretized reflection function is $r[n] = r(n/F_s)$. The Fourier transform of $r$ is then transformed into an impedance, which is decomposed into modes. The number of retained samples for the reflection function is 101, at sample rate $F_s = 44100$ Hz. The comparison below is obtained with a relatively high value of $\zeta = 0.95$ which highlights the differences between the formalisms and the effect of the modal truncation.

## 3.2 Result for the reflection function formalism

Bifurcation diagrams are constructed by continuation, starting with a time-domain integration for $\gamma = 0.9$, which was found to consistently lead to a first register regime – i.e., whose frequency is around the first resonance peak. The comparison is made on only the first register branch, although other branches of solution exist, notably corresponding to second register regimes with fundamental frequency close to that of the second mode. For the sake of simplicity in this first approach, the stability of the solutions is not taken into account.

The bifurcation diagram of the reflection function formalism is represented in Figure 1. The amplitude of the pressure in the mouthpiece is represented as a function of the blowing pressure parameter $\gamma$. The branch starts at around $\gamma = 0.4$ with a direct Hopf bifurcation and rises steadily until a saturation point around $\gamma = 2.5$ quickly followed by a turning point, beyond which the branch no longer exists. The branch then rejoins equilibrium via an inverse Hopf bifurcation around $\gamma = 1$, which corresponds to the point where the reed becomes clamped against the mouthpiece in the static regime. This does not occur at precisely $\gamma = 1$ because of the regularization. This structure of bifurcation diagram for wind instruments is classical and found for instance in [15].
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Bifurcation diagrams for the reflection function formalism.</figcaption></figure>

### 3.3 Results for the modal formalism

The bifurcation diagrams of Figure 2 all correspond to the modal formalism, with a different number of modes $N_m$. All diagrams share common features: the location and nature of the Hopf bifurcation at the appearance and disappearance of the oscillations appear consistently at $\gamma \approx 0.4$ and $\gamma \approx 1$ respectively, and the first part of the branch after the apparition of the oscillation are very similar. However, things change around saturation for $\gamma > 2$. Here, instead of the single turning point that is observed on Figure 1, multiple turning points now appear. Their number is directly dependent on the number of modes taken into account, with two turning points being added by each addition of a mode. These double turning points become closer and closer together as more modes are added. For example, the addition of the sixth mode between the two lowest graphs results in two very close turning points being added just above $\gamma = 2$.

In terms of dynamic behavior, these turning points are quite important, because they mean the system is subject to multistability and therefore, hysteresis. Indeed, beyond $\gamma = 2$, several solutions coexist and their emergence depends on the systems' history. The more turning points, the more solutions potentially coexist, which increases the possibilities for hysteretic phenomena. However, the more modes, the closer turning points become, which shortens the range of control parameter $\gamma$ on which the hysteresis is felt, and reduces the difference between the coexisting regimes along the branch.

<figure><img src="image_3.png" /><figcaption>**Figure 2.** Bifurcation diagrams for the modal system.</figcaption></figure>

### 4. CONCLUSION

For a large value of the reed opening parameter $\zeta$, modal and waveguide based models of clarinet-like instruments exhibit fundamental differences in the nature of their first register solutions depending on the blowing pressure parameter. In particular, the modal formalism yields a bi-
furcation diagram which features multiple turning points, absent from its reflection function counterpart. There is a certain paradoxical evolution in the comparison of the bifurcation diagrams as one takes into account more and more modes in the modal formalism: on the one hand, more modes lead to more spurious turning points appearing on the diagram. On the other hand, these numerous turning points become much closer together, and the global shape of the diagram appears to be resembling its reflection function counterpart more and more. This can lead to a surprising conclusion: depending on the criterion for comparing dynamical responses, one could deem that a better approximation of the reflection function behavior is given by a low number of modes (when simply counting the bifurcations and comparing their nature) rather than a high number of modes.

## 5. ACKNOWLEDGMENTS

This work was supported by the French ANR LabCom Liamfi (ANR-16-LCV2-007-01)

## 6. REFERENCES

[1] J.-M. Adrien, "The missing link: modal synthesis," in *Representations of musical signals*, The MIT Press, 1991.

[2] P.-A. Taillard, F. Silva, P. Guillemain, and J. Kergomard, "Modal analysis of the input impedance of wind instruments. application to the sound synthesis of a clarinet," *Applied Acoustics*, vol. 141, pp. 271–280, 2018.

[3] S. Bilbao, A. Torin, and V. Chatziioannou, "Numerical modeling of collisions in musical instruments," *Acta Acustica united with Acustica*, vol. 101, no. 1, pp. 155–173, 2015.

[4] M. E. McIntyre, R. T. Schumacher, and J. Woodhouse, "On the oscillations of musical instruments," *The Journal of the Acoustical Society of America*, vol. 74, no. 5, pp. 1325–1345, 1983.

[5] J. Kergomard, P. Guillemain, F. Silva, and S. Karkar, "Idealized digital models for conical reed instruments, with focus on the internal pressure waveform," *The Journal of the Acoustical Society of America*, vol. 139, no. 2, pp. 927–937, 2016.

[6] A. Muñoz Arancón, B. Gazengel, J.-P. Dalmont, and E. Conan, "Estimation of saxophone reed parameters

during playing," *The Journal of the Acoustical Society of America*, vol. 139, no. 5, pp. 2754–2765, 2016.

[7] P. Guillemain, C. Vergez, D. Ferrand, and A. Farcy, "An instrumented saxophone mouthpiece and its use to understand how an experienced musician plays," *Acta Acustica united with Acustica*, vol. 96, no. 4, pp. 622–634, 2010.

[8] V. Chatziioannou, S. Schmutzhard, M. Pàmies-Vilà, and A. Hofmann, "Investigating clarinet articulation using a physical model and an artificial blowing machine," *Acta Acustica United with Acustica*, vol. 105, no. 4, pp. 682–694, 2019.

[9] T. Colinot, L. Guillot, C. Vergez, P. Guillemain, J.-B. Doc, and B. Cochelin, "Influence of the "ghost reed" simplification on the bifurcation diagram of a saxophone model," *Acta Acustica united with Acustica*, vol. [In press], no. -, pp. -, 2019.

[10] T. A. Wilson and G. S. Beavers, "Operating modes of the clarinet," *The Journal of the Acoustical Society of America*, vol. 56, no. 2, pp. 653–658, 1974.

[11] A. Chaigne and J. Kergomard, *Acoustique des instruments de musique*. Belin, 2008.

[12] P. Guillemain, J. Kergomard, and T. Voinier, "Real-time synthesis of clarinet-like instruments using digital impedance models," *The Journal of the Acoustical Society of America*, vol. 118, no. 1, pp. 483–494, 2005.

[13] A. H. Nayfeh and D. T. Mook, *Nonlinear oscillations*. John Wiley & Sons, 2008.

[14] M. Peeters, R. Viguié, G. Sérandour, G. Kerschen, and J.-C. Golinval, "Nonlinear normal modes, part ii: Toward a practical computation using numerical continuation techniques," *Mechanical systems and signal processing*, vol. 23, no. 1, pp. 195–216, 2009.

[15] J.-P. Dalmont and C. Frappe, "Oscillation and extinction thresholds of the clarinet: Comparison of analytical results and experiments," *The Journal of the Acoustical Society of America*, vol. 122, no. 2, pp. 1173–1179, 2007.