uj-mastering-master/web/dsp.js

// SPDX-License-Identifier: Apache-2.0
//
// Metric computation, ported from `metrics.py` / `master_core.py` to run in the
// browser.
//
// Pure functions over a mono Float32Array at the file's native sample rate — no
// DOM, no Web Audio — so the same module runs inside the analysis Web Worker and
// under Node for the parity tests (tests/parity), which check every metric
// against the Python implementation on identical samples.
//
// The numerics follow the Python where it matters:
//  * running sums of squares are sequential float64 sums, like np.cumsum, but
//    sampled only at the window/block edges that are read — the same values
//    without an N-length array;
//  * K-weighting is pyloudnorm's biquad pair run through scipy.signal.lfilter's
//    direct-form-II-transposed recurrence in float64, streamed sample by sample;
//  * true peak reproduces scipy.signal.resample_poly(y, 4, 1) — the same 81-tap
//    Kaiser(beta=5) firwin filter and zero-padded alignment — then an exact
//    sliding maximum (what maximum_filter1d sampled at window centres yields);
//  * the spectrogram follows librosa.stft (centred, zero padded, periodic Hann)
//    and amplitude_to_db(ref=np.max, top_db=80), dropping the DC row.
// What remains is float32 rounding where numpy works in float32 (resampling,
// RMS, spectrogram dB): far below anything a plot or readout can show.

const EPS = 1e-12;             // keeps 20*log10 finite on perfect silence
const LUFS_OFFSET = -0.691;    // BS.1770
const ABS_GATE = -70.0;        // BS.1770 absolute gate
const SILENCE_FLOOR = -70.0;

const toDbfs = (v) => 20 * Math.log10(Math.max(v, EPS));

/** Number of full windows of `windowN` stepped by `hopN` over `n` samples. */
function windowCount(n, windowN, hopN) {
  return 1 + Math.floor((n - windowN) / hopN);
}

function sampleStream(x) {
  let i = 0;
  return () => x[i++];
}

function absStream(x) {
  let i = 0;
  return () => Math.abs(x[i++]);
}

function lowerBound(sorted, value) {
  let lo = 0;
  let hi = sorted.length;
  while (lo < hi) {
    const mid = (lo + hi) >>> 1;
    if (sorted[mid] < value) lo = mid + 1;
    else hi = mid;
  }
  return lo;
}

/**
 * csq(i) = sum of x[k]^2 for k < i — np.concatenate(([0], np.cumsum(x * x))) —
 * at the requested indices only.
 *
 * `next()` yields x[0], x[1], ... in order, so a filtered signal can be streamed
 * straight in. Indices past `n` read csq(n), i.e. a zero-padded tail.
 */
function prefixSquares(n, next, indexArrays) {
  let total = 0;
  for (const a of indexArrays) total += a.length;
  const marks = new Float64Array(total);
  let offset = 0;
  for (const a of indexArrays) {
    marks.set(a, offset);
    offset += a.length;
  }
  marks.sort();

  const values = new Float64Array(marks.length);
  let m = 0;
  while (m < marks.length && marks[m] <= 0) values[m++] = 0;
  let acc = 0;
  for (let k = 0; k < n && m < marks.length; k++) {
    const v = next();
    acc += v * v;
    while (m < marks.length && marks[m] === k + 1) values[m++] = acc;
  }
  while (m < marks.length) values[m++] = acc;
  return (index) => values[lowerBound(marks, index)];
}

/**
 * max(x[s .. s + windowN)) for s = 0, hopN, 2*hopN, ... (`count` windows).
 *
 * A monotonic deque: one O(N) pass, O(windowN + hopN) memory, exact maxima.
 */
function windowMaxima(count, windowN, hopN, next) {
  const out = new Float64Array(Math.max(count, 0));
  if (count <= 0) return out;
  const cap = windowN + hopN + 1;
  const idx = new Float64Array(cap);
  const val = new Float64Array(cap);
  let head = 0;
  let size = 0;
  let w = 0;
  let end = windowN;                         // exclusive end of window w
  const last = (count - 1) * hopN + windowN;

  for (let k = 0; k < last; k++) {
    const v = next();
    while (size > 0) {
      let back = head + size - 1;
      if (back >= cap) back -= cap;
      if (val[back] <= v) size--;
      else break;
    }
    let pos = head + size;
    if (pos >= cap) pos -= cap;
    idx[pos] = k;
    val[pos] = v;
    size++;

    if (k + 1 === end) {
      const start = end - windowN;
      while (idx[head] < start) {
        head++;
        if (head === cap) head = 0;
        size--;
      }
      out[w++] = val[head];
      if (w === count) break;
      end += hopN;
    }
  }
  return out;
}

// ---- K-weighting and loudness ------------------------------------------------

/** pyloudnorm IIRfilter.generate_coefficients, normalised by a0. */
function rbjBiquad(type, G, Q, fc, rate) {
  const A = 10 ** (G / 40.0);
  const w0 = 2.0 * Math.PI * (fc / rate);
  const alpha = Math.sin(w0) / (2.0 * Q);
  let b0, b1, b2, a0, a1, a2;
  if (type === 'high_shelf') {
    b0 = A * (A + 1 + (A - 1) * Math.cos(w0) + 2 * Math.sqrt(A) * alpha);
    b1 = -2 * A * (A - 1 + (A + 1) * Math.cos(w0));
    b2 = A * (A + 1 + (A - 1) * Math.cos(w0) - 2 * Math.sqrt(A) * alpha);
    a0 = A + 1 - (A - 1) * Math.cos(w0) + 2 * Math.sqrt(A) * alpha;
    a1 = 2 * (A - 1 - (A + 1) * Math.cos(w0));
    a2 = A + 1 - (A - 1) * Math.cos(w0) - 2 * Math.sqrt(A) * alpha;
  } else if (type === 'high_pass') {
    b0 = (1 + Math.cos(w0)) / 2;
    b1 = -(1 + Math.cos(w0));
    b2 = (1 + Math.cos(w0)) / 2;
    a0 = 1 + alpha;
    a1 = -2 * Math.cos(w0);
    a2 = 1 - alpha;
  } else {
    throw new Error(`unsupported filter type ${type}`);
  }
  return { b0: b0 / a0, b1: b1 / a0, b2: b2 / a0, a1: a1 / a0, a2: a2 / a0 };
}

/** pyloudnorm Meter(rate) 'K-weighting': high shelf, then high pass. */
export function kWeightingCoefficients(rate) {
  return [
    rbjBiquad('high_shelf', 4.0, 1 / Math.sqrt(2), 1500.0, rate),
    rbjBiquad('high_pass', 0.0, 0.5, 38.0, rate),
  ];
}

/**
 * The K-weighted signal as a stream (metrics._kweight).
 *
 * Each stage is scipy's lfilter recurrence — y = z0 + b0*x;
 * z0 = z1 + x*b1 - y*a1; z1 = x*b2 - y*a2 — in the same operation order, so the
 * filtered samples match pyloudnorm's whole-array filtering.
 */
function kWeightedStream(mono, rate) {
  const [s, p] = kWeightingCoefficients(rate);
  let i = 0;
  let sz0 = 0, sz1 = 0, pz0 = 0, pz1 = 0;
  return () => {
    const x = mono[i++];
    const y1 = sz0 + s.b0 * x;
    sz0 = sz1 + x * s.b1 - y1 * s.a1;
    sz1 = x * s.b2 - y1 * s.a2;
    const y2 = pz0 + p.b0 * y1;
    pz0 = pz1 + y1 * p.b1 - y2 * p.a1;
    pz1 = y1 * p.b2 - y2 * p.a2;
    return y2;
  };
}

/** np.round: round half to even. */
function roundHalfEven(x) {
  const r = Math.round(x);
  return Math.abs(x % 1) === 0.5 && r % 2 !== 0 ? r - 1 : r;
}

/** Block edges as metrics._block_loudness lays them out, or null when none fit. */
function blockEdges(nSignal, sr, blockS, stepPct) {
  const T = nSignal / sr;
  const nBlocks = Math.trunc(roundHalfEven((T - blockS) / (blockS * stepPct)) + 1);
  if (!(nBlocks >= 1)) return null;
  const lo = new Float64Array(nBlocks);
  const up = new Float64Array(nBlocks);
  for (let j = 0; j < nBlocks; j++) {
    lo[j] = Math.trunc(blockS * (j * stepPct) * sr);
    up[j] = Math.min(Math.trunc(blockS * (j * stepPct + 1) * sr), nSignal);
  }
  return { lo, up, n: nBlocks, divisor: blockS * sr };
}

function blockLoudness(blocks, csq) {
  const z = new Float64Array(blocks.n);
  const l = new Float64Array(blocks.n);
  for (let j = 0; j < blocks.n; j++) {
    z[j] = (csq(blocks.up[j]) - csq(blocks.lo[j])) / blocks.divisor;
    l[j] = LUFS_OFFSET + 10.0 * Math.log10(z[j]);
  }
  return { z, l };
}

/** metrics._integrated_lufs: two-stage gated mean of the 400 ms blocks. */
function gatedIntegrated(z, l) {
  let sum = 0;
  let count = 0;
  for (let j = 0; j < l.length; j++) {
    if (l[j] >= ABS_GATE) { sum += z[j]; count++; }
  }
  if (!count) return -Infinity;
  const gammaR = LUFS_OFFSET + 10.0 * Math.log10(sum / count) - 10.0;
  sum = 0;
  count = 0;
  for (let j = 0; j < l.length; j++) {
    if (l[j] > gammaR && l[j] > ABS_GATE) { sum += z[j]; count++; }
  }
  if (!count) return -Infinity;
  return LUFS_OFFSET + 10.0 * Math.log10(sum / count);
}

/** np.percentile(values, q), default 'linear' method, on a sorted array. */
function percentileSorted(sorted, q) {
  const index = (q / 100) * (sorted.length - 1);
  const lo = Math.floor(index);
  const hi = Math.min(lo + 1, sorted.length - 1);
  const t = index - lo;
  const a = sorted[lo];
  const b = sorted[hi];
  return t >= 0.5 ? b - (b - a) * (1 - t) : a + (b - a) * t;
}

/** metrics._loudness_range: EBU Tech 3342 LRA from 3 s block loudness. */
function loudnessRange(l) {
  const absGated = [];
  for (let j = 0; j < l.length; j++) if (l[j] >= ABS_GATE) absGated.push(l[j]);
  if (!absGated.length) return NaN;
  let energy = 0;
  for (const v of absGated) energy += 10.0 ** (v / 10.0);
  const stlIntegrated = 10.0 * Math.log10(energy / absGated.length);
  const rel = Float64Array.from(absGated.filter((v) => v >= stlIntegrated - 20.0));
  if (!rel.length) return NaN;
  rel.sort();
  return percentileSorted(rel, 95) - percentileSorted(rel, 10);
}

/**
 * Every loudness quantity from one K-weighting pass: the ungated short-term
 * (3 s / 0.5 s) series, integrated loudness and LRA.
 *
 * Memoised on `memo` (the AudioFile attribute caching in Python) so LUFS and PSR
 * share it. Returned arrays are owned by the memo — callers copy before handing
 * them out.
 */
function loudness(mono, sr, memo) {
  if (memo && memo.loudness) return memo.loudness;
  const n = mono.length;
  const windowN = Math.max(Math.trunc(3.0 * sr), 1);
  const hopN = Math.max(Math.trunc(0.5 * sr), 1);

  const stCount = n >= windowN ? windowCount(n, windowN, hopN) : 0;
  const stStarts = new Float64Array(stCount);
  const stEnds = new Float64Array(stCount);
  for (let i = 0; i < stCount; i++) {
    stStarts[i] = i * hopN;
    stEnds[i] = i * hopN + windowN;
  }
  const integ = blockEdges(n, sr, 0.4, 0.25);
  // LRA appends 1.5 s of silence; indices past n read csq(n), which is exactly that.
  const lraBlocks = blockEdges(n + Math.trunc(1.5 * sr), sr, 3.0, 0.03);

  const wanted = [stStarts, stEnds, Float64Array.of(0, n)];
  if (integ) wanted.push(integ.lo, integ.up);
  if (lraBlocks) wanted.push(lraBlocks.lo, lraBlocks.up);
  const csq = prefixSquares(n, kWeightedStream(mono, sr), wanted);

  let times;
  let lufs;
  if (n < windowN) {
    const ms = n ? csq(n) / n : 0;
    times = Float64Array.of(n / (2.0 * sr));
    lufs = Float64Array.of(LUFS_OFFSET + 10.0 * Math.log10(Math.max(ms, EPS)));
  } else {
    times = new Float64Array(stCount);
    lufs = new Float64Array(stCount);
    for (let i = 0; i < stCount; i++) {
      const ms = (csq(stEnds[i]) - csq(stStarts[i])) / windowN;
      lufs[i] = LUFS_OFFSET + 10.0 * Math.log10(Math.max(ms, EPS));
      times[i] = (stStarts[i] + windowN / 2.0) / sr;
    }
  }

  let integrated = -Infinity;
  if (integ) {
    const { z, l } = blockLoudness(integ, csq);
    integrated = gatedIntegrated(z, l);
  }

  let lra = NaN;
  if (n >= Math.trunc(3.0 * sr) && lraBlocks) {
    lra = loudnessRange(blockLoudness(lraBlocks, csq).l);
  }

  const result = { times, lufs, integrated, lra };
  if (memo) memo.loudness = result;
  return result;
}

// ---- metrics -------------------------------------------------------------------

/** RMSPowerMetric: librosa.feature.rms, 10 s frames / 2 s hop, centred. */
export function computeRmsPower(mono, sr, window = 10, hop = 2) {
  const n = mono.length;
  const frameLength = window * sr;
  const hopLength = hop * sr;
  const half = Math.floor(frameLength / 2);
  const padded = n + 2 * half;
  const count = padded >= frameLength ? 1 + Math.floor((padded - frameLength) / hopLength) : 0;

  // Frame i covers padded[i*hop, i*hop + frameLength): the zero padding adds
  // nothing to the sum but still counts in the mean, as in librosa.
  const starts = new Float64Array(count);
  const ends = new Float64Array(count);
  for (let i = 0; i < count; i++) {
    const a = i * hopLength - half;
    starts[i] = Math.min(Math.max(a, 0), n);
    ends[i] = Math.min(Math.max(a + frameLength, 0), n);
  }
  const csq = prefixSquares(n, sampleStream(mono), [starts, ends]);

  const times = new Float64Array(count);
  const rms = new Float64Array(count);
  for (let i = 0; i < count; i++) {
    times[i] = (i * hopLength) / sr;
    rms[i] = Math.sqrt((csq(ends[i]) - csq(starts[i])) / frameLength);
  }
  return { times, rms };
}

/** WaveformMetric: min/max envelope over ~4000 columns. */
export function computeWaveform(mono, sr, targetColumns = 4000) {
  const n = mono.length;
  if (n <= targetColumns) {
    const times = new Float64Array(n);
    for (let i = 0; i < n; i++) times[i] = i / sr;
    return { times, lo: Float32Array.from(mono), hi: Float32Array.from(mono) };
  }
  const chunk = Math.floor(n / targetColumns);
  const times = new Float64Array(targetColumns);
  const lo = new Float32Array(targetColumns);
  const hi = new Float32Array(targetColumns);
  for (let c = 0; c < targetColumns; c++) {
    let mn = Infinity;
    let mx = -Infinity;
    const base = c * chunk;
    for (let k = base; k < base + chunk; k++) {
      const v = mono[k];
      if (v < mn) mn = v;
      if (v > mx) mx = v;
    }
    lo[c] = mn;
    hi[c] = mx;
    times[c] = (c * chunk + chunk / 2) / sr;
  }
  return { times, lo, hi };
}

/** LUFSMetric: short-term series clipped to [-70, 0], plus integrated and LRA. */
export function computeLufs(mono, sr, memo) {
  const L = loudness(mono, sr, memo);
  const lufs = new Float64Array(L.lufs.length);
  for (let i = 0; i < lufs.length; i++) {
    const v = Number.isFinite(L.lufs[i]) ? L.lufs[i] : SILENCE_FLOOR;
    lufs[i] = Math.min(Math.max(v, SILENCE_FLOOR), 0.0);
  }
  return { times: L.times.slice(), lufs, integrated: L.integrated, lra: L.lra };
}

/** CrestFactorMetric: 20*log10(peak / RMS) per 1 s window, 0.25 s hop. */
export function computeCrestFactor(mono, sr) {
  const n = mono.length;
  const windowN = Math.trunc(1.0 * sr);
  const hopN = Math.trunc(0.25 * sr);

  if (n < windowN) {
    let peak = 0;
    let sumSq = 0;
    for (let i = 0; i < n; i++) {
      const v = mono[i];
      if (Math.abs(v) > peak) peak = Math.abs(v);
      sumSq += v * v;
    }
    const rms = n ? Math.sqrt(sumSq / n) : 0;
    return {
      times: Float64Array.of(n / (2.0 * sr)),
      crest_db: Float64Array.of(20.0 * Math.log10(Math.max(peak, EPS) / Math.max(rms, EPS))),
    };
  }

  const count = windowCount(n, windowN, hopN);
  const starts = new Float64Array(count);
  const ends = new Float64Array(count);
  for (let i = 0; i < count; i++) {
    starts[i] = i * hopN;
    ends[i] = i * hopN + windowN;
  }
  const csq = prefixSquares(n, sampleStream(mono), [starts, ends]);
  const peaks = windowMaxima(count, windowN, hopN, absStream(mono));

  const times = new Float64Array(count);
  const crest = new Float64Array(count);
  for (let i = 0; i < count; i++) {
    const meanSq = (csq(ends[i]) - csq(starts[i])) / windowN;
    const rms = Math.sqrt(Math.max(meanSq, EPS));
    crest[i] = 20.0 * Math.log10(Math.max(peaks[i], EPS) / rms);
    times[i] = (starts[i] + windowN / 2.0) / sr;
  }
  return { times, crest_db: crest };
}

/** PSRMetric: sample peak minus short-term LUFS (3 s), NaN below the gate. */
export function computePsr(mono, sr, memo) {
  const n = mono.length;
  const windowN = Math.max(Math.trunc(3.0 * sr), 1);
  const hopN = Math.max(Math.trunc(0.5 * sr), 1);
  const L = loudness(mono, sr, memo);
  const count = L.times.length;

  let peaksDb;
  if (n < windowN) {
    let peak = 0;
    for (let i = 0; i < n; i++) peak = Math.max(peak, Math.abs(mono[i]));
    peaksDb = Float64Array.of(n ? toDbfs(peak) : SILENCE_FLOOR);
  } else {
    peaksDb = windowMaxima(count, windowN, hopN, absStream(mono)).map(toDbfs);
  }

  const psr = new Float64Array(count);
  for (let i = 0; i < count; i++) {
    const l = L.lufs[i];
    psr[i] = Number.isFinite(l) && l > SILENCE_FLOOR ? peaksDb[i] - l : NaN;
  }
  return { times: L.times.slice(), psr };
}

// ---- true peak -------------------------------------------------------------------

const OVERSAMPLE = 4;
let truePeakTaps = null;

function besselI0(x) {
  const q = (x * x) / 4;
  let term = 1;
  let sum = 1;
  for (let k = 1; k < 500; k++) {
    term *= q / (k * k);
    sum += term;
    if (term < 1e-17 * sum) break;
  }
  return sum;
}

/**
 * The anti-imaging filter resample_poly(x, 4, 1) builds for float32 input:
 * firwin(81, 1/4, window=('kaiser', 5.0)), cast to float32, times 4.
 */
function truePeakFilter() {
  if (truePeakTaps) return truePeakTaps;
  const numtaps = 2 * 10 * OVERSAMPLE + 1;
  const cutoff = 1.0 / OVERSAMPLE;
  const beta = 5.0;
  const alpha = 0.5 * (numtaps - 1);
  const i0beta = besselI0(beta);
  const h = new Float64Array(numtaps);
  let sum = 0;
  for (let k = 0; k < numtaps; k++) {
    const m = k - alpha;
    const x = Math.PI * (m === 0 ? 1e-20 : cutoff * m);
    const sinc = Math.sin(x) / x;
    const r = (k - alpha) / alpha;
    h[k] = cutoff * sinc * (besselI0(beta * Math.sqrt(1 - r * r)) / i0beta);
    sum += h[k];
  }
  truePeakTaps = new Float32Array(numtaps);
  for (let k = 0; k < numtaps; k++) truePeakTaps[k] = Math.fround(h[k] / sum) * OVERSAMPLE;
  return truePeakTaps;
}

/**
 * |resample_poly(x, 4, 1)| as a stream, one output sample per call.
 *
 * resample_poly zero-stuffs x, convolves with the taps and keeps outputs
 * [41, 41 + 4n); output k therefore sums x[s] * h[40 + k - 4s] over the input
 * samples within the filter's reach — 20 or 21 taps per output.
 */
function upsampledAbsStream(x, h) {
  const n = x.length;
  const half = 10 * OVERSAMPLE;
  let k = 0;
  return () => {
    const q = Math.floor(k / OVERSAMPLE);
    const r = k - OVERSAMPLE * q;
    const sLo = Math.max(0, r === 0 ? q - 10 : q - 9);
    const sHi = Math.min(n - 1, q + 10);
    let acc = 0;
    for (let s = sLo; s <= sHi; s++) acc += x[s] * h[half + r + OVERSAMPLE * (q - s)];
    k++;
    return Math.abs(Math.fround(acc));
  };
}

/** TruePeakMetric: 4x oversampled dBTP per 250 ms window, 100 ms hop. */
export function computeTruePeak(mono, sr) {
  const n = mono.length;
  const windowN = Math.trunc(0.25 * sr);
  const hopN = Math.trunc(0.1 * sr);
  const h = truePeakFilter();

  if (n < windowN) {
    let peak = 0;
    const next = upsampledAbsStream(mono, h);
    for (let k = 0; k < n * OVERSAMPLE; k++) peak = Math.max(peak, next());
    const db = toDbfs(peak);
    return { times: Float64Array.of(n / (2.0 * sr)), tp_db: Float64Array.of(db), integrated_tp_db: db };
  }

  const count = windowCount(n, windowN, hopN);
  const peaks = windowMaxima(count, windowN * OVERSAMPLE, hopN * OVERSAMPLE, upsampledAbsStream(mono, h));
  const times = new Float64Array(count);
  const tpDb = new Float64Array(count);
  let max = -Infinity;
  for (let i = 0; i < count; i++) {
    tpDb[i] = toDbfs(peaks[i]);
    times[i] = (i * hopN + windowN / 2.0) / sr;
    if (tpDb[i] > max) max = tpDb[i];
  }
  return { times, tp_db: tpDb, integrated_tp_db: max };
}

// ---- spectrogram -----------------------------------------------------------------

const N_FFT = 4096;
const SPEC_TARGET_COLUMNS = 4000;
const SPEC_DB_FLOOR = -80.0;
const fftPlans = new Map();

/** In-place radix-2 complex FFT for `size` (a power of two). */
function fftPlan(size) {
  if (fftPlans.has(size)) return fftPlans.get(size);
  const levels = Math.round(Math.log2(size));
  const rev = new Uint32Array(size);
  for (let i = 0; i < size; i++) {
    let r = 0;
    for (let b = 0; b < levels; b++) r |= ((i >> b) & 1) << (levels - 1 - b);
    rev[i] = r;
  }
  const cos = new Float64Array(size / 2);
  const sin = new Float64Array(size / 2);
  for (let i = 0; i < size / 2; i++) {
    cos[i] = Math.cos((2 * Math.PI * i) / size);
    sin[i] = Math.sin((2 * Math.PI * i) / size);
  }
  const plan = (re, im) => {
    for (let i = 0; i < size; i++) {
      const j = rev[i];
      if (j > i) {
        let t = re[i]; re[i] = re[j]; re[j] = t;
        t = im[i]; im[i] = im[j]; im[j] = t;
      }
    }
    for (let len = 2; len <= size; len <<= 1) {
      const halfLen = len >> 1;
      const step = size / len;
      for (let i = 0; i < size; i += len) {
        for (let j = 0; j < halfLen; j++) {
          const c = cos[j * step];
          const s = sin[j * step];
          const u = i + j;
          const v = u + halfLen;
          const tre = re[v] * c + im[v] * s;
          const tim = im[v] * c - re[v] * s;
          re[v] = re[u] - tre;
          im[v] = im[u] - tim;
          re[u] += tre;
          im[u] += tim;
        }
      }
    }
  };
  fftPlans.set(size, plan);
  return plan;
}

/**
 * SpectrogramMetric: |STFT| in dB relative to the loudest bin, floored at -80.
 *
 * `s_db` is row-major (rows x cols), row 0 the lowest non-DC bin — the same
 * layout as the Python `s_db[1:, :]`.
 */
export function computeSpectrogram(mono, sr) {
  const n = mono.length;
  const half = N_FFT / 2;
  const hop = Math.max(N_FFT / 4, Math.floor(n / SPEC_TARGET_COLUMNS));
  const cols = 1 + Math.floor(n / hop);
  const rows = half;                           // bins 1..N_FFT/2

  const window = new Float64Array(N_FFT);      // periodic Hann (fftbins=True)
  for (let j = 0; j < N_FFT; j++) window[j] = 0.5 - 0.5 * Math.cos((2 * Math.PI * j) / N_FFT);

  const fft = fftPlan(N_FFT);
  const re = new Float64Array(N_FFT);
  const im = new Float64Array(N_FFT);
  const sDb = new Float32Array(rows * cols);
  let maxPower = 0;

  for (let c = 0; c < cols; c++) {
    // Centred framing: frame c is padded[c*hop, c*hop + N_FFT) with N_FFT/2
    // zeros either side of the signal.
    const offset = c * hop - half;
    for (let j = 0; j < N_FFT; j++) {
      const t = offset + j;
      re[j] = t >= 0 && t < n ? mono[t] * window[j] : 0;
      im[j] = 0;
    }
    fft(re, im);
    const p0 = re[0] * re[0] + im[0] * im[0];
    if (p0 > maxPower) maxPower = p0;
    for (let b = 1; b <= half; b++) {
      const p = re[b] * re[b] + im[b] * im[b];
      if (p > maxPower) maxPower = p;
      sDb[(b - 1) * cols + c] = p;
    }
  }

  // amplitude_to_db(ref=np.max, amin=1e-5, top_db=80): in power terms
  // 10*log10(max(1e-10, p)) - 10*log10(max(1e-10, max p)), then clipped to
  // (maximum - 80). The maximum of that is 0 at the reference bin, so the clip
  // is a fixed -80 floor.
  const amin = 1e-5 ** 2;
  const refDb = 10.0 * Math.log10(Math.max(amin, maxPower));
  for (let i = 0; i < sDb.length; i++) {
    sDb[i] = Math.max(10.0 * Math.log10(Math.max(amin, sDb[i])) - refDb, SPEC_DB_FLOOR);
  }

  const binHz = 1.0 / (N_FFT * (1.0 / sr));    // np.fft.rfftfreq spacing
  const freqs = new Float64Array(rows);
  for (let b = 0; b < rows; b++) freqs[b] = (b + 1) * binHz;
  const times = new Float64Array(cols);
  for (let c = 0; c < cols; c++) times[c] = (c * hop + half) / sr;

  return { freqs, times, s_db: sDb, rows, cols, nyquist: sr / 2.0 };
}

// ---- registry ----------------------------------------------------------------------

/** metric id -> compute(mono, sampleRate, memo). `memo` is per track. */
export const COMPUTE = {
  rms_power: (mono, sr) => computeRmsPower(mono, sr),
  waveform: (mono, sr) => computeWaveform(mono, sr),
  lufs: (mono, sr, memo) => computeLufs(mono, sr, memo),
  crest_factor: (mono, sr) => computeCrestFactor(mono, sr),
  psr: (mono, sr, memo) => computePsr(mono, sr, memo),
  true_peak: (mono, sr) => computeTruePeak(mono, sr),
  spectrogram: (mono, sr) => computeSpectrogram(mono, sr),
};

/** The ArrayBuffers inside a computed result, for a zero-copy postMessage. */
export function transferables(data) {
  const buffers = new Set();
  for (const value of Object.values(data)) {
    if (ArrayBuffer.isView(value)) buffers.add(value.buffer);
  }
  return [...buffers];
}