Module 1 · Weeks 2–3 · Undergraduate → Graduate

Acoustics & Fourier Analysis

Sound is vibration; vibration is governed by differential equations; and differential equations are solved by decomposing functions into sinusoidal components. This is the Fourier programme, and it connects the physics of a vibrating string directly to the mathematics of Hilbert spaces.

1. The Vibrating String

Consider a string of length \(L\), fixed at both ends, with linear mass density \(\mu\) under tension \(T\). The transverse displacement \(u(x,t)\) satisfies the one-dimensional wave equation:

\[\frac{\partial^2 u}{\partial t^2} = c^2 \frac{\partial^2 u}{\partial x^2}, \qquad c = \sqrt{\frac{T}{\mu}}\]

with boundary conditions \(u(0,t) = u(L,t) = 0\) (fixed ends) and initial conditions \(u(x,0) = f(x)\), \(\dot{u}(x,0) = g(x)\).

D'Alembert's Solution

On the infinite line, the general solution takes the elegant form:

\[u(x,t) = F(x - ct) + G(x + ct)\]

This reveals that solutions are superpositions of left-moving and right-moving waves travelling at speed \(c\). On the finite string with fixed ends, reflections create standing waves — which leads us to the eigenvalue problem.

2. The Eigenvalue Problem

Separation of variables \(u(x,t) = X(x)T(t)\) yields the Sturm–Liouville eigenvalue problem:

\[X''(x) + \lambda X(x) = 0, \qquad X(0) = X(L) = 0\]

\[\text{Eigenvalues: } \lambda_n = \left(\frac{n\pi}{L}\right)^2, \qquad n = 1, 2, 3, \ldots\]

\[\text{Eigenfunctions: } X_n(x) = \sin\left(\frac{n\pi x}{L}\right)\]

Normal Mode Frequencies

Each eigenfunction oscillates at its own natural frequency:

\[f_n = \frac{n}{2L}\sqrt{\frac{T}{\mu}} = n \cdot f_1\]

The frequencies form an arithmetic progression: \(f_1, 2f_1, 3f_1, \ldots\)This is precisely the harmonic series from Module 0. The eigenvalue problem thus provides the mathematical explanation for why vibrating strings produce harmonics.

Physical interpretation

The \(n\)-th mode has \(n-1\) nodes (stationary points) along the string. The fundamental (\(n=1\)) has no internal nodes; the second harmonic (\(n=2\)) has one node at the midpoint, etc. A real string vibrates as a superposition of all these modes, with amplitudes determined by the initial pluck shape.

3. Fourier Series

The normal modes \(\{\sin(n\pi x/L)\}\) form an orthogonal basis for \(L^2([0,L])\). Any square-integrable function on \([0,L]\)can be expanded:

\[f(x) = \sum_{n=1}^{\infty} b_n \sin\left(\frac{n\pi x}{L}\right)\]

\[b_n = \frac{2}{L}\int_0^L f(x)\sin\left(\frac{n\pi x}{L}\right)\,dx\]

More generally, Fourier series on \([-\pi, \pi]\) use both sines and cosines:

\[f(x) = \frac{a_0}{2} + \sum_{n=1}^{\infty}\left[a_n \cos(nx) + b_n \sin(nx)\right]\]

Parseval's Identity

Energy is conserved in the Fourier domain:

\[\frac{1}{\pi}\int_{-\pi}^{\pi}|f(x)|^2\,dx = \frac{a_0^2}{2} + \sum_{n=1}^{\infty}(a_n^2 + b_n^2)\]

In acoustic terms: the total power of a sound equals the sum of the powers of its individual harmonics. This is why the amplitudes \(a_n, b_n\) (or equivalently, the power spectrum \(|c_n|^2\)) completely characterize the timbre.

Gibbs Phenomenon

At jump discontinuities, partial Fourier sums overshoot by approximately 9% of the jump — no matter how many terms are included. This Gibbs phenomenon is audible as a “ringing” artifact when synthesizing square waves with finitely many harmonics. The overshoot converges to \(\frac{1}{\pi}\int_0^{\pi}\frac{\sin t}{t}\,dt - \frac{1}{2} \approx 0.0895\), or about 8.95% of the jump height.

Simulation: Fourier Series Convergence & Gibbs Phenomenon

The simulation below demonstrates square-wave and sawtooth-wave partial-sum convergence, the Gibbs overshoot at discontinuities, Parseval energy convergence, and the power-spectrum envelope for a square wave.

Python

script.py151 lines

import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec

plt.rcParams['font.family'] = 'monospace'
BG = '#0a0a1a'
COLORS = ['#60a5fa', '#34d399', '#f59e0b', '#a78bfa', '#ef4444']

fig = plt.figure(figsize=(14, 12), facecolor=BG)
fig.suptitle('Fourier Series: Convergence & Gibbs Phenomenon', color='white', fontsize=16, fontweight='bold', y=0.98)

gs = gridspec.GridSpec(2, 2, figure=fig, hspace=0.42, wspace=0.35)

x = np.linspace(-np.pi, np.pi, 2000, endpoint=False)

# ── Plot 1: Square wave convergence ──────────────────────────────────────
ax1 = fig.add_subplot(gs[0, 0], facecolor=BG)

harmonics_sq = [1, 3, 5, 10, 50]
# True square wave
sq_true = np.sign(np.sin(x))
sq_true[sq_true == 0] = 1.0
ax1.plot(x, sq_true, color='#1e3a5f', lw=1.5, ls='--', label='True square', zorder=1)

for idx, N in enumerate(harmonics_sq):
    y = np.zeros_like(x)
    for k in range(1, N + 1, 2):     # odd harmonics only
        y += (4 / (np.pi * k)) * np.sin(k * x)
    ax1.plot(x, y, color=COLORS[idx], lw=1.2 if N < 50 else 1.8,
             label=f'N={N}', alpha=0.85, zorder=2 + idx)

# Annotate Gibbs overshoot on N=50 curve
y50 = np.zeros_like(x)
for k in range(1, 51, 2):
    y50 += (4 / (np.pi * k)) * np.sin(k * x)
peak_idx = np.argmax(y50[:1000])
ax1.annotate('~9% overshoot', xy=(x[peak_idx], y50[peak_idx]),
             xytext=(x[peak_idx] + 0.4, y50[peak_idx] + 0.08),
             arrowprops=dict(arrowstyle='->', color='#ef4444', lw=1),
             color='#ef4444', fontsize=8)

ax1.set_xlim(-np.pi, np.pi)
ax1.set_ylim(-1.4, 1.6)
ax1.set_title('Square Wave: Partial Sums', color='white', fontsize=10)
ax1.set_xlabel('x', color='#94a3b8', fontsize=9)
ax1.legend(fontsize=7, facecolor='#0f172a', labelcolor='white', edgecolor='#1e3a5f', ncol=2)
ax1.set_facecolor(BG)
ax1.tick_params(colors='#94a3b8', labelsize=8)
for sp in ax1.spines.values():
    sp.set_color('#1e3a5f')

gibbs_overshoot = (y50.max() - 1.0) / 2.0 * 100
print(f'Gibbs overshoot (N=50 square wave): {gibbs_overshoot:.2f}% of jump height')

# ── Plot 2: Sawtooth wave convergence ────────────────────────────────────
ax2 = fig.add_subplot(gs[0, 1], facecolor=BG)

harmonics_saw = [1, 5, 20, 100]
saw_true = x / np.pi    # true sawtooth on [-pi, pi]
ax2.plot(x, saw_true, color='#1e3a5f', lw=1.5, ls='--', label='True sawtooth', zorder=1)

for idx, N in enumerate(harmonics_saw):
    y = np.zeros_like(x)
    for k in range(1, N + 1):
        y += (2 * (-1)**(k + 1) / (np.pi * k)) * np.sin(k * x)
    ax2.plot(x, y, color=COLORS[idx], lw=1.2 if N < 100 else 1.8,
             label=f'N={N}', alpha=0.85, zorder=2 + idx)

ax2.set_xlim(-np.pi, np.pi)
ax2.set_ylim(-1.35, 1.35)
ax2.set_title('Sawtooth Wave: Partial Sums', color='white', fontsize=10)
ax2.set_xlabel('x', color='#94a3b8', fontsize=9)
ax2.legend(fontsize=7, facecolor='#0f172a', labelcolor='white', edgecolor='#1e3a5f', ncol=2)
ax2.set_facecolor(BG)
ax2.tick_params(colors='#94a3b8', labelsize=8)
for sp in ax2.spines.values():
    sp.set_color('#1e3a5f')

print('Sawtooth wave plotted with N =', harmonics_saw)

# ── Plot 3: Parseval convergence ─────────────────────────────────────────
ax3 = fig.add_subplot(gs[1, 0], facecolor=BG)

# Square wave: ||f||^2 = (1/pi) * integral of 1 dx over [-pi,pi] = 2/pi * pi = 2...
# exact: (1/pi) int_{-pi}^{pi} (sq)^2 dx = (1/pi)*2*pi = 2, but convention varies.
# Using L2 norm: (1/(2*pi)) * integral |f|^2 dx = 1 for unit square wave.
# Parseval: sum |b_n|^2 = ||f||^2, where b_n = 4/(pi*n) for odd n.
# True value: sum_{k odd} (4/(pi*k))^2 / 2 = (8/pi^2) * (pi^2/8) = 1. Correct.
true_norm_sq = 1.0   # (1/2pi)*int |sq|^2 dx for unit amplitude = 1

N_max = 60
partial_parseval = []
running = 0.0
for k in range(1, 2 * N_max + 1, 2):
    b_k = 4 / (np.pi * k)
    running += 0.5 * b_k**2
    partial_parseval.append(running)

ns = list(range(1, len(partial_parseval) + 1))
ax3.axhline(true_norm_sq, color='#ef4444', lw=1.5, ls='--', label='||f||^2 = 1.0')
ax3.plot(ns, partial_parseval, color='#60a5fa', lw=2, label='Partial sum of |b_n|^2')
ax3.fill_between(ns, partial_parseval, true_norm_sq, alpha=0.15, color='#ef4444')

ax3.set_xlabel('Number of odd harmonics included', color='#94a3b8', fontsize=9)
ax3.set_ylabel('Cumulative energy', color='#94a3b8', fontsize=9)
ax3.set_title("Parseval's Theorem: Energy Convergence", color='white', fontsize=10)
ax3.legend(fontsize=8, facecolor='#0f172a', labelcolor='white', edgecolor='#1e3a5f')
ax3.set_facecolor(BG)
ax3.tick_params(colors='#94a3b8', labelsize=8)
for sp in ax3.spines.values():
    sp.set_color('#1e3a5f')

pct_at_5 = partial_parseval[2] / true_norm_sq * 100
pct_at_10 = partial_parseval[4] / true_norm_sq * 100
print(f'Parseval convergence (square wave):')
print(f'  After 5 harmonics:  {pct_at_5:.1f}% of total energy')
print(f'  After 11 harmonics: {pct_at_10:.1f}% of total energy')
print(f'  After {N_max} harmonics: {partial_parseval[-1]/true_norm_sq*100:.2f}% of total energy')

# ── Plot 4: Power spectrum of square wave ────────────────────────────────
ax4 = fig.add_subplot(gs[1, 1], facecolor=BG)

odd_ns = np.arange(1, 52, 2)
power = (4 / (np.pi * odd_ns))**2

ax4.bar(odd_ns, power, width=0.8, color='#60a5fa', alpha=0.8, label='|b_n|^2 (square wave)')

# 1/n^2 envelope
n_env = np.linspace(1, 51, 200)
env = (4 / np.pi)**2 / n_env**2
ax4.plot(n_env, env, color='#f59e0b', lw=2, ls='--', label='Envelope: (4/pi)^2 / n^2')

ax4.set_xlabel('Harmonic n (odd only)', color='#94a3b8', fontsize=9)
ax4.set_ylabel('|b_n|^2', color='#94a3b8', fontsize=9)
ax4.set_title('Power Spectrum: Square Wave', color='white', fontsize=10)
ax4.legend(fontsize=8, facecolor='#0f172a', labelcolor='white', edgecolor='#1e3a5f')
ax4.set_facecolor(BG)
ax4.tick_params(colors='#94a3b8', labelsize=8)
ax4.set_xlim(0, 52)
for sp in ax4.spines.values():
    sp.set_color('#1e3a5f')

print('Power spectrum |b_n|^2 for n = 1,3,5,...,51:')
for n, p in zip(odd_ns[:6], power[:6]):
    print(f'  n={n:2d}: |b_n|^2 = {p:.5f}')

plt.savefig('output.png', dpi=120, bbox_inches='tight', facecolor=BG)
print('Saved output.png')

Click Run to execute the Python code

Code will be executed with Python 3 on the server

4. Fourier Transform & the Spectrogram

For non-periodic signals, the Fourier series generalizes to the Fourier transform:

\[\hat{f}(\xi) = \int_{-\infty}^{\infty} f(t)\,e^{-2\pi i \xi t}\,dt\]

\[f(t) = \int_{-\infty}^{\infty} \hat{f}(\xi)\,e^{2\pi i \xi t}\,d\xi\]

Short-Time Fourier Transform (STFT)

Music changes in time, so we need time-frequency analysis. The STFT applies the Fourier transform through a sliding window function \(w(t)\):

\[\text{STFT}\{f\}(\tau, \xi) = \int_{-\infty}^{\infty} f(t)\,w(t-\tau)\,e^{-2\pi i \xi t}\,dt\]

The squared modulus \(|\text{STFT}(\tau,\xi)|^2\) is the spectrogram — the standard visual representation of music.

Heisenberg–Gabor Uncertainty

\[\Delta t \cdot \Delta \xi \geq \frac{1}{4\pi}\]

You cannot simultaneously localize a signal precisely in time and in frequency. A short window gives good time resolution but poor frequency resolution, and vice versa. This is the time-frequency uncertainty principle, the signal-processing analogue of Heisenberg's uncertainty in quantum mechanics. The Gaussian window (Gabor atom) achieves the theoretical minimum product.

5. Interactive Fourier Synthesizer

Build a timbre from its harmonic components. Each slider controls the amplitude of the\(n\)-th harmonic (\(n \cdot f_0\) Hz). The waveform is \(y(t) = \sum_{n=1}^{8} a_n \sin(2\pi n f_0 t)\). Try the presets to hear how different amplitude profiles create familiar timbres.

Interactive Fourier Synthesizer

a₁ · 1f₀1.00

a₂ · 2f₀0.00

a₃ · 3f₀0.00

a₄ · 4f₀0.00

a₅ · 5f₀0.00

a₆ · 6f₀0.00

a₇ · 7f₀0.00

a₈ · 8f₀0.00

Adjust harmonic amplitudes and click Play. The waveform updates in real time via the Web Audio API analyser node.

Why these presets?

Square wave: only odd harmonics, amplitudes \(1/n\)
Sawtooth: all harmonics, amplitudes \(1/n\)
Triangle: odd harmonics only, amplitudes \(1/n^2\)
Clarinet: predominantly odd harmonics (cylindrical bore resonance)
Cello: rich spectrum with all harmonics, slowly decaying

6. Psychoacoustics

The ear is not a perfect Fourier analyser. The cochlea performs a mechanical frequency decomposition via the basilar membrane, but with finite resolution determined by critical bands.

Place Theory

Different frequencies excite different positions along the basilar membrane. High frequencies resonate near the base (stiff end); low frequencies near the apex (flexible end). Each position responds to a band of frequencies approximately 1/3 octave wide — these are the critical bands (Bark scale).

Helmholtz Consonance Theory

Hermann von Helmholtz (1863) proposed that dissonance arises when partials of two tones fall within the same critical band, producing audible beating. Consonant intervals have partials that either coincide or are well-separated on the basilar membrane.

The Plomp–Levelt curve

Plomp and Levelt (1965) measured the perceived dissonance of two pure tones as a function of their frequency separation. Maximum dissonance occurs at about 25% of the critical bandwidth; negligible dissonance occurs beyond one critical bandwidth. For complex tones, total dissonance is the sum over all partial pairs — reproducing the traditional ranking of intervals.

Modern theories (Terhardt, Parncutt) add virtual pitch: the auditory system infers a missing fundamental from its harmonics, explaining why a phone speaker with no bass output still conveys the pitch of a low cello note. This pattern-matching mechanism is distinct from place coding and involves central neural processing.

The Math Behind Music and Sound Synthesis

← Module 0: Prologue Module 2: Tuning →