d8098edAn interactive tutorial on the wave equation in acoustics: strings, drums and rooms, measured by earJeremy Magland 1% Driver for Chapters 4 and 5: the square membrane.
2%
3% The editable script supplies
4%
5% [y, U] = drum2d(u, uold, lam, M, scheme, ir, jr, K)
6%
7% which advances a scheme M steps on the unit square, rim held at zero.
8%
9% The exact solution is a sum of the modes sin(m pi x) sin(n pi y), each with
10% frequency (c/2) sqrt(m^2 + n^2). Those same sine vectors, sampled on the
11% grid, are exact eigenvectors of every stencil in these chapters, so a
12% scheme started from a single mode keeps it and changes only its frequency,
13% which the driver measures from the recording and nothing else.
14%
15% mode 'strike' a struck drum, with any scheme: the membrane at rest, given a
16% Gaussian blow of velocity at one point, simulated and heard
17% at another. The reference is the modal sum, evaluated
18% exactly at the receiver. The time step is lam h / c, so the
19% sound comes out at its own sample rate, 1/dt.
20% mode 'cents' the pitch error of modes in every direction, at three
21% wavelengths, for several Courant numbers (c = 1 here).
22% mode 'orders' the pitch error of three fixed modes against N, for each
23% scheme at its own Courant number, and the time each run took.
24p = wv_params();
25out = struct();
27if strcmp(p.mode, 'strike')
28 c = p.c;
29 N = p.N;
30 h = 1 / N;
31 lam = p.lam;
32 dt = lam * h / c;
33 fs = 1 / dt;
34 M = round(p.T / dt);
35 K = p.kmax;
36 [nn, mm] = meshgrid(1:K, 1:K);
37 k2 = (mm .^ 2 + nn .^ 2) * pi ^ 2;
38 w = c * sqrt(k2);
39 sw = p.width;
40 % sine coefficients of the initial velocity exp(-r^2 / (2 sw^2)); the
41 % Gaussian is far enough from the rim that its truncation there is below
42 % rounding, so these are exact
43 A = 4 * 2 * pi * sw ^ 2 * sin(mm * pi * p.xs) .* sin(nn * pi * p.ys) .* exp(-k2 * sw ^ 2 / 2);
44 x = (0:N)' * h;
45 S = sin(pi * x * (1:K));
46 u0 = zeros(N + 1, N + 1);
47 um = S * (A .* sin(-w * dt) ./ w) * S';
48 ir = round(p.xr * N) + 1;
49 jr = round(p.yr * N) + 1;
50 t0 = tic;
51 [y, ~] = drum2d(u0, um, lam, M, p.scheme, ir, jr, 0);
52 out.ms = toc(t0) * 1000;
53 % the exact signal at the same node
54 amp = A .* sin(mm * pi * x(ir)) .* sin(nn * pi * x(jr)) ./ w;
55 keep = find(abs(amp) > 1e-10 * max(abs(amp(:))));
56 t = (1:M)' * dt;
57 ye = zeros(M, 1);
58 for q = keep'
59 ye = ye + amp(q) * sin(w(q) * t);
60 end
61 e = y - ye;
62 out.N = N;
63 out.lam = lam;
64 out.fs = fs;
65 out.rel = sqrt(sum(e .^ 2) / sum(ye .^ 2));
66 out.nmodes = numel(keep);
67 sc = max(abs(ye));
68 out.clipExact = wv_audio(ye, fs, sc);
69 out.clipFd = wv_audio(y, fs, sc);
70 out.clipDiff = wv_audio(e, fs, sc);
71 nw = min(M, round(0.4 * fs));
72 [fr, dbe] = wv_spectrum(ye(1:nw), fs, p.fmax, 800);
73 [~, dbf] = wv_spectrum(y(1:nw), fs, p.fmax, 800);
74 out.fr = fr;
75 out.dbe = dbe;
76 out.dbf = dbf;
77 % the lowest few mode frequencies that the blow excites, for the plot
78 fm = w(:) / (2 * pi);
79 am = abs(amp(:));
80 [fm, ord] = sort(fm);
81 am = am(ord);
82 sel = am > 1e-3 * max(am) & fm < p.fmax;
83 fsel = fm(sel);
84 out.fmodes = fsel(1:min(end, 40))';
85 out.f11 = c * sqrt(2) / 2;
86 w0 = find(t >= 0.1 & t < 0.1 + 4 / out.f11);
87 out.tw = t(w0)' * 1000;
88 out.yw = y(w0)';
89 out.yew = ye(w0)';
90elseif strcmp(p.mode, 'cents')
91 N = p.N;
92 h = 1 / N;
93 x = (0:N)' * h;
94 lams = p.lams(:)';
95 Rs = p.rs(:)';
96 th = linspace(4, 86, p.nang) * pi / 180;
97 MN = zeros(0, 3);
98 for r = 1:numel(Rs)
99 R = Rs(r);
100 for q = 1:numel(th)
101 m = round(R * cos(th(q)));
102 n = round(R * sin(th(q)));
103 if m >= 1 && n >= 1 && m < N && n < N
104 MN(end + 1, :) = [m, n, r];
105 end
106 end
107 end
108 % no mode twice
109 [~, iu] = unique(MN(:, 1) * 10000 + MN(:, 2));
110 MN = MN(sort(iu), :);
111 nq = size(MN, 1);
112 Phi = zeros((N + 1) ^ 2, nq);
113 for q = 1:nq
114 f = sin(MN(q, 1) * pi * x) * sin(MN(q, 2) * pi * x)';
115 Phi(:, q) = f(:);
116 end
117 wq = pi * sqrt(MN(:, 1) .^ 2 + MN(:, 2) .^ 2); % exact, with c = 1
118 fex = wq / (2 * pi);
119 MEAS = NaN(numel(lams), nq);
120 PRED = NaN(numel(lams), nq);
121 for k = 1:numel(lams)
122 lam = lams(k);
123 dt = lam * h;
124 Kst = ceil(12 / (min(fex) * dt));
125 u0 = reshape(sum(Phi, 2), N + 1, N + 1);
126 um = reshape(Phi * cos(wq * dt), N + 1, N + 1);
127 [~, U] = drum2d(u0, um, lam, Kst, 'fd2', 1, 1, Kst);
128 A = Phi' * U;
129 for q = 1:nq
130 fmq = wv_peakfreq(A(q, :), dt, 0, 0.5 / dt);
131 MEAS(k, q) = wv_cents(fmq, fex(q));
132 end
133 s2 = sin(MN(:, 1) * pi * h / 2) .^ 2 + sin(MN(:, 2) * pi * h / 2) .^ 2;
134 wp = 2 * asin(lam * sqrt(s2)) / dt;
135 PRED(k, :) = wv_cents(wp / (2 * pi), fex)';
136 end
137 out.N = N;
138 out.lams = wv_vec(lams);
139 out.m = MN(:, 1)';
140 out.n = MN(:, 2)';
141 out.ring = MN(:, 3)';
142 out.rs = Rs;
143 out.angle = atan2(MN(:, 2), MN(:, 1))' * 180 / pi;
144 out.ppw = (2 * N ./ sqrt(MN(:, 1) .^ 2 + MN(:, 2) .^ 2))';
145 out.meas = wv_rows(MEAS);
146 out.pred = wv_rows(PRED);
147 out.agree = max(abs(MEAS(:) - PRED(:)));
148else
149 Ns = p.Ns(:)';
150 schemes = p.schemes;
151 lams = p.lams(:)';
152 modes = reshape(p.modes, 2, [])';
153 nmode = size(modes, 1);
154 ns = numel(schemes);
155 CE = NaN(ns * nmode, numel(Ns));
156 TM = NaN(ns, numel(Ns));
157 WORST = NaN(ns, numel(Ns));
158 for s = 1:ns
159 drum2d(zeros(12), zeros(12), 0.5, 20, schemes{s}, 3, 3, 0); % warm-up
160 for j = 1:numel(Ns)
161 N = Ns(j);
162 h = 1 / N;
163 x = (0:N)' * h;
164 lam = lams(s);
165 dt = lam * h;
166 % the first interior node: sin(m pi / N) is never zero for m < N,
167 % so no mode has a nodal line through it
168 ir = 2;
169 jr = 2;
170 wst = 0;
171 for q = 1:nmode
172 m = modes(q, 1);
173 n = modes(q, 2);
174 if m >= N || n >= N
175 continue;
176 end
177 wq = pi * sqrt(m ^ 2 + n ^ 2);
178 f = sin(m * pi * x) * sin(n * pi * x)';
179 M = ceil(12 * 2 * pi / (wq * dt));
180 if q == 1
181 % the lowest mode's run is also the one that is timed:
182 % twelve of its periods, whatever the scheme, fastest of
183 % at least three runs while they are cheap
184 best = Inf;
185 tot = 0;
186 rep = 0;
187 while rep == 0 || (rep < 20 && tot < 0.3 && (rep < 3 || tot < 0.05))
188 t0 = tic;
189 y = drum2d(f, f * cos(wq * dt), lam, M, schemes{s}, ir, jr, 0);
190 el = toc(t0);
191 best = min(best, el);
192 tot = tot + el;
193 rep = rep + 1;
194 end
195 TM(s, j) = best * 1000;
196 else
197 y = drum2d(f, f * cos(wq * dt), lam, M, schemes{s}, ir, jr, 0);
198 end
199 fm = wv_peakfreq(y, dt, 0, 0.5 / dt);
200 cq = wv_cents(fm, wq / (2 * pi));
201 CE((s - 1) * nmode + q, j) = cq;
202 wst = max(wst, abs(cq));
203 end
204 WORST(s, j) = wst;
205 end
206 end
207 out.Ns = Ns;
208 out.schemes = schemes;
209 out.lams = wv_vec(lams);
210 out.modes = wv_rows(modes);
211 out.cents = wv_rows(CE);
212 out.worst = wv_rows(WORST);
213 out.ms = wv_rows(TM);
214 Q = NaN(ns, numel(Ns));
215 for s = 1:ns
216 Q(s, :) = wv_orders(WORST(s, :), Ns);
217 end
218 out.orders = wv_rows(Q);
219end
221wv_write(out);