1% Driver for Chapter 8: a room, by image sources.
2%
3% The editable script supplies
4%
5% h = images(L, xs, xr, beta, fs, tmax, taps)
6%
7% the impulse response of a rectangular room whose walls all reflect a
8% fraction beta of the pressure, from source xs to receiver xr.
9%
10% What the driver measures is the decay of the response, and what it compares
11% it with are three predictions that share no algebra with the script:
12%
13% Sabine (1900) RT = 24 ln(10) V / (c S alpha)
14% Eyring (1930) RT = 24 ln(10) V / (c S (-ln(1 - alpha)))
15% the lattice the decay that the image lattice itself implies once
16% the images are too many to count: an image at distance
17% ct in direction u has been reflected about
18% ct (|ux|/Lx + |uy|/Ly + |uz|/Lz) times, so the energy
19% arriving at time t is the average over directions of
20% beta^(2 ct g(u)), computed here by quadrature on the
21% sphere.
22%
23% with alpha = 1 - beta^2 the absorption coefficient. The reverberation time
24% is read off the Schroeder curve, the backward integral of h^2, by a straight
25% line through its -5 and -25 dB points (T20, scaled to 60 dB), in the band
26% from 0.7 to 5.7 kHz: the three octaves around 1, 2 and 4 kHz together,
27% which averages away much of what one octave at one receiver fluctuates by.
28%
29% The band matters. Every image has the same sign, so arrivals that land
30% within a sample or two of each other add coherently at low frequency, and
31% the broadband response carries a slowly decaying low-frequency build-up
32% (for rigid walls, a pressure offset that never decays at all: a point
33% source has pushed air into a closed box). The lattice prediction is the
34% incoherent, high-frequency limit, which is what an octave band well above
35% the room's lowest modes measures, and it is also what the two formulas are
36% for. The broadband curve is returned too, so that the difference can be
37% seen.
38%
39% mode 'room' one room: the response, its decay, and what it does to a sound.
40% mode 'decay' reverberation time against absorption, by all four routes.
41p = wv_params();
42out = struct();
43L = p.L(:)';
44V = prod(L);
45S = 2 * (L(1) * L(2) + L(2) * L(3) + L(1) * L(3));
46c = 343;
47fs = p.fs;
48K = 24 * log(10) / c;
50if strcmp(p.mode, 'room')
51 alpha = p.alpha;
52 beta = sqrt(1 - alpha);
53 te = K * V / (S * (-log(1 - alpha)));
54 tmax = min(p.tcap, 0.75 * te + 0.1);
55 t0 = tic;
56 h = images(L, p.xs, p.xr, beta, fs, tmax, p.taps);
57 out.ms = toc(t0) * 1000;
58 t = (0:numel(h) - 1)' / fs;
59 [sch, rtb] = drv_schroeder(h, fs);
60 [schB, rt] = drv_schroeder(drv_band(h, fs, 700, 5600), fs);
61 out.rt = rt;
62 out.rtbroad = rtb;
63 out.sabine = K * V / (S * alpha);
64 out.eyring = te;
65 [tp, sp, rtp] = drv_lattice(L, beta, tmax, 300);
66 out.lattice = rtp;
67 step = max(1, floor(numel(t) / 1200));
68 out.ts = t(1:step:end)' * 1000;
69 out.sch = sch(1:step:end)';
70 out.schB = schB(1:step:end)';
71 out.tp = tp * 1000;
72 out.sp = sp;
73 out.tmax = tmax;
74 % the first 60 ms, finely
75 ne = min(numel(h), round(0.06 * fs));
76 out.te = t(1:ne)' * 1000;
77 out.he = (h(1:ne) / max(abs(h)))';
78 out.direct = norm(p.xs - p.xr) / c * 1000;
79 % the sound: a phrase of plucked notes, dry and in the room
80 dry = drv_phrase(p.notes, p.step, fs);
81 wet = wv_fftconv(dry, h);
82 out.clipDry = wv_audio(dry, fs);
83 out.clipWet = wv_audio(wet, fs);
84 out.clipIr = wv_audio(h, fs);
85 out.alpha = alpha;
86else
87 alphas = p.alphas(:)';
88 na = numel(alphas);
89 RT = NaN(1, na);
90 LAT = NaN(1, na);
91 MS = NaN(1, na);
92 for k = 1:na
93 beta = sqrt(1 - alphas(k));
94 te = K * V / (S * (-log(1 - alphas(k))));
95 tmax = min(p.tcap, 0.75 * te + 0.1);
96 t0 = tic;
97 h = images(L, p.xs, p.xr, beta, fs, tmax, p.taps);
98 MS(k) = toc(t0) * 1000;
99 [~, RT(k)] = drv_schroeder(drv_band(h, fs, 700, 5600), fs);
100 [~, ~, LAT(k)] = drv_lattice(L, beta, tmax, 300);
101 end
102 out.alphas = alphas;
103 out.rt = RT;
104 out.lattice = LAT;
105 out.sabine = K * V ./ (S * alphas);
106 out.eyring = K * V ./ (S * (-log(1 - alphas)));
107 out.ms = MS;
108end
109out.V = V;
110out.S = S;
111wv_write(out);
113function [db, rt] = drv_schroeder(h, fs)
114% Schroeder's backward integral of the squared response, in dB, and the
115% reverberation time from a line through its -5 and -25 dB crossings
116e = flipud(cumsum(flipud(h .^ 2)));
117db = 10 * log10(max(e / e(1), 1e-30));
118t = (0:numel(h) - 1)' / fs;
119i1 = find(db <= -5, 1);
120i2 = find(db <= -25, 1);
121rt = NaN;
122if ~isempty(i1) && ~isempty(i2) && i2 > i1 + 2
123 cf = polyfit(t(i1:i2), db(i1:i2), 1);
124 rt = -60 / cf(1);
125end
126end
128function y = drv_band(h, fs, flo, fhi)
129% the band from flo to fhi, by FFT, with a raised-cosine edge a sixth of an
130% octave wide at each end so that the filter itself rings only briefly
131n = numel(h);
132nf = 2 ^ nextpow2(2 * n);
133H = fft(h, nf);
134f = (0:nf - 1)' * fs / nf;
135f(f > fs / 2) = fs - f(f > fs / 2);
136hw = 0.5 * log2(fhi / flo); % half the width, in octaves
137o = abs(log2(max(f, 1e-9) / sqrt(flo * fhi))); % octaves from the centre
138B = zeros(nf, 1);
139B(o <= hw - 1 / 12) = 1;
140edge = o > hw - 1 / 12 & o < hw + 1 / 12;
141B(edge) = 0.5 + 0.5 * cos(pi * (o(edge) - (hw - 1 / 12)) / (1 / 6));
142y = real(ifft(H .* B));
143y = y(1:n);
144end
146function [t, db, rt] = drv_lattice(L, beta, tmax, nt)
147% The decay the image lattice implies, averaged over directions: the energy
148% arriving at time t is proportional to mean_u exp(-a g(u) t), with
149% a = -2 ln(beta) c and g(u) = sum |u_i| / L_i, so its backward integral is
150% mean_u exp(-a g t) / (a g). Midpoint quadrature in (cos theta, phi) over
151% one octant, which is enough by symmetry.
152c = 343;
153a = -2 * log(beta) * c;
154nq = 120;
155mu = ((1:nq) - 0.5) / nq; % cos(theta), uniform in solid angle
156ph = ((1:nq) - 0.5) / nq * pi / 2;
157[M, P] = meshgrid(mu, ph);
158st = real(sqrt(1 - M .^ 2));
159g = abs(st .* cos(P)) / L(1) + abs(st .* sin(P)) / L(2) + abs(M) / L(3);
160g = g(:);
161t = linspace(0, tmax, nt);
162e = zeros(1, nt);
163for k = 1:nt
164 e(k) = mean(exp(-a * g * t(k)) ./ (a * g));
165end
166db = 10 * log10(e / e(1));
167i1 = find(db <= -5, 1);
168i2 = find(db <= -25, 1);
169rt = NaN;
170if ~isempty(i1) && ~isempty(i2)
171 cf = polyfit(t(i1:i2), db(i1:i2), 1);
172 rt = -60 / cf(1);
173end
174end
176function y = drv_phrase(notes, step, fs)
177% plucked notes by modal synthesis, as in Chapter 1, one every `step` seconds
178ns = numel(notes);
179ring = 1.2;
180n = round((step * (ns - 1) + ring) * fs);
181y = zeros(n, 1);
182t = (0:round(ring * fs) - 1)' / fs;
183for i = 1:ns
184 f0 = wv_note(notes(i));
185 v = zeros(size(t));
186 for m = 1:floor(0.45 * fs / f0)
187 f = m * f0;
188 s = 6.908 * (1 + (f / 2000) ^ 2) / 1.5;
189 v = v + sin(m * pi * 0.2) / m * exp(-s * t) .* cos(2 * pi * f * t);
190 end
191 i0 = round((i - 1) * step * fs);
192 y(i0 + (1:numel(t))) = y(i0 + (1:numel(t))) + v;
193end
194end