A PHYSICAL MODEL OF THE TROMBONE USING DYNAMIC GRIDS FOR FINITE-DIFFERENCE SCHEMES

In this paper, a complete simulation of a trombone using finite-difference time-domain (FDTD) methods is proposed. In particular, we propose the use of a novel method to dynamically vary the number of grid points associated to the FDTD method, to simulate the fact that the physical dimension of the trombone’s resonator dynamically varies over time. We describe the different elements of the model and present the results of a real-time simulation


INTRODUCTION
The trombone is a musical instrument that presents distinct challenges from the perspective of physical modelling synthesis.In particular, the excitation mechanism between the lips and the player has been extensively studied, and simulated mostly using a simple mass-spring damper system [1].Because the majority of the bore is cylindrical, nonlinear effects can appear at high blowing pressures [2], leading to changes in timbre, or brassiness; such effects have been investigated and simulated [1,3,4].However, the defining characteristic of the trombone is that the physical dimensions of the resonator vary during playing.Synthesis techniques such as digital waveguides allow an approach to dynamic resonator changes in a simple and computationally efficient way, simply by varying the length of the corresponding delay line.This feature has been used in real-time sound synthesis [5], for simplified bore profiles suitable for modelling in terms of travelling waves.
However, when attempting more fine-grained modelling of the trombone resonator using finite-difference time-domain (FDTD) methods, the issue of the change in the tube length is not trivial.Previous implementations of brass instruments using these methods focus on the trumpet [6] and various brass instruments (including the trombone bore) under static conditions [7].To our knowledge, the simulation of a trombone varying the shape of the resonator in real time using FDTD methods has not been approached.We can tackle this problem by having a grid that dynamically changes while the simulation is running as presented in a companion paper [8].Briefly described, we modify the grid configurations of the FDTD method by adding and removing grid points based on parameters describing the system.
In this paper, we propose a full simulation of a trombone, describing all its elements in detail with a specific focus on the dynamic grid simulation.Section 2 presents the models for the tube and lip reed interaction in continuous time.Section 3 briefly introduces FDTD methods and the discretisation of the aforementioned continuous equations.Section 4 presents the dynamic grid used to simulate the trombone slide and details on the implementation are provided in Section 5. Section 6 presents simulation results, and some concluding remarks appear in Section 7.

CONTINUOUS SYSTEM
Wave propagation in an acoustic tube can be approximated using a 1-dimensional (1D) model, for wavelengths that are long relative to the largest lateral dimension of the tube.Consider a tube of time-varying length L = L(t) (in m) defined over spatial domain x ∈ [0, L] and time t ≥ 0. Using operators ∂t and ∂x denoting partial derivatives with respect to time t and spatial coordinate x, respectively, a system of first-order partial differential equations (PDEs) describing the wave propagation in an acoustic tube can then be written as: with acoustic pressure p = p(x, t) (in N/m 2 ), particle velocity v = v(x, t) (in m/s) and (circular) cross-sectional area S(x) (in m 2 ).Furthermore, ρ0 is the density of air (in kg/m 3 ) and c is the speed of sound in air (in m/s).System (1) can be condensed into a second-order equation in p alone, often referred to as Webster's equation [9].For simplicity, effects of viscothermal losses have been neglected in (1).For a full time domain model of such effects in an acoustic tube, see, e.g.[10].System (1) requires two boundary conditions, one at either end of the domain.The left boundary condition, at x = 0, will be set according to an excitation model to be described in Section 2.1.The right boundary, at x = L, is set according to a radiation condition.The radiation model used here, is the one for the unflanged cylindrical pipe proposed by Levine and Schwinger in [11] and discretised by Silva et al. in [12].As this model is not important for the contribution of this work it will not be detailed here in full.The interested reader is instead referred to [7,13] for a comprehensive explanation.

Coupling to a Lip Reed
To excite the system, a lip reed can be modelled as a mass-springdamper system including two nonlinearities due to flow, and the collision of the lip against the mouthpiece.In the following, y can be seen as the moving upper lip where the lower lip is left static and rigid.A diagram of the full lip-reed model is shown in Figure 1.Using dots to indicate time-derivatives, the lip reed is modelled as with displacement from the equilibrium y = y(t), lip mass Mr (in kg), externally supplied (angular) frequency of oscillation ωr = ωr(t) = Kr/Mr (in rad/s) and stiffness Kr = Kr(t) (in N/m).We extend the existing models of lip reeds [1] by introducing a nonlinear collision between the lips based on potential quadratisation proposed by [14].The collision potential is defined as with collision stiffness Kc > 0 and dimensionless nonlinear collision coefficient αc ≥ 1, The inverted distance between the lips η = η(t) ≜ −y − H0 (in m), for static equilibrium separation H0 (in m).
[η]+ = 0.5(η + |η|) indicates the "positive part of η".Notice, that if η ≥ 0, the lips are closed and the collision potential will be non-zero.This quadratic form of a collision potential allows for a non-iterative implementation [14].This will be explained further in Section 3.
Finally, Sr (in m 2 ) is the effective surface area and is the difference between the externally supplied pressure in the mouth Pm = Pm(t) and the pressure in the mouthpiece p(0, t) (all in Pa).This pressure difference causes a volume flow velocity following the Bernoulli equation (in m 3 /s) with effective lip-reed width wr (m).Another volume flow is generated by the lip reed itself according to (in m 3 /s).Assuming that the volume flow velocity is conserved, the total air volume entering the system is defined as This condition serves as a boundary condition at x = 0 for system (1).

DISCRETISATION
The continuous system described in the previous section is discretised using FDTD methods, through an approximation over a grid in space and time.Before presenting this discretisation, we briefly summarize the operation of FDTD methods.

Numerical Methods
Consider a 1D system of (static) length L described by state variable u = u(x, t) with spatial domain x ∈ [0, L] and time t ≥ 0.
The spatial domain can be disctretised according to x = lh with spatial index l ∈ {0, . . ., N }, number of intervals between the grid points N , grid spacing h (in m) and time as t = nk with tem- poral index n ∈ Z 0+ and time step k (in s).The grid function u n l represents an approximation to u(x, t) at x = lh and t = nk.Shift operators can then be applied to grid function u n l .Temporal and spatial shift operators are from which more complex operators can be derived.First-order derivatives can be approximated using forward, backward and centred difference operators in time (all approximating ∂t) and space (all approximating ∂x) where 1 is the identity operator.Furthermore, forward, backward and centred averaging operators can be defined in time and space Finally, an approximation δtt to a second time derivative may be defined as
In addition, a discrete cross-sectional area S l ≈ S(x = lh) with l ∈ {0, . . ., N } is assumed known.System (1) can then be discretised as Sl where S l+1/2 = µx+S l and Sl = µx−S l+1/2 are approximations to the continuous cross-sectional area S(x).The values for Sl at the boundaries, i.e., S0 and SN , are set equal to S(0) and S(L).
Expanding the operators, we obtain the following recursion where λ = ck/h is referred to as the Courant number and in order for the scheme to be stable [15].In implementation, the following steps are taken to calculate λ: where ⌊•⌋ denotes the flooring operation and is necessary because N is an integer.This causes (16) to not be satisfied with equality for all choices of L.
Equations ( 15a) and (15b) hold for l ∈ {0, . . ., N } and l ∈ {0, . . ., N − 1} respectively, and thus, in analogy with the continuous case, two numerical boundary conditions are required in order to update p n+1 0 and p n+1 N .These are provided by numerical equivalents of the excitation condition (see Section 3.3 below) and the radiation condition (in [13]).

Lip reed
As the lip reed interacts with the particle velocity of the tube via Eq.( 7), it is discretised to the interleaved temporal grid, but to the regular spatial grid as it interacts with the boundary at x = 0. Equations ( 2) -( 7) are then discretised as follows: Here, following [14], where κ = 1 if ψ n ≥ 0, otherwise κ = −1.It should be noted that condition (19c) has been added to the definition of g from [14] to prevent a division by 0 in (19b).Finally, η ⋆ = −y ⋆ − H0 where y ⋆ is the value of y n+3/2 calculated using system (18) (after expansion) without the collision potential.This means that system (18) needs to be calculated twice every iteration, once without the collision term and once with.The process of calculating the pressure difference ∆p n+1/2 in (18) will not be given here, but the interested reader is referred to [13, Ch. 5] for a derivation.

DYNAMIC GRID
The defining feature of the trombone is its slide that alters the length of the tube, changing the resonant frequencies.In a companion paper [8], we present a method to dynamically change grid configurations of FD schemes by inserting and deleting grid points based on an instantaneous value of the time-varying wave speed c(t).Although here, the tube length L(t) is varied, the method still applies.Note that this method only works for slow (sub-audio rate) parameter changes.
We can split a tube with time-varying length L n into two smaller sections with lengths L n p and L n q (in m) such that L n = L n p + L n q .Splitting the schemes in (14) in this way yields two sets of firstorder systems.The pressure and particle velocity of the first (left) system p n lp and v n+1/2 lp+1/2 are both defined over discrete domain lp ∈ {0, . . ., M n }, and those of the second (right) system q n lq and w n+1/2 lq −1/2 are defined over discrete domain lq ∈ {0, . . ., M n q }, with M n = ⌈L n p /h⌉, and where ⌈•⌉ denotes the ceiling operation.Note, that the domains for v and w have an extra grid point when compared to the regular case in (14) and that w is indexed with lq − 1/2 rather than lq + 1/2.
The resulting system of FD schemes then becomes Sl Sl Here, due to the different indexing for w, the spatial derivatives for the right system are flipped (δx+ became δx− and vice versa).Also note, that l is still used for the spatial indices of S and S which now approximate S(x) according to M n and q n+1 0 are shown.The most important difference with the usual case is that the virtual grid points p n M n +1 and q n −1 are the result of the interpolation of known pressure values at n using Eq.(27).
The conditions for the outer boundaries of this system, i.e., at lp = 0 and lq = M n q , are the same as for the full system.The inner boundaries, lp = M n and lq = 0 are connected according to the method described in [8] to be explained shortly.To be able to calculate p n+1 M n and q n+1 0 , the domains of v and w have been extended at the inner boundaries to include v n+1/2 M n +1/2 and w n+1/2 −1/2 .These, however, require points outside of the domains of p n lp and q n lq , i.e., p n M n +1 and q n −1 .In [8] we propose to calculate these virtual grid points based on known values of the system.Despite the fact that [8] presents the method using a second-order system, it can still be applied here.The process of how p n+1 M n and q n+1 0 are calculated is visualised in in Figure 2. Notice that all time steps use the same value of M n and M n q .In other words, the expansion of the temporal operators in (9) do not affect the temporal indices n in M n and M n q .

Changing the Tube Length
In the following, the location of a grid point u l along the grid (in m from the left boundary) at time index n is denoted as x n u l .The two pairs of first order systems in (21) are placed on the same domain x with x n p lp = lph, and describing the locations of the left system and right system respectively.Here, it can be observed that as the tube length L n changes, the locations of the grid points of the right system will change.More specifically, as the trombone-slide is extended and L n increases, all grid points of the right system move to the right, and to the left for a contracting slide.If L n is changed in a smooth fashion, the continuous domain x ∈ [0, L n ] will not necessarily be subdivided into an integer amount of intervals N n (of size h = ck).This is where a fractional number of intervals is introduced and is defined as which is essentially the calculation of N in Eq. ( 17) without the flooring operation, and N n = ⌊N n ⌋.The fractional part of N n can then be calculated using which describes the distance between the inner boundaries along the grid in terms of how many times h would fit in-between (which is always less than once).If N n = N n and α = 0, the inner boundary locations perfectly overlap, and x n p M n = x n q 0 .This also means that the domain x can be exactly divided into N n equal intervals of size h = ck.As the virtual grid points p n M n +1 and q n −1 perfectly overlap with q n 1 and p n M n −1 respectively, these values can be used directly to calculate the grid points at the inner boundaries.This situation effectively acts as a rigid connection between the grid points at the inner boundaries defined as If α ̸ = 0, some other definition for p n M n +1 and q n −1 needs to be found.We use quadratic Lagrangian interpolation according to which can then be used to calculate v n+1/2 M n +1/2 and w n+1/2 −1/2 and consequently p n+1 M n and q n+1 0 (see Figure 2).This process is repeated every sample.It can be shown through the rigid connection in (26), that if α = 0, the definitions in (27) reduce to p n M n +1 = q n 1 and q n −1 = p n M n −1 as stated before.

Adding and removing grid points
As the tube length L n changes, L n p and L n q also change according to where DAFx.4 Proceedings of the 24 th International Conference on Digital Audio Effects (DAFx20in21), Vienna, Austria, September 8-10, 2021 which causes the number of intervals between grid points M n and M n q to change as well, according to Eq. ( 20).
The following state vectors are introduced for the pressure, defined for n + 1 and n and for the velocity, defined for n + 1/2 and n − 1/2 and contain the different states over the discrete domains defined at the beginning of this section.Here, T denotes the transpose operation.
If N n > N n−1 , points are added to the left and right system in an alternating fashion: where and cubic Lagrangian interpolator Here, adds an offset to half of the elements in the z vectors depending on the difference between v n−1/2 M n +1/2 and w n−1/2 −1/2 .Why this is necessary will be further explained in Section 4.3.Finally, I ← 3 and η ← are flipped versions of (34) and (35) respectively.
If N n < N n−1 , points are simply removed from the vectors according to Notice that the even and odd conditions in Eqs. ( 32) and (36) can be swapped.To stay as close to the desired location of adding and removing grid points as possible, this requires the ceiling and flooring operations in (20) to be swapped as well.

Drift of w
The inner boundaries of the pressure states p and q are connected by ( 27), but no such connection exists for the velocity states v and w.As the radiating boundary is implemented on the pressure grid, this leaves w without any boundary condition; it is only "held in place" by the pressure values of q, or more specifically, by derivatives (both spatial and temporal).As FD schemes are an approximation, it does not give a perfect solution and w tends to 'drift' during the simulation, especially when L n is changed.
Luckily, as the pressure values are also calculated from derivatives of the velocity, the absolute state of w does not matter.The difference in values at the connection point is also irrelevant as there is no spatial derivative taken between v and w (refer to Figure 2).Finally, the pressure values are used for the output audio of the simulation, so the drift does not affect the audio.
The absolute states of the velocity vectors do, however, need to be accounted for when adding points to the v and w using (32).The current drift can be approximated by observing the difference between ) when a grid point is added.This is then used in a drift-correction vector η n−1/2 presented in (35).When a point is added to v, the values of w in zv are offset by the aforementioned difference and when a point is added to w the same happens (inverted) for the values of v in zw.This way, the drift is allowed, but does not affect the state of the newly added grid points.Notice that the drift does not affect the operations of point removal in (36).

State Correction
As L n , and consequently the number of grid points, is decreased, it might occur that the grid points at the inner boundaries p n M n and q n 0 have a very different value when α ≳ 0, i.e., right before a point is removed.This violates the rigid connection in Eq. ( 26).
We propose in [8] to add an artificial spring-like connection between the grid points at the inner boundaries that "corrects" the state of these points.Applying this to system (21) extends Eqs.(21a) and (21c) according to Sl where the spreading operators are defined as Furthermore, the correction effect is defined as with spring damping σsc, pressure difference and scaling coefficient Table 1: Geometry of a measured trombone taken from [16].Numbers correspond to Figure 3. Here, ε ≪ 1 to prevent division by 0. Just like in [8], the implementation of the correction effect allows for an infinite β when α = ε = 0 acting like a rigid connection between Eqs. (37a) and (37b).

IMPLEMENTATION
The implementation has been done in C++ using the JUCE framework 1 , and is available online 2 as well as a demo showcasing it. 3 The audio output of the system can be retrieved by selecting a grid point on the pressure grid and listening to this at the given sample rate fs.Here, the radiating boundary q n M n q is chosen, as this is where the sound enters the listening space in the real world.To mimic low-pass filtering happening due to a distributed radiating area, a 4 th -order low-passing Butterworth filter with a cutoff frequency of fc = c 2 π/S(L) ≈ 3245 Hz is used.This equation is retrieved by choosing the listening point to be at the bell surface 1 https://juce.com/ 2https://github.com/SilvinWillemsen/cppBrass/releases/ 3https://youtu.be/Ht5gVNrshYoNumbers correspond to the parts of the tube found in Table 1 and dashed lines highlight where the different parts are separated.The tube is split in the middle of the slide crook with the colours corresponding to those in Figure 2.
and integrating over the bell area.

Parameters
For the most part, the parameters used in the simulation have been obtained from [13,16,17].The lengths and radii of different parts of the tube can be found in Table 1 and a diagram showing this geometry is shown in Figure 3.The system is split in the middle of the slide crook such that the ranges for the lengths of the two tubes are L n p ∈ [0.797, 1.327] and L n q ∈ [1.796, 2.326].Other parameters used in the simulation can be found in Table 2.Not included here is λ, which has been set slightly lower than the stability condition in (16), i.e., λ = 0.999.Although the implementation works when λ = 1, this is done to tolerate (much) higher speeds of change in L n before instability occurs (see Section 5.2).Not satisfying condition (16) causes bandlimiting and dispersive effects [15], but such a small deviation from the condition has no perceptual influence on the output sound and outweighs the problems caused by instability.
As the tube acts mainly as an amplifier for specific resonant frequencies it is important to match the frequency of the lip reed to a resonating mode of the tube.This frequency depends on L n in the following way where L n+1/2 = L n and scalar multiplier F = 2.4 was heuristically found to best match the 4 th resonating mode of the tube and generates a recognisable brass sound.

Limit on speed of change
To reduce audible artifacts and instability issues from adding and removing points, and to stay in the sub-audio rate regime, a limit can be placed on (29) as where Nmaxdiff is the maximum change in N per sample and has been set to Nmaxdiff = 1/20.This means that a grid point can be added or removed every 20 samples and allows the entire range of L to be traversed in ca.0.06 s at a sample rate of fs = 44100 Hz.

DAFx.6
Figure 4: Screenshot of the graphical user interface (GUI).The geometry (in orange) as well as the states of the pressure (in blue) and velocity scaled by S (in green) are shown.For clarity, the start and end of the outer slide are denoted by dashed lines.The drift of w as explained in Section 4.3 is visible from the "kink" in the green line exactly in the middle of the outer slide.

State correction
The introduction of system states at n + 1 through the centred operators in Eq. (39) seem to make the scheme implicit.It is, however, possible to calculate Fsc explicitly [15,18].The same operators also introduce the need for values at n − 1, i.e., p n−1 M n and q n−1 0 .Therefore, the vectors p n−1 and q n−1 will need to be stored, and the operations to add and remove grid points as described in 4.2 need to be applied to these as well.One could argue that only two points at the inner boundaries are needed for the calculation and to create r in (33) at n − 1.For generality, we continue with the entire vectors defined over the same domains as p n and q n respectively.

Graphical User Interface and Control Mapping
A screenshot of the graphical user interface (GUI) is shown in Figure 4.The geometry of the tube is plotted along with paths showing the pressure states in blue and the velocity (scaled by the geometry S) in green.The audio thread of the application runs at 44100 Hz whereas the graphics are updated at a rate of 15 Hz.
The real-time application is controlled by interacting with the bottom panel using the mouse.The x-axis is mapped to tubelength L n and also modifies the lip-reed frequency ωr according to Eq. ( 42).The y-axis changes the multiplier F in Eq. ( 42) and the black line in the vertical middle of the control panel is mapped to F = 2.4.The pressure is modulated by a slider at the bottom of the control panel.As of now, no focus has been put on intuitive parameter mapping; it has only been implemented for simple parameter exploration.

Order of Calculation
Algorithm 1 shows the order in which the different parts of the system presented in this paper are calculated.the realism of the simulation sound could be to add viscothermal losses [19] or nonlinear effects [3].Furthermore, for lower values of the lip frequency ωr, the sound exhibits extra oscillatory behaviour making the output "non-smooth".This might be due a higher average displacement of y for lower ωr and the nonlinear collision present in the lip model will have a greater effect on its displacement.Variable collision stiffness might solve this issue but is left for future work.
Informal listening by the authors shows that the method used to implement the dynamic grid does not introduce perceivable audible artifacts, even when L n is changed very rapidly.Naturally, this needs to be confirmed by formal listening tests.Despite the limit placed on the speed of change of L n in (43) the control of the application does not exhibit a noticeable delay and changes in L n feel immediate.
The main difference between the method in [8] and the version used here, is that the method is applied to a system of first-order equations rather than the second-order 1D wave equation.Because the connection between the inner boundaries is only applied to the grid functions describing pressure, a drift occurs in w as it is left without boundary conditions.Although this drift does not have an effect on the output sound, as discussed in Section 4.3, too high or low values might cause rounding errors in the simulation.As it is expected that this only happens at extremely high or low values after a long simulation length, the drift is not considered an issue at this point.

CONCLUSION
In this paper, we have presented a full implementation of the trombone including a lip reed, radiation and a tube, discretised using FDTD methods on a dynamic grid.Informal evaluation by the authors shows that the implementation exhibits no audible artifacts when grid points are added and removed, even under relatively fast variation in tube length.Naturally, this needs to be confirmed by formal listening tests.Moreover, the simulation easily runs in real-time allowing it to be used as an audio plugin.
Future work will include extending the tube model to include more realistic viscothermal and nonlinear effects and variable collision stiffness in the lip model.Furthermore, the investigation of more intuitive control parameter mappings is a necessary step towards a real-time instrument.

Figure 1 :
Figure 1: Diagram of the lip-reed system with the equilibrium at 0 and the distance from the lower lip H0.The various symbols relate to those used in Eq. (2).

) DAFx. 3 Figure 2 :
Figure 2: Schematic showing data flow of how different grid points at time index n + 1 are calculated with α = 0.25 in Eq. (25).To prevent cluttering, arrows going straight up (indicating that the state of a grid point at time step n is needed to calculate the state of that grid point at n + 1) are suppressed.As an example of the usual case (refer to Eq. (15)), the points required to calculate p n+1 2 are shown.Furthermore, the points needed to calculate p n+1 M n and q n+1

Figure 3 :
Figure 3: Diagram showing the trombone geometry (not to scale).Numbers correspond to the parts of the tube found in Table1and dashed lines highlight where the different parts are separated.The tube is split in the middle of the slide crook with the colours corresponding to those in Figure2.