Tutorial
Introduction to mass-interaction physical modelling

In this first of a series of online MI tutorials, we’ll be looking into the fundamentals of mass-interaction physical modelling for sound synthesis. First off, we’ll start with a short introduction placing its history and place within the wider scope of physical modelling techniques, and then we’ll get down into the grizzly algorithms and implementation details. For this, we will be working with the online Mass Interaction Model Scripting tool (MIMS), as well as the Faust environment (specifically the online FaustIDE which will allow us to run all of our models directly from a browser).

By the end of this tutorial you should have a decent understanding of mass-interaction principles, know enough to start creating your own physical models with MIMS, be able to test these models using Faust, and even use the great Faust utilities to compile these models into code for almost any target (VST plugin, embedded platforms, Unity, Max, Pure Data, etc) !

Without further ado, let’s get started.

Background in classical mechanics

One of the most elementary approaches to solve problems in classical mechanics is to model systems as constructions of basic elements (such as punctual masses) which act upon each other through various interaction forces, such as gravitational forces, springs obeying Hookes Law, and so forth. The system is then solved by combining the partial differential equations (PDEs) governing each element and interaction into a global PDE system.

Two masses linked by springs, between two fixed points

In simple cases like the above, systems may be solved exactly to obtain the equations of motion for a given mass of the system. However, for more complex scenarios, systems will generally be solved using numerical methods.

Numerical methods for discrete-time physics simulation

There are a panoply of ways in which computers can be used to model and simulate physical systems: finite difference schemes, finite element analysis, lumped methods (such as mass-interaction systems)... in the realm of sound synthesis, we can also mention methods more grounded in signal processing, such as digital waveguides, wave digital filters, or modal synthesis techniques.

In musical acoustics, we are often interested in simulating physical structures more complex than the simple mechanical example presented above: we might want to approximate a string, a plate, or even a 3D acoustical propagation space. This brings us to some key questions:

  • Do we still want to describe our physical system as a construction of basic physical components, or is it preferable to mathematically express a more complex “integrated” system?
  • A consequence and follow-up question stemming from the above is : where do we put the line between real (i.e. continuous time) physics and virtual (i.e. discrete time) simulation: in other words, do we formulate a physical problem formally, then solve it numerically, or do we create a discrete-time paradigm in which we can directly express physical constructions ?

Basic approaches to physically-based sound synthesis

Different paradigms for physically-based sound synthesis approach these considerations from different angles. For instance:

  • In modal synthesis, we formulate mathematical models of (generally linear) physical systems in order to obtain their modal decomposition into spectral components (an assembly of sinusoidal vibration modes with specific frequencies and damping times. Once we have this, we can simply synthesize the sound of the model by means of additive synthesis.
  • In waveguide methods, we decompose physical motion into travelling waves in opposite directions, linked together by scattering junctions (such as reflexions at the bridge of a string, contact with an obstacle in which some waves are reflected backwards and others are transmitted forwards, etc.).
  • In finite difference numerical methods, we formulate the mathematics for a given problem. Once we can express the maths of the continuous time system, we can figure out ways to discretise the system so that it can be solved numerically.
  • In mass-interaction modelling (often referred to as mass-spring), the idea is to build systems using discretised approximations of elementary physical elements that can be assembled to form complex networks: the physical behaviour emerges from the designed physical network.

Illustration by example: an ideal string

To ground this in a simple concrete example, let’s say we wanted to model and simulated a simple ideal string:

  • Using modal synthesis we would calculate the modes then use a resonator bank to synthesise the sound.
  • Using waveguides we would implement some delay lines with appropriate lengths, boundaries and junctions for input.
  • Using the finite difference approach, we would express equations for the continuous time system using the 1D wave equation with given boundary conditions, then discretise it using a spatial grid (decomposing continuous matter into discrete set of grid points) and a temporal update scheme (allowing to calculate the dynamics of the system from a given initial state).
  • Using the mass-interaction approach, we would consider the string as an assembly of punctual masses connected by springs, and then use the discrete equations for each of them to update the system over time.

Remarks on finite difference schemes and mass-interaction models:

Finite Difference methods rely on mathematical formulations of (possibly complex) physical systems. This integration gives leverage for a variety of numerical resolution schemes that may be explicit (future values may be calculated only based on current (and possibly previous) values) or implicit (requiring iterative solvers to calculate future values). The choice of a discretisation scheme affects both numerical stability and the amount of discrepancy (such as numerical dispersion) between the discrete-time system and the continuous time model.

However, these methods are generally only semi-modular in that they allow to connect various macro structures together (such as strings, plates, beams) through interactions. There is no simple way to create a structure that has not been mathematically formalised beforehand, and for some irregular or inhomogeneous physical constructions this might well be impossible.


Mass-interaction schemes are more limited in terms of discretisation schemes, due to their modularity: since all systems are broken down into discrete-time mass algorithms and discrete-time interaction algorithms, discretisation schemes are rather limited. On the other hand, this modularity means that any kind of model can be assembled from the ground up.

Note: there are several cases (notably ideal linear strings and meshes) in which lumped mass-interaction methods and finite difference schemes yield exactly the same result, with grid points being equivalent to punctual masses. For more on this, we refer readers to the first chapters of Stefan Bilbao’s book Numerical Sound Synthesis.

History of mass-interaction physics

  • CORDIS-ANIMA can be considered as the original mass-interaction physical modelling formalism, coming into existence in its prototypical form at ACROE as soon as the early eighties, including pioneering views as to the potential of coupling with force-feedback technologies. It forms the basis for Mimesis, an environment for 3D physical modelling destined for animation, and Genesis, an environment for physical modelling sound synthesis based on 1D mass-interaction networks - both of which are off-line modelling and simulation tools providing advanced user interfaces for designing complex mass-interaction physical models.
  • Following years saw the emergence of several direct variations on Cordis-Anima, providing open implementations for sound synthesis in the form of Tao, Pmpd's integration into Pure Data, or Cymatic, a tool allowing for model design and real-time force-feedback interaction.
  • A third wave of mass-interaction tools have appeared in the last decade, driven by open-source initiatives: HSP (haptic signal processing) provided a first means for audio rate simulation in Max/MSP, whereas Synth-a-modeler provides a first Faust-based engine allowing compilation for a variety of targets and platforms. It has since been extended with a modelling user interface and bridges allowing for interconnection between mass-interaction, waveguide and modal synthesis elements. And in the last two years, we can add mi-faust (a different take on mass-interaction models in Faust, and the basis of this tutorial material), miPhysics (Java library for visual audio and haptic 3D mass-interaction models) and mi-gen~ (an efficient mass-interaction framework allowing to compile models into gen~ code).

A Basic Mass-Interaction Scheme

So, we’ve established that mass-interaction schemes form discrete time algorithms for elementary physical components of two types: masses and interactions. As you probably know, many different numerical approximations can be used to integrate PDEs, from simple Euler methods, to higher order integrators such as 4th order Runge-Kutta, etc.

For an extensive review of numerical methods for sound synthesis, we refer you to Stefan Bilbao’s book, and Valimäki et al.’s state of the art paper (provide links).

A possible discrete-time scheme for punctual masses

The motion equation for a continuous time mass is given by Newton’s second law:

\begin{equation} \label{eq:fma} \begin{aligned} f = ma = m \frac{d^2 x}{dt^2} \end{aligned} \end{equation}

Where f is the force applied to the mass, m is its inertia a its acceleration, and x its position. Applying the second-order central difference scheme, with the sampling interval noted ΔT, a discrete equation of the mass can be formulated as follows:

\begin{equation} \label{eq:centered} \begin{aligned} f(t) = m.\frac{x(t\!+\!\Delta{T}) - 2 x (t) + x(t\!-\!\Delta{T}) }{\Delta{T}^2} \end{aligned} \end{equation}

Equation (2) can be normalized to unity, and rearranged in order to express the mass’ position update scheme (discrete-time positions and forces are noted X and F) :

\begin{equation} \label{eq:discretetimemass} \begin{aligned} X_{(n+1)} = 2 X_{(n)}- X_{(n-1)} + \frac{F_{(n)}}{M} \end{aligned} \end{equation}

With M, the discrete time inertial parameter defined as :

\begin{equation} \label{eq:Mm} \begin{aligned} M =\frac{m}{\Delta{T}^2} \end{aligned} \end{equation}

Hence, the basic discrete-time mass module produces new position data based on its current position, previous position, the “discrete-time” mass parameter M, and the sum of forces applied to the mass from the previous interaction computation step.

The initial position X(0), delayed initial position X( − 1) (which infers initial velocity) and initial force F(0) must be supplied at the start of the computation.

Implementation in Faust

FAUST (Functional Audio Streams) is, as the name suggests, a functional programming language for audio or any kind of digital signal processing. Its syntax allows expressing signal processing algorithms efficiently, which are transformed into functional block diagrams.
Atop of the language, FAUST allows exporting algorithms to virtually any target language or platform, so it’s a great tool for highly cross-platform audio development. And to finish things off, the FAUST team are at the forefront of audio technologies on the web: you can write, compile and test FAUST code directly on your browser with a dedicated IDE !

For more information: https://faust.grame.fr/
Online IDE: https://faustide.grame.fr/
Faust documentation: https://faustdoc.grame.fr/
Faust mass-interaction library documentation: https://faustlibraries.grame.fr/libs/mi/ 


mass(m, grav, x0, xr0) = equation
with{
  A = 2;
  B = -1;
  C = 1/m;
  equation = x
    letrec{
        'x = A*(x : initState(x0))
           + B*(x' : initState((xr0,x0)))
           + (_-grav)*(C);
    };
};

Block diagram of the mass algorithm

A possible implementation of a damped spring

The elastic force applied by a linear spring with a stiffness k and a resting length of l0 = 0 connecting a mass m2 at the position x2 to a mass m1 at the position x1 is given by Hookes law:

\begin{equation}\label{eq:HookeSpring} \begin{aligned} f_{s{1\rightarrow2}} = -k.(x_2 - x_1) \end{aligned} \end{equation}

The exact equivalent of this equation in discrete time is :

\begin{equation}\label{eq:discreteSpring2} \begin{aligned} F_{s{1\rightarrow2}(n)} = -K.(X_{2(n)} - X_{1(n)}) \end{aligned} \end{equation}

Where the discrete-time stiffness parameter K = k. The friction force applied by a linear damper with a damping parameter z connecting the same two masses is :

\begin{equation}\label{eq:Damper} \begin{aligned} f_{d{1\rightarrow2}} = -z.\frac{d (x_2\! -\! x_1)}{dt} \end{aligned} \end{equation}

Using the Backward Euler difference scheme, the frictional force can be formulated in discrete-time as :

\begin{equation}\label{eq:DiscreteDamper} \begin{aligned} f_{d{1\rightarrow2}}(t) = -z .\frac{(x_2(t)\! -\! x_1(t)) - (x_2(t\! -\! \Delta{T}) - x_1(t\! -\! \Delta{T}))}{\Delta{T}} \end{aligned} \end{equation}

Which after normalization becomes :

\begin{equation}\label{eq:normDamper} \begin{aligned} F_{d{1\rightarrow2}(n)} = -Z .((X_{2(n)}\!-\!X_{2(n-1)})\!-\!(X_{1(n)} - X_{1(n-1)})) \end{aligned} \end{equation}

With the discrete time inertial parameter Z defined as :

\begin{equation}\label{eq:discreteDampParam} \begin{aligned} Z =\frac{z}{\Delta{T}} \end{aligned} \end{equation}

The global equation of the force applied by the dampened spring is composed of Fs and Fd :

\begin{equation}\label{eq:springDampAll} \begin{aligned} F_{(n)} =& -K.(X_{2(n)}\!-\!X_{1(n)})\\ & -Z .((X_{2(n)}\!-\!X_{2(n-1)}) - (X_{1(n)}\!-\!X_{1(n-1)})) \end{aligned} \end{equation}

It is applied symmetrically to each mass (Newton’s third law):

\begin{equation}\label{eq:springFrcApply} \begin{aligned} F_{{2\rightarrow1}(n)} &= - F_{(n)} \\ F_{{1\rightarrow2}(n)} &= + F_{(n)} \end{aligned} \end{equation}

Implementation in Faust 


springDamper(k, z, x1r0, x2r0, x1, x2) = k*deltapos + z*deltavel <: *(-1),_
with{
    deltapos = x1-x2;
    deltavel = (x1 - x1' : initState(x1r0)) - (x2 - x2' : initState(x2r0));
};

Faust implementation of the spring damper algorithm

Block diagram of the spring damper algorithm

Composing elements: oscillator

A linear harmonic oscillator is obtained by combining equations (3) and (11), in the case where X1 is a fixed point set to X1(n) = 0, n ∈ Z. This results in :

\begin{equation}\label{eq:discreteOsc} \begin{aligned} X_{(n+1)} = \Big(2\!-\!\frac{K\!+\!Z}{M}\Big).X_{(n)} + \Big(\frac{Z}{M}\!-\!1\Big).X_{(n-1)} + \frac{F_{(n)}}{M} \end{aligned} \end{equation}

Put it all together and voilà. We have a discrete time oscillator.

We invite you to check SMC 2019 paper about stability considerations.

How this looks in Faust

Block diagram of the coupled ground, mass and spring damper system.

Is this a good scheme?

The discretisation scheme developed here is explicit, not unconditionally stable, presents some numerical dispersion and frequency warping… but it is simple and efficient, and it is energy conserving, which is a very important property for simulated acoustical systems, which are generally stiff and rather undamped.

As a comparison point: many computer graphics mass-spring approaches aim to obtain numerical stability at the lowest possible sampling rate, and are not too concerned with overdamping often induced by implicit integration schemes (Keyser and House’s book Foundations of Physically-based Modeling and Animation is a great resource on this topic).

If you would like to go deeper into mass-interaction discretisation schemes and their properties, dive into these papers by Dan Morgan and Sanzheng Qiao: Accuracy and stability in mass-spring systems for sound synthesis and Analysis of damped mass-spring systems for sound synthesis. Alexandros K’s paper CORDIS-ANIMA System Analysis is also a good entry point.

Model Design and Implementation

Ground Rules

To recap briefly, mass-interactions models are topological networks build with very simple connexion rules:

  • Masses are nodes of the network. They receive forces from connected interactions and output their position to the interactions.
  • An interaction connects a mass to another mass through a certain interaction force. Interactions are calculated by calculating a force according to the positions and/or velocities of the two connected masses, and output equal and opposite forces to both masses.
  • Any number of interactions may connect to a mass. However, an interaction may only connect two physical elements together.

The most extensive theoretical description of modular mass-interaction formalism is given in the founding CORDIS-ANIMA article. Although implementations have varied after this, core concepts remain very similar.

A more general block diagram for the harmonic oscillator

The above oscillator is correct in terms of signal routing, but we can’t do much with it as it is totally hermetic and has no inputs and no outputs. To solve this, we can create a slightly more complex and general routing for the block diagram. It contains routing functions that pass mass positions to interactions, and interaction forces back to masses, and passes any input and output functions through. 


in1 = button("Frc Input 1"): ba.impulsify * 0.25;
OutGain = 1;
model = (
    mass(1., 0, 0., 0.),
    ground(0.),
    par(i, nbFrcIn,_):
    RoutingMassToLink ,
    par(i, nbFrcIn,_):
    springDamper(0.1, 0.0003, 0., 0.),
    par(i, nbOut+nbFrcIn, _):
    RoutingLinkToMass
)~par(i, nbMass, _):
par(i, nbMass, !), par(i, nbOut , _)
with{
    RoutingMassToLink(m0, m1) = /* routed positions */ m1, m0, /* outputs */ m0;
    RoutingLinkToMass(l0_f1, l0_f2, p_out1, f_in1) = /* routed forces  */ f_in1 + l0_f2, l0_f1, /* pass-through */ p_out1;
    nbMass = 2;
    nbFrcIn = 1;
    nbOut = 1;
};
process = in1 : model:*(OutGain);

A slightly more general block diagram for the harmonic oscillator.

Building Physical Networks: MIMS

MIMS (Mass Interaction Model Scripter) is a scripting tool that allows expressing mass-interaction models in the simplest possible manner, by describing the topological network abstracted from any implementation specificities. It allows to:

  • Create and connect physical elements
  • Define labelled parameters that can be used for physical parameters.
  • Create position and/or force inputs and outputs allowing to interact with the physical model.

An online version of MIMS is available here:
MIMS Online
It allows generating model code for Faust and gen~ (a low level DSP engine in Max/MSP).

A guided tour based on the oscillator

Now that we have a relatively clear picture of how basic mass-interaction models are built and how elements can be expressed in Faust, let’s get to some practical work:

Below, you will find three codes for the harmonic oscillator. Each one is expressed in the form of a MIMS model script (.mdl file) and a Faust implementation, partially generated from the MIMS description.

Open the three variations of the harmonic oscillator

All three models are equivalent, and only differ by use they make of integrated elements (such as the oscillator, or spring-damper). Explore the code and block diagrams to understand what differs between them.







If we want some real time control over the models that we create, we can use labelled parameters to specify module properties. This way, we can modify the parameters of several modules at the same time, hook a parameter up to some user interface controls, etc.

In MIMS scripts this is done via the param and audioParam primitives, which allow defining parameters at a control rate or at audio signal rate. There is no real actual difference between control rate and audio rate parameters in Faust, as all signals run at audio rate.

Open the Faust codes below to familiarize yourself with the way parameters are defined and can be controlled or hooked up to UI elements.





So, now we have an oscillator that we can control by sending force impulses to it and changing its physical properties, but where can we go from here? Let’s stick with the oscillator as the simplest vibrating object, and explore other ways of interacting with it.

Note: listening to a position of a mass is like placing a contact microphone on a table: although you will get sound produced by the motion of the table, you it will not sound the same as the table’s vibration creating acoustic pressure waves through a medium (air), bouncing of surfaces (reflexions), some combination of which makes its way to your ears… but more complex listening situations are a topic for another day.

Hammer Time

If we want some more interesting excitation than just sending bolts of force to the oscillator, we’ll have to build some things around it, for instance a hammer.

A hammer can be as simple as creating another mass in the system, initially placed above or below the oscillator and configured with an initial velocity, and adding a contact interaction between the hammer and the oscillator (a viscoelastic interaction that is active when the modules are interpenetrating).

If we want to be able to hit the oscillator several times, the above won’t quite do: the hammer will strike once, then be sent flying, never to come back. We can solve this by attaching the hammer to a fixed point via a slack spring damper, configured near critical damping. Sending force impulses to the hammer will launch it towards the oscillator, and the critically damped spring system will bring the hammer gently back to its resting position.

Open the hammer example code below.

Using the Faust IDE plotting tools, we can observe the contact between the hammer and the oscillator.

Try changing the contact parameters (stiffness and damping) to see how the behaviour changes.

You can also change the properties of the hammer spring system to experiment with rebounds.



Using gravity

Another way that we could use the hammer is to set it up with some gravity, so that it falls down onto the oscillator, bounces for a while until it stops:

Open the gravity example code below.

Experiment with some parameters: the gravity value, stiffness and damping of the contact between hammer and oscillator, etc.



Now, instead of inputting forces into our physical system, we are going to create a special type of mass, one whose position can be controlled from the outside world.

Plucking the Oscillator

A “pluck” interaction may be constructed as a piecewise linear elastic interaction between two masses: upon contact, the masses push away from each other, until reaching a tipping point at which the plucking mass “moves through” the other mass. Note that in order to be able to pluck either upwards or downwards, the non linear function must be symmetrical.

Open the plucked oscillator example below.

Move the position controlled mass around with the slider: can you see when the mass presses into the oscillator, and when it is released? Use the plot tools in Faust IDE to visualise both the picking mass and the oscillator



Bowing the Oscillator

Similarly, non-linear viscous interactions can be designed to model friction between physical elements (this is covered extensively in theoretical musical acoustics, such as in Rossing and Fletcher’s book). Up to a given threshold in relative velocity, a sticking force tends to make the masses drag each other along. Past the tipping point, the masses slide until reaching the sticky section of the interaction again. This is called stick-slip motion.

Open the bowed oscillator example below.

Move the position controlled mass around with the slider: isn’t this the coolest squeaky sound that you’ve ever heard?

Play around with the bowing speed. The resulting pitch varies dynamically, since the bow/oscillator combination is a nonlinear system



Non Linear Oscillator

Of course, the oscillator itself may display some non-linear properties, for instance by implementing springs that are slightly more complex and whose stiffness depends on their elongation.

Open the nonlinear oscillator example below.

Use the slider to change the scale of the nonlinear term in the oscillator’s spring.

Can you hear the pitch glide when striking the oscillator? Try hitting it harder, softer, changing the coefficients of the nonlinear spring, etc.



Building More Complex Models

Two Mass Chain

Moving on to bigger things, let's build a model slightly more complex than the oscillator: a chain with two masses and a fixed point:

Open the 2 mass chain example below.

Listen to the sound, have a glimpse at the spectrogram: how does the second mass contribute to the system?

Open up the MIMS online scripter and try out some other constructions: what happens with three chained masses? Or when inhomogeneous parameters are used?



Tiny String

Continuing on this path, we can build a short string, composed of 8 masses, two fixed points and spring-dampers.



If we look at the block diagram and the routing functions: it’s beginning to get quite hairy… this is why generation tools such as MIMS are vital for designing more complex models.

Some Larger Examples

The examples below show some larger constructions that are possible with tools like MIMS, starting to sound like plausible physical vibrating structures. Have a look through all of them, and test them out in the Faust IDE:









Leveraging physical models within Faust

There are some pretty neat things that Faust can help us with when designing physically-based sound synthesis systems. One of them is inherent handling of polyphony. While we won’t go too deep into specifics here (head over to the Faust documentation for that), we will quickly demonstrate how easy it is to create a polyphonic instrument with dynamically allocated mass-interaction models.

Open the Polyphonic example below in the Faust IDE

Activate Faust’s poly voices mode (located on the left panel) to a number of voices (e.g 16), and activate the computer keyboard as MIDI input (on the right panel).

Tap some keys, mess around with model parameters: there you have it, a polyphonic instrument based on several voices of a very simple physical model!



Another possibility is using physical models not for their output as synthesised sound, but as any other sound or control unit in a modular setup. For instance, the model below uses a very slack short string as an amplitude LFO to modulate a white noise signal:

Open the PhysicalLFO example below in the Faust IDE.

Can you explain why the modulation pattern starts quite complex and gradually becomes more regular?

Try changing physical parameters and observe the effect on the sound output.