d8098edAn interactive tutorial on the wave equation in acoustics: strings, drums and rooms, measured by earJeremy Magland 1This is a tutorial about the wave equation as it is used in acoustics: for strings, for drums, and for the air in a room. Every method in it is written as a MATLAB-syntax script that you can edit and run in this page, and almost every result can be played as sound. There is no server, nothing to install, and nothing is sent anywhere.
3The organising idea is the same as in the two sibling tutorials on the Dirichlet problem and on particle filters. Every problem here has an answer that is known exactly, so the error of every computation is known exactly too, and what we measure is not whether a method works but what accuracy it delivers for what cost. The difference is in the units. A numerical error in acoustics is heard, and the ear is a peculiar instrument: it is very sensitive to some errors and almost blind to others. We therefore measure errors the way a listener would notice them, as a pitch in cents or a level in decibels, and we let the reader listen to the error as well as plot it.
5### The equation
7Sound in air is a small disturbance of pressure $p$ and velocity $\mathbf{v}$ about a state of rest. Linearising the conservation of mass and of momentum about that state gives two first-order equations,
9$$\frac{\partial p}{\partial t} = -K\, \nabla\cdot\mathbf{v}, \qquad \rho\, \frac{\partial \mathbf{v}}{\partial t} = -\nabla p,$$
11where $\rho$ is the density of the air and $K$ its bulk modulus. Eliminating $\mathbf{v}$ gives the wave equation
13$$\frac{\partial^2 p}{\partial t^2} = c^2\, \nabla^2 p, \qquad c = \sqrt{K/\rho},$$
15with $c \approx 343$ m/s in air at room temperature. The same equation, with a different $c$, governs the small transverse motion of a stretched string ($c^2$ is the tension divided by the mass per unit length) and of a stretched membrane ($c^2$ is the tension per unit length divided by the mass per unit area). We use one dimension for the string, two for the drum, and three for the room, and the pressure-velocity form above reappears in Chapter 7, where it is the natural setting for an absorbing layer.
17<!--fig:settings-->
19What distinguishes the three settings is mostly the boundary. A string is fixed at both ends, so its displacement vanishes there; this is a **Dirichlet** condition. A drumhead is clamped at its rim, which is Dirichlet again. A hard wall stops the air from moving through it, so the normal velocity vanishes, and by the momentum equation so does the normal derivative of the pressure; this is a **Neumann** condition, and a wall that absorbs some of the sound is modelled by a **Robin** (impedance) condition that relates the two. A computation that is supposed to model open space needs a boundary that pretends not to be there, which is harder than any of these, and Chapter 7 is about it.
21### Where the exact answers come from
23Four sources, and each chapter says which it is using.
25**Modes.** On a string, a rectangle or a disk, the wave equation separates, and its solutions are sums of standing waves whose shapes and frequencies are known in closed form: sines for the string and the rectangle, Bessel functions for the disk. The frequencies of a disk are zeros of Bessel functions, which the shared library computes to rounding error.
27**d'Alembert's solution.** On a string, the solution is also a sum of two travelling waves, $u = \tfrac12\big(\phi(x - ct) + \phi(x + ct)\big)$, which shares no algebra with the modal sum. Chapter 1 computes the same plucked string both ways.
29**Images.** A rigid wall acts as a mirror, so the sound in a rectangular room with rigid walls is exactly the sound in free space of the source and an infinite lattice of mirror images of it. Chapters 8 and 9 use this.
31**Causality.** An explicit scheme moves information by at most one grid node per step. The same scheme on a grid that is larger by $M$ nodes on every side therefore cannot tell, during $M$ steps, that it has edges at all. Chapter 7 uses such a grid as an exact reference for the scheme's own free-space solution, which is not a formula but is correct to rounding error.
33### What is measured
35Most of what a numerical method gets wrong in a wave computation is the speed at which waves of different lengths travel, which is called **numerical dispersion**. For a standing wave, a wrong speed is a wrong frequency, and a wrong frequency is heard as a wrong pitch. We measure it in **cents**: a cent is a hundredth of an equal-tempered semitone, so an octave is 1200 cents and a frequency ratio $r$ is $1200\log_2 r$ cents. A relative frequency error $\varepsilon$ is about $1731\,\varepsilon$ cents when it is small. Every frequency in the tutorial is read off a computed signal by one routine, `wv_peakfreq`, which finds the peak of a tapered spectrum and polishes it with Newton's method; on a clean partial it is good to about ten digits, so the error it reports is the method's and not its own.
37How large a pitch error matters is a question about listeners rather than about numerical analysis, and it has no single answer. Frequency discrimination thresholds for pure tones at moderate levels in the middle of the audible range are commonly reported at around 0.2 to 0.5 per cent, which is 3 to 9 cents, and they are larger at low and high frequencies and for untrained listeners. We draw a line at 5 cents on the pitch plots as a rough guide, and nothing in the tutorial depends on its exact position.
39Two other measures appear. Levels and energies are in **decibels**, $10\log_{10}$ of an energy ratio. And a waveform error, the root-mean-square difference between a computed signal and the exact one relative to the size of the exact one, is reported where it is informative. It is often not, and the reason is worth stating at the start: a note held for one second at 440 Hz lasts 440 periods, so a pitch error of one cent, which few listeners can hear, puts the computed wave a quarter of a period out of phase by the end of it, and the waveform error is then of order one. Much of what looks like failure on a waveform plot is inaudible. Chapter 1 has an example of the converse: an error of a few per cent in the waveform that can be heard as tones which do not belong to the note.
41Cost is wall-clock time, measured inside the MATLAB layer around the method and nothing else. Computing the exact answer, measuring frequencies, and encoding the sound are not charged to any method.
43### The chapters
45**Strings**
471. **Sound from an exact solution**: the plucked string by d'Alembert and by Fourier, the harmonic series and the gaps a pluck leaves in it, and why an exact solution sampled directly aliases.
482. **Leapfrog, and the magic time step**: the standard explicit scheme, its stability limit, the one Courant number at which it is exact, the pitch of every partial measured everywhere else, and the sound of numerical dispersion.
493. **The digital waveguide**: the exact scheme without its grid, a delay loop that costs one operation per sample, and the tuning problem that comes with it.
51**Drums**
534. **A square drum**: the same scheme in two dimensions, an inharmonic spectrum, and a pitch error that depends on the direction of travel.
545. **Fewer points per wavelength**: fourth-order stencils in space, in time and in both, and why raising the order in space alone buys little.
556. **A round drum on a square grid**: Bessel modes and the sound of an inharmonic membrane, and what a curved rim costs a grid that cannot follow it.
57**Rooms**
597. **Letting sound out**: four ways of ending a grid, from a rigid wall to a perfectly matched layer, measured against a reference that is exact by causality.
608. **A room, by image sources**: the room's impulse response from a lattice of mirror images, and the reverberation formulas of Sabine and Eyring measured against it.
619. **A room, by finite differences**: the leapfrog scheme in three dimensions against the exact image sum for rigid walls, and the arithmetic that confines it to low frequencies.
63A **shared MATLAB library** page lists the function files every lab boots with.
65### Using it
67Each lab has controls above its figures and, behind the "show script" button, the method itself. Edit the script and press run: the session reboots with your version compiled in, and the figures and sounds are whatever your code produces. Ctrl-Enter (Cmd-Enter) runs from inside the editor, and "reset" restores the original. Some labs run as soon as the chapter opens; the heavier ones wait for the run button, and their captions say so.
69The listen buttons play through your browser's audio output. Start with the volume low: some of the clips are deliberately harsh, and within one comparison they share a scale, so that a louder clip means a larger signal.