1% Driver for Chapter 1: an exact string, and the sound of sampling it.
2%
3% The editable script supplies
4%
5% y = pluck(t, f0, beta, where, xp, fmax, T60)
6%
7% the modal sum for a plucked string, restricted to partials below fmax.
8%
9% The driver's own ground truth is d'Alembert's solution, which shares no
10% algebra with the modal sum: the string's shape is the average of two copies
11% of the initial triangle travelling in opposite directions, so the bridge
12% force is a square-cornered pulse train and the pickup displacement is
13% piecewise linear in time. Both are evaluated exactly, in closed form, at
14% any time, and neither is band-limited.
15%
16% mode 'pluck' one note: the first three periods drawn finely, the harmonic
17% amplitudes, and the sound with a listening loss applied.
18% mode 'alias' the same note sampled two ways: d'Alembert's signal read off
19% at the sampling instants (every partial, most of them above
20% the Nyquist frequency and folded back below it), and the
21% modal sum stopped at the Nyquist frequency. Their difference
22% is exactly the folded part, and its energy is also known in
23% closed form, which the driver prints alongside the measurement.
24p = wv_params();
25fs = p.fs;
26f0 = p.f0;
27beta = p.beta;
28out = struct();
29out.fs = fs;
30out.f0 = f0;
31out.nyq = floor(fs / 2 / f0);
33if strcmp(p.mode, 'pluck')
34 % three periods, finely, for the waveform panel
35 tf = linspace(0, 3 / f0, 900)';
36 ts = (0:floor(3 * fs / f0))' / fs;
37 out.tf = tf' * 1000;
38 out.exact = drv_dalembert(tf, f0, beta, p.where, p.xp)';
39 out.ts = ts' * 1000;
40 out.synth = pluck(ts, f0, beta, p.where, p.xp, fs / 2, Inf)';
42 % harmonic amplitudes, from the exact coefficients
43 nh = min(out.nyq, 60);
44 n = 1:nh;
45 c = drv_coef(n, beta, p.where, p.xp);
46 out.n = n;
47 out.cdb = drv_db(c, max(abs(c)));
48 out.missing = sum(abs(sin(n * pi * beta)) < 1e-9);
50 % the sound: band-limited, with the listening loss
51 t = (0:round(p.T * fs) - 1)' / fs;
52 y = pluck(t, f0, beta, p.where, p.xp, fs / 2, p.T60);
53 [fr, db] = wv_spectrum(y(1:min(end, round(0.5 * fs))), fs, p.fmax, 700);
54 out.fr = fr;
55 out.db = db;
56 out.clip = wv_audio(y, fs);
58 % how far the band-limited sum is from d'Alembert at the sample instants
59 yb = pluck(t(1:round(0.25 * fs)), f0, beta, p.where, p.xp, fs / 2, Inf);
60 ye = drv_dalembert(t(1:round(0.25 * fs)), f0, beta, p.where, p.xp);
61 out.reldiff = sqrt(sum((yb - ye) .^ 2) / sum(yb .^ 2));
62else
63 t = (0:round(p.T * fs) - 1)' / fs;
64 env = exp(-6.908 * t / p.T60); % the same envelope on both clips
65 ye = drv_dalembert(t, f0, beta, p.where, p.xp);
66 yb = pluck(t, f0, beta, p.where, p.xp, fs / 2, Inf);
67 d = ye - yb;
68 out.measured = 10 * log10(sum(d .^ 2) / sum(yb .^ 2));
69 out.predicted = drv_predicted(f0, fs, beta, p.where, p.xp);
70 nw = round(0.5 * fs);
71 [fr, dbe] = wv_spectrum(ye(1:nw), fs, fs / 2, 900);
72 [~, dbb] = wv_spectrum(yb(1:nw), fs, fs / 2, 900);
73 out.fr = fr;
74 out.dbnaive = dbe;
75 out.dbband = dbb;
76 % where the first few folded partials land
77 nmax = 3 * out.nyq;
78 nn = out.nyq + 1:nmax;
79 fa = mod(nn * f0, fs);
80 fa(fa > fs / 2) = fs - fa(fa > fs / 2);
81 cc = drv_coef(nn, beta, p.where, p.xp);
82 keep = abs(cc) > 1e-9 * max(abs(drv_coef(1:out.nyq, beta, p.where, p.xp)));
83 out.falias = fa(keep);
84 out.harm = (1:out.nyq) * f0;
85 sc = max(abs(ye));
86 out.clipNaive = wv_audio(ye .* env, fs, sc);
87 out.clipBand = wv_audio(yb .* env, fs, sc);
88 out.clipDiff = wv_audio(d .* env, fs, sc);
89 out.ratio = fs / f0;
90end
92wv_write(out);
94function y = drv_dalembert(t, f0, beta, where, xp)
95% d'Alembert: u(x, t) = (phi(x - ct) + phi(x + ct)) / 2, with phi the odd,
96% period-2 extension of the initial triangle and c = 2 f0 on a string of
97% length 1.
98s = 2 * f0 * t;
99if strcmp(where, 'bridge')
100 % u_x(0, t) = phi'(ct), and phi' is even: 1/beta on (0, beta),
101 % -1/(1 - beta) on (beta, 2 - beta), 1/beta again up to 2. A sample that
102 % lands exactly on a jump takes the midpoint, as a Fourier series would.
103 r = mod(s, 2);
104 hi = 1 / beta;
105 lo = -1 / (1 - beta);
106 y = lo * ones(size(r));
107 y(r < beta | r > 2 - beta) = hi;
108 y(abs(r - beta) < 1e-12 | abs(r - (2 - beta)) < 1e-12) = 0.5 * (hi + lo);
109else
110 y = 0.5 * (drv_phi(xp - s, beta) + drv_phi(xp + s, beta));
111end
112end
114function v = drv_phi(x, beta)
115% the initial triangle, extended oddly about 0 and periodically with period 2
116r = mod(x, 2);
117sg = ones(size(r));
118sg(r > 1) = -1;
119r(r > 1) = 2 - r(r > 1);
120v = min(r / beta, (1 - r) / (1 - beta));
121v = sg .* v;
122end
124function c = drv_coef(n, beta, where, xp)
125A = 2 * sin(n * pi * beta) ./ (n .^ 2 * pi ^ 2 * beta * (1 - beta));
126if strcmp(where, 'bridge')
127 c = n * pi .* A;
128else
129 c = A .* sin(n * pi * xp);
130end
131end
133function d = drv_db(c, ref)
134d = 20 * log10(max(abs(c) / ref, 1e-8));
135end
137function r = drv_predicted(f0, fs, beta, where, xp)
138% Energy above the Nyquist frequency over energy below it, from the exact
139% coefficients. For the bridge force the total is known in closed form by
140% Parseval, mean(F^2) = 1 / (beta (1 - beta)), so the tail is the total less
141% the head; for the pickup the coefficients fall like 1/n^2 and the tail is
142% summed directly to n = 2e5, beyond which it is below rounding.
143K = floor(fs / 2 / f0);
144head = drv_coef(1:K, beta, where, xp);
145eh = sum(head .^ 2) / 2;
146if strcmp(where, 'bridge')
147 et = 1 / (beta * (1 - beta)) - eh;
148else
149 tail = drv_coef(K + 1:200000, beta, where, xp);
150 et = sum(tail .^ 2) / 2;
151end
152r = 10 * log10(et / eh);
153end