1% Driver for Chapter 3: the digital waveguide and the Karplus-Strong string.
2%
3% The editable script supplies
4%
5% y = ks(x0, D, tuning, M)
6%
7% a delay loop of length D samples, started from x0.
8%
9% mode 'equiv' the lossless loop against the leapfrog scheme of Chapter 2 at
10% lam = 1, on grids from 25 to 3200 intervals. The leapfrog
11% scheme is run by the driver's own copy; its first 2N output
12% samples seed the loop, and the two are then compared for the
13% rest of the run, sample for sample, and timed.
14% mode 'keys' one note on each of 25 keys from C2 to C8, with each tuning.
15% The fundamental's frequency is measured by wv_peakfreq and
16% compared with the equal-tempered target; its decay time is
17% measured from the envelope. A two-octave scale is also played
18% with each tuning, and once by exact modal synthesis at the
19% target pitches, for the ear to compare.
20p = wv_params();
21fs = p.fs;
22out = struct();
24if strcmp(p.mode, 'equiv')
25 Ns = p.Ns(:)';
26 M = round(p.T * fs);
27 TL = NaN(1, numel(Ns));
28 TW = NaN(1, numel(Ns));
29 DIFF = NaN(1, numel(Ns));
30 for k = 1:numel(Ns)
31 N = Ns(k);
32 x = (0:N)' / N;
33 u0 = min(x / p.beta, (1 - x) / (1 - p.beta));
34 um = 0.5 * ([u0(2:end); 0] + [0; u0(1:end-1)]); % exact at lam = 1
35 um(1) = 0;
36 um(end) = 0;
37 drv_lf(u0, um, 50);
38 t0 = tic;
39 yl = drv_lf(u0, um, M);
40 TL(k) = toc(t0) * 1000;
41 seed = yl(1:2 * N);
42 ks(seed, 2 * N, 'lossless', 4 * N);
43 t0 = tic;
44 yw = ks(seed, 2 * N, 'lossless', M);
45 TW(k) = toc(t0) * 1000;
46 DIFF(k) = max(abs(yw - yl)) / max(abs(yl));
47 end
48 out.Ns = Ns;
49 out.tl = TL;
50 out.tw = TW;
51 out.diff = DIFF;
52 out.M = M;
53else
54 notes = p.notes(:)';
55 tunings = {'integer', 'linear', 'allpass'};
56 nt = numel(tunings);
57 Mn = round(p.Tn * fs);
58 CE = NaN(nt, numel(notes));
59 T60 = NaN(nt, numel(notes));
60 for i = 1:numel(notes)
61 f = wv_note(notes(i));
62 D = fs / f;
63 x0 = drv_seed(f, fs, p.beta);
64 for k = 1:nt
65 y = ks(x0, D, tunings{k}, Mn);
66 fm = wv_peakfreq(y, 1 / fs, f / 1.1, f * 1.1);
67 CE(k, i) = wv_cents(fm, f);
68 T60(k, i) = drv_t60(y, fm, fs);
69 end
70 end
71 out.notes = notes;
72 out.f = wv_note(notes);
73 out.tunings = tunings;
74 out.cents = wv_rows(CE);
75 out.t60 = wv_rows(T60);
77 % a scale, played with each tuning and exactly
78 sc = p.scale(:)';
79 step = round(p.step * fs);
80 ring = round(p.ring * fs);
81 Mt = step * (numel(sc) - 1) + ring;
82 Y = zeros(Mt, nt + 1);
83 for i = 1:numel(sc)
84 f = wv_note(sc(i));
85 x0 = drv_seed(f, fs, p.beta);
86 idx = (i - 1) * step + (1:ring);
87 for k = 1:nt
88 Y(idx, k) = Y(idx, k) + ks(x0, fs / f, tunings{k}, ring);
89 end
90 % the reference: the same pluck, modally, at exactly the target pitch,
91 % with each partial given the decay the integer loop gives it
92 t = (0:ring - 1)' / fs;
93 Y(idx, nt + 1) = Y(idx, nt + 1) + drv_modal(t, f, fs, p.beta);
94 end
95 s = max(abs(Y(:)));
96 out.clipInt = wv_audio(Y(:, 1), fs, s);
97 out.clipLin = wv_audio(Y(:, 2), fs, s);
98 out.clipAp = wv_audio(Y(:, 3), fs, s);
99 out.clipExact = wv_audio(Y(:, 4), fs, s);
100 out.scale = sc;
101end
103wv_write(out);
105function y = drv_lf(u, uold, M)
106% leapfrog at lam = 1 on a fixed string, recording the first node off the end
107n1 = numel(u);
108y = zeros(M, 1);
109unew = zeros(n1, 1);
110for n = 1:M
111 for j = 2:n1 - 1
112 unew(j) = u(j + 1) + u(j - 1) - uold(j);
113 end
114 uold = u;
115 u = unew;
116 y(n) = u(2);
117end
118end
120function x0 = drv_seed(f, fs, beta)
121% The loop's initial content, the same for every tuning: the band-limited
122% bridge force of a string plucked at beta, for a little over one period,
123% because the loop reads up to two samples further back than its length.
124P = ceil(fs / f) + 2;
125t = (0:P - 1)' / fs;
126x0 = zeros(P, 1);
127for n = 1:floor(fs / 2 / f)
128 x0 = x0 + 2 * sin(n * pi * beta) / (n * pi * beta * (1 - beta)) * cos(2 * pi * n * f * t);
129end
130% A loop with a two-point average in it passes a constant unchanged for
131% ever, so a seed with a nonzero mean leaves a DC offset that never decays.
132% It is inaudible, but it would sit under the fundamental's decay and spoil
133% the measurement of it, so the seed is made to have none.
134x0 = x0 - mean(x0);
135end
137function y = drv_modal(t, f, fs, beta)
138% the same pluck by modal synthesis, each partial decaying as the two-point
139% average makes it decay in a loop of exactly the right length
140y = zeros(size(t));
141for n = 1:floor(fs / 2 / f)
142 g = abs(cos(pi * n * f / fs));
143 s = -log(max(g, 1e-300)) * f; % amplitude lost per second
144 y = y + 2 * sin(n * pi * beta) / (n * pi * beta * (1 - beta)) * exp(-s * t) .* cos(2 * pi * n * f * t);
145end
146end
148function T = drv_t60(y, f, fs)
149% decay time of the partial at f: its amplitude in successive tapered frames,
150% a straight line fitted to the logarithm, extrapolated to -60 dB
151nw = round(8 * fs / f);
152nw = max(nw, 256);
153hop = nw;
154nfr = floor((numel(y) - nw) / hop);
155w = wv_window(nw);
156tt = (0:nw - 1)' / fs;
157e = exp(-2i * pi * f * tt);
158a = zeros(nfr, 1);
159tc = zeros(nfr, 1);
160for k = 1:nfr
161 seg = y((k - 1) * hop + (1:nw));
162 a(k) = abs(sum(w .* seg .* e));
163 tc(k) = ((k - 1) * hop + nw / 2) / fs;
164end
165% only the frames within 100 dB of the loudest, where the partial is still
166% far above rounding error and the window's leakage from anything else
167use = find(a > 1e-5 * max(a));
168use = use(use > 1);
169if numel(use) < 3
170 use = 2:min(4, nfr);
171end
172c = polyfit(tc(use), log(a(use)), 1);
173T = -log(1000) / c(1);
174if c(1) >= 0
175 T = Inf;
176end
177end