Resources Numerical accuracy Calculation policy
Works out how sound energy arrives at a listener over time in a room, by image sources or by ray tracing, and from that the reverberation and clarity measures per octave band.
The echogram, the energy arriving in each time bin, is built in Float64Array. Before the measures are taken, the square root of each bin is stored in a Float32Array, so the amplitudes the measures read carry about 7 digits. The methods themselves set a coarser limit: arrival times are binned at the sample rate, ray tracing is a sample of 20,000 rays by default, and the speed of sound is the fixed value 343 m/s.
Open item. Build the amplitude array in echogramMetrics as a Float64Array instead of a Float32Array, so the measures are taken from 64-bit values.
6 calculations in 1 source file. The code is read from the released app when this page loads. It is the calculation people run, not a copy kept for this page.
Tests None found by the review
The speed of sound used for every delay.
c = 343 m/s
A fixed, rounded value. The env settings that are passed on to the air absorption do not change it, so every delay is distance ÷ 343 whatever those settings say.
export const SOUND_SPEED = 343;The energy arriving at the listener from the source and every mirror image of it within the time limit.
image = ((1 − 2p)·s + 2·i·L); e[round(d ÷ c × sr)] += Π (1 − α)^hits × exp(−m·d) ÷ (4π·d²)
Distances by Math.hypot and energies added in 64-bit into a Float64Array. Each arrival is put in the nearest sample bin, so its time is known to half a sample (62.5 µs at the default 8000 Hz). The distance in the divisor is floored at 0.05 m.
/** Image sources of a box lx × ly × lz. `walls` = { x0, x1, y0, y1, z0, z1 } each an array of 6 band absorptions. */
export function imageSourceBox(box, source, listener, { walls, sr = 8000, tmax = 0.8, env = {} } = {}) {
const c = SOUND_SPEED, R = c * tmax, N = Math.ceil(tmax * sr);
const L = [box.lx, box.ly, box.lz];
const bands = BANDS.map((f) => ({ f, e: new Float64Array(N) })), m = BANDS.map((f) => (env.air === false ? 0 : airM(f, env)));
const rf = BANDS.map((_, b) => ({ x0: 1 - walls.x0[b], x1: 1 - walls.x1[b], y0: 1 - walls.y0[b], y1: 1 - walls.y1[b], z0: 1 - walls.z0[b], z1: 1 - walls.z1[b] }));
const n = L.map((l) => Math.ceil(R / (2 * l)) + 1);
let count = 0;
for (let i = -n[0]; i <= n[0]; i++) for (let j = -n[1]; j <= n[1]; j++) for (let k = -n[2]; k <= n[2]; k++) for (let p = 0; p < 8; p++) {
const px = p & 1, py = p >> 1 & 1, pz = p >> 2 & 1;
const ix = (1 - 2 * px) * source[0] + 2 * i * L[0] + 2 * px * 0 + (px ? 0 : 0), iy = (1 - 2 * py) * source[1] + 2 * j * L[1], iz = (1 - 2 * pz) * source[2] + 2 * k * L[2];
const d = Math.hypot(ix - listener[0], iy - listener[1], iz - listener[2]);
if (d > R) continue;
// reflections on each wall pair: image index (i, px) → |2i − px| hits of the x1 side and |2i| … of x0; standard counts
const hx0 = Math.abs(i - px), hx1 = Math.abs(i), hy0 = Math.abs(j - py), hy1 = Math.abs(j), hz0 = Math.abs(k - pz), hz1 = Math.abs(k);
const idx = Math.round(d / c * sr);
if (idx >= N) continue;
count++;
for (let b = 0; b < bands.length; b++) {
const r = rf[b], g = Math.pow(r.x0, hx0) * Math.pow(r.x1, hx1) * Math.pow(r.y0, hy0) * Math.pow(r.y1, hy1) * Math.pow(r.z0, hz0) * Math.pow(r.z1, hz1);
bands[b].e[idx] += g * Math.exp(-m[b] * d) / (4 * Math.PI * Math.max(d, 0.05) ** 2);
}
}
return { sr, bands, method: 'image-source', images: count, direct: Math.hypot(source[0] - listener[0], source[1] - listener[1], source[2] - listener[2]) };
}The unit normal and plane distance of each room surface, from its first three corners.
n = (b − a) × (c − a) ÷ |…|; d = n · a
The normal comes from the first three corners only; the rest are assumed to lie in the same plane. The division by the cross product's length is not guarded against three corners in a line.
function prep(surfaces) {
return surfaces.map((s) => {
const [a, b, c] = s.pts, u = sub(b, a), v = sub(c, a); let n = cross(u, v); const l = Math.hypot(...n); n = n.map((x) => x / l);
return { ...s, n, d: dot(n, a) };
});
}Whether a point on a surface's plane lies inside its convex outline.
the sign of ((b − a) × (p − a)) · n is the same for every edge
Ordinary 64-bit cross and dot products. A value within 1e-12 of zero (absolute) is treated as on the edge and accepted. A ray through a corner or an edge shared by two surfaces can therefore be accepted by both; the nearer one wins.
/** Is point p (on the plane) inside the convex polygon? */
function inside(s, p) {
const pts = s.pts; let sign = 0;
for (let i = 0; i < pts.length; i++) {
const a = pts[i], b = pts[(i + 1) % pts.length], k = dot(cross(sub(b, a), sub(p, a)), s.n);
if (Math.abs(k) < 1e-12) continue;
const sg = k > 0 ? 1 : -1;
if (sign === 0) sign = sg; else if (sg !== sign) return false;
}
return true;
}The energy arriving at a spherical receiver from rays sent out evenly from the source and reflected round the room.
e[⌊t × sr⌋] += E × chord ÷ (4⁄3·π·r³) × exp(−m·c·t); after a hit E = E × (1 − α)
The result is a sample: 20,000 rays by default on a golden-angle spiral, with a fixed seed so a run repeats exactly. The chord through the receiver comes from a quadratic whose discriminant b² − c subtracts nearly equal numbers when the ray just grazes the sphere. A specular reflection subtracts 2(d·n)n from the direction, which is not renormalized between reflections.
/** Ray tracing. room = { surfaces: [{ pts (convex planar polygon), alpha[6], scatter }] }. */
export function traceRoom(room, source, listener, { rays = 20000, sr = 8000, tmax = 0.8, radius = 0.5, seed = 1, env = {} } = {}) {
const c = SOUND_SPEED, N = Math.ceil(tmax * sr), surf = prep(room.surfaces), rand = rng(seed);
const nb = BANDS.length, bands = BANDS.map((f) => ({ f, e: new Float64Array(N) })), m = BANDS.map((f) => (env.air === false ? 0 : airM(f, env)));
const vs = 4 / 3 * Math.PI * radius ** 3, golden = Math.PI * (3 - Math.sqrt(5));
const en = new Float64Array(nb);
for (let r = 0; r < rays; r++) {
const z = 1 - 2 * (r + 0.5) / rays, rr = Math.sqrt(1 - z * z);
let dir = [rr * Math.cos(golden * r), rr * Math.sin(golden * r), z], pos = source.slice(), t = 0;
en.fill(1 / rays);
for (let bounce = 0; bounce < 400 && t < tmax; bounce++) {
let best = Infinity, hit = -1;
for (let i = 0; i < surf.length; i++) {
const s = surf[i], den = dot(s.n, dir);
if (Math.abs(den) < 1e-12) continue;
const tt = (s.d - dot(s.n, pos)) / den;
if (tt > 1e-7 && tt < best && inside(s, [pos[0] + dir[0] * tt, pos[1] + dir[1] * tt, pos[2] + dir[2] * tt])) { best = tt; hit = i; }
}
if (hit < 0) break;
// receiver: the chord of this segment inside the sphere
const w = sub(pos, listener), bq = dot(w, dir), cq = dot(w, w) - radius * radius, disc = bq * bq - cq;
if (disc > 0) {
const sq = Math.sqrt(disc), t0 = Math.max(0, -bq - sq), t1 = Math.min(best, -bq + sq);
if (t1 > t0) {
const chord = t1 - t0, tm = t + (t0 + t1) / 2 / c, idx = Math.floor(tm * sr);
if (idx >= 0 && idx < N) for (let b = 0; b < nb; b++) bands[b].e[idx] += en[b] * chord / vs * Math.exp(-m[b] * (tm * c));
}
}
pos = [pos[0] + dir[0] * best, pos[1] + dir[1] * best, pos[2] + dir[2] * best]; t += best / c;
const s = surf[hit];
let any = 0;
for (let b = 0; b < nb; b++) { en[b] *= 1 - s.alpha[b]; any += en[b]; }
if (any < 1e-12 / rays) break;
let n = s.n; if (dot(n, dir) > 0) n = n.map((x) => -x);
if (rand() < s.scatter) { // diffuse: Lambert about the normal
const u1 = rand(), u2 = rand(), rho = Math.sqrt(u1), ph = 2 * Math.PI * u2;
const a = Math.abs(n[0]) < 0.9 ? [1, 0, 0] : [0, 1, 0], t1v = cross(n, a), l1 = Math.hypot(...t1v), tt1 = t1v.map((x) => x / l1), tt2 = cross(n, tt1);
dir = [0, 1, 2].map((q) => tt1[q] * rho * Math.cos(ph) + tt2[q] * rho * Math.sin(ph) + n[q] * Math.sqrt(1 - u1));
} else { const k = 2 * dot(dir, n); dir = [dir[0] - k * n[0], dir[1] - k * n[1], dir[2] - k * n[2]]; }
}
}
return { sr, bands, method: 'ray-tracing', rays, radius, direct: Math.hypot(source[0] - listener[0], source[1] - listener[1], source[2] - listener[2]) };
}EDT, T20, T30, C50, C80, D50 and Ts per octave band, from the energy in each time bin.
amplitude[i] = √e[i]; the measures are taken from the amplitudes
The square root of each 64-bit energy is stored in a Float32Array before roomMetrics reads it. Every measure this returns is therefore computed from amplitudes that carry about 7 digits, and an amplitude below the smallest 32-bit number is stored as zero. roomMetrics copies its input into a Float64Array, so a 64-bit array could be passed as it is.
/** EDT, T20, T30, C50, C80, D50, Ts per band from an echogram (the square root of its energy is the amplitude the measures read). */
export function echogramMetrics(eg) {
return eg.bands.map((b) => {
const amp = new Float32Array(b.e.length);
for (let i = 0; i < amp.length; i++) amp[i] = Math.sqrt(b.e[i]);
return { f: b.f, ...roomMetrics(amp, eg.sr, null) };
});
}One step stores or exchanges 32-bit values that reach a number a person reads. It is listed here with its open item, and the step is either moved to 64-bit or the app will warn where the number is read.