"""TopEdu project 7: Why a plucked note changes as it fades A string of length L vibrates as u_tt + a u_t = c^2 u_xx, fixed at both ends. The term a u_t is a simple drag. In this model every harmonic fades as e^(-a t / 2), all at the same rate. A real guitar loses its high harmonics sooner, so the sound gets mellower as it fades. Measuring by how much is the research. Stage 1: run this file. It plucks a low E string (82.4 Hz), steps it for four periods, and compares the whole string with the exact answer. Then it measures each harmonic's fade from the simulated sound. The CHECK line should say PASS. Stage 2: record real strings, measure each harmonic's fade the same way, and show where one damping rate fails. Run: python 07_damped_string.py [--show] or python 07_damped_string.py my_pluck.wav Needs: numpy and matplotlib (reading .wav uses Python's own wave module) """ import os import sys import wave import numpy as np import matplotlib if "--show" not in sys.argv: matplotlib.use("Agg") import matplotlib.pyplot as plt OUT = os.path.join(os.path.dirname(os.path.abspath(__file__)), "output") os.makedirs(OUT, exist_ok=True) L = 0.65 # string length, m (a common guitar scale) F = 82.4 # lowest note, Hz (low E) A = 1.0 # damping, 1/s C = 2 * L * F # wave speed that gives that note N = 129 def exact(x, t, modes=((1, 1.0), (2, 0.3))): """Each mode n fades as e^(-a t/2) and swings at sqrt(w_n^2 - a^2/4), with w_n = n pi c / L, starting at rest.""" total = 0.0 for n, size in modes: w = n * np.pi * C / L omega = np.sqrt(w ** 2 - A ** 2 / 4) total = total + size * np.exp(-A * t / 2) * (np.cos(omega * t) + A / (2 * omega) * np.sin(omega * t)) * np.sin(n * np.pi * x / L) return total def simulate(final, n=N, pickup=0.12, start=None): """Centred steps of the damped wave equation. Returns x, the string at the end, the time step and the motion at the pickup point (a fraction of the length from the bridge) at every step.""" x = np.linspace(0, L, n) h = x[1] - x[0] dt = 0.5 * h / C steps = int(round(final / dt)); dt = final / steps u_old = exact(x, 0.0) if start is None else start(x) r2 = (C * dt / h) ** 2 lap = lambda u: np.concatenate(([0.0], u[2:] - 2 * u[1:-1] + u[:-2], [0.0])) u = u_old + 0.5 * r2 * lap(u_old) * (1 - A * dt / 2) # first step from rest probe = int(round(pickup * (n - 1))) heard = [u_old[probe], u[probe]] for _ in range(steps - 1): # (u_new - 2u + u_old)/dt^2 + a (u_new - u_old)/(2 dt) = c^2 u_xx, solved for u_new u, u_old = (2 * u - (1 - A * dt / 2) * u_old + r2 * lap(u)) / (1 + A * dt / 2), u heard.append(u[probe]) return x, u, dt, np.array(heard) def harmonic_decay(signal, dt, f0, harmonics=5, window=0.1): """Fade rate of each harmonic: its strength in short windows, then the slope of log(strength) over time. Returns a list of (harmonic number, decay rate in 1/s); a rate of a means strength ~ e^(-a t).""" size = int(window / dt) starts = range(0, len(signal) - size, size // 2) times = np.array([(s + size / 2) * dt for s in starts]) t_local = np.arange(size) * dt taper = np.hanning(size) tracks = [] for k in range(1, harmonics + 1): probe = np.exp(-2j * np.pi * k * f0 * t_local) * taper tracks.append(np.array([abs(np.dot(signal[s:s + size], probe)) for s in starts])) loudest = max(track[0] for track in tracks) rates = [] for k, strength in enumerate(tracks, start=1): if strength[0] < 0.01 * loudest: continue # this harmonic is barely there: its rate would be noise good = strength > 1e-3 * strength[0] # stop once it has faded into the noise if good.sum() > 3: rates.append((k, -np.polyfit(times[good], np.log(strength[good]), 1)[0])) return rates def read_wav(path): with wave.open(path) as handle: rate, channels = handle.getframerate(), handle.getnchannels() data = np.frombuffer(handle.readframes(handle.getnframes()), dtype=np.int16).astype(float) return 1.0 / rate, data.reshape(-1, channels).mean(axis=1) def main(): final = 4 / F x, u, _, _ = simulate(final) worst = np.abs(u - exact(x, final)).max() print(f"after four periods: largest error along the string {worst:.2e} (pluck size 1.3)") print(f"CHECK string against the exact damped answer -> {'PASS' if worst < 0.013 else 'FAIL'} (limit 1% of the pluck)") # A pluck at one point excites many harmonics. Listen for 2 s and measure how each fades. pluck = lambda x: np.where(x < 0.2 * L, x / (0.2 * L), (L - x) / (0.8 * L)) * 0.01 _, _, dt, heard = simulate(2.0, start=pluck) print(f"each harmonic's fade in the model (strength ~ e^(-rate t)); the model says every rate is a/2 = {A / 2}") for k, rate in harmonic_decay(heard, dt, F): print(f" harmonic {k} ({k * F:6.1f} Hz): rate {rate:.3f} per second") fig, ax = plt.subplots(figsize=(6, 3.4)) t = np.arange(len(heard)) * dt ax.plot(t, heard, lw=0.4) ax.set_xlabel("time (s)"); ax.set_ylabel("string motion at the pickup"); ax.set_title("Simulated pluck, 2 seconds") fig.tight_layout(); fig.savefig(os.path.join(OUT, "07_pluck.png"), dpi=150) for path in [arg for arg in sys.argv[1:] if arg.lower().endswith(".wav")]: rdt, samples = read_wav(path) start = int(np.argmax(np.abs(samples))) # measure from the loudest moment onward print(f"{os.path.basename(path)} (assuming the lowest note is {F} Hz; change F to your note):") for k, rate in harmonic_decay(samples[start:], rdt, F, harmonics=8): print(f" harmonic {k}: rate {rate:.2f} per second") print(f"plots saved in {OUT}") # ---- YOUR TURN ---- # 1. Record one plucked string on a phone (a quiet room, 3 seconds, saved as .wav). Set F to its note and run: # python 07_damped_string.py my_pluck.wav # 2. Plot rate against harmonic number. One a cannot fit it. Try rate = b0 + b2 * k^2. Does that fit better? # 3. Compare a wound and a plain string, or nylon and steel. Which loses its high harmonics fastest? # 4. Then fix the model: what term would make high harmonics fade faster? (Hint: a drag on u_xxt.) if "--show" in sys.argv: plt.show() if __name__ == "__main__": main()