d8098edAn interactive tutorial on the wave equation in acoustics: strings, drums and rooms, measured by earJeremy Magland 1% Driver for Chapter 9: a room, by finite differences.
2%
3% The editable script supplies
4%
5% y = room3d(n, lam, M, is, ir, src)
6%
7% leapfrog in a rigid box on a grid of n(1) by n(2) by n(3) nodes, driven at
8% node is and heard at node ir.
9%
10% The source is a point source of strength f(t), the derivative of a
11% Gaussian, which puts no net volume of air into the box. For rigid walls
12% the exact pressure is the image sum
13%
14% u(t) = sum over images f(t - d/c) / (4 pi c^2 d),
15%
16% which the driver evaluates directly; Chapter 8's image method with beta = 1
17% is the same formula, and here it is not a model but the solution.
18%
19% mode 'converge' the error at the receiver against the grid, and the time
20% per node per step, from which the driver extrapolates what
21% a full-band simulation would cost.
22% mode 'hear' a longer run on two grids, against the exact response.
23p = wv_params();
24out = struct();
25c = 343;
26L = p.L(:)';
27lam = p.lam;
28sg = p.sigma;
29t0p = 5 * sg;
31if strcmp(p.mode, 'converge')
32 Ns = p.Ns(:)';
33 E = NaN(1, numel(Ns));
34 TM = NaN(1, numel(Ns));
35 NS = NaN(1, numel(Ns));
36 cells = NaN(1, numel(Ns));
37 for k = 1:numel(Ns)
38 [y, ye, dt, ms, n] = drv_run(Ns(k), L, lam, p.T, p.xs, p.xr, sg, t0p, c);
39 E(k) = sqrt(sum((y - ye) .^ 2) / sum(ye .^ 2));
40 TM(k) = ms;
41 cells(k) = prod(n);
42 NS(k) = ms * 1e6 / (prod(n) * numel(y)); % nanoseconds per node per step
43 if k == numel(Ns)
44 out.t = (1:numel(y)) * dt * 1000;
45 out.y = y';
46 out.ye = ye';
47 end
48 if k == 1
49 out.t1 = (1:numel(y)) * dt * 1000;
50 out.y1 = y';
51 end
52 end
53 out.Ns = Ns;
54 out.err = E;
55 out.ms = TM;
56 out.cells = cells;
57 out.ns = NS;
58 out.orders = wv_orders(E, Ns);
59 % points per wavelength at the pulse's upper frequency, taken as the
60 % point where its spectrum is 40 dB below its peak
61 f = linspace(1, 20000, 20000);
62 S = f .* exp(-0.5 * (2 * pi * f * sg) .^ 2);
63 S = S / max(S);
64 fhi = f(find(S > 0.01, 1, 'last'));
65 out.fhi = fhi;
66 out.ppw = c / fhi ./ (L(1) ./ Ns);
67 % what a full-band run of Chapter 8's room would cost at the measured
68 % rate: 10 points per shortest wavelength at 20 kHz, one second of sound.
69 % The rate is the best of the grids from N = 32 up, where the fixed costs
70 % no longer matter and a single slow run (a garbage collection, another
71 % tab) does not decide the answer.
72 rate = min(NS(3:end));
73 out.rate = rate;
74 V8 = 81;
75 h8 = c / 20000 / 10;
76 dt8 = lam * h8 / c;
77 upd = V8 / h8 ^ 3 / dt8;
78 out.fullband = upd * rate * 1e-9; % seconds
79 out.fullcells = V8 / h8 ^ 3;
80 h2 = c / 2000 / 10; % the same up to 2 kHz only
81 out.lowband = V8 / h2 ^ 3 / (lam * h2 / c) * rate * 1e-9;
82else
83 Ns = p.Ns(:)';
84 fs = p.fs;
85 te = (1:round(p.T * fs))' / fs;
86 ye = drv_exact(te, L, p.xs, p.xr, sg, t0p, c);
87 sc = max(abs(ye));
88 out.clipExact = wv_audio(ye, fs, sc);
89 clips = cell(1, numel(Ns));
90 for k = 1:numel(Ns)
91 [y, ~, dt, ~, ~] = drv_run(Ns(k), L, lam, p.T, p.xs, p.xr, sg, t0p, c, true);
92 clips{k} = wv_audio(y, 1 / dt, sc);
93 fsk(k) = 1 / dt;
94 end
95 out.clips = clips;
96 out.Ns = Ns;
97 out.fsk = fsk;
98 nw = round(0.25 * fs);
99 [fr, db] = wv_spectrum(ye(1:nw), fs, 3000, 600);
100 out.fr = fr;
101 out.dbe = db;
102end
103wv_write(out);
105function [y, ye, dt, ms, n] = drv_run(N, L, lam, T, xs, xr, sg, t0p, c, noexact)
106h = L(1) / N;
107n = round(L / h) + 1;
108dt = lam * h / c;
109M = round(T / dt);
110t = (1:M)' * dt;
111% step k computes u at t(k) from the two levels before it, so the source term
112% it adds is the one at the middle level, f(t(k) - dt)
113ts = t - dt;
114f = drv_pulse(ts, sg, t0p);
115src = dt ^ 2 * f / h ^ 3;
116is = round(xs / h) + 1;
117ir = round(xr / h) + 1;
118room3d(n, lam, 3, is, ir, src(1:3)); % warm-up
119tic0 = tic;
120y = room3d(n, lam, M, is, ir, src);
121ms = toc(tic0) * 1000;
122if nargin > 9 && noexact
123 ye = [];
124else
125 ye = drv_exact(t, L, xs, xr, sg, t0p, c);
126end
127end
129function f = drv_pulse(t, sg, t0p)
130% the derivative of a Gaussian: no DC, so no permanent pressure in the box
131s = (t - t0p) / sg;
132f = -s .* exp(-0.5 * s .^ 2);
133end
135function y = drv_exact(t, L, xs, xr, sg, t0p, c)
136% the image sum for rigid walls, every image close enough to be heard by
137% the end of t
138y = zeros(size(t));
139tend = t(end);
140rmax = c * (tend + 0.001);
141nl = ceil(rmax ./ (2 * L)) + 1;
142dtt = t(2) - t(1);
143w = 7 * sg; % the pulse is below 1e-10 beyond this
144for p = 0:1
145 for q = 0:1
146 for r = 0:1
147 for l = -nl(1):nl(1)
148 dx = (1 - 2 * p) * xs(1) + 2 * l * L(1) - xr(1);
149 for m = -nl(2):nl(2)
150 dy = (1 - 2 * q) * xs(2) + 2 * m * L(2) - xr(2);
151 for k = -nl(3):nl(3)
152 dz = (1 - 2 * r) * xs(3) + 2 * k * L(3) - xr(3);
153 d = real(sqrt(dx * dx + dy * dy + dz * dz));
154 if d > rmax
155 continue;
156 end
157 tc = d / c + t0p;
158 i0 = max(1, floor((tc - w - t(1)) / dtt) + 1);
159 i1 = min(numel(t), ceil((tc + w - t(1)) / dtt) + 1);
160 a = 1 / (4 * pi * c * c * d);
161 for i = i0:i1
162 s = (t(i) - tc) / sg;
163 y(i) = y(i) - a * s * exp(-0.5 * s * s);
164 end
165 end
166 end
167 end
168 end
169 end
170end
171end