Reads a recorded or computed impulse response and works out reverberation times, early decay time, clarity, definition and center time.
How the precision is kept
A WAV file's samples are read into a Float32Array. Samples of 8, 16 and 24 bits and 32-bit float samples fit a 32-bit float exactly; 32-bit integer samples are cut to a 24-bit significand. Everything after that is 64-bit: the squared response, the backward integral and its decibel curve are Float64Array, and decay times come from a least-squares line.
Open item. Read WAV samples into a Float64Array so that 32-bit integer recordings keep all their bits.
Calculations under this policy
4 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.
metrics.jsBimCraft 0.8.1 32-bit in part4 calculationsReviewed
TestsNone found by the review
Reading WAV samples
32-bit in part
The first channel of a WAV file as numbers between −1 and 1.
x = sample ÷ 2^(bits − 1)
Samples are divided by 128, 32768, 8388608 or 2147483648, all powers of two, and stored in a Float32Array. That is exact for 8-, 16- and 24-bit samples and for float files; a 32-bit integer sample is rounded to 24 significant bits.
metrics.js
lines 25 to 49
BimCraft 0.8.1
export function parseWav(bytes) {
const b = bytes instanceof Uint8Array ? bytes : new Uint8Array(bytes), dv = new DataView(b.buffer, b.byteOffset, b.byteLength);
const tag = (o) => String.fromCharCode(b[o], b[o + 1], b[o + 2], b[o + 3]);
if (b.length < 44 || tag(0) !== 'RIFF' || tag(8) !== 'WAVE') throw new Error('not a WAV file');
let fmt = null, data = null;
for (let o = 12; o + 8 <= b.length;) {
const id = tag(o), size = dv.getUint32(o + 4, true), body = o + 8;
if (id === 'fmt ') fmt = { code: dv.getUint16(body, true) === 0xFFFE && size >= 40 ? dv.getUint16(body + 24, true) : dv.getUint16(body, true), channels: dv.getUint16(body + 2, true), sr: dv.getUint32(body + 4, true), bits: dv.getUint16(body + 14, true) };
else if (id === 'data') { data = { at: body, size: Math.min(size, b.length - body) }; break; }
o = body + size + (size & 1);
}
if (!fmt || !data) throw new Error('the WAV file has no format or data');
const bytesPer = fmt.bits / 8, frame = bytesPer * fmt.channels, n = Math.floor(data.size / frame);
if (![8, 16, 24, 32].includes(fmt.bits) || !fmt.channels || !fmt.sr || n < 2) throw new Error('this WAV layout is not supported');
const x = new Float32Array(n);
for (let i = 0; i < n; i++) {
const o = data.at + i * frame;
if (fmt.code === 3 && fmt.bits === 32) x[i] = dv.getFloat32(o, true);
else if (fmt.bits === 8) x[i] = (b[o] - 128) / 128;
else if (fmt.bits === 16) x[i] = dv.getInt16(o, true) / 32768;
else if (fmt.bits === 24) x[i] = (((b[o + 2] << 24) | (b[o + 1] << 16) | (b[o] << 8)) >> 8) / 8388608;
else x[i] = dv.getInt32(o, true) / 2147483648;
}
return { sr: fmt.sr, x, channels: fmt.channels };
}
Decay time from a line fit
64-bit
A decay time from the slope of the least-squares line through part of the decay curve.
slope = (m Σ t d − Σ t Σ d) ÷ (m Σ t² − (Σ t)²); T = 60 ÷ −slope
The slope is (m Σtd − Σt Σd) ÷ (m Σt² − (Σt)²) from running sums, a form that subtracts two products of similar size. Times are counted from the direct sound, which keeps them small. Fewer than 8 points or a rising line gives no value.
metrics.js
lines 51 to 61
BimCraft 0.8.1
function fitTime(db, sr, from, to, scale) {
let sx = 0, sy = 0, sxx = 0, sxy = 0, m = 0;
for (let i = 0; i < db.length; i++) {
const d = db[i];
if (d > from || d < to) continue;
const t = i / sr; sx += t; sy += d; sxx += t * t; sxy += t * d; m++;
}
if (m < 8) return null;
const slope = (m * sxy - sx * sy) / (m * sxx - sx * sx);
return slope < 0 ? scale / -slope : null;
}
Decay curve, clarity, definition and center time
64-bit
The energy decay curve and the measures read from it.
curve = 10 lg(∫ₜ p² ÷ ∫₀ p²); C = 10 lg(early ÷ late); D50 = early ÷ total; Ts = Σ t p² ÷ Σ p²
The backward integral adds the tail first, so small late values are summed before large early ones. Clarity takes the late energy as the total minus the early energy, a subtraction that loses digits when almost all the energy is early. The 50 and 80 ms limits are rounded to whole samples.
metrics.js
lines 63 to 95
BimCraft 0.8.1
export function roomMetrics(ir, sr, band = null, { noise = 'none' } = {}) {
const x = band === null ? Float64Array.from(ir) : bandFilter(ir, sr, band, true), n = x.length;
const e = new Float64Array(n);
let peak = 0, at = 0;
for (let i = 0; i < n; i++) { e[i] = x[i] * x[i]; if (e[i] > peak) { peak = e[i]; at = i; } }
if (!(peak > 0)) return null;
let t0 = 0;
while (t0 < at && e[t0] < peak * 0.01) t0++; // within 20 dB of the peak
let end = n, floor = 0;
if (noise === 'auto') {
const k = Math.max(8, Math.floor(n * 0.1));
for (let i = n - k; i < n; i++) floor += e[i];
floor /= k;
const w = Math.max(8, Math.round(0.01 * sr));
let acc = 0;
for (let i = at; i < n; i++) { acc += e[i]; if (i - w >= at) acc -= e[i - w]; if (i - at >= w && acc / w <= 2 * floor) { end = i - w; break; } }
end = Math.max(end, at + w);
}
const p = new Float64Array(end - t0);
for (let i = t0; i < end; i++) p[i - t0] = Math.max(0, e[i] - floor);
const m = p.length, cum = new Float64Array(m);
let s = 0;
for (let i = m - 1; i >= 0; i--) { s += p[i]; cum[i] = s; }
if (!(cum[0] > 0)) return null;
const db = new Float64Array(m);
for (let i = 0; i < m; i++) db[i] = 10 * Math.log10(cum[i] / cum[0] + 1e-30);
const upTo = (ms) => Math.min(m, Math.round((ms / 1000) * sr));
const early = (ms) => { let a = 0; const k = upTo(ms); for (let i = 0; i < k; i++) a += p[i]; return a; };
const tot = cum[0], c = (ms) => { const a = early(ms), l = tot - a; return l > 0 ? 10 * Math.log10(a / l) : null; };
let st = 0;
for (let i = 0; i < m; i++) st += (i / sr) * p[i];
return { edt: fitTime(db, sr, 0, -10, 60), t20: fitTime(db, sr, -5, -25, 60), t30: fitTime(db, sr, -5, -35, 60), c50: c(50), c80: c(80), d50: early(50) / tot, ts: st / tot, t0: t0 / sr };
}
Predicted against measured
64-bit
The difference between a predicted and a measured value in each band, against the just-noticeable difference.
within when |predicted − measured| ≤ JND
One subtraction per value. The threshold is 5 % of the measured value for the decay times and a fixed 1 dB, 0.05 or 0.010 s for the others, compared with no tolerance.
metrics.js
lines 101 to 111
BimCraft 0.8.1
/** Predicted against measured, per octave band and measure; `within` = the difference is under one just-noticeable difference. */
export function compare(pred, meas) {
const out = [];
for (let i = 0; i < pred.bands.length; i++) for (const key of Object.keys(JND)) {
const a = pred.bands[i]?.[key], b = meas.bands[i]?.[key];
if (a == null || b == null) continue;
const j = JND[key], jnd = j.rel ? j.rel * b : j.abs, diff = a - b;
out.push({ f: pred.bands[i].f, key, predicted: a, measured: b, diff, jnd, within: Math.abs(diff) <= jnd });
}
return out;
}
In detail
Method
Schroeder backward integration; decay times from a least-squares line through the curve between 0 and −10 dB (× 6), −5 and −25 dB (× 3), −5 and −35 dB (× 2), as the file's header describes from ISO 3382-1; clarity, definition and center time from energy sums after the direct sound.
Number types
Float32Array for the samples read from a WAV file; Float64Array and 64-bit numbers for the rest.
Tolerances
Direct sound: the first sample whose energy is at least 1 % of the peak. 1e-30 is added inside the logarithm. A fit needs at least 8 points. With noise handling on, the floor is the mean energy of the last 10 % and the integral stops where a 10 ms average meets twice the floor.
Limits
The 50 and 80 ms limits are rounded to the nearest sample. The speech transmission index is not computed.
Approach: 32-bit step, flagged
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.