Plucked string harmonics: why a pluck at a third has no third harmonic

import json
from pathlib import Path

import numpy as np
import matplotlib.pyplot as plt
from scipy.fft import dst

plt.rcParams.update({
    "figure.figsize": (7, 3.6), "figure.dpi": 110,
    "axes.spines.top": False, "axes.spines.right": False,
    "axes.grid": True, "grid.alpha": 0.25,
    "font.size": 11, "lines.linewidth": 1.8,
})
INK, ACCENT, SECOND, MUTED = "#1f2a44", "#c8553d", "#2a7f9e", "#8a8f98"

L, F1 = 0.65, 110.0                 # the A string of a classical guitar: scale length / m, fundamental / Hz
C = 2 * L * F1                      # wave speed / m/s, from the round trip
T = 1 / F1                          # period / s
N_MODES = 50                        # the modes the game carries
SLOW = 1000                         # the picture runs 1000 times slower than the sound
A_RANGE = (0.05, 0.95)              # where the game lets you pluck, as a fraction of L from the bridge
PLUCKS = [1 / 2, 1 / 3, 1 / 10]
SNAP, SNAP_WITHIN = [1 / 2, 1 / 3, 2 / 3, 1 / 4, 3 / 4, 1 / 5, 1 / 10, 9 / 10], 0.01
FS, SECONDS, TAPER = 44100, 3.0, 0.25   # sample rate / Hz, length of a pluck's sound / s, its fade at the end / s
RMS = 0.15                          # every pluck's sound at the same rms, so the ear compares tone, not volume
n = np.arange(1, N_MODES + 1)
TAU_N = 3.0 / (1 + n / 8)           # decay time of harmonic n / s, a model: the high ones fade first
# No random numbers anywhere: every number on the page follows from the formulas.

def triangle(a, x):
    """The released string, height 1 at x = aL; x is measured from the bridge."""
    return np.where(x <= a * L, x / (a * L), (L - x) / ((1 - a) * L))

def coefficients(a, modes=N_MODES):
    """b_n / h for n = 1 .. modes: the triangle's Fourier coefficients."""
    k = np.arange(1, modes + 1)
    return 2 * np.sin(k * np.pi * a) / (k**2 * np.pi**2 * a * (1 - a))

def mode_sum(a, x, t):
    """y(x, t) / h as the game draws it: the sum of 50 modes, mode n swinging at n·f1."""
    b = coefficients(a)
    return (b * np.cos(2 * np.pi * n * F1 * t)) @ np.sin(np.outer(n, np.pi * x / L))

def extension(a, u):
    """d'Alembert's F: the triangle on [0, L], odd about 0, period 2L. Reduce u modulo 2L into [0, 2L);
    on the second half, oddness and the period give F(u) = F(u − 2L) = −F(2L − u)."""
    u = np.mod(u, 2 * L)
    return np.where(u <= L, triangle(a, np.minimum(u, L)), -triangle(a, np.clip(2 * L - u, 0, L)))

def dalembert(a, x, t):
    """y(x, t) / h exactly: half of F running right, half running left."""
    return 0.5 * (extension(a, x - C * t) + extension(a, x + C * t))

def kinks(a, t):
    """Where the kinks are at time t, in m: the corners of F, at ±aL modulo 2L, carried along by ±ct."""
    s = np.array([a * L - C * t, -a * L - C * t, a * L + C * t, -a * L + C * t])
    x = np.mod(s, 2 * L)
    return np.sort(x[(x > 1e-9) & (x < L - 1e-9)])

def bars(a, modes=16):
    """|n·b_n| / b_1: how loud harmonic n sounds against the fundamental, what the game's bars show."""
    k = np.arange(1, modes + 1)
    return np.abs(k * coefficients(a, modes) / coefficients(a, 1)[0])

def energy_shares(a, modes=N_MODES):
    """E_n / E, exact: 2 sin²(nπa) / (n²π² a(1 − a))."""
    k = np.arange(1, modes + 1)
    return 2 * np.sin(k * np.pi * a) ** 2 / (k**2 * np.pi**2 * a * (1 - a))

def pluck_sound(a, damped=True, seconds=SECONDS, fs=FS):
    """The game's synthesis in NumPy: the force on the bridge as 50 cosines at n·f1 with amplitudes n·b_n,
    each fading as exp(−t/τ_n), scaled to the rms RMS; the last TAPER seconds fall to zero."""
    t = np.arange(int(round(seconds * fs))) / fs
    amp = n * coefficients(a)
    amp = amp * RMS / np.sqrt(np.sum(amp**2) / 2)
    y = np.zeros_like(t)
    for k in range(N_MODES):
        wave = amp[k] * np.cos(2 * np.pi * n[k] * F1 * t)
        y += wave * np.exp(-t / TAU_N[k]) if damped else wave
    if damped:
        y *= np.clip((seconds - t) / TAPER, 0, 1)
    return y

def draw_scene(ax_string, ax_bars, a, readout=True):
    """The game's scene at the moment of release: the held string above, the bars |n·b_n| / b_1 below."""
    ax_string.set_axis_off()
    ax_string.set(xlim=(-4, 69), ylim=(-0.55, 1.45))
    ax_string.axhline(0, color=MUTED, lw=1, xmin=4 / 73, xmax=69 / 73)
    for frac, label in [(1 / 10, "1/10"), (1 / 3, "1/3"), (1 / 2, "1/2")]:
        ax_string.plot([100 * frac * L] * 2, [-0.18, -0.28], color=MUTED, lw=1)
        ax_string.text(100 * frac * L, -0.32, label, color=MUTED, ha="center", va="top")
    for x0, name, ha in [(0, "bridge", "left"), (100 * L, "nut", "right")]:
        ax_string.plot([x0, x0], [-0.12, 0.12], color=INK, lw=4, solid_capstyle="butt")
        ax_string.text(x0, 1.2, name, color=MUTED, ha=ha, va="bottom")
    x = np.linspace(0, L, 201)
    ax_string.plot(100 * x, triangle(a, x), color=INK, lw=2)
    ax_string.plot(100 * a * L, 1, "o", color=INK, ms=6)
    if readout:                                                     # the game's readout, top right, above the tallest bar
        ax_bars.text(0.99, 0.98, f"fundamental: {100 * energy_shares(a, 1)[0]:.0f} % of the energy",
                     transform=ax_bars.transAxes, color=ACCENT, ha="right", va="top", weight="bold")
    k = np.arange(1, 17)
    v = bars(a)
    ax_bars.bar(k, v, width=0.62, color=ACCENT)
    for kk in k[v < 0.005]:
        ax_bars.text(kk, 0.02, "0", color=MUTED, ha="center", va="bottom")
    ax_bars.set(xlim=(0.4, 16.6), ylim=(0, 1.3), xticks=[1, 4, 8, 12, 16], yticks=[0, 0.5, 1], xlabel="harmonic n")

print(f"c = {C:.0f} m/s, period {1e3 * T:.2f} ms, a kink crosses the string in {1e3 * L / C:.2f} ms")

x = np.linspace(0, L, 400)
fig, axes = plt.subplots(5, 1, figsize=(7, 5), sharex=True, sharey=True, gridspec_kw=dict(hspace=0.25))
for ax, frac in zip(axes, [0, 1 / 12, 1 / 6, 1 / 3, 1 / 2]):
    t = frac * T
    ax.axhline(0, color=MUTED, lw=1)
    ax.plot(100 * x, triangle(1 / 3, x), color=MUTED, lw=1, alpha=0.5)
    ax.plot(100 * x, mode_sum(1 / 3, x, t), color=ACCENT)
    xk = kinks(1 / 3, t)
    ax.plot(100 * xk, dalembert(1 / 3, xk, t), "o", color=INK, ms=6)
    ax.text(0.99, 0.95, f"t = {1e3 * t:.2f} ms", transform=ax.transAxes, ha="right", va="top", color=INK)
    ax.set(ylim=(-1.25, 1.25), yticks=[-1, 0, 1])
axes[-1].set(xlim=(0, 100 * L), xlabel="x / cm from the bridge")
fig.supylabel("y / h", fontsize=11)
plt.show()

fig, axes = plt.subplots(2, 2, figsize=(7.5, 4), sharey="row", gridspec_kw=dict(height_ratios=[1, 1.3], hspace=0.15, wspace=0.08))
for col, a in enumerate([1 / 2, 1 / 3]):
    draw_scene(axes[0, col], axes[1, col], a)
axes[1, 0].set_ylabel("loudness / fundamental")
plt.show()

fig, (ax_s, ax_b) = plt.subplots(2, 1, figsize=(7, 3.6), gridspec_kw=dict(height_ratios=[1, 1.3], hspace=0.15))
draw_scene(ax_s, ax_b, 1 / 10)
ax_b.set_ylabel("loudness / fundamental")
plt.show()

a_grid = np.linspace(0.05, 0.5, 400)
fund = np.array([energy_shares(a, 1)[0] for a in a_grid])
high = np.array([1 - energy_shares(a, 4).sum() for a in a_grid])
fig, ax = plt.subplots(figsize=(7, 3.4))
ax.plot(a_grid, 100 * fund, color=ACCENT)
ax.plot(a_grid, 100 * high, color=SECOND)
for a in PLUCKS:
    f1, hi = 100 * energy_shares(a, 1)[0], 100 * (1 - energy_shares(a, 4).sum())
    ax.plot([a, a], [f1, hi], "o", color=INK, ms=6)
    if hi > f1:                                                      # at 1/10 the curves have crossed: one label below, one above right
        ax.annotate(f"{f1:.0f} %", (a, f1), xytext=(0, -10), textcoords="offset points", color=ACCENT, ha="center", va="top")
        ax.annotate(f"{hi:.0f} %", (a, hi), xytext=(8, 4), textcoords="offset points", color=SECOND, ha="left", va="bottom")
    else:                                                            # the rising red curve: its label above left, the flat teal one's below
        ax.annotate(f"{f1:.0f} %", (a, f1), xytext=(-6, 6), textcoords="offset points", color=ACCENT, ha="right", va="bottom")
        ax.annotate(f"{hi:.0f} %", (a, hi), xytext=(0, -8), textcoords="offset points", color=SECOND, ha="center", va="top")
ax.text(0.22, 52, "fundamental", color=ACCENT, ha="left")
ax.text(0.22, 20, "above the 4th harmonic", color=SECOND, ha="left")
ax.set(xlim=(0, 0.52), ylim=(0, 90), xticks=[0, 0.1, 1 / 3, 0.5], xticklabels=["0", "1/10", "1/3", "1/2"],
       xlabel="pluck point a (distance from the bridge / length)", ylabel="share of the energy / %")
plt.show()

fig, axes = plt.subplots(3, 1, figsize=(7, 5.4), sharex=True, sharey=True, gridspec_kw=dict(hspace=0.2))
worst = 0.0
for ax, a in zip(axes, PLUCKS):
    y = pluck_sound(a, damped=False, seconds=1.0)
    power = np.abs(np.fft.rfft(y)) ** 2
    f = np.fft.rfftfreq(len(y), 1 / FS)
    rel = 10 * np.log10(np.maximum(power / power[110], 1e-30))
    shares = energy_shares(a)
    present = shares > 1e-12 * shares[0]
    formula = 10 * np.log10(shares[present] / shares[0])
    worst = max(worst, np.abs(rel[110 * n[present]] - formula).max())
    ax.vlines(110 * n[~present & (n <= 20)], -40, 0, color=MUTED, ls="--", lw=1)   # the missing harmonics, under the label row
    shown = f <= 2300
    ax.plot(f[shown], np.maximum(rel[shown], -60), color=ACCENT, lw=1.4)
    keep = n[present] <= 20
    ax.plot(110 * n[present][keep], formula[keep], "o", mfc="none", color=INK, ms=6)
    ax.text(0.99, 0.98, f"a = 1/{1 / a:.0f}", transform=ax.transAxes, ha="right", va="top", color=INK)
    ax.set(ylim=(-40, 10), yticks=[-40, -20, 0])
axes[-1].set(xlim=(0, 2300), xlabel="f / Hz")
fig.supylabel("power / dB rel. fundamental", fontsize=11)
plt.show()
print(f"largest difference between the peaks and the formula: {worst:.0e} dB")

# The page around the game: the site's light and dark colors, the scene, the readouts, and the controls.
# game.js fills in every number; the notebook fills in DATA and SCRIPT.
TEMPLATE = """<!doctype html>
<html lang="en">
<head>
<meta charset="utf-8">
<meta name="viewport" content="width=device-width, initial-scale=1">
<title>Pluck a guitar string</title>
<style>
:root { --bg: #ffffff; --fg: #1f2a44; --ink: #1f2a44; --accent: #c8553d; --second: #2a7f9e; --muted: #8a8f98; }
@media (prefers-color-scheme: dark) {
  :root { --bg: #141a2a; --fg: #e9ecf3; --ink: #c9d0e3; --accent: #e3826b; --second: #6ec1de; --muted: #8e97ab; }
}
* { box-sizing: border-box; }
html, body { margin: 0; height: 100%; overflow: hidden; background: var(--bg); color: var(--fg);
  font: 15px/1.3 system-ui, sans-serif; -webkit-text-size-adjust: 100%; }
body { display: grid; gap: 4px; padding: 4px 6px 6px; grid-template-columns: minmax(0, 1fr);
  grid-template-rows: minmax(0, 1fr) auto auto; grid-template-areas: "scene" "readouts" "controls"; }
#view { grid-area: scene; position: relative; min-height: 0; }
#scene { position: absolute; inset: 0; width: 100%; height: 100%; touch-action: none; }
#readouts { grid-area: readouts; display: flex; flex-wrap: wrap; align-items: baseline; gap: 0 12px;
  font-variant-numeric: tabular-nums; }
#share { color: var(--accent); font-weight: 600; font-size: 18px; }
#where, #pitch { color: var(--muted); font-size: 14px; }
#controls { grid-area: controls; display: flex; flex-wrap: wrap; align-items: center; gap: 6px; }
.label { color: var(--muted); }
.seg { display: flex; } .seg button + button { margin-left: -1.5px; }
button { font: inherit; color: var(--second); background: none; border: 1.5px solid var(--second); border-radius: 4px;
  padding: 0.2em 0.7em; touch-action: manipulation; }
button[aria-pressed="true"], #play { color: var(--bg); background: var(--second); }
@media (pointer: coarse) { button { min-height: 44px; min-width: 44px; } }
@media (min-width: 600px) { body { gap: 6px 16px; padding: 10px; } #share { font-size: 22px; } }
@media (min-width: 600px) and (max-height: 399.98px) {   /* a phone held sideways: controls in a panel beside */
  body { grid-template-columns: minmax(0, 1fr) auto; grid-template-rows: auto minmax(0, 1fr);
    grid-template-areas: "scene readouts" "scene controls"; }
  #readouts { display: grid; align-content: start; max-width: 15em; }
  #controls { align-content: start; max-width: 15em; }
}
</style>
</head>
<body>
<div id="view"><canvas id="scene" aria-label="A guitar string between bridge and nut, with bars for the loudness of its harmonics"></canvas></div>
<div id="readouts"><span id="share"></span><span id="where"></span><span id="pitch"></span></div>
<div id="controls">
  <span class="label">pluck at, from the bridge:</span>
  <div class="seg" role="group" aria-label="Pluck at a fraction of the length from the bridge">
    <button id="half" title="1">1/2</button><button id="third" title="2">1/3</button><button id="tenth" title="3">1/10</button></div>
  <button id="play" title="Space"></button><button id="reset" title="R">Reset</button><button id="sound" title="S"></button>
</div>
<script type="application/json" id="data">{{DATA}}</script>
<script>{{SCRIPT}}</script>
</body>
</html>
"""