STIPA analyser

src/App.jsx

@@ -20,6 +20,9 @@ import {
   stiRatingLabel,
   stiRatingColor,
   parseRewRT60,
+  generateStipaSignal,
+  encodeWAV,
+  computeSTIPAfromRecording,
 } from './sti';
 
 // ─── UI PRIMITIVES ────────────────────────────────────────────────────────────
@@ -490,6 +493,12 @@ function AnalysisTab({ hall, speakers, settings, addLog }) {
   const [streamText,     setStreamText]     = useState('');
   const [analysis,       setAnalysis]       = useState(null);
   const [stipaMeasured,  setStipaMeasured]  = useState(null);
+  const [stipaGenStatus, setStipaGenStatus] = useState('idle');
+  const [stipaDuration,  setStipaDuration]  = useState(18);
+  const [stipaAnalysing, setStipaAnalysing] = useState(false);
+  const [stipaResult,    setStipaResult]    = useState(null);
+  const [stipaFile,      setStipaFile]      = useState(null);
+  const stipaFileRef = useRef(null);
 
   const rewBase  = `${settings.rewHost}:${settings.rewPort}`;
   const rowData  = computeHallGeometry(hall);
@@ -578,6 +587,49 @@ function AnalysisTab({ hall, speakers, settings, addLog }) {
     }
   }
 
+  async function generateAndDownload() {
+    setStipaGenStatus('generating');
+    addLog(`Generating ${stipaDuration}s STIPA test signal (IEC 60268-16)…`);
+    try {
+      const buf = await generateStipaSignal(stipaDuration, 48000);
+      const wav = encodeWAV(buf);
+      const url = URL.createObjectURL(wav);
+      const a   = document.createElement('a');
+      a.href     = url;
+      a.download = `stipa-signal-${stipaDuration}s.wav`;
+      document.body.appendChild(a);
+      a.click();
+      document.body.removeChild(a);
+      URL.revokeObjectURL(url);
+      setStipaGenStatus('done');
+      addLog('STIPA signal downloaded ✓ — play at full PA level and record', 'good');
+    } catch (e) {
+      setStipaGenStatus('error');
+      addLog(`STIPA generator: ${e.message}`, 'error');
+    }
+  }
+
+  async function analyseRecording(file) {
+    setStipaAnalysing(true);
+    setStipaResult(null);
+    addLog(`Analysing STIPA recording: ${file.name}…`);
+    try {
+      const arrayBuf = await file.arrayBuffer();
+      const audioCtx = new AudioContext();
+      const decoded  = await audioCtx.decodeAudioData(arrayBuf);
+      await audioCtx.close();
+      const result = await computeSTIPAfromRecording(decoded);
+      setStipaResult(result);
+      setStipaMeasured(result.sti);
+      addLog(`STIPA result: ${result.sti.toFixed(2)} — ${result.rating} ✓`, 'good');
+      if (result.warning) addLog(`⚠ ${result.warning}`, 'warn');
+    } catch (e) {
+      addLog(`STIPA analyser: ${e.message}`, 'error');
+    } finally {
+      setStipaAnalysing(false);
+    }
+  }
+
   function buildPrompt() {
     const flutter  = computeFlutter(hall);
     const row8Clr  = hall.ceilFlat - hall.dishRise;
@@ -862,6 +914,76 @@ ARTA window size per row (limited by ceiling reflection time quoted above).`;
             </div>
           ))}
         </div>
+
+        <Divider />
+        <div style={{ display: 'grid', gridTemplateColumns: '1fr 1fr', gap: 16 }}>
+          <div>
+            <div style={{ fontSize: 10, color: C.muted, letterSpacing: 1, textTransform: 'uppercase', marginBottom: 8 }}>Generate STIPA Test Signal</div>
+            <div style={{ fontSize: 11, color: C.muted, marginBottom: 10, lineHeight: 1.6 }}>Play through the PA at full level, record on any phone, upload below.</div>
+            <div style={{ display: 'flex', gap: 8, alignItems: 'center', marginBottom: 10 }}>
+              <span style={{ fontSize: 11, color: C.muted, whiteSpace: 'nowrap' }}>Duration</span>
+              <NumIn value={stipaDuration} onChange={setStipaDuration} min={15} max={60} step={1} unit="s" style={{ width: 70 }} />
+            </div>
+            <Btn onClick={generateAndDownload} disabled={stipaGenStatus === 'generating'} color={C.accentDim}>
+              {stipaGenStatus === 'generating' ? 'Generating…' : '⬇ Download STIPA WAV'}
+            </Btn>
+            {stipaGenStatus === 'done'  && <div style={{ fontSize: 11, color: C.good, marginTop: 6 }}>✓ Downloaded — play through PA and record</div>}
+            {stipaGenStatus === 'error' && <div style={{ fontSize: 11, color: C.danger, marginTop: 6 }}>⚠ Generation failed — check browser console</div>}
+          </div>
+          <div>
+            <div style={{ fontSize: 10, color: C.muted, letterSpacing: 1, textTransform: 'uppercase', marginBottom: 8 }}>Analyse Recording</div>
+            <div style={{ fontSize: 11, color: C.muted, marginBottom: 10, lineHeight: 1.6 }}>Upload the WAV recorded through the hall. Result auto-populates STIPA value above.</div>
+            <input ref={stipaFileRef} type="file" accept=".wav,audio/*" style={{ display: 'none' }}
+              onChange={e => { const f = e.target.files?.[0]; if (f) setStipaFile(f); }} />
+            <div style={{ display: 'flex', gap: 8, flexWrap: 'wrap', alignItems: 'center' }}>
+              <Btn onClick={() => stipaFileRef.current?.click()} color={C.accentDim}>Choose WAV…</Btn>
+              {stipaFile && (
+                <Btn onClick={() => analyseRecording(stipaFile)} disabled={stipaAnalysing} color={C.good}>
+                  {stipaAnalysing ? 'Analysing…' : '▶ Analyse'}
+                </Btn>
+              )}
+            </div>
+            {stipaFile && <div style={{ fontSize: 11, color: C.muted, marginTop: 6 }}>{stipaFile.name}</div>}
+          </div>
+        </div>
+
+        {stipaResult && (
+          <div style={{ marginTop: 14, borderTop: `1px solid ${C.border}`, paddingTop: 14 }}>
+            {stipaResult.warning && (
+              <div style={{ color: C.gold, fontSize: 11, marginBottom: 8 }}>⚠ {stipaResult.warning}</div>
+            )}
+            <div style={{ display: 'flex', alignItems: 'center', gap: 10, marginBottom: 10 }}>
+              <div style={{ width: 12, height: 12, borderRadius: '50%', background: stiRatingColor(stipaResult.sti) }} />
+              <span style={{ color: stiRatingColor(stipaResult.sti), fontWeight: 700, fontSize: 15 }}>{stipaResult.rating}</span>
+              <span style={{ color: C.heading, fontSize: 18, fontFamily: 'monospace' }}>{stipaResult.sti.toFixed(3)}</span>
+              <span style={{ color: C.muted, fontSize: 11 }}>STIPA measured (IEC 60268-16)</span>
+              {stipaResult.sti < 0.75 && (
+                <span style={{ color: C.danger, fontSize: 11 }}>▲ {(0.75 - stipaResult.sti).toFixed(3)} below Excellent</span>
+              )}
+            </div>
+            <div style={{ fontSize: 10, color: C.muted, marginBottom: 6, textTransform: 'uppercase', letterSpacing: 1 }}>
+              Modulation depth mk per octave band
+            </div>
+            <div style={{ display: 'grid', gridTemplateColumns: 'repeat(7, 1fr)', gap: 4 }}>
+              {STI_BANDS_HZ.map((f, k) => (
+                <div key={f} style={{ background: C.bgAlt, borderRadius: 6, padding: '6px 4px', textAlign: 'center' }}>
+                  <div style={{ fontSize: 10, color: C.muted, marginBottom: 3 }}>{f >= 1000 ? `${f / 1000}k` : f} Hz</div>
+                  {stipaResult.mk[k].map((m, j) => (
+                    <div key={j} style={{ fontSize: 11, fontFamily: 'monospace', color: m > 1.3 ? C.danger : m > 0.6 ? C.good : m > 0.3 ? C.gold : C.danger }}>
+                      {m.toFixed(3)}
+                    </div>
+                  ))}
+                  <div style={{ fontSize: 10, color: stiRatingColor(stipaResult.tiPerBand[k]), marginTop: 3, fontWeight: 600 }}>
+                    TI {stipaResult.tiPerBand[k].toFixed(2)}
+                  </div>
+                </div>
+              ))}
+            </div>
+            <div style={{ fontSize: 10, color: C.muted, marginTop: 4 }}>
+              Two mk rows per band (fm1 / fm2). mk &gt; 1.3 indicates impulsive noise ⚠. TI = transmission index per IEC 60268-16.
+            </div>
+          </div>
+        )}
       </Panel>
 
       <Btn onClick={runAnalysis} disabled={analysisState === 'loading' || analysisState === 'streaming'}

src/sti.js

@@ -27,14 +27,16 @@ const MOD_FREQS_FULL = [
 // STIPA: 2 modulation frequencies per octave band (IEC 60268-16:2020)
 // Row = octave band [125, 250, 500, 1k, 2k, 4k, 8k Hz]
 // Source: zawi01/stipa, confirmed against IEC standard Table C.1
+// Source: zawi01/stipa stipa.m fm matrix, read column-by-column (k=1:7 bands)
+// fm = [1.6, 1, 0.63, 2, 1.25, 0.8, 2.5; 8, 5, 3.15, 10, 6.25, 4, 12.5]
 export const STIPA_MOD_FREQS = [
-  [1.00, 5.00],   // 125 Hz
-  [0.63, 3.15],   // 250 Hz
-  [2.00, 10.00],  // 500 Hz
-  [1.25, 8.00],   // 1 kHz
-  [0.80, 4.00],   // 2 kHz
-  [2.50, 12.50],  // 4 kHz
-  [6.30, 1.60],   // 8 kHz
+  [1.60, 8.00],   // 125 Hz
+  [1.00, 5.00],   // 250 Hz
+  [0.63, 3.15],   // 500 Hz
+  [2.00, 10.00],  // 1 kHz
+  [1.25, 6.25],   // 2 kHz  (6.25, not 6.3)
+  [0.80, 4.00],   // 4 kHz
+  [2.50, 12.50],  // 8 kHz
 ];
 
 // STIPA test-signal band levels relative to 0 dBFS (Revision 5, IEC 60268-16)
@@ -142,6 +144,330 @@ export function stiRatingColor(v) {
   return '#f87171';
 }
 
+// ─── STIPA SIGNAL GENERATOR (IEC 60268-16:2020) ──────────────────────────────
+//
+// Port of zawi01/stipa generateStipaSignal.m
+// Generates the standard STIPA test signal as an AudioBuffer.
+// Pink noise carrier, amplitude-modulated per octave band at 2 modulation
+// frequencies, with band-level scaling (Revision 5).
+//
+// NOTE: Uses the Web Audio API — call only in browser context.
+
+/**
+ * Generate the STIPA test signal and return an AudioBuffer.
+ * Caller can play via AudioContext or export as WAV via encodeWAV().
+ *
+ * @param {number} duration   Signal length in seconds (15–25 recommended, default 18)
+ * @param {number} sampleRate Audio sample rate in Hz (min 22050, default 48000)
+ * @returns {Promise<AudioBuffer>}
+ */
+export async function generateStipaSignal(duration = 18, sampleRate = 48000) {
+  const N   = Math.round(duration * sampleRate);
+  const ctx = new OfflineAudioContext(1, N, sampleRate);
+
+  // ── Pink noise via Paul Kellett's approximation (voss-mcartney filtered white noise)
+  const whiteBuffer = ctx.createBuffer(1, N, sampleRate);
+  const whiteData   = whiteBuffer.getChannelData(0);
+  // Approximate pink noise using sum of white noise sources at octave intervals
+  let b0 = 0, b1 = 0, b2 = 0, b3 = 0, b4 = 0, b5 = 0, b6 = 0;
+  for (let i = 0; i < N; i++) {
+    const white = Math.random() * 2 - 1;
+    b0 = 0.99886 * b0 + white * 0.0555179;
+    b1 = 0.99332 * b1 + white * 0.0750759;
+    b2 = 0.96900 * b2 + white * 0.1538520;
+    b3 = 0.86650 * b3 + white * 0.3104856;
+    b4 = 0.55000 * b4 + white * 0.5329522;
+    b5 = -0.7616 * b5 - white * 0.0168980;
+    whiteData[i] = (b0 + b1 + b2 + b3 + b4 + b5 + b6 + white * 0.5362) / 7;
+    b6 = white * 0.115926;
+  }
+
+  // ── Per-band: octave bandpass filter → amplitude modulate → apply gain
+  const BAND_GAINS = STIPA_BAND_LEVELS_DB.map(db => Math.pow(10, db / 20));
+  const REF_DEPTH  = 0.55; // IEC standard reference modulation depth
+
+  let outputData = new Float32Array(N);
+
+  for (let k = 0; k < 7; k++) {
+    const fc    = STI_BANDS_HZ[k];
+    const [fm1, fm2] = STIPA_MOD_FREQS[k];
+
+    // Single-pole bandpass approximation around each octave centre
+    // Q = √2 for 1-octave bandwidth (2nd-order Butterworth bandpass)
+    const bandSrc   = ctx.createBufferSource();
+    bandSrc.buffer  = whiteBuffer;
+    const bpFilter  = ctx.createBiquadFilter();
+    bpFilter.type      = 'bandpass';
+    bpFilter.frequency.value = fc;
+    bpFilter.Q.value   = Math.SQRT2;
+
+    const gainNode  = ctx.createGain();
+    gainNode.gain.value = BAND_GAINS[k];
+
+    bandSrc.connect(bpFilter);
+    bpFilter.connect(gainNode);
+    gainNode.connect(ctx.destination);
+    bandSrc.start(0);
+
+    // We need the filtered samples to apply time-varying AM modulation.
+    // Web Audio API doesn't expose per-sample access in OfflineAudioContext
+    // mid-graph, so we render each band separately then modulate in JS.
+    const bandCtx = new OfflineAudioContext(1, N, sampleRate);
+    const bandSrc2  = bandCtx.createBufferSource();
+    bandSrc2.buffer = whiteBuffer;
+    const bpFilt2   = bandCtx.createBiquadFilter();
+    bpFilt2.type       = 'bandpass';
+    bpFilt2.frequency.value = fc;
+    bpFilt2.Q.value    = Math.SQRT2;
+    bandSrc2.connect(bpFilt2);
+    bpFilt2.connect(bandCtx.destination);
+    bandSrc2.start(0);
+
+    const rendered = await bandCtx.startRendering();
+    const bandData = rendered.getChannelData(0);
+
+    // Amplitude modulation envelope:
+    // envelope = √(0.5 × (1 + 0.55 × (sin(2π·fm1·t) − sin(2π·fm2·t))))
+    for (let i = 0; i < N; i++) {
+      const t   = i / sampleRate;
+      const env = Math.sqrt(Math.max(0, 0.5 * (1 + REF_DEPTH * (
+        Math.sin(2 * Math.PI * fm1 * t) - Math.sin(2 * Math.PI * fm2 * t)
+      ))));
+      outputData[i] += bandData[i] * env * BAND_GAINS[k];
+    }
+  }
+
+  // Normalise to target RMS = 0.1
+  const rms = Math.sqrt(outputData.reduce((s, v) => s + v * v, 0) / N);
+  if (rms > 0) {
+    const scale = 0.1 / rms;
+    for (let i = 0; i < N; i++) outputData[i] *= scale;
+  }
+
+  const outBuf = ctx.createBuffer(1, N, sampleRate);
+  outBuf.copyToChannel(outputData, 0);
+  return outBuf;
+}
+
+/**
+ * Encode an AudioBuffer (mono, Float32) as a 16-bit PCM WAV Blob.
+ * @param {AudioBuffer} audioBuffer
+ * @returns {Blob}
+ */
+export function encodeWAV(audioBuffer) {
+  const samples    = audioBuffer.getChannelData(0);
+  const sampleRate = audioBuffer.sampleRate;
+  const numSamples = samples.length;
+  const buf        = new ArrayBuffer(44 + numSamples * 2);
+  const view       = new DataView(buf);
+  const write      = (off, str) => { for (let i = 0; i < str.length; i++) view.setUint8(off + i, str.charCodeAt(i)); };
+
+  write(0, 'RIFF');
+  view.setUint32(4,  36 + numSamples * 2, true);
+  write(8, 'WAVE');
+  write(12, 'fmt ');
+  view.setUint32(16, 16, true);
+  view.setUint16(20, 1,  true);                 // PCM
+  view.setUint16(22, 1,  true);                 // mono
+  view.setUint32(24, sampleRate, true);
+  view.setUint32(28, sampleRate * 2, true);     // byte rate
+  view.setUint16(32, 2,  true);                 // block align
+  view.setUint16(34, 16, true);                 // bits per sample
+  write(36, 'data');
+  view.setUint32(40, numSamples * 2, true);
+
+  let offset = 44;
+  for (let i = 0; i < numSamples; i++, offset += 2) {
+    const s = Math.max(-1, Math.min(1, samples[i]));
+    view.setInt16(offset, s < 0 ? s * 0x8000 : s * 0x7fff, true);
+  }
+  return new Blob([buf], { type: 'audio/wav' });
+}
+
+// ─── STIPA ANALYSER FROM RECORDED WAV ────────────────────────────────────────
+//
+// Port of zawi01/stipa stipa.m — direct STIPA measurement from recorded audio.
+// Feed the result of playing generateStipaSignal() through the PA system and
+// recording it (any smartphone recorder is sufficient).
+//
+// Pipeline: band-filter → discard 200ms → envelope detect → single-bin DFT
+//           per modulation freq → mk matrix → SNR → TI → STI
+
+/**
+ * Compute STIPA from a recorded AudioBuffer.
+ *
+ * @param {AudioBuffer} recordedBuffer   Recording of the STIPA signal through the system
+ * @param {{ Lsk?: number[], Lnk?: number[] }} opts
+ *   Lsk — per-band signal levels in dB (for auditory masking correction)
+ *   Lnk — per-band ambient noise levels in dB (for noise correction)
+ * @returns {Promise<{ sti: number, mk: number[][], tiPerBand: number[], rating: string, warning?: string }>}
+ */
+export async function computeSTIPAfromRecording(recordedBuffer, { Lsk, Lnk } = {}) {
+  const fs      = recordedBuffer.sampleRate;
+  const raw     = recordedBuffer.getChannelData(0);
+  const trimMs  = 200;
+  const trimSmp = Math.round((trimMs / 1000) * fs);
+
+  // ── Per-band octave filtering (cascaded biquad, 3 sections ≈ 6th-order)
+  const bandSignals = await Promise.all(STI_BANDS_HZ.map(fc => filterOctave(raw, fs, fc)));
+
+  // ── Discard first 200ms (suppress IIR transients), detect envelope
+  const envelopes = bandSignals.map(band => {
+    const trimmed = band.slice(trimSmp);
+    return envelopeDetect(trimmed, fs);
+  });
+
+  // ── MTF: single-bin DFT at each modulation frequency (from zawi01/stipa MTF fn)
+  const mk = STIPA_MOD_FREQS.map(([fm1, fm2], k) => {
+    const env     = envelopes[k];
+    const seconds = env.length / fs;
+    return [fm1, fm2].map(Fm => singleBinDFT(env, fs, Fm, seconds) / 0.55);
+  });
+
+  // ── Validity check (IEC: mk > 1.3 indicates impulsive noise)
+  const maxMk  = Math.max(...mk.flat());
+  const warning = maxMk > 1.3
+    ? `One or more m-values exceed 1.3 (max ${maxMk.toFixed(2)}) — impulsive noise or invalid recording`
+    : undefined;
+
+  // ── Optional corrections
+  let mkCorrected = mk.map(row => [...row]);
+  if (Lsk && Lnk) {
+    mkCorrected = adjustAmbientNoise(mkCorrected, Lsk, Lnk);
+    mkCorrected = adjustAuditoryMasking(mkCorrected, Lsk, Lnk);
+  } else if (Lsk) {
+    mkCorrected = adjustAuditoryMasking(mkCorrected, Lsk, Array(7).fill(0));
+  }
+
+  // ── Clamp mk to [0, 1]
+  mkCorrected = mkCorrected.map(row => row.map(v => Math.min(1, Math.max(0, v))));
+
+  // ── TI per band (mean of 2 modulation-frequency TIs)
+  const tiPerBand = mkCorrected.map(row => {
+    const avgSnr = row.reduce((s, m) => s + mtfToSnrDb(m), 0) / row.length;
+    return (avgSnr + 15) / 30;
+  });
+
+  // ── STI
+  const { alpha, beta } = WEIGHTS.male;
+  let sti = 0;
+  for (let k = 0; k < 7; k++) {
+    sti += alpha[k] * tiPerBand[k];
+    if (k < 6) sti -= beta[k] * Math.sqrt(tiPerBand[k] * tiPerBand[k + 1]);
+  }
+  sti = Math.max(0, Math.min(1, sti));
+
+  return { sti, mk, tiPerBand, rating: stiRatingLabel(sti), warning };
+}
+
+// ─── SIGNAL PROCESSING HELPERS ───────────────────────────────────────────────
+
+/**
+ * Apply a 1-octave bandpass filter centred at fc using cascaded biquad sections.
+ * 3 cascaded 2nd-order Butterworth bandpass sections ≈ 6th-order response.
+ * Not as sharp as the 18th-order filter in zawi01/stipa but adequate for STIPA.
+ */
+async function filterOctave(samples, fs, fc) {
+  const offlineCtx = new OfflineAudioContext(1, samples.length, fs);
+  const buf        = offlineCtx.createBuffer(1, samples.length, fs);
+  buf.copyToChannel(samples instanceof Float32Array ? samples : Float32Array.from(samples), 0);
+
+  const src = offlineCtx.createBufferSource();
+  src.buffer = buf;
+
+  // 3 cascaded biquad bandpass sections
+  const Q = Math.SQRT2; // Q = √2 → 1-octave −3 dB bandwidth
+  let node = src;
+  for (let i = 0; i < 3; i++) {
+    const bp = offlineCtx.createBiquadFilter();
+    bp.type            = 'bandpass';
+    bp.frequency.value = fc;
+    bp.Q.value         = Q;
+    node.connect(bp);
+    node = bp;
+  }
+  node.connect(offlineCtx.destination);
+  src.start(0);
+  const rendered = await offlineCtx.startRendering();
+  return rendered.getChannelData(0);
+}
+
+/**
+ * Square-law envelope detection followed by 100 Hz low-pass filter.
+ * Mirrors stipa.m: envelope = lowpass(x.*x, 100, fs)
+ */
+function envelopeDetect(bandSamples, fs) {
+  // Single-pole IIR low-pass at 100 Hz
+  const wc  = 2 * Math.PI * 100 / fs;
+  const a1  = Math.exp(-wc);
+  const b0  = 1 - a1;
+  const out = new Float32Array(bandSamples.length);
+  let prev  = 0;
+  for (let i = 0; i < bandSamples.length; i++) {
+    const sq = bandSamples[i] * bandSamples[i];
+    prev     = b0 * sq + a1 * prev;
+    out[i]   = prev;
+  }
+  return out;
+}
+
+/**
+ * Single-bin DFT at modulation frequency Fm — port of zawi01/stipa MTF computation.
+ * Uses only complete periods of Fm within the signal to avoid spectral leakage.
+ *
+ * mk(Fm) = 2 × |Σ I(t)·exp(j·2π·Fm·t)| / Σ I(t)
+ */
+function singleBinDFT(envelope, fs, Fm, totalSeconds) {
+  const wholePeriodSeconds = Math.floor(Fm * totalSeconds) / Fm;
+  const N = Math.round(wholePeriodSeconds * fs);
+  if (N < 1) return 0;
+
+  let re = 0, im = 0, dc = 0;
+  for (let i = 0; i < N; i++) {
+    const t   = (i / fs);
+    const val = envelope[i] ?? 0;
+    re  += val * Math.cos(2 * Math.PI * Fm * t);
+    im  += val * Math.sin(2 * Math.PI * Fm * t);
+    dc  += val;
+  }
+  if (dc < 1e-12) return 0;
+  return 2 * Math.sqrt(re * re + im * im) / dc;
+}
+
+/** Adjust mk for ambient noise: mk_corrected = mk × Isk / (Isk + Ink) */
+function adjustAmbientNoise(mk, Lsk, Lnk) {
+  return mk.map((row, k) => {
+    const Isk = Math.pow(10, Lsk[k] / 10);
+    const Ink = Math.pow(10, Lnk[k] / 10);
+    const ratio = Isk / (Isk + Ink);
+    return row.map(v => v * ratio);
+  });
+}
+
+/** Adjust mk for auditory masking and reception threshold (IEC 60268-16:2020). */
+function adjustAuditoryMasking(mk, Lsk, Lnk) {
+  const Ak  = AUDITORY_MASKING_THRESHOLDS;
+  const Isk = Lsk.map(l => Math.pow(10, l / 10));
+  const Ink = Lnk.map(l => Math.pow(10, l / 10));
+  const Ik  = Isk.map((s, k) => s + Ink[k]);
+
+  // Level-dependent upward masking from adjacent lower band
+  const La = Lsk.slice(0, 6).map(L => {
+    if (L < 63)            return 0.5 * L - 65;
+    if (L < 67)            return 1.8 * L - 146.9;
+    if (L < 100)           return 0.5 * L - 59.8;
+    return -10;
+  });
+  const a    = La.map(la => Math.pow(10, la / 10));
+  const Iamk = [0, ...Isk.slice(0, 6).map((is, k) => is * a[k])];
+  const Irtk = Ak.map(ak => Math.pow(10, ak / 10));
+
+  return mk.map((row, k) => {
+    const ratio = Ik[k] / (Ik[k] + Iamk[k] + Irtk[k]);
+    return row.map(v => v * ratio);
+  });
+}
+
 // ─── REW RT60 RESPONSE PARSER ─────────────────────────────────────────────────
 
 /**