d8098edAn interactive tutorial on the wave equation in acoustics: strings, drums and rooms, measured by earJeremy Magland 1% Driver for Chapter 2: the leapfrog scheme on a string and in a tube.
2%
3% The editable script supplies
4%
5% [y, U] = leapfrog(u, uold, lam, M, bc, jr, K)
6%
7% which advances the scheme M steps from two given time levels.
8%
9% The driver supplies exact data at both starting levels, so that the only
10% error left is the scheme's, and it measures that error against solutions it
11% computes independently:
12%
13% mode 'run' a plucked string at the audio rate. The time step is one
14% sample, 1/fs, and the grid is N intervals, so the Courant number
15% is lam = 2 f0 N / fs and choosing the grid chooses lam. The
16% reference is d'Alembert's solution at the receiver.
17% mode 'cents' every partial's frequency, measured. The field is projected
18% onto the discrete sine vectors sin(n pi j / N), which are exact
19% eigenvectors of the three-point difference, so each projection
20% is a single sinusoid whose frequency wv_peakfreq reads off.
21% The exact frequency of partial n is n f0.
22% mode 'click' a short pulse circulating in a periodic tube, recorded at one
23% point, against the exact pulse train.
24p = wv_params();
25out = struct();
27if strcmp(p.mode, 'run')
28 fs = p.fs;
29 f0 = p.f0;
30 N = p.N;
31 c = 2 * f0; % a string of length 1
32 h = 1 / N;
33 dt = 1 / fs;
34 lam = c * dt / h;
35 M = round(p.T * fs);
36 x = (0:N)' * h;
37 u0 = drv_phi(x, p.beta);
38 um = 0.5 * (drv_phi(x + c * dt, p.beta) + drv_phi(x - c * dt, p.beta));
39 if strcmp(p.where, 'bridge')
40 jr = 2; % the first node off the bridge
41 else
42 jr = round(p.xp * N) + 1;
43 end
44 t0 = tic;
45 [y, ~] = leapfrog(u0, um, lam, M, 'fixed', jr, 0);
46 out.ms = toc(t0) * 1000;
47 t = (1:M)' * dt;
48 ye = 0.5 * (drv_phi(x(jr) - c * t, p.beta) + drv_phi(x(jr) + c * t, p.beta));
50 sc = max(abs(ye));
51 bad = find(~isfinite(y) | abs(y) > 50 * sc, 1);
52 out.unstable = ~isempty(bad);
53 if out.unstable
54 out.blowup = bad * dt * 1000;
55 y(bad:end) = 0;
56 else
57 out.blowup = NaN;
58 end
59 out.lam = lam;
60 out.N = N;
61 e = y - ye;
62 % relative rms error, period by period
63 P = fs / f0;
64 np = floor(M / P);
65 er = NaN(1, np);
66 for k = 1:np
67 i0 = floor((k - 1) * P) + 1;
68 i1 = floor(k * P);
69 er(k) = sqrt(sum(e(i0:i1) .^ 2) / sum(ye(i0:i1) .^ 2));
70 end
71 out.tp = ((1:np) - 0.5) / f0 * 1000;
72 out.err = er;
73 last = t > p.T - 1;
74 out.relerr = sqrt(sum(e(last) .^ 2) / sum(ye(last) .^ 2));
75 % three periods starting at one second, finely enough to see the drift
76 w = find(t >= 1 & t < 1 + 3 / f0);
77 out.tw = t(w)' * 1000;
78 out.yw = y(w)';
79 out.yew = ye(w)';
80 % predicted pitch of partials 1 and 10, for the status strip
81 out.pred1 = drv_pred(1, N, lam, dt, f0);
82 out.pred10 = drv_pred(10, N, lam, dt, f0);
83 out.clipExact = wv_audio(ye, fs, sc);
84 out.clipFd = wv_audio(max(min(y, 2 * sc), -2 * sc), fs, sc);
85 out.clipDiff = wv_audio(max(min(e, 2 * sc), -2 * sc), fs, sc);
86elseif strcmp(p.mode, 'cents')
87 f0 = p.f0;
88 N = p.N;
89 c = 2 * f0;
90 h = 1 / N;
91 x = (0:N)' * h;
92 lams = p.lams(:)';
93 nl = numel(lams);
94 % the partials to measure: every one of the first ten, and twenty more
95 % spread over the rest of the grid's range
96 sel = unique([1:min(10, N - 1), round(linspace(1, N - 1, 20))]);
97 S = sin(pi * sel(:) * (1:N - 1) / N); % discrete sine vectors
98 nall = 1:N - 1;
99 MEAS = NaN(nl, numel(sel));
100 PRED = NaN(nl, N - 1);
101 MS = NaN(1, nl);
102 for k = 1:nl
103 lam = lams(k);
104 dt = lam * h / c;
105 K = ceil(12 / (f0 * dt)); % twelve periods of the fundamental
106 % every partial present, with amplitude 1/n, released from rest
107 u0 = zeros(N + 1, 1);
108 um = zeros(N + 1, 1);
109 for n = nall
110 u0 = u0 + sin(n * pi * x) / n;
111 um = um + sin(n * pi * x) * cos(2 * pi * n * f0 * dt) / n;
112 end
113 t0 = tic;
114 [~, U] = leapfrog(u0, um, lam, K, 'fixed', 1, K);
115 MS(k) = toc(t0) * 1000;
116 A = S * U(2:N, :);
117 for i = 1:numel(sel)
118 fm = wv_peakfreq(A(i, :), dt, 0, 0.5 / dt);
119 MEAS(k, i) = wv_cents(fm, sel(i) * f0);
120 end
121 PRED(k, :) = drv_pred(nall, N, lam, dt, f0);
122 end
123 out.N = N;
124 out.f0 = f0;
125 out.lams = wv_vec(lams);
126 out.sel = sel;
127 out.nall = nall;
128 out.meas = wv_rows(MEAS);
129 out.pred = wv_rows(PRED);
130 out.ms = wv_vec(MS);
131 % the largest disagreement between measurement and prediction, in cents
132 D = NaN(nl, numel(sel));
133 for k = 1:nl
134 D(k, :) = MEAS(k, :) - PRED(k, sel);
135 end
136 out.agree = max(abs(D(:)));
137else
138 % the click in a periodic tube
139 c = p.c;
140 fs = p.fs;
141 dt = 1 / fs;
142 L = p.L;
143 N = round(L * p.lam / (c * dt));
144 h = L / N;
145 lam = c * dt / h;
146 M = round(p.T * fs);
147 x = (0:N - 1)' * h;
148 sx = c * p.sigt;
149 u0 = drv_train(x, 0, L, sx);
150 um = drv_train(x + c * dt, 0, L, sx); % moving right: u(x, -dt) = g(x + c dt)
151 jr = round(N / 2) + 1;
152 t0 = tic;
153 [y, ~] = leapfrog(u0, um, lam, M, 'periodic', jr, 0);
154 out.ms = toc(t0) * 1000;
155 t = (1:M)' * dt;
156 ye = drv_train(x(jr) - c * t, 0, L, sx);
157 out.N = N;
158 out.lam = lam;
159 out.fcut = asin(min(lam, 1)) / (pi * dt);
160 sc = max(abs(ye));
161 out.clipExact = wv_audio(ye, fs, sc);
162 out.clipFd = wv_audio(y, fs, sc);
163 % windows around the first and the twentieth arrival
164 tp = L / c;
165 ta = x(jr) / c;
166 for k = [1, 20]
167 tk = ta + (k - 1) * tp;
168 w = find(t >= tk - 0.002 & t < tk + 0.020);
169 out.(sprintf('tw%d', k)) = (t(w)' - tk) * 1000;
170 out.(sprintf('y%d', k)) = y(w)' / sc;
171 out.(sprintf('ye%d', k)) = ye(w)' / sc;
172 end
173 out.rel = sqrt(sum((y - ye) .^ 2) / sum(ye .^ 2));
174 % spectrogram of the simulated recording
175 nw = 512;
176 hop = 256;
177 nfr = floor((M - nw) / hop) + 1;
178 win = wv_window(nw);
179 nb = nw / 2 + 1;
180 G = zeros(nb, nfr);
181 for k = 1:nfr
182 seg = y((k - 1) * hop + (1:nw));
183 F = abs(fft(win .* seg));
184 G(:, k) = F(1:nb);
185 end
186 G = max(20 * log10(max(G / max(G(:)), 1e-6)), -80);
187 keep = 1:2:nb; % every other bin is plenty to draw
188 out.spec = wv_rows(G(keep, :));
189 out.specF = [0, (keep(end) - 1) * fs / nw];
190 out.specT = [0, ((nfr - 1) * hop + nw / 2) / fs];
191end
193wv_write(out);
195function v = drv_phi(x, beta)
196% the plucked shape, extended oddly about 0 and periodically with period 2
197r = mod(x, 2);
198sg = ones(size(r));
199sg(r > 1) = -1;
200r(r > 1) = 2 - r(r > 1);
201v = min(r / beta, (1 - r) / (1 - beta));
202v = sg .* v;
203end
205function g = drv_train(x, x0, L, sx)
206% a Gaussian pulse of width sx at x0, repeated with period L
207d = mod(x - x0 + L / 2, L) - L / 2;
208g = exp(-0.5 * (d / sx) .^ 2);
209end
211function c = drv_pred(n, N, lam, dt, f0)
212% the leapfrog dispersion relation sin(w dt / 2) = lam sin(n pi / (2 N)),
213% as a pitch error in cents against the exact n f0
214w = 2 * asin(lam * sin(n * pi / (2 * N))) / dt;
215c = wv_cents(w / (2 * pi), n * f0);
216end