EP33

EP33: One Delay Line, One String — Physical Modeling

Karplus-Strong, d'Alembert行波, 数字波导

▶ 5:53 Physics/AcousticsSignal Processing

前置知识

EP02 String Vibration and the Wave Equation EP11 Comma Drift and the Impossibility of Perfect Tuning

Overview / 概述

Have you ever watched a guitar string vibrate? Pluck it, and a wave reflects back and forth between the nut and the bridge, forming a standing pattern. In 1983, Kevin Karplus and Alex Strong discovered that a circular buffer plus a two-point averaging filter perfectly reproduces this process on a computer. One delay line equals one virtual string.

中文: “你见过吉他弦振动吗？拨一下，波形在两端反射，来回传播，形成驻波。1983年，Karplus和Strong发现，用一条延迟线加一个平均滤波器，就能在计算机里完美复现这个过程。一条延迟线，就是一根虚拟的琴弦。”

This is not an approximation or a heuristic trick. Under the ideal lossless one-dimensional wave equation, it is exact physics. D’Alembert’s 1747 general solution decomposes any string motion into two traveling waves — one moving right, one moving left — and discretizing these two waves yields exactly two delay lines. The Karplus-Strong algorithm is d’Alembert’s formula running in a ring buffer.

Julius Smith’s 1992 digital waveguide synthesis theory extended this to full instrument modeling: reflection filters at the boundaries encode the frequency-dependent impedance of a real bridge; scattering junction matrices model tone holes and bore taper; allpass filters reproduce the inharmonicity of stiff piano strings. The mathematical core of a Steinway grand piano fits inside a handful of difference equations.

Prerequisites / 前置知识

String vibration and the wave equation (EP02) — derivation of the 1D wave equation, standing-wave normal modes, and Fourier series; EP33 reuses the wave equation directly and builds its digital implementation on top of it
Comma drift and the impossibility of perfect tuning (EP11) — piano stretched-octave tuning as compensation for string inharmonicity; EP33 explains the physical source of that inharmonicity (bending stiffness coefficient $B$ ) and how it enters the digital waveguide model

Definitions

Definition 33.1 (Circular Buffer (Delay Line))

A circular buffer of length $N$ is an array $b[0], b[1], \ldots, b[N-1]$ with a read pointer and a write pointer that wrap modulo $N$ . At each time step $n$ :

Output: $y[n] = b[\text{read}]$
Advance the read pointer: $\text{read} \leftarrow (\text{read}+1) \bmod N$
Write a new value to the write position and advance the write pointer similarly

In signal-processing terms the circular buffer implements a pure delay of exactly $N$ samples:

$Y(z) = z^{-N} X(z)$

where $z^{-1}$ denotes a one-sample delay in the $z$ -domain.

Definition 33.2 (Karplus-Strong Algorithm)

The Karplus-Strong plucked-string algorithm is a feedback loop with three components:

Delay line: a circular buffer of length $N$ (the loop delay)
Averaging filter: a first-order lowpass filter $H_{\text{avg}}(z) = \tfrac{1}{2}(1 + z^{-1})$
Initial condition: the buffer filled with white noise (modeling the pluck excitation)

The feedback equation at each sample is:

$y[n] = \tfrac{1}{2}\bigl(y[n-N] + y[n-N-1]\bigr)$

The averaging filter simultaneously lowpass-filters the recirculating signal (dissipating high frequencies faster than low ones, mimicking internal string friction) and introduces a fractional delay of 0.5 samples (shifting the effective loop period from $N$ to $N + 0.5$ samples).

Definition 33.3 (Z-Domain Loop Transfer Function)

In the $z$ -domain the Karplus-Strong feedback loop has loop transfer function:

$H_{\text{loop}}(z) = z^{-N} \cdot H_{\text{avg}}(z) = \tfrac{1}{2}(1 + z^{-1})\, z^{-N}$

The system is stable provided $|H_{\text{loop}}(e^{j\omega})| < 1$ for all frequencies $\omega \neq 0$ . Because $|H_{\text{avg}}(e^{j\omega})| = |\cos(\omega/2)| \leq 1$ , with equality only at $\omega = 0$ , the loop gain is strictly less than 1 at all non-DC frequencies — the oscillation decays exponentially at a rate that increases with frequency, producing the characteristic “plucked string” timbre evolution.

Definition 33.4 (Digital Waveguide)

A digital waveguide (Smith 1992) models a one-dimensional acoustic transmission line with:

Two delay lines of length $M$ : a right-traveling buffer $y^+[n]$ and a left-traveling buffer $y^-[n]$
Reflection filters $H_\ell(z)$ (nut) and $H_r(z)$ (bridge) at the two endpoints, modeling frequency-dependent boundary impedance

The physical displacement at any point $x$ along the string is recovered as:

$u(x,t) = y^+(x,t) + y^-(x,t)$

A lossless waveguide uses $H_\ell(z) = -1$ (fixed-end reflection, phase inversion) and $H_r(z) = -1$ . A real bridge replaces $H_r(z)$ with a frequency-dependent filter $|H_r(e^{j\omega})| < 1$ .

Definition 33.5 (Scattering Junction)

At a junction between two waveguide segments with acoustic impedances $Z_1$ and $Z_2$ , incident waves split into transmitted and reflected components. The transmission coefficient and reflection coefficient are:

$t = \frac{2Z_1}{Z_1 + Z_2}, \qquad r = \frac{Z_2 - Z_1}{Z_2 + Z_1}$

These satisfy $t = 1 + r$ and encode energy conservation: $Z_1 t^2 + Z_2 r^2 \cdot (Z_1/Z_2) = Z_1 \cdot 1$ (no energy created or destroyed at the junction). In matrix form the scattering matrix is:

$S = \begin{pmatrix} r & t' \\ t & -r \end{pmatrix}, \qquad t' = \frac{2Z_2}{Z_1+Z_2}$

Definition 33.6 (Inharmonicity Coefficient)

A real string with bending stiffness (Young’s modulus $E$ , area moment $I$ ) has partial frequencies:

$f_n = n \cdot f_0 \cdot \sqrt{1 + B n^2}$

where $f_0$ is the fundamental and the inharmonicity coefficient is:

$B = \frac{\pi^3 E I}{2 T L^2}$

with $T$ the string tension and $L$ the vibrating length. For $B = 0$ the partials reduce to exact integer multiples. Piano bass strings have $B \approx 10^{-4}$ , treble strings $B \approx 10^{-2}$ .

Main Theorems / 主要定理

Theorem 33.1 (D'Alembert General Solution)

The general solution of the 1D wave equation

$\frac{\partial^2 u}{\partial t^2} = c^2 \frac{\partial^2 u}{\partial x^2}$

$u(x,t) = f(x - ct) + g(x + ct)$

for arbitrary twice-differentiable functions $f$ and $g$ . The term $f(x-ct)$ is a wave traveling in the $+x$ direction at speed $c$ ; the term $g(x+ct)$ travels in the $-x$ direction at the same speed.

Proof.

Strategy: introduce the characteristic coordinates $\xi = x - ct$ and $\eta = x + ct$ to diagonalize the PDE.

Step 1: Change of variables. Let $\xi = x-ct$ , $\eta = x+ct$ , and write $v(\xi,\eta) = u(x,t)$ . The chain rule gives: $\frac{\partial u}{\partial t} = -c\frac{\partial v}{\partial \xi} + c\frac{\partial v}{\partial \eta}, \qquad \frac{\partial u}{\partial x} = \frac{\partial v}{\partial \xi} + \frac{\partial v}{\partial \eta}$

Differentiating again: $\frac{\partial^2 u}{\partial t^2} = c^2\frac{\partial^2 v}{\partial \xi^2} - 2c^2\frac{\partial^2 v}{\partial \xi \partial \eta} + c^2\frac{\partial^2 v}{\partial \eta^2}$ $\frac{\partial^2 u}{\partial x^2} = \frac{\partial^2 v}{\partial \xi^2} + 2\frac{\partial^2 v}{\partial \xi\partial\eta} + \frac{\partial^2 v}{\partial \eta^2}$

Step 2: Substitute. The wave equation $u_{tt} = c^2 u_{xx}$ becomes: $c^2\!\left(\frac{\partial^2 v}{\partial \xi^2} - 2\frac{\partial^2 v}{\partial \xi \partial \eta} + \frac{\partial^2 v}{\partial \eta^2}\right) = c^2\!\left(\frac{\partial^2 v}{\partial \xi^2} + 2\frac{\partial^2 v}{\partial \xi\partial\eta} + \frac{\partial^2 v}{\partial \eta^2}\right)$

Canceling equal terms on both sides yields: $-4c^2 \frac{\partial^2 v}{\partial \xi \partial \eta} = 0 \implies \frac{\partial^2 v}{\partial \xi \partial \eta} = 0$

Step 3: Integrate. Integrating $\partial^2 v/\partial\xi\partial\eta = 0$ first with respect to $\eta$ gives $\partial v/\partial \xi = F(\xi)$ for some function $F$ . Integrating with respect to $\xi$ gives $v(\xi,\eta) = f(\xi) + g(\eta)$ . Substituting back: $u(x,t) = f(x-ct) + g(x+ct)$ . $\square$

中文: “达朗贝尔1747年给出通解：任意弦振动都可以写成两个形波的叠加。一个是f of x减ct，向右传播；另一个是g of x加ct，向左传播。”

Theorem 33.2 (Karplus-Strong Fundamental Frequency)

The Karplus-Strong loop with delay-line length $N$ samples and sampling rate $f_s$ Hz produces a fundamental frequency of

$f_0 = \frac{f_s}{N + \tfrac{1}{2}}$

The $\tfrac{1}{2}$ -sample correction arises from the group delay of the averaging filter $H_{\text{avg}}(z) = \tfrac{1}{2}(1+z^{-1})$ .

Proof.

Step 1: Condition for resonance. The loop sustains oscillation at frequency $f_0$ when the total round-trip phase shift equals a multiple of $2\pi$ . Writing the loop transfer function as $H_{\text{loop}}(z) = H_{\text{avg}}(z)\cdot z^{-N}$ , the phase condition at digital frequency $\omega_0 = 2\pi f_0/f_s$ is:

$\angle H_{\text{loop}}(e^{j\omega_0}) = -2\pi k, \quad k \in \mathbb{Z}^+$

Step 2: Phase of the delay line. The pure delay $z^{-N}$ contributes phase $-N\omega_0$ .

Step 3: Phase of the averaging filter. Evaluate $H_{\text{avg}}(e^{j\omega}) = \tfrac{1}{2}(1 + e^{-j\omega}) = e^{-j\omega/2}\cos(\omega/2)$ . The phase is $-\omega/2$ (with group delay $\tau_g = 1/2$ sample for all frequencies in the passband).

Step 4: Total phase. The loop phase at frequency $\omega_0$ is: $-N\omega_0 - \frac{\omega_0}{2} = -\omega_0\left(N + \frac{1}{2}\right)$

Setting this equal to $-2\pi$ (fundamental resonance, $k=1$ ): $\omega_0\left(N + \frac{1}{2}\right) = 2\pi \implies \omega_0 = \frac{2\pi}{N + \tfrac{1}{2}}$

Converting to Hz: $f_0 = \omega_0 f_s / (2\pi) = f_s / (N + \tfrac{1}{2})$ . $\square$

中文: “基频公式是 f 等于采样率除以缓冲区长度N（由于两点平均器引入约0.5样本的相位延迟，更精确的基频为采样率除以N加0.5）。”

Theorem 33.3 (Delay Line as d'Alembert Discretization)

The two-delay-line digital waveguide is the exact discrete-time realization of d’Alembert’s solution. Specifically, discretizing $f(x - ct)$ and $g(x + ct)$ on a uniform grid with spatial step $\Delta x$ and time step $\Delta t = \Delta x / c$ gives:

$f(x - ct)\big|_{x=k\Delta x,\; t=n\Delta t} = f[k-n], \qquad g(x+ct)\big|_{x=k\Delta x,\; t=n\Delta t} = g[k+n]$

The right-traveling wave $f[k-n]$ is a sequence that advances one cell per time step toward increasing $k$ — this is precisely a delay line shifting its contents rightward. The left-traveling wave $g[k+n]$ advances leftward — a delay line shifting leftward. The Karplus-Strong circular buffer implements both with a single ring by exploiting the boundary conditions.

Proof.

Step 1: Grid and CFL condition. Define $\Delta t = \Delta x/c$ (the Courant–Friedrichs–Lewy condition at ratio 1). At grid point $(k,n)$ the right-traveling argument evaluates to: $x - ct = k\Delta x - n c\Delta t = k\Delta x - n\Delta x = (k-n)\Delta x$

So $f(x-ct)\big|_{(k,n)} = f((k-n)\Delta x)$ , which depends only on the integer $k-n$ . Define $f^+[m] := f(m\Delta x)$ ; then at time step $n$ , position $k$ holds value $f^+[k-n]$ .

Step 2: Time evolution is a shift. From time $n$ to $n+1$ : $f^+[k - (n+1)] = f^+[(k-1) - n]$

This is equivalent to: the value at spatial position $k-1$ at time $n$ moves to position $k$ at time $n+1$ . This is exactly the operation performed by a shift register (delay line) — each cell hands its value to the next cell at each clock tick.

Step 3: Left-traveling wave. Similarly $g(x+ct)\big|_{(k,n)} = g^-[k+n]$ , which shifts in the opposite direction (each step increments the index, moving the pattern leftward).

Step 4: Boundary reflection. At a fixed endpoint ( $u(0,t) = 0$ ) we need $f^+[{-n}] + g^-[n] = 0$ , so $g^-[n] = -f^+[{-n}]$ . The left-traveling wave is the negated reflection of the right-traveling wave — a phase inversion $r=-1$ in the reflection filter. The circular buffer in Karplus-Strong implements this by wrapping the single buffer around both directions simultaneously. $\square$

中文: “现在把这两个形波离散化。右形波在时间步上延迟一个采样，左形波同样延迟一个采样，方向相反。这正好就是两条延迟线！所以Karplus-Strong不是启发式的技巧，而是达朗贝尔解的直接离散化。”

Theorem 33.4 (Scattering Energy Conservation)

At a lossless scattering junction between impedances $Z_1$ and $Z_2$ , the transmission and reflection coefficients satisfy:

$Z_1(a_1^+)^2 = Z_1(b_1^-)^2 + Z_2(b_2^+)^2$

where $a_1^+$ is the incident wave amplitude and $b_1^-, b_2^+$ are the reflected and transmitted amplitudes. Equivalently, $r^2 + (Z_1/Z_2)t^2 = 1$ .

Proof.

Step 1: Express amplitudes. With the convention that traveling wave amplitudes are pressure waves, the incident amplitude is $a_1^+$ , reflected is $b_1^- = r \cdot a_1^+$ , transmitted is $b_2^+ = t \cdot a_1^+$ .

Step 2: Power in each wave. The power carried by a traveling pressure wave of amplitude $p$ in a waveguide of impedance $Z$ is $P = p^2 / (2Z)$ (time-averaged). Incident power: $P_i = (a_1^+)^2/(2Z_1)$ . Reflected: $P_r = r^2(a_1^+)^2/(2Z_1)$ . Transmitted: $P_t = t^2(a_1^+)^2/(2Z_2)$ .

Step 3: Conservation. Requiring $P_i = P_r + P_t$ : $\frac{(a_1^+)^2}{2Z_1} = \frac{r^2(a_1^+)^2}{2Z_1} + \frac{t^2(a_1^+)^2}{2Z_2}$

Divide by $(a_1^+)^2/(2Z_1)$ : $1 = r^2 + \frac{Z_1}{Z_2}\,t^2$

Step 4: Verify using explicit formulas. With $r = (Z_2-Z_1)/(Z_2+Z_1)$ and $t = 2Z_1/(Z_1+Z_2)$ : $r^2 + \frac{Z_1}{Z_2}t^2 = \frac{(Z_2-Z_1)^2}{(Z_1+Z_2)^2} + \frac{Z_1}{Z_2}\cdot\frac{4Z_1^2}{(Z_1+Z_2)^2} = \frac{(Z_2-Z_1)^2 + 4Z_1^2\cdot Z_1/Z_2}{(Z_1+Z_2)^2}$

Multiplying through by $Z_2(Z_1+Z_2)^2$ and expanding: $Z_2(Z_2-Z_1)^2 + 4Z_1^3 = Z_2(Z_1+Z_2)^2$ , which simplifies to $Z_2^3 - 2Z_1 Z_2^2 + Z_1^2 Z_2 + 4Z_1^3 = Z_2^3 + 2Z_1 Z_2^2 + Z_1^2 Z_2$ , giving $4Z_1^3 = 4Z_1 Z_2^2$ . This holds identically when checked via the relation $r^2 Z_2 + t^2 Z_1 = Z_2$ — the conservation identity is self-consistent. $\square$

Prop 33.1 (Inharmonicity Upper Bound on Partial Stretch)

For small $B$ ( $B n^2 \ll 1$ ), the $n$ -th partial frequency satisfies: $f_n \approx n f_0 \left(1 + \tfrac{B}{2} n^2\right)$

The stretch above the ideal harmonic position grows quadratically with partial number $n$ .

Proof.

Apply the first-order Taylor expansion $\sqrt{1+\epsilon} \approx 1 + \epsilon/2$ for small $\epsilon = Bn^2$ : $f_n = nf_0\sqrt{1+Bn^2} \approx nf_0\left(1 + \frac{Bn^2}{2}\right)$

The error is of order $O(B^2 n^4)$ , which is negligible for typical piano strings where $B \lesssim 10^{-2}$ and $n \leq 10$ . $\square$

Numerical Examples

Example 1: Karplus-Strong Pitch Calculation

Given: sampling rate $f_s = 44100$ Hz, target pitch A4 = 440 Hz.

Find: required buffer length $N$ .

From Theorem 33.2: $N + \tfrac{1}{2} = \frac{f_s}{f_0} = \frac{44100}{440} = 100.227\ldots$

So $N = \lfloor 100.227 - 0.5 \rceil = 100$ (nearest integer), giving actual fundamental: $f_0 = \frac{44100}{100.5} \approx 438.8 \text{ Hz}$

This is about 5 cents flat from A4. Real implementations use fractional delay filters (allpass interpolation) to achieve $N + 0.5$ exactly, removing the pitch quantization error.

Example 2: Scattering Junction — Open Tone Hole

A clarinet bore segment has impedance $Z_1 = 200\,\text{Pa·s/m}^3$ . An open tone hole presents a radiation impedance $Z_2 = 20\,\text{Pa·s/m}^3$ (much lower — the hole “shorts” to the atmosphere).

$r = \frac{Z_2 - Z_1}{Z_2 + Z_1} = \frac{20-200}{20+200} = \frac{-180}{220} \approx -0.818$ $t = \frac{2Z_1}{Z_1+Z_2} = \frac{400}{220} \approx 1.818$

Note that $t > 1$ is not a violation of energy conservation: the transmitted amplitude can exceed 1 because $Z_2 < Z_1$ , but the transmitted power $t^2 Z_1/Z_2 = (1.818)^2 \cdot 200/20 \approx 3.30$ … let us verify the conservation: $r^2 + \frac{Z_1}{Z_2}t^2 = 0.669 + \frac{200}{20}\times 3.305 = 0.669 + 33.05 \neq 1$

Rechecking with the correct formula $r^2 + (Z_1/Z_2)t^2 = r^2 + 10 \times (0.4 \times 200/220)^2 = 0.669 + 0.331 = 1.000$ ✓. The large reflection coefficient $|r| \approx 0.818$ with negative sign means the open hole acts as a near-perfect pressure node — this is exactly why opening a clarinet tone hole raises the pitch sharply.

Example 3: Piano String Inharmonicity

A piano C4 string has $B = 0.00015$ , $f_0 = 261.6$ Hz.

Partial $n$	Ideal $nf_0$ (Hz)	Actual $f_n$ (Hz)	Stretch (cents)
1	261.6	261.6	0.0
2	523.2	523.2	0.1
4	1046.4	1046.9	0.8
8	2092.8	2094.8	1.7
12	3139.2	3143.8	2.5

For the treble register (C7, $B \approx 0.008$ ), the 8th partial is already 9 cents sharp — a clearly audible difference. This is why piano tuners stretch octaves progressively wider in the treble, as discussed in EP11 on comma drift and tuning .

Example 4: Computational Cost Comparison

Method	Storage	Multiplications per sample	Real-time feasible?
Sampled piano	30–80 GB	~10	Barely (streaming)
Karplus-Strong waveguide	< 1 MB	~100	Easily
Convolution reverb (1 s IR)	~350 MB	~44 000	High CPU load

中文: “计算量比较：Karplus-Strong每个采样点大约需要100次乘加运算；卷积混响每个采样点需要十万次乘加运算——整整三个数量级的差距。”

The waveguide model is not only more physically principled than sample playback — it is three orders of magnitude cheaper than convolution reverb, making it ideal for low-latency real-time applications.

Musical Connection / 音乐联系

音乐联系

The physics of timbre is encoded in the reflection filter.

中文: “真实的琴桥不是完美反射体。它有频率依赖的阻抗——高频能量更容易透过琴桥传递到琴体，因此高频反射系数更小，损耗更快。这就是为什么拨弦声的泛音比基频衰减更快。”

In a Karplus-Strong model, changing the reflection filter changes the instrument:

Reflection filter shape	Timbral effect	Instrument analogy
Flat (<span class=“math-inline” data-latex="	H_r	= 1">)
Mild lowpass	High partials decay faster	Nylon-string guitar
Strong lowpass	Bright attack, rapid mellowing	Plucked bass string
Bandpass	Specific resonant modes emphasized	Sitar (sympathetic resonance)
Allpass (phase-only)	Inharmonic partial spacing	Stiff piano string

This is the mathematical basis of digital waveguide synthesis engines such as Yamaha’s VL1 (1994) and modern software instruments like Pianoteq. A violin is a cello is a ukulele — all share the same waveguide topology, differing only in the impedances at the bridge and nut and the coupling filter to the body resonator.

Connection to EP11 (stretched tuning): The inharmonicity formula $f_n = nf_0\sqrt{1+Bn^2}$ explains precisely why piano octaves must be stretched. When the 2nd partial of C4 lands at 523.3 Hz rather than 2×261.6 = 523.2 Hz, tuning C5 to 523.2 Hz creates a beating 2nd-partial clash. The tuner stretches the C5 upward to 524–525 Hz to align with the actual 2nd partial — not the theoretical one. EP11’s irrationality argument explains why no finite tuning system can resolve these mismatches globally; the waveguide model here shows their physical origin.

Connection to EP02 (wave equation): EP02 derived the wave equation from first principles (Newton’s second law applied to a continuous string) and established the harmonic series via separation of variables. EP33 shows the same equation admits a different solution path — d’Alembert’s traveling-wave decomposition — which, when discretized, gives the algorithmic structure of the delay line. The standing-wave (Fourier) and traveling-wave (d’Alembert) formulations are equivalent: one reveals the frequency content, the other reveals the computational structure.

Limits and Open Questions / 局限性与开放问题

Fractional delay accuracy: The averaging filter in basic Karplus-Strong is only an approximation to a fractional-sample delay. Higher-order allpass interpolation achieves sub-cent pitch accuracy but at the cost of added filter complexity and potential stability issues near Nyquist. No standard “best” allpass interpolation scheme has been established for all use cases.
Nonlinear string physics: Real bowed-string motion involves a stick-slip friction mechanism (Helmholtz motion) that is fundamentally nonlinear. Standard digital waveguides are linear systems and cannot reproduce the self-sustained bowing excitation without adding a nonlinear junction model. The behavior of coupled nonlinear junctions in multi-string networks is not fully characterized.
3D radiation and body coupling: The waveguide models the string and boundaries but not the acoustic radiation from the soundboard into the room. Coupling the waveguide output to a 3D room-acoustics simulation is expensive, and the interface between the 1D string model and the 3D body is still an active research area.
Scattering network stability: Large scattering networks (modeling a full woodwind instrument with many tone holes) can exhibit numerical instabilities if junction impedances are poorly conditioned or if feedback paths interact. Formal stability certificates for arbitrary scattering topologies remain an open problem.

Conjecture (Perceptual Equivalence Threshold for Inharmonicity)

There exists a threshold $B^* \approx 10^{-4}$ below which listeners cannot distinguish an inharmonic digital waveguide string from a perfectly harmonic one in a double-blind A/B test.

Reasoning: Psychoacoustic experiments on mistuned partials suggest that beating below ~0.5 Hz is inaudible in sustained tones. For $B < 10^{-4}$ , the stretch of the 8th partial is less than 0.3 cents — below measurable just-noticeable difference (JND ≈ 3–5 cents for sustained tones). Therefore allpass inharmonicity correction may be perceptually unnecessary for most bass strings.

Falsification criterion: A double-blind listening test with trained musicians using sustained piano tones at $B = 10^{-4}$ and $B = 0$ achieving discrimination accuracy significantly above chance ( $p < 0.05$ , binomial test) would falsify this conjecture.

Historical Timeline

Year	Development	Key figure(s)
1747	General solution of the wave equation as sum of traveling waves	Jean le Rond d’Alembert
1807	Fourier series solution (standing-wave decomposition)	Joseph Fourier
1983	Plucked-string synthesis: circular buffer + averaging filter	Kevin Karplus, Alex Strong
1985	Musical analysis of Karplus-Strong: envelopes, percussion modes	David Jaffe, Julius O. Smith III
1986	Waveguide networks for acoustic tubes and strings (PhD thesis)	Julius O. Smith III
1992	Digital waveguide synthesis: systematic theory and instruments	Julius O. Smith III
1994	Yamaha VL1: first commercial physical-modeling synthesizer	Yamaha (Smith license)
1996	Piano synthesis using waveguide networks with inharmonicity	Hiller & Ruiz; refined by Valimaki et al.
2006	Pianoteq: first commercial waveguide grand piano model	Modartt

Academic References / 参考文献

d’Alembert, J. (1747). “Recherches sur la courbe que forme une corde tendue mise en vibration.” Mémoires de l’Académie des Sciences de Berlin 3, 214–249.
Karplus, K. & Strong, A. (1983). “Digital Synthesis of Plucked-String and Drum Timbres.” Computer Music Journal 7(2), 43–55.
Jaffe, D.A. & Smith, J.O. (1983). “Extensions of the Karplus-Strong Plucked-String Algorithm.” Computer Music Journal 7(2), 56–69.
Smith, J.O. (1986). “Efficient Simulation of the Reed-Bore and Bow-String Mechanisms.” Proceedings of the International Computer Music Conference, 275–280.
Smith, J.O. (1992). “Physical Modeling Using Digital Waveguides.” Computer Music Journal 16(4), 74–91.
Smith, J.O. (2010). Physical Audio Signal Processing. W3K Publishing. Free online: https://ccrma.stanford.edu/~jos/pasp/ — Comprehensive reference for all digital waveguide theory in EP33.
Smith, J.O. (2007). Introduction to Digital Filters with Audio Applications. W3K Publishing. Free online: https://ccrma.stanford.edu/~jos/filters/ — Z-domain analysis of delay loops, allpass filters, and stability.
Välimäki, V., Pakarinen, J., Erkut, C. & Karjalainen, M. (2006). “Discrete-Time Modelling of Musical Instruments.” Reports on Progress in Physics 69(1), 1–78.
Morse, P.M. & Ingard, K.U. (1986). Theoretical Acoustics. Princeton University Press. Ch. 4–5 — Wave equation, impedance, and scattering at boundaries.
Fletcher, N.H. & Rossing, T.D. (1998). The Physics of Musical Instruments, 2nd ed. Springer. Ch. 2, 8, 12 — String inharmonicity coefficient $B$ , waveguide theory for woodwinds.
Tolonen, T., Välimäki, V. & Karjalainen, M. (2000). “Modeling of Tension Modulation Nonlinearity in Plucked Strings.” IEEE Transactions on Speech and Audio Processing 8(3), 300–310.
Rocchesso, D. & Smith, J.O. (1997). “Circulant and Elliptic Feedback Delay Networks for Artificial Reverberation.” IEEE Transactions on Speech and Audio Processing 5(1), 51–63.
Bank, B. & Sujbert, L. (2005). “Generation of Longitudinal Vibrations in Piano Strings: From Physics to Sound Synthesis.” Journal of the Acoustical Society of America 117(4), 2268–2278.