Real-time implementation of a friction drum inspired instrument using finite difference schemes

Physical modelling sound synthesis is a powerful method for constructing virtual instruments aiming to mimic the sound of real-world counterparts, while allowing for the possibility of engaging with these instruments in ways which may be impossible in person. Such a case is explored in this paper: particularly the simulation of a friction drum inspired instrument. It is an instrument played by causing the membrane of a drum head to vibrate via friction. This involves rubbing the membrane via a stick or a cord attached to its center, with the induced vibrations being transferred to the air inside a sound box. This paper describes the development of a real-time audio application which models such an instrument as a bowed membrane connected to an acoustic tube. This is done by means of a numerical simulation using finite-difference time-domain (FDTD) meth-ods in which the excitation, whose position is free to change in real-time, is modelled by a highly non-linear elasto-plastic friction model. Additionally, the virtual instrument allows for dynamically modifying physical parameters of the model, thereby allowing the user to generate new and interesting sounds that go beyond a real-world friction drum.


INTRODUCTION
The friction drum has been described as a peculiar musical instrument or even a noisy toy [1].In Scandinavian tradition, the friction drum was used in the Middle Ages as a rhythmic instrument.Over time, children used it in the 19 th century to play when they went door to door during the Christmas holidays and sang [2]. Figure 1 shows the friction drum present at the Danish Music Museum.As can be seen, the drum is a combination of a stick inserted in the middle of a cylindrical drum.Similar drums are found in other cultures, for example the cuica in Brazil, putipù in Italy, zambomba in Spain, buhai in Romania and variations on rummelpot in Germanic countries like Denmark, Germany or the Netherlands.The sound is produced by rubbing the stick connected to the membrane, generating a frictional excitation, hence the name friction drum.This excitation is transferred to the membrane and then further to the air inside the acoustic tube, with the output radiating at the open end of the tube.As the stick is of limited length the sound produced by these means cannot be sustained indefinitely and these instruments are typically used for percussive sounds, characterised by a  sharp attack.The pitch produced by these instruments can be also modulated by applying pressure to the drum head and therefore changing its tension.
For our virtual model, we envision the excitation mechanism somewhat differently: we consider the stick as an infinitely long bow that excites the membrane directly via frictional interaction.This allows for an output sound sustained indefinitely.Therefore our virtual instrument does not necessarily sound identical straight "out of the box" compared to either of the real friction drums in particular, but could hopefully be tuned to achieve a particular sound.In addition, as the virtual instrument is based on a physical model, there comes a learning curve for any new user with regards to learning how to play it, as there is for any real instrument.
Friction has been extensively investigated in the sound synthesis literature, being the sound excitation mechanism of several musical instruments such as the violin and the musical saw, but also everyday sounds such as squeaking doors and rubbed wine glasses [3].The literature has examined different ways of simulating friction, and elasto-plastic friction models have proven to be an accurate way to simulate dry interactions between rubbed surfaces [4].In the elasto-plastic friction model, friction does not depend only on the relative velocity between bodies in contact, but also on the relative displacement of the micro-interactions between such DAFx.1 bodies [5].Such models have been recently used in combination with finite difference schemes [6].In this paper, we combine an elasto-plastic friction model together with a membrane simulation based on FDTD methods to model the interaction between a stick, acting as a bow, and the drum head of a friction drum inspired virtual instrument.This is furthermore coupled with a 1D wave model with radiating boundary conditions at the open end used to describe the drum's sound box.
We mathematically describe the different elements of the drum, and present a real-time implementation which uses the Sensel Morph, a pressure sensitive touch pad controller [7], in order to play the virtual instrument.

FRICTION DRUM MODEL
We propose to model the friction drum as two main components and an excitation mechanism: a membrane connected to an acoustic tube, where the membrane can be bowed via a non-linear elastoplastic friction model.This section describes the partial differential equations (PDEs) describing these components in isolation.

Membrane
First off, for a membrane defined over a domain Dm = [0, Lx] × [0, Ly], with Lx and Ly being the lengths of the membrane [m] in the Cartesian coordinates (x, y), the transverse displacement at a time t [s]: u(x, y, t) [m] can be described by the following PDE, as per [8]: where the 2D Laplacian operator is defined as: and the parameter c = T /ρmH is a measure of wave speed resulting from the membrane's tension per meter T [N/m], its density ρm [kg/m 3 ] and its thickness H [m]. Furthermore, σ0 [s −1 ] and σ1 [m 2 /s] are parameters controlling the frequency-dependent and frequency-independent loss respectively.
Dirichlet boundary conditions are considered for the membrane, as it is a reasonable assumption that no energy is lost at the boundaries.Thus:

Acoustic Tube
The longitudinal vibration of an air column ζ(χ, t) [m] in a tube of uniform cross-section and length Lχ can be described by the following equation: χ ∈ [0, Lχ] is a spatial coordinate along the length of the tube.This 1D wave approximation for a tube holds true if the length scale in the longitudinal direction is significantly greater than in the others.For the case of the friction drum this is not true and this approximation will actually not produce the desired "drum" type sound for the wave speed γ resulting from the bulk modulus and density of air.However, this value can be tuned to produce the desired sound.The choice for this value is given in Section 4. A more accurate model would be the 3D wave equation or propagating 2D cross-sections but would be too computationally demanding for a real-time implementation.
A Neumann boundary condition is imposed at the side of the tube connected to the membrane, while at the open end of the tube a radiating boundary condition is chosen: with the constants α1 and α2 modelling the inertia and loss at the open end of the tube.

Connection
A connection between the two components is considered, which means that a connection force is added to the PDEs describing the two separate components in the following way: where fc,m [N] is the connection force acting on the membrane, while fc,t [N] acts on the tube.A [m 2 ] is the area of the tube and the terms Em [m −2 ] and Et [m −1 ] represent some distributions over which the connection force is applied on.Notice that the distributions are given over the appropriate domain, i.e. over a surface (m 2 ) for the membrane and over a length (m) for the tube.This essentially means that the connection forces applied to the zero-input PDEs given in Equations ( 1) and ( 4) are scaled with the mass per unit length or area of each considered component.For the current model, a rigid connection is assumed meaning that: with ⟨•, •⟩ being an L 2 inner product over the appropriate domain (1D for the tube and 2D for the membrane).Thus η represents the relative displacement of the two components over the connection.A more detailed discussion on the choice of connection distributions is given in the Section 3.

Excitation -Bowing Model
An elasto-plastic friction model, first proposed in [9], is used for the friction drum excitation mechanism.Such a model, applied to the bowing of a stiff string, was shown by both Serafin et.al [10] and Willemsen et.al [6] to produce a hysteresis loop in the bowing force versus relative velocity plane, a detail which was experimentally observed by Smith and Woodhouse in [11].This was arrived at by both a digital waveguide model, implemented by [10] as well as a FDTD method used by [6], who presented a working real-time implementation of a bowed stiff string.
The elasto-plastic friction model assumes the contact between two interacting elements is highly irregular at the microscopic level, i.e. not all the overlapping surface is actually in contact.Instead, the contact can be modelled via a large group of bristles each contributing to the total friction force.These bristles are modelled as damped stiff springs and therefore each generates increasing contact force with increasing displacement, describing an elastic DAFx.2 regime (stick).However, each bristle can only displace so far before it "breaks" and not all bristles "break" at the same time.This represents the elasto-plastic regime of the friction model, where only some bristles have reached the breaking point (slide).Once all bristles "break", a completely plastic regime is entered (slip).In the case of the bowed membrane, after the slip a "new" portion of the bow gets in contact with the membrane and the stick-slide-slip cycle restarts.
Adding the bowing force to the membrane can be done by introducing an extra term to (6a), the PDE which governs the membrane connected to the tube: Again the force is applied over some distribution which is in fact a single point on the Cartesian grid of the membrane given by the bowing position at some time where v can be computed as the difference between the velocity of the membrane at the bowing location and the externally supplied velocity of the bow, vB(t) [m/s]: Furthermore, s0 is the bristle stiffness [N/m], s1 is the damping coefficient of the bristles [kg/s] and s2 is the viscous friction [kg/s].s3 [N] is a force coefficient proportional to the normal bowing force fN (t) (which is an external input and can vary over time) scaled with a pseudorandom function w(t) ∈ [−1, 1] and is used to add noise to the total bowing force, as per [3].The time derivative of the average bristle displacement, ż [m/s] is given by: Perhaps the most important function in this elasto-plastic friction model is introduced above: the adhesion map α(z, v) which controls the transition between the various regimes of friction.The function is defined as follows and is illustrated in Figure 2a: ) where z ba [m] is the breakaway bristle displacement, below which the friction regime is purely elastic.Indeed, when α = 0, it follows that ż = v.Then, when z ba is surpassed, the elasto-plastic regime is entered where the value of z will be a proportion, governed by αm(v, z), of the steady-state bristle displacement zss(v), with: At steady-state, when slipping occurs and therefore ż = 0, α = 1 and together with Equation ( 11) it follows that z = zss(v), with zss(v) [m] being defined as: where the Coulomb force FC = µC fN [N] and stiction force FS = µSfN [N] are given as a proportion of the normal force, fN (t) of the bow acting on the membrane.These proportions are controlled by the dimensionless dynamic and static friction coefficients respectively, µC and µS.Looking at Figure 2b which illustrates zss(v) for a fixed fN one can see that the Stribeck effect, i.e. the dip of force at low velocities is captured, as the magnitude of zss at values of v close to zero is larger than for higher relative velocities.This allows for a larger total friction force to be obtained in this region before the plastic regime is reached.After slipping occurs, the "grip" of the bow on the membrane is briefly lost and the membrane displaces in the opposite direction, hence sgn(v) ̸ = sgn(z) and α becomes again zero, meaning that the bow again sticks to the membrane and a new stick-slip cycle begins.

Complete System
The complete system for the friction drum can be therefore written in continuous time as:

DISCRETIZATION
The system given in ( 15) is discretized using FDTD methods, which subdivides the continuous model into grid points in space Using these discrete definitions for space and time, the continuous state variables presented in the previous section can then be approximated by grid functions as u(x, y, t) ≈ u n l,m for the membrane and ζ(χ, t) ≈ ζ n p for the tube.Furthermore, approximations to the derivatives can be described in the following way: Note that the same continuous operation can be approximated in different ways.With these definitions in place we can move on to discretize the individual components of the friction drum model.An important thing to take into account when it comes to numerical models is the issue of stability, from which limitations arise on the possible size of the grid spacing hm and ht.Stability conditions are available for each individual component and will be presented in the upcoming subsections.Working with grid spacings that satisfy the stability conditions as close to equality as possible ensures a more accurate numerical scheme.

Membrane
The complete membrane including the bowing force and connection to the tube can be discretized as Notice in Equation ( 17) the use of the δt− operator in the mixed time/space derivative term, which is used in order to keep the numerical scheme explicit.
otherwise, (18) with lB = floor(xB/hm), mB = floor(yB/hm), αx B = xB/hm − lB and αy B = yB/hm − mB.This spreading function, necessary for exciting a discretized grid has a dual: the interpolation function IB which is of interest when obtaining the state of a discrete grid between grid points and will be of use further along for the discretization of the complete system.See [8] for more details on this.They are termed dual functions because they serve inverse purposes: one is used for adding input to a distributed object at a specific location and the other is used to extract a state at a specific location on the same object.Furthermore, they are related as such: Drum heads are typically of circular shape and although the membrane is defined in Cartesian coordinates over some rectangle of length Lx × Ly, one can "sculpt" a circular grid using a staircase approximation, as done in [12], as long as boundary conditions are satisfied.Since Dirichlet conditions are assumed, the only requirement is that points on the rows and columns at the edge of the square grid need to be fixed to zero.Regarding the connection with the tube, it is clear that the entire membrane contributes to movement of the air column inside the tube.However, there is a factor which points towards skewing the weight of the membrane displacements towards its center.Due to the boundary layer effect, the air at the edges of the tube will be semi-stationary.Therefore a 2D Hann distribution over 72.25% of the area of the grid is used, centered at the middle of the membrane.This is illustrated in Figure 3 together with the grid points of the circular membrane approximated from the initial rectangular grid.
This connection distribution, named Im, is normalized such that its integral is equal to 1 and can be seen to act as an interpolation function acting on u n l,m .Therefore it's dual spreading function Jm will be defined in the same way as Equation ( 19).
Using von Neumann analysis [8], a stability condition can be derived and is given by the following inequality: (20)

Acoustic Tube
A discretized version of the acoustic tube and its connection is given by: with Jt being the spreading operator for the connection force acting on the tube, essentially the discretized version of Et.This is DAFx.4 related to its dual interpolant function It in the following way: with It taken as a normalized half-Hann window such that its integral is 1, spread over 4% of the length of the tube, with its peak at the connection point (the top of the acoustic tube).This was preferred due to a Dirac type connection in order to dampen out some of the high frequencies which would result from an excitation of the tube at a single point and thus produce a more realistic friction drum sound.Notice in Equation ( 22) that ht is not squared, as was the case for hm in Equation ( 19).This is due to the different spatial dimensions of the components.The boundary conditions of the tube presented in Equation 5 are discretized in the following way: and a stability condition on the grid size ht is given by [8]: ht ≥ ht,min = γk. (24)

Connection
The rigid connection given in Equation (7b) can be discretized as with Dm and Dt being the domains of the membrane and of the tube respectively.This means that the connection location is described by the two spreading functions Jm and Jt for the membrane and the tube respectively.If the equality in Equation ( 25) is true at sample n then it follows that it will be true at sample n + 1 as well.This together with with the following identity will provide valuable information for solving the complete discretized system: where f is a grid function in some domain D and I and J must be dual interpolation and spreading functions.This results in the following equality:

Excitation -Bowing Model
For the bowing force, the discrete counterpart of Equation ( 9) is taken as: Additionally, the relative velocity between the bow and the membrane described in Equation ( 10) will be: (30)

Solving the System
In order to calculate the update values for the grid functions: u n+1 l,m and ζ n+1 p , three unknown variables must first be determined: v n , z n and f n c and for this we need a system of three equations dependent on these variables at each sample n, which can then be solved using a multivariate Newton-Raphson method.An interesting observation is that the reaction of the air inside the acoustic tube, i.e. the connection force, will instantaneously affect the bowing force and vice-versa.This would not be the case in a simpler model where bowing would not occur at the connection point.
The first equation is g1(v n , z n , f n c ) = 0 and can be found by making use of the following identity: and introducing it together with Equation (30) in Equation ( 17), which results in: with The second equation needed is, as per [6] and [8]: where the operators applied to z n describe the trapezoid rule.Finally, the third equation comes from the rigid connection condition in Equation ( 27).The displacements of the membrane u n+1 l,m and for the tube ζ n+1 p can be extracted and expressed only in terms of values at current or previous samples by expanding the operators in Equations ( 17) and ( 21).This results in DAFx.5 with The following iteration is then used to calculate the unknown values v n , z n and f n c : where i is the iteration number.The threshold for convergence is set at 10 −7 , with a maximum number of iterations of 99.
Once the three values at at the sample n are known, update values for the grid points u n+1 l,m and ζ n+1 p can be found by expanding the operators in Equation (17) and Equation (21).

IMPLEMENTATION
The implementation of the finite difference scheme presented in Section 3 has been carried out in C++ using the JUCE framework [13] and a demonstration video is available at [14].The parameters used can be found in Table 1, and have been chosen starting from the work of Serafin [3] and Willemsen et al. [6], but tuning them to achieve the desired sound for the instrument.The number of grid intervals in the discretization of both elements is limited as compared to the stability condition in order for the model to be able to run in real-time without audio dropouts.

Prototype Model Results
Before implementing the audio application in JUCE, tests were done in MATLAB [15] to identify that the model is stable and that results are in line with expectations.Figure 4 shows a snapshot of the circular membrane being bowed with fN = 12 [N] and vB = 0.1 [m/s] coupled with the acoustic tube at some time step in the middle of a simulation.Looking at the tube, one can see the free and the radiating boundaries at its endpoints.
The next step was to test whether the vibrations of the membrane exhibit the Helmholtz motion, which is typical for bowed instruments and tends to produce triangular-shaped wave forms.Finally, the presence of a hysteresis loop in the force vs. relative velocity plane is investigated, which is an expected behavior as per experimental observations of bowed strings by Woodhouse and Smith [11].This is illustrated in Figure 5b.tice the orange square, highlighting the bowing position, which has different opacities in the two snapshots, their meaning being further described in the following.An important part in designing the real-time application was to have a natural type of interaction.Since there are 4 dimensions of input to the model, i.e., the bowing position (xB(t), yB(t)), bowing force fN (t) and velocity vB(t), it was desired to find a way to somehow control all these inputs simultaneously.An ideal match for this task was the Sensel Morph which is a tablet-sized pressure sensitive controller which is very fast and extremely sensitive [7].The work of [16] provided an open source library for allowing easy communication between the Sensel and JUCE.This Other parameters which can be modulated via sliders are the tuning of the membrane, i.e. the wave speed c ∈ [15, 150] [m/s], which is named "Tuning" thus allowing for a more intuitive understanding of the parameter by a non-technical user.Similarly the damping parameters σ0 ∈ [0, 6] [s −1 ] and σ1 ∈ [0, 0.00266] [m 2 /s] are combined into one value called "Damping".Note, that grid spacing in Eq. ( 20) is initialised using the highest values for c and σ1 so that the stability condition is not violated.Also note, that even when the damping parameters are set to 0 the radiation damping parameters for the tube α1 and α2 are fixed.Hence, even with zero damping, there will still be decay present.A third slider at the top of the graphical user interface (GUI) window controls the s3 ∈ [0fN , 0.04fN ] [N] term in the bowing force, and is called "Noise" as it adds some white noise to the friction force proportional to the normal force fN .

Real-Time Application
Furthermore, a vibrato effect is added where one can modulate via a sine wave the tuning of the membrane by a chosen frequency and with a chosen amount.This is introduced in the GUI as the sliders, named "Variation" which adds an oscillation between [0, 3] [m/s] to the wave speed c and "Rate", which controls oscillation frequency of the sine wave and is in the range [0, 10] [Hz].
All the ranges mentioned above are mapped in the GUI to be in a [0, 10] non-dimensional scale as to not confuse the user with different scales.
To add to the natural feel of the interaction another important addition is included in the GUI: the vibration of the membrane is plotted in real-time in a gray scale, inspired by [16], together with the bowing position, plotted with an orange color and an opacity given by the amount of pressure one applies to the Sensel.This is in fact, in the first author's view, one of the paramount features of this digital instrument, the fact that one can hear and see what the membrane is doing and change the bowing accordingly to find "sweet spots" for the sound.
The output sound is retrieved from the model by following the state of the open end of the tube (ζ n Nχ ) and amplified to the usual range of amplitudes [-1,1].Since the amplitudes of the model states are higher when using a lower value for c, an adjustable gain is used.

EVALUATION
The real-time simulation presented in the previous section was demoed by the first author during a Zoom session with 17 students enrolled in a physical modelling for sound synthesis class.
After a demo where the different parameters of the interface were explored, a qualitative interview and discussion took place.To the question regarding which instrument it was, the answers were varied.One student said the sounds were inspired by the Theremin, another mentioned a gong, hand drum, metallic drum, bowed bar, low frequency saw, a cymbal that is "contact-miced" and bowed or even "a chair being dragged across the floor".Particularly, one student mentioned that it sounded like one was "inside" the instrument or that it resembles the sounds a contact microphone might pick up.This was encouraging to hear as the sound is being picked up right at the end of the tube so in some way the listener is inside the drum itself.Overall, the answers gave some indications on how the sonorities of the physical model remind of a bowed inharmonic resonator and references to friction were abundant in the DAFx.students' responses.Even if it was not possible for the viewers to play with the interface, they found the use of the Sensel intuitive and the sound produced felt natural.This informal evaluation is obviously not ideal; the feeling of the quality of the instrument is better experienced from the viewpoint of the player.Nonetheless, it provided some indications for further development of the interaction and the GUI associated with the instrument.An important addition to the evaluation would be an objective comparison with the sound produced by various real friction drums.To this end, however, the virtual instrument would need to be tuned and played accordingly, as the excitation mechanisms are not entirely analogous.

CONCLUSION
In this paper, the development of a real-time audio implementation of a virtual friction drum inspired instrument using physical modelling has been presented.FDTD methods are used for simulating the friction drum as a bowed membrane connected to an acoustic tube.Furthermore, an advanced elasto-plastic friction model is used for the excitation mechanism which is shown to exhibit physically consistent behavior observed in experiments on real music instruments such as the presence of a hysteresis loop in the resulting bowing force-velocity plane.
Future work may involve the investigation of other possible mappings to the various parameters or the use of different type of controllers.Another important direction should be the optimization of the C++ implementation with the aim of reducing the grid size intervals in the numerical model and therefore work closer to the stability conditions.This would produce more accurate results with a broader bandwidth and allow for the application to be implemented on a micro-controller and developed as a stand-alone digital instrument.
Copyright: © 2021 Marius George Onofrei 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.

Figure 1 :
Figure 1: A picture of the friction drum present at the Danish music museum.
with the wave speed γ = B/ρt [m/s] resulting from the bulk modulus B [Pa] and density ρt[kg/m 3 ] of the air inside the tube.
t: (xB(t), yB(t)).This is achieved by the use of a 2D Dirac delta function δ(x − xB, y − yB) [m −2 ].Furthermore, the force needs to be scaled with the mass per unit area of the membrane of the component.As for the bowing force itself, fB [N], it is a function of the relative velocity between the bow and the membrane, v [m/s], and the average bristle displacement, z [m]:

Figure 2 :
Figure 2: (a) A plot of the adhesion map α(v, z) plotted against z when the signs of v and z are the same.(b) Steady-state bristle displacement zss(v) for a constant normal force fN .

DAFx. 3 and
samples in time.The (x, y)-plane of the membrane is discretised as x = lhm and y = mhm, with l ∈ [0, ..., Nx] and m ∈ [0, ..., Ny].Here, Nx = Lx/hm and Ny = Ly/hm are the horizontal and vertical number of grid intervals the membrane is divided in with grid spacing hm [m].For simplicity the same spacing is used in both directions.Similarly for the tube, χ = pht, where p ∈ [0, ..., Nχ] and Nχ = Lχ/ht is the total number of grid intervals along the tube's length with a grid spacing ht [m].Time t is discretized as t = nk where k = 1/fS, sampling frequency fS [Hz] and temporal index n ∈ N.

Furthermore
Jm and JB(xB, yB) are spreading operators.The former is a discretized version of the connection distribution Em and the latter, a discrete 2D Dirac delta function which defines the bowing position in the continuous model.Here we use a first order 2D spreading function defined as

Figure 3 :
Figure 3: Circular grid approximation from a rectangular grid and the normalized Hann distribution used for the connection to the tube, Im.The green crosses are the original grid points from the square Lx × Ly grid, while the red circles are the points used in the calculation.

Figure 4 :
Figure 4: Snapshot showing the displacements of the friction drum's components at a time step in the middle of a bowing simulation, (a) being the longitudinal displacements of the air column in the acoustic tube ζ and (b) being the transverse displacements of the membrane u.The magenta cross highlights the bowing position.

Figure
Figure5ashows the displacement of the membrane at the bowing location during a simulation, as well as the relative velocity and the resulting displacements at the open end of the tube.The membrane indeed shows Helmholtz motion, while the relative velocity exhibits the stick-slip behavior with values hovering around zero followed by an abrupt drop after which a new portion of the bow sticks again and the cycle restarts.The wave form of the displacement at the open end of the tube is somewhat more complex and highlights the effect of the interaction of the elements, with the entire membrane contributing to the motion of the air inside the tube, which then feeds back into the membrane.Finally, the presence of a hysteresis loop in the force vs. relative velocity plane is investigated, which is an expected behavior as per experimental observations of bowed strings by Woodhouse and Smith[11].This is illustrated in Figure5b.

Figure 6 Figure 5 :
Figure 6 shows snapshots of the friction drum audio application during use, where due to variation of the bowing position and force/velocity different modes of vibration are in resonance.No- allowed to map the (x, y) touch position to the bowing position while the pressure was mapped to the bowing force and velocity, linearly coupled.The normal force is limited in the range of fN ∈ [0, 20] [N] while the bowing velocity is mapped in the range vB ∈ [0, 0.2] [m/s].Naturally, the bowing position is limited in the range of [0, Lx] and [0, Ly].

Figure 6 :
Figure 6: A screenshot of the real-time audio application where a resonance occurs with (a) mode 2 of vibration and (b) mode 3.

Table 1 :
Parameter values used for the friction drum simulation.