1% Driver for Chapter 6: the round drum.
2%
3% The editable script supplies
4%
5% [C, E, W, Nn, S, in] = diskstencil(N, method)
6%
7% the five-point Laplacian (times h^2) on the unknowns of an N-by-N grid over
8% [-1, 1]^2 that lie inside the unit disk, with some treatment of the rim.
9%
10% The exact modes of the unit disk are J_m(j_mk r) cos(m theta), with
11% frequency j_mk / (2 pi) when c = 1, where j_mk is the k-th zero of the
12% Bessel function J_m; wv_bjzeros computes them to rounding error.
13%
14% mode 'hear' the exact round drum, struck at radius r0, heard at a fixed
15% point, by modal synthesis; and a square drum with the same
16% fundamental, also exact, for comparison.
17% mode 'grid' the pitch error of the lowest modes on grids of increasing
18% size, for each rim treatment. The driver runs leapfrog itself,
19% at a Courant number of 0.95 times the largest the stencil
20% allows (found by power iteration on the stencil, and never
21% more than 1/sqrt(2)); the grid starts in the exact mode, and
22% the frequency is read off the projection of the field onto it.
23p = wv_params();
24out = struct();
26if strcmp(p.mode, 'hear')
27 fs = p.fs;
28 f01 = p.f01;
29 T = p.T;
30 t = (0:round(T * fs) - 1)' / fs;
31 c = 2 * pi * f01 / 2.404825557695773; % so that the lowest mode is f01, R = 1
32 sw = p.width;
33 y = zeros(size(t));
34 F = [];
35 Amp = [];
36 Lab = zeros(0, 2);
37 for m = 0:p.mmax
38 z = wv_bjzeros(m, p.kmax);
39 for k = 1:numel(z)
40 j = z(k);
41 f = c * j / (2 * pi);
42 if f > p.fmax
43 continue;
44 end
45 % mode shape normalised: int J_m(j r)^2 r dr = J_{m+1}(j)^2 / 2,
46 % times pi (m > 0) or 2 pi (m = 0) for the angle
47 nrm = besselj(m + 1, j) ^ 2 / 2 * pi * (1 + (m == 0));
48 ps = besselj(m, j * p.r0); % strike at angle 0
49 pr = besselj(m, j * p.rr) * cos(m * p.thr); % receiver
50 a = ps * pr / nrm * exp(-0.5 * (j * sw) ^ 2) / (2 * pi * f);
51 if abs(a) < 1e-12
52 continue;
53 end
54 y = y + a * sin(2 * pi * f * t);
55 F(end + 1) = f;
56 Amp(end + 1) = abs(a);
57 Lab(end + 1, :) = [m, k];
58 end
59 end
60 % the square drum with the same fundamental: unit side, c chosen so that
61 % f11 = f01, struck and heard at comparable points. Partials more than
62 % 80 dB below the loudest are left out; they cannot be heard over it.
63 cs = 2 * f01 / sqrt(2);
64 ysq = zeros(size(t));
65 nq = floor(2 * p.fmax / cs);
66 [nn, mm] = meshgrid(1:nq, 1:nq);
67 fq = cs / 2 * sqrt(mm .^ 2 + nn .^ 2);
68 aq = sin(mm * pi * 0.5 * (1 + p.r0 / 1.2)) .* sin(nn * pi * 0.5) ...
69 .* sin(mm * pi * 0.77) .* sin(nn * pi * 0.63) ...
70 .* exp(-0.5 * (pi * sqrt(mm .^ 2 + nn .^ 2) * sw / 2) .^ 2) ./ (2 * pi * fq);
71 aq(fq > p.fmax) = 0;
72 for q = find(abs(aq(:)) > 1e-4 * max(abs(aq(:))))'
73 ysq = ysq + aq(q) * sin(2 * pi * fq(q) * t);
74 end
75 % one decay for every partial of both drums, so that neither rings forever
76 env = exp(-t / p.tau);
77 y = y .* env;
78 ysq = ysq .* env;
79 [F, ord] = sort(F);
80 Amp = Amp(ord);
81 Lab = Lab(ord, :);
82 out.f = F;
83 out.adb = 20 * log10(Amp / max(Amp));
84 out.m = Lab(:, 1)';
85 out.k = Lab(:, 2)';
86 out.ratio = F / f01;
87 [fr, db] = wv_spectrum(y, fs, p.fmax, 800);
88 out.fr = fr;
89 out.db = db;
90 out.clipRound = wv_audio(y, fs);
91 out.clipSquare = wv_audio(ysq, fs);
92 % the first dozen partials, whatever the strike, for the table
93 tbl = zeros(0, 3);
94 for m = 0:5
95 z = wv_bjzeros(m, 3);
96 for k = 1:3
97 tbl(end + 1, :) = [m, k, z(k)];
98 end
99 end
100 [~, o] = sort(tbl(:, 3));
101 tbl = tbl(o(1:12), :);
102 out.tm = tbl(:, 1)';
103 out.tk = tbl(:, 2)';
104 out.tratio = (tbl(:, 3) / tbl(1, 3))';
105else
106 Ns = p.Ns(:)';
107 methods = p.methods;
108 modes = reshape(p.modes, 2, [])';
109 nm = numel(methods);
110 nq = size(modes, 1);
111 CE = NaN(nm * nq, numel(Ns));
112 LAM = NaN(nm, numel(Ns));
113 TM = NaN(nm, numel(Ns));
114 WORST = NaN(nm, numel(Ns));
115 for s = 1:nm
116 for j = 1:numel(Ns)
117 N = Ns(j);
118 h = 2 / N;
119 [C, E, W, Nn, S, in] = diskstencil(N, methods{s});
120 rho = drv_radius(C, E, W, Nn, S, in);
121 lam = 0.95 * min(2 / sqrt(rho), 1 / sqrt(2));
122 LAM(s, j) = lam;
123 dt = lam * h;
124 x = -1 + (0:N) * h;
125 [X, Y] = meshgrid(x, x);
126 R = sqrt(X .^ 2 + Y .^ 2);
127 TH = atan2(Y, X);
128 wst = 0;
129 el = 0;
130 for q = 1:nq
131 m = modes(q, 1);
132 z = wv_bjzeros(m, modes(q, 2));
133 jm = z(end);
134 phi = besselj(m, jm * R) .* cos(m * TH) .* in;
135 w = jm; % c = 1
136 M = ceil(p.periods * 2 * pi / (w * dt));
137 t0 = tic;
138 a = drv_leapfrog(phi, phi * cos(w * dt), C, E, W, Nn, S, lam, M, phi);
139 el = el + toc(t0);
140 fm = wv_peakfreq(a, dt, 0, 0.5 / dt);
141 cq = wv_cents(fm, w / (2 * pi));
142 CE((s - 1) * nq + q, j) = cq;
143 wst = max(wst, abs(cq));
144 end
145 TM(s, j) = el * 1000;
146 WORST(s, j) = wst;
147 end
148 end
149 out.Ns = Ns;
150 out.methods = methods;
151 out.modes = wv_rows(modes);
152 out.cents = wv_rows(CE);
153 out.worst = wv_rows(WORST);
154 out.lam = wv_rows(LAM);
155 out.ms = wv_rows(TM);
156 Q = NaN(nm, numel(Ns));
157 for s = 1:nm
158 Q(s, :) = wv_orders(WORST(s, :), Ns);
159 end
160 out.orders = wv_rows(Q);
161end
163wv_write(out);
165function a = drv_leapfrog(u, uold, C, E, W, Nn, S, lam, M, phi)
166% leapfrog with the script's stencil; returns the projection of the field
167% onto phi after each step
168n1 = size(u, 1);
169l2 = lam ^ 2;
170unew = zeros(n1, n1);
171a = zeros(M, 1);
172for n = 1:M
173 s = 0;
174 for j = 2:n1 - 1
175 for i = 2:n1 - 1
176 v = 2 * u(i, j) - uold(i, j) + l2 * (C(i, j) * u(i, j) + E(i, j) * u(i, j + 1) ...
177 + W(i, j) * u(i, j - 1) + Nn(i, j) * u(i + 1, j) + S(i, j) * u(i - 1, j));
178 unew(i, j) = v;
179 s = s + phi(i, j) * v;
180 end
181 end
182 uold = u;
183 u = unew;
184 a(n) = s;
185end
186end
188function r = drv_radius(C, E, W, Nn, S, in)
189% the spectral radius of the stencil, by power iteration: leapfrog with this
190% stencil is stable when lam^2 r <= 4
191n1 = size(C, 1);
192v = double(in) .* (mod((1:n1)' + (1:n1), 2) * 2 - 1); % a checkerboard: the roughest vector
193r = 0;
194w = zeros(n1, n1);
195for it = 1:300
196 for j = 2:n1 - 1
197 for i = 2:n1 - 1
198 w(i, j) = C(i, j) * v(i, j) + E(i, j) * v(i, j + 1) + W(i, j) * v(i, j - 1) ...
199 + Nn(i, j) * v(i + 1, j) + S(i, j) * v(i - 1, j);
200 end
201 end
202 nv = sqrt(sum(w(:) .^ 2));
203 rn = nv / sqrt(sum(v(:) .^ 2));
204 v = w / nv;
205 if abs(rn - r) < 1e-6 * rn
206 r = rn;
207 break;
208 end
209 r = rn;
210end
211end