d8098edAn interactive tutorial on the wave equation in acoustics: strings, drums and rooms, measured by earJeremy Magland 1% The image-source method for a rectangular room (Allen and Berkley 1979).
2%
3% A wall that reflects sound specularly behaves, for the listener, like a
4% window onto a mirror-image room with a mirror-image source in it. Mirror
5% the mirror rooms in their turn and the whole of space fills with a lattice
6% of image sources, eight to each cell of size 2 Lx by 2 Ly by 2 Lz:
7%
8% x = (1 - 2p) xs + 2 l Lx, p in {0, 1}, l any integer,
9%
10% and the same in y and z. The path from image to listener crosses the
11% walls x = 0 and x = Lx |l - p| and |l| times, so the image's pressure is
12%
13% beta^(number of reflections) / (4 pi d), arriving d / c after the source,
14%
15% with beta the pressure reflection factor of the walls. For rigid walls,
16% beta = 1, this is not a model: it is the exact solution of the wave
17% equation in the room, by the method of images. For beta < 1 it is a model,
18% one in which every wall reflects the same fraction at every angle and every
19% frequency, which no real wall does.
20%
21% Each arrival is a delta at a time that falls between samples. With taps = 0
22% it is put on the nearest sample, as Allen and Berkley did; otherwise it is
23% spread over `taps` samples as a Hann-windowed sinc, which is the delta
24% band-limited to the Nyquist frequency, cut off at the window.
25%
26% L = [Lx Ly Lz], xs and xr the source and the receiver; the result is the
27% impulse response, sampled at fs, up to tmax seconds.
29function h = images(L, xs, xr, beta, fs, tmax, taps)
30% the two ways of placing an arrival are two separate loops, because numbl
31% will not compile a loop in which a variable is assigned on only one branch
32% of an `if`, and a loop nest with the choice inside it is such a loop
33if taps == 0
34 h = img_nearest(L, xs, xr, beta, fs, tmax);
35else
36 h = img_sinc(L, xs, xr, beta, fs, tmax, taps);
37end
38end
40function h = img_nearest(L, xs, xr, beta, fs, tmax)
41c = 343;
42rmax = c * tmax;
43h = zeros(ceil(tmax * fs) + 2, 1);
44nl = ceil(rmax / (2 * L(1))) + 1;
45nm = ceil(rmax / (2 * L(2))) + 1;
46nk = ceil(rmax / (2 * L(3))) + 1;
47% beta^n from a table: numbl compiles an indexed read, but not a power whose
48% exponent it cannot prove is a whole number
49bp = beta .^ (0:2 * (nl + nm + nk) + 6);
50for p = 0:1
51 for q = 0:1
52 for r = 0:1
53 for l = -nl:nl
54 dx = (1 - 2 * p) * xs(1) + 2 * l * L(1) - xr(1);
55 ax = abs(l - p) + abs(l);
56 for m = -nm:nm
57 dy = (1 - 2 * q) * xs(2) + 2 * m * L(2) - xr(2);
58 ay = abs(m - q) + abs(m);
59 for k = -nk:nk
60 dz = (1 - 2 * r) * xs(3) + 2 * k * L(3) - xr(3);
61 d = real(sqrt(dx * dx + dy * dy + dz * dz));
62 if d <= rmax
63 a = bp(ax + ay + abs(k - r) + abs(k) + 1) / (4 * pi * d);
64 jn = floor(d / c * fs + 0.5) + 1;
65 h(jn) = h(jn) + a;
66 end
67 end
68 end
69 end
70 end
71 end
72end
73end
75function h = img_sinc(L, xs, xr, beta, fs, tmax, taps)
76c = 343;
77rmax = c * tmax;
78h = zeros(ceil(tmax * fs) + taps + 2, 1);
79nl = ceil(rmax / (2 * L(1))) + 1;
80nm = ceil(rmax / (2 * L(2))) + 1;
81nk = ceil(rmax / (2 * L(3))) + 1;
82bp = beta .^ (0:2 * (nl + nm + nk) + 6);
83half = taps / 2;
84for p = 0:1
85 for q = 0:1
86 for r = 0:1
87 for l = -nl:nl
88 dx = (1 - 2 * p) * xs(1) + 2 * l * L(1) - xr(1);
89 ax = abs(l - p) + abs(l);
90 for m = -nm:nm
91 dy = (1 - 2 * q) * xs(2) + 2 * m * L(2) - xr(2);
92 ay = abs(m - q) + abs(m);
93 for k = -nk:nk
94 dz = (1 - 2 * r) * xs(3) + 2 * k * L(3) - xr(3);
95 d = real(sqrt(dx * dx + dy * dy + dz * dz));
96 if d > rmax
97 continue;
98 end
99 a = bp(ax + ay + abs(k - r) + abs(k) + 1) / (4 * pi * d);
100 tau = d / c * fs;
101 j0 = floor(tau);
102 for j = max(0, j0 - half + 1):j0 + half
103 x = j - tau;
104 if abs(x) < 1e-12
105 s = 1;
106 else
107 s = sin(pi * x) / (pi * x);
108 end
109 h(j + 1) = h(j + 1) + a * s * (0.5 + 0.5 * cos(2 * pi * x / taps));
110 end
111 end
112 end
113 end
114 end
115 end
116end
117end