% Driver for Chapter 8: a room, by image sources.
%
% The editable script supplies
%
%     h = images(L, xs, xr, beta, fs, tmax, taps)
%
% the impulse response of a rectangular room whose walls all reflect a
% fraction beta of the pressure, from source xs to receiver xr.
%
% What the driver measures is the decay of the response, and what it compares
% it with are three predictions that share no algebra with the script:
%
%   Sabine (1900)    RT = 24 ln(10) V / (c S alpha)
%   Eyring (1930)    RT = 24 ln(10) V / (c S (-ln(1 - alpha)))
%   the lattice      the decay that the image lattice itself implies once
%                    the images are too many to count: an image at distance
%                    ct in direction u has been reflected about
%                    ct (|ux|/Lx + |uy|/Ly + |uz|/Lz) times, so the energy
%                    arriving at time t is the average over directions of
%                    beta^(2 ct g(u)), computed here by quadrature on the
%                    sphere.
%
% with alpha = 1 - beta^2 the absorption coefficient.  The reverberation time
% is read off the Schroeder curve, the backward integral of h^2, by a straight
% line through its -5 and -25 dB points (T20, scaled to 60 dB), in the band
% from 0.7 to 5.7 kHz: the three octaves around 1, 2 and 4 kHz together,
% which averages away much of what one octave at one receiver fluctuates by.
%
% The band matters.  Every image has the same sign, so arrivals that land
% within a sample or two of each other add coherently at low frequency, and
% the broadband response carries a slowly decaying low-frequency build-up
% (for rigid walls, a pressure offset that never decays at all: a point
% source has pushed air into a closed box).  The lattice prediction is the
% incoherent, high-frequency limit, which is what an octave band well above
% the room's lowest modes measures, and it is also what the two formulas are
% for.  The broadband curve is returned too, so that the difference can be
% seen.
%
% mode 'room'   one room: the response, its decay, and what it does to a sound.
% mode 'decay'  reverberation time against absorption, by all four routes.
p = wv_params();
out = struct();
L = p.L(:)';
V = prod(L);
S = 2 * (L(1) * L(2) + L(2) * L(3) + L(1) * L(3));
c = 343;
fs = p.fs;
K = 24 * log(10) / c;

if strcmp(p.mode, 'room')
    alpha = p.alpha;
    beta = sqrt(1 - alpha);
    te = K * V / (S * (-log(1 - alpha)));
    tmax = min(p.tcap, 0.75 * te + 0.1);
    t0 = tic;
    h = images(L, p.xs, p.xr, beta, fs, tmax, p.taps);
    out.ms = toc(t0) * 1000;
    t = (0:numel(h) - 1)' / fs;
    [sch, rtb] = drv_schroeder(h, fs);
    [schB, rt] = drv_schroeder(drv_band(h, fs, 700, 5600), fs);
    out.rt = rt;
    out.rtbroad = rtb;
    out.sabine = K * V / (S * alpha);
    out.eyring = te;
    [tp, sp, rtp] = drv_lattice(L, beta, tmax, 300);
    out.lattice = rtp;
    step = max(1, floor(numel(t) / 1200));
    out.ts = t(1:step:end)' * 1000;
    out.sch = sch(1:step:end)';
    out.schB = schB(1:step:end)';
    out.tp = tp * 1000;
    out.sp = sp;
    out.tmax = tmax;
    % the first 60 ms, finely
    ne = min(numel(h), round(0.06 * fs));
    out.te = t(1:ne)' * 1000;
    out.he = (h(1:ne) / max(abs(h)))';
    out.direct = norm(p.xs - p.xr) / c * 1000;
    % the sound: a phrase of plucked notes, dry and in the room
    dry = drv_phrase(p.notes, p.step, fs);
    wet = wv_fftconv(dry, h);
    out.clipDry = wv_audio(dry, fs);
    out.clipWet = wv_audio(wet, fs);
    out.clipIr = wv_audio(h, fs);
    out.alpha = alpha;
else
    alphas = p.alphas(:)';
    na = numel(alphas);
    RT = NaN(1, na);
    LAT = NaN(1, na);
    MS = NaN(1, na);
    for k = 1:na
        beta = sqrt(1 - alphas(k));
        te = K * V / (S * (-log(1 - alphas(k))));
        tmax = min(p.tcap, 0.75 * te + 0.1);
        t0 = tic;
        h = images(L, p.xs, p.xr, beta, fs, tmax, p.taps);
        MS(k) = toc(t0) * 1000;
        [~, RT(k)] = drv_schroeder(drv_band(h, fs, 700, 5600), fs);
        [~, ~, LAT(k)] = drv_lattice(L, beta, tmax, 300);
    end
    out.alphas = alphas;
    out.rt = RT;
    out.lattice = LAT;
    out.sabine = K * V ./ (S * alphas);
    out.eyring = K * V ./ (S * (-log(1 - alphas)));
    out.ms = MS;
end
out.V = V;
out.S = S;
wv_write(out);

function [db, rt] = drv_schroeder(h, fs)
% Schroeder's backward integral of the squared response, in dB, and the
% reverberation time from a line through its -5 and -25 dB crossings
e = flipud(cumsum(flipud(h .^ 2)));
db = 10 * log10(max(e / e(1), 1e-30));
t = (0:numel(h) - 1)' / fs;
i1 = find(db <= -5, 1);
i2 = find(db <= -25, 1);
rt = NaN;
if ~isempty(i1) && ~isempty(i2) && i2 > i1 + 2
    cf = polyfit(t(i1:i2), db(i1:i2), 1);
    rt = -60 / cf(1);
end
end

function y = drv_band(h, fs, flo, fhi)
% the band from flo to fhi, by FFT, with a raised-cosine edge a sixth of an
% octave wide at each end so that the filter itself rings only briefly
n = numel(h);
nf = 2 ^ nextpow2(2 * n);
H = fft(h, nf);
f = (0:nf - 1)' * fs / nf;
f(f > fs / 2) = fs - f(f > fs / 2);
hw = 0.5 * log2(fhi / flo);                   % half the width, in octaves
o = abs(log2(max(f, 1e-9) / sqrt(flo * fhi)));  % octaves from the centre
B = zeros(nf, 1);
B(o <= hw - 1 / 12) = 1;
edge = o > hw - 1 / 12 & o < hw + 1 / 12;
B(edge) = 0.5 + 0.5 * cos(pi * (o(edge) - (hw - 1 / 12)) / (1 / 6));
y = real(ifft(H .* B));
y = y(1:n);
end

function [t, db, rt] = drv_lattice(L, beta, tmax, nt)
% The decay the image lattice implies, averaged over directions: the energy
% arriving at time t is proportional to mean_u exp(-a g(u) t), with
% a = -2 ln(beta) c and g(u) = sum |u_i| / L_i, so its backward integral is
% mean_u exp(-a g t) / (a g).  Midpoint quadrature in (cos theta, phi) over
% one octant, which is enough by symmetry.
c = 343;
a = -2 * log(beta) * c;
nq = 120;
mu = ((1:nq) - 0.5) / nq;               % cos(theta), uniform in solid angle
ph = ((1:nq) - 0.5) / nq * pi / 2;
[M, P] = meshgrid(mu, ph);
st = real(sqrt(1 - M .^ 2));
g = abs(st .* cos(P)) / L(1) + abs(st .* sin(P)) / L(2) + abs(M) / L(3);
g = g(:);
t = linspace(0, tmax, nt);
e = zeros(1, nt);
for k = 1:nt
    e(k) = mean(exp(-a * g * t(k)) ./ (a * g));
end
db = 10 * log10(e / e(1));
i1 = find(db <= -5, 1);
i2 = find(db <= -25, 1);
rt = NaN;
if ~isempty(i1) && ~isempty(i2)
    cf = polyfit(t(i1:i2), db(i1:i2), 1);
    rt = -60 / cf(1);
end
end

function y = drv_phrase(notes, step, fs)
% plucked notes by modal synthesis, as in Chapter 1, one every `step` seconds
ns = numel(notes);
ring = 1.2;
n = round((step * (ns - 1) + ring) * fs);
y = zeros(n, 1);
t = (0:round(ring * fs) - 1)' / fs;
for i = 1:ns
    f0 = wv_note(notes(i));
    v = zeros(size(t));
    for m = 1:floor(0.45 * fs / f0)
        f = m * f0;
        s = 6.908 * (1 + (f / 2000) ^ 2) / 1.5;
        v = v + sin(m * pi * 0.2) / m * exp(-s * t) .* cos(2 * pi * f * t);
    end
    i0 = round((i - 1) * step * fs);
    y(i0 + (1:numel(t))) = y(i0 + (1:numel(t))) + v;
end
end