Room impulse response and listening

const aMid = Math.min(0.95, Math.max(0.02, (res.bands[2].alpha + res.bands[3].alpha) / 2)), R = res.area * aMid / (1 - aMid);
const ir = new Float64Array(N);
ir[0] = 1;
// six first reflections of a box (a room of this volume and floor area when no box is given), source and listener `distance` apart
const V = res.volume, H = box?.h || Math.min(4, Math.cbrt(V)), base = V / H, lx = box?.lx || Math.sqrt(base * 1.4), ly = box?.ly || base / lx;
// 0.8 EG.49 (bim-sim_PART_3): a source and a listener the person placed (metres from the room's low corner, inside the box) replace the default pair
const inBox = (p) => (p ? [Math.min(lx - 0.05, Math.max(0.05, p[0])), Math.min(ly - 0.05, Math.max(0.05, p[1])), Math.min(H - 0.05, Math.max(0.05, p[2] ?? 1.2))] : null);
let src = inBox(source) || [0.3 * lx, 0.45 * ly, Math.min(1.5, H * 0.5)], lis = inBox(listener);
if (!lis) {
  lis = [src[0] + distance, src[1], src[2]];
  if (lis[0] > lx - 0.1) lis = [Math.max(0.1, src[0] - distance), src[1], src[2]];
}
const g = Math.sqrt(1 - aMid), d0 = Math.max(0.3, Math.hypot(lis[0] - src[0], lis[1] - src[1], lis[2] - src[2]) || distance);
let want = 16 * Math.PI * d0 * d0 / R;   // tail energy with the direct sound's energy as 1
let er = 0;
const images = [[-src[0], src[1], src[2]], [2 * lx - src[0], src[1], src[2]], [src[0], -src[1], src[2]], [src[0], 2 * ly - src[1], src[2]], [src[0], src[1], -src[2]], [src[0], src[1], 2 * H - src[2]]];
for (const im of images) {
  const d = Math.hypot(lis[0] - im[0], lis[1] - im[1], lis[2] - im[2]), idx = Math.round(((d - d0) / SOUND) * sr);
  if (idx > 0 && idx < N) { const a = g * d0 / d; ir[idx] += a; er += a * a; }
}
want = Math.max(want - er, want * 0.2);
let e = 0;
for (let n = 0; n < N; n++) e += tail[n] * tail[n];
const scale = e > 0 ? Math.sqrt(want / e) : 0;
for (let n = 0; n < N; n++) ir[n] += tail[n] * scale;
const out = new Float32Array(N);
for (let n = 0; n < N; n++) out[n] = ir[n];
return { ir: out, sr, rt, distance: d0, source: src, listener: lis, criticalDistance: Math.sqrt(R / (16 * Math.PI)) };

export function convolve(x, h) {
  const n = pow2(x.length + h.length - 1), ar = new Float64Array(n), ai = new Float64Array(n), br = new Float64Array(n), bi = new Float64Array(n);
  ar.set(x); br.set(h);
  fft(ar, ai); fft(br, bi);
  for (let i = 0; i < n; i++) { const r = ar[i] * br[i] - ai[i] * bi[i]; ai[i] = ar[i] * bi[i] + ai[i] * br[i]; ar[i] = r; }
  fft(ar, ai, true);
  const out = new Float32Array(x.length + h.length - 1);
  for (let i = 0; i < out.length; i++) out[i] = ar[i] / n;
  return out;
}

/** The dry sound in the room, its loudness matched to the dry sound's (over the dry sound's own length) so a switch between them changes the room, not the volume. */
export function listen(dry, ir, { tail = 0.5 } = {}) {
  const y = convolve(dry, ir), keep = Math.min(y.length, dry.length + Math.round(tail * 44100 * 4));
  const out = y.slice(0, keep), a = rms(dry), b = rms(out.subarray(0, dry.length));
  const k = b > 0 ? a / b : 1;
  for (let i = 0; i < out.length; i++) out[i] *= k;
  return out;
}

export function wavBytes(x, sr = 44100) {
  const n = x.length, buf = new ArrayBuffer(44 + n * 2), v = new DataView(buf);
  const str = (o, s) => { for (let i = 0; i < s.length; i++) v.setUint8(o + i, s.charCodeAt(i)); };
  str(0, 'RIFF'); v.setUint32(4, 36 + n * 2, true); str(8, 'WAVE'); str(12, 'fmt '); v.setUint32(16, 16, true); v.setUint16(20, 1, true); v.setUint16(22, 1, true);
  v.setUint32(24, sr, true); v.setUint32(28, sr * 2, true); v.setUint16(32, 2, true); v.setUint16(34, 16, true); str(36, 'data'); v.setUint32(40, n * 2, true);
  let peak = 0;
  for (let i = 0; i < n; i++) peak = Math.max(peak, Math.abs(x[i]));
  const k = peak > 1 ? 1 / peak : 1;   // never clip: scale down when the response is hot
  for (let i = 0; i < n; i++) v.setInt16(44 + i * 2, Math.round(Math.max(-1, Math.min(1, x[i] * k)) * 32767), true);
  return new Uint8Array(buf);
}