Browse Source

STIPA modulation

Benjamin Harris 1 month ago
parent
commit
1405d22685
4 changed files with 156 additions and 59 deletions
  1. 9 1
      .claude/settings.local.json
  2. 1 0
      .gitignore
  3. 39 10
      CLAUDE.md
  4. 107 48
      src/sti.js

+ 9 - 1
.claude/settings.local.json

@@ -3,7 +3,15 @@
     "allow": [
       "PowerShell(Get-Command node, npm, npx -ErrorAction SilentlyContinue)",
       "PowerShell(docker --version 2>&1)",
-      "WebFetch(domain:www.roomeqwizard.com)"
+      "WebFetch(domain:www.roomeqwizard.com)",
+      "WebFetch(domain:github.com)",
+      "WebFetch(domain:raw.githubusercontent.com)",
+      "WebSearch",
+      "WebFetch(domain:www.nti-audio.com)",
+      "WebFetch(domain:www.utko.fekt.vut.cz)",
+      "WebFetch(domain:cdn.standards.iteh.ai)",
+      "WebFetch(domain:acousplan.com)",
+      "WebFetch(domain:en.wikipedia.org)"
     ]
   }
 }

+ 1 - 0
.gitignore

@@ -4,3 +4,4 @@ coverage
 .vite
 documents
 CLAUDE_CHAT.md
+.claude

+ 39 - 10
CLAUDE.md

@@ -195,7 +195,23 @@ Mono audio source
 
 ### STIPA measurement
 
-**STIPA** (Speech Transmission Index for PA systems) is the simplified single-measurement STI method defined in IEC 60268-16. It uses a pre-recorded test signal with specific modulation at frequencies across seven octave bands (125 Hz–8 kHz), played through the live PA system and measured with a handheld analyser.
+**STIPA** (Speech Transmission Index for PA systems) is the simplified single-measurement STI method defined in IEC 60268-16. It uses a pre-recorded test signal with specific amplitude modulation at 2 frequencies per octave band (14 combinations total across 7 bands), played through the live PA system and measured with a handheld analyser.
+
+**STIPA test signal — modulation frequencies per octave band (IEC 60268-16:2020, confirmed: zawi01/stipa):**
+
+| Band | Mod freq 1 | Mod freq 2 | Band level |
+| ---- | ---------- | ---------- | ---------- |
+| 125 Hz | 1.00 Hz | 5.00 Hz | −2.5 dB |
+| 250 Hz | 0.63 Hz | 3.15 Hz | +0.5 dB |
+| 500 Hz | 2.00 Hz | 10.00 Hz | 0 dB |
+| 1 kHz | 1.25 Hz | 8.00 Hz | −6 dB |
+| 2 kHz | 0.80 Hz | 4.00 Hz | −12 dB |
+| 4 kHz | 2.50 Hz | 12.50 Hz | −18 dB |
+| 8 kHz | 6.30 Hz | 1.60 Hz | −24 dB |
+
+Signal: pink noise carrier, RMS = 0.1, duration ≈ 18 s (15–25 s valid). Reference implementation: [zawi01/stipa](https://github.com/zawi01/stipa) — validated at MAD = 0.0033 vs NTI XL2 across 12 measurement points.
+
+**Full STI vs STIPA:** Full STI uses all 14 modulation frequencies per band (98 combinations); STIPA uses only 2 per band (14 combinations). The app computes STI from REW RT60 using the full 14-frequency method. Handheld measurement with the STIPA test signal is faster (~18 s) and accurate to ±0.02–0.03.
 
 **IEC 60268-16 rating scale:**
 
@@ -376,17 +392,30 @@ REW can calculate STI internally (via its STI panel, IEC 60268-16:2020 indirect
 3. Click **Compute STI** — the app calls `GET /measurements/:id/rt60?octaveFrac=1` and runs the IEC 60268-16 Schroeder/Houtgast–Steeneken calculation
 4. The result shows overall STI, per-band Transmission Index (TI), and per-band T60
 
-**Calculation method** (implemented in `src/sti.js`):
+**Calculation method** (implemented in `src/sti.js`, reference: [zawi01/stipa](https://github.com/zawi01/stipa)):
+
+For each of the 7 octave bands and all 14 modulation frequencies:
+
+```text
+MTF(Fm, T60, SNR) = [1 / √(1 + (2πFm·T60/13.8)²)] × [1 / (1 + 10^(−SNR/10))]
+SNR_apparent      = clamp(10·log10(MTF / (1−MTF)), −15, +15) dB
+TI(band)          = mean[(SNR_apparent + 15) / 30] across 14 mod frequencies
+STI               = Σ(α_k·TI_k) − Σ(β_k·√(TI_k·TI_{k+1}))
+```
+
+**Speech weighting coefficients (IEC 60268-16:2020 Table B.1):**
+
+| Band | 125 | 250 | 500 | 1k | 2k | 4k | 8k |
+| ---- | --- | --- | --- | -- | -- | -- | -- |
+| α male | 0.085 | 0.127 | 0.230 | 0.233 | 0.309 | 0.224 | 0.173 |
+| β male | 0.085 | 0.078 | 0.065 | 0.011 | 0.047 | 0.095 | — |
+| α female | 0.000 | 0.117 | 0.223 | 0.216 | 0.328 | 0.250 | 0.194 |
+| β female | 0.000 | 0.099 | 0.066 | 0.062 | 0.025 | 0.076 | — |
 
-- For each of the 7 STI octave bands (125 Hz–8 kHz) and 14 modulation frequencies (0.63–12.5 Hz):
-  - MTF from T60: `1 / √(1 + (2π × Fm × T60 / 13.8)²)` (Schroeder formula)
-  - Optional noise correction: `MTF_eff = MTF × 10^(SNR/10) / (1 + 10^(SNR/10))`
-  - Apparent SNR: `10 × log10(MTF_eff / (1 − MTF_eff))`, clamped to [−15, +15] dB
-  - Band TI: `(SNR_apparent + 15) / 30`, averaged over 14 modulation frequencies
-- STI = `Σ(α_k × TI_k) − Σ(β_k × √(TI_k × TI_{k+1}))` using male speech weights (IEC Table B.1)
-- No background noise correction is applied unless per-band SNR is provided
+**Auditory masking reception thresholds** (not yet applied — future enhancement):
+125 Hz: 46 dB | 250 Hz: 27 dB | 500 Hz: 12 dB | 1 kHz: 6.5 dB | 2 kHz: 7.5 dB | 4 kHz: 8 dB | 8 kHz: 12 dB
 
-**Important:** The RT60 method gives a theoretically exact STI only for a diffuse reverberant field. In this hall with anti-phasing, each listener hears primarily one speaker, so the effective RT60 driving that listener's intelligibility is different from the room RT60. Treat the REW-derived STI as an upper bound; compare with measured STIPA at representative seat positions.
+**Important:** The RT60 method assumes a diffuse reverberant field. In this hall with anti-phasing, each listener hears primarily one speaker segment, so the effective acoustic environment differs from a simple diffuse field. Treat the REW-derived STI as an upper-bound estimate; compare with measured STIPA at representative seat positions to calibrate.
 
 ---
 

+ 107 - 48
src/sti.js

@@ -1,66 +1,127 @@
-// ─── STI CALCULATION — IEC 60268-16:2020 INDIRECT METHOD ────────────────────
-// Computes Speech Transmission Index from per-octave RT60 and SNR values
-// (Schroeder/Houtgast–Steeneken approximation).
+// ─── STI / STIPA — IEC 60268-16:2020 ────────────────────────────────────────
 //
-// Input:  rt60PerBand[7]  — T60 in seconds for 125, 250, 500, 1k, 2k, 4k, 8k Hz
-//         snrPerBand[7]   — signal-to-noise ratio in dB per octave band
-//                           (use 30 dB where noise is unknown — effectively no degradation)
-// Output: { sti, tiPerBand, rating }
+// Two distinct computations live here:
+//
+//   computeSTIfromRT60()  — full STI using all 14 modulation frequencies per
+//                           octave band (Schroeder/Houtgast–Steeneken indirect
+//                           method). Used when REW RT60 data is available.
+//                           Validated: zawi01/stipa MAD = 0.0033 vs NTI XL2.
+//
+//   STIPA_MOD_FREQS       — the 2-per-band subset used by the STIPA test signal
+//                           (IEC 60268-16:2020). Exported for reference and for
+//                           signal-based analysis if raw audio is ever added.
+//
+// Reference: https://github.com/zawi01/stipa
+//            IEC 60268-16:2020
 
-// ─── CONSTANTS ────────────────────────────────────────────────────────────────
+// ─── BAND DEFINITIONS ────────────────────────────────────────────────────────
 
 export const STI_BANDS_HZ = [125, 250, 500, 1000, 2000, 4000, 8000];
 
-// 14 modulation frequencies (Hz) defined in IEC 60268-16:2020
-const MOD_FREQS = [0.63, 0.80, 1.00, 1.25, 1.60, 2.00, 2.50, 3.15, 4.00, 5.00, 6.30, 8.00, 10.00, 12.50];
-
-// Male speech weighting (IEC 60268-16:2020 Table B.1, auditory masking excluded)
-// α weights per band; β inter-band redundancy (0 for last band)
-const ALPHA = [0.085, 0.127, 0.230, 0.233, 0.309, 0.224, 0.173];
-const BETA  = [0.085, 0.078, 0.065, 0.011, 0.047, 0.095, 0.000];
-// Σα − Σβ = 1.381 − 0.381 = 1.000 (identity check)
-
-// ─── CORE CALCULATION ─────────────────────────────────────────────────────────
+// Full STI: 14 modulation frequencies per octave band (IEC 60268-16:2020)
+const MOD_FREQS_FULL = [
+  0.63, 0.80, 1.00, 1.25, 1.60, 2.00,
+  2.50, 3.15, 4.00, 5.00, 6.30, 8.00, 10.00, 12.50,
+];
+
+// 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
+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
+];
+
+// STIPA test-signal band levels relative to 0 dBFS (Revision 5, IEC 60268-16)
+// [125, 250, 500, 1k, 2k, 4k, 8k Hz]
+export const STIPA_BAND_LEVELS_DB = [-2.5, 0.5, 0, -6, -12, -18, -24];
+
+// ─── SPEECH WEIGHTING COEFFICIENTS (IEC 60268-16:2020 Table B.1) ─────────────
+// α: per-band importance weight; β: inter-band redundancy correction
+// Identity: Σα − Σβ = 1.000 for both speaker types
+
+const WEIGHTS = {
+  male: {
+    alpha: [0.085, 0.127, 0.230, 0.233, 0.309, 0.224, 0.173],
+    beta:  [0.085, 0.078, 0.065, 0.011, 0.047, 0.095, 0.000],
+  },
+  female: {
+    alpha: [0.000, 0.117, 0.223, 0.216, 0.328, 0.250, 0.194],
+    beta:  [0.000, 0.099, 0.066, 0.062, 0.025, 0.076, 0.000],
+  },
+};
+
+// Auditory masking reception thresholds per band (dB) — IEC 60268-16:2020 Table B.4
+// [125, 250, 500, 1k, 2k, 4k, 8k Hz]
+export const AUDITORY_MASKING_THRESHOLDS = [46, 27, 12, 6.5, 7.5, 8, 12];
+
+// ─── CORE MTF ────────────────────────────────────────────────────────────────
 
 /**
- * Compute STI from RT60 per octave band using the Schroeder formula.
+ * Modulation Transfer Function for a single modulation frequency and band.
+ * Combines exponential reverberation decay (Schroeder) with additive noise.
  *
- * @param {number[]} rt60PerBand  T60 in seconds for each of the 7 STI octave bands
- * @param {number[]} snrPerBand   SNR in dB per band (default 30 dB if not available)
- * @returns {{ sti: number, tiPerBand: number[], rating: string }}
+ *   m(Fm, T60, snr) = [1 / √(1 + (2πFm·T60/13.8)²)] × [1 / (1 + 10^(−snr/10))]
+ *
+ * Constant 13.8 ≈ −60 dB / (ln(10) × 2) — the RT60 decay constant.
  */
-export function computeSTIfromRT60(rt60PerBand, snrPerBand = Array(7).fill(30)) {
-  const tiPerBand = STI_BANDS_HZ.map((_, k) => {
-    const T  = Math.max(0.01, rt60PerBand[k]);
-    const snr = snrPerBand[k];
+function mtf(Fm, T60, snrDb) {
+  const reverbTerm = 1 / Math.sqrt(1 + Math.pow((2 * Math.PI * Fm * T60) / 13.8, 2));
+  const noiseTerm  = 1 / (1 + Math.pow(10, -snrDb / 10));
+  return reverbTerm * noiseTerm;
+}
 
-    const modTIs = MOD_FREQS.map(Fm => {
-      // Schroeder: MTF from exponential decay with T60
-      const mtfT60 = 1 / Math.sqrt(1 + Math.pow((2 * Math.PI * Fm * T) / 13.8, 2));
+/** Convert MTF value to apparent SNR in dB, clamped to [−15, +15] per IEC. */
+function mtfToSnrDb(m) {
+  const clamped = Math.max(1e-9, Math.min(1 - 1e-9, m));
+  return Math.max(-15, Math.min(15, 10 * Math.log10(clamped / (1 - clamped))));
+}
 
-      // Noise correction (additive noise reduces MTF)
-      const linearSNR = Math.pow(10, snr / 10);
-      const mtfEff    = mtfT60 * linearSNR / (1 + linearSNR);
+/** Transmission Index for a single octave band from an array of MTF values. */
+function tiFromMtfArray(mtfValues) {
+  const avgSnr = mtfValues.reduce((s, m) => s + mtfToSnrDb(m), 0) / mtfValues.length;
+  return (avgSnr + 15) / 30;
+}
 
-      // Apparent SNR, clamped to [−15, +15] dB
-      const ratio    = Math.max(1e-6, Math.min(1 - 1e-6, mtfEff));
-      const snrApp   = Math.max(-15, Math.min(15, 10 * Math.log10(ratio / (1 - ratio))));
+// ─── MAIN EXPORT ─────────────────────────────────────────────────────────────
 
-      return (snrApp + 15) / 30;
-    });
+/**
+ * Compute STI from per-octave-band RT60 values (IEC 60268-16:2020 indirect method).
+ * Uses all 14 modulation frequencies — this is full STI, not the STIPA subset.
+ *
+ * @param {number[]} rt60PerBand  T60 in seconds — [125, 250, 500, 1k, 2k, 4k, 8k Hz]
+ * @param {number[]} snrPerBand   SNR in dB per band. Default 30 dB (no noise degradation).
+ * @param {'male'|'female'} speaker  Speech weighting. Default 'male'.
+ * @returns {{ sti: number, tiPerBand: number[], rating: string }}
+ */
+export function computeSTIfromRT60(
+  rt60PerBand,
+  snrPerBand  = Array(7).fill(30),
+  speaker     = 'male',
+) {
+  const { alpha, beta } = WEIGHTS[speaker] ?? WEIGHTS.male;
 
-    return modTIs.reduce((s, v) => s + v, 0) / modTIs.length;
+  const tiPerBand = STI_BANDS_HZ.map((_, k) => {
+    const T   = Math.max(0.01, rt60PerBand[k]);
+    const snr = snrPerBand[k];
+    const mtfValues = MOD_FREQS_FULL.map(Fm => mtf(Fm, T, snr));
+    return tiFromMtfArray(mtfValues);
   });
 
-  // Weighted combination across bands with inter-band redundancy correction
+  // Weighted combination with inter-band redundancy correction
   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 += 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, tiPerBand, rating: stiRatingLabel(sti) };
+  const result = Math.max(0, Math.min(1, sti));
+  return { sti: result, tiPerBand, rating: stiRatingLabel(result) };
 }
 
 // ─── RATING HELPERS (IEC 60268-16:2020) ──────────────────────────────────────
@@ -84,12 +145,11 @@ export function stiRatingColor(v) {
 // ─── REW RT60 RESPONSE PARSER ─────────────────────────────────────────────────
 
 /**
- * Parse REW's RT60 API response into per-octave-band T60 values.
- * REW returns Base64-encoded big-endian float32 pairs [frequency, T60_seconds].
- * We pick the 7 STI octave band centre frequencies.
+ * Parse REW's RT60 API response into per-octave T60 values for STI_BANDS_HZ.
+ * REW encodes data as Base64 big-endian float32 pairs [frequency_Hz, T60_seconds].
  *
  * @param {object} rewRT60Response  Raw JSON from GET /measurements/:id/rt60
- * @returns {number[]}  T60 in seconds for STI_BANDS_HZ, or null if parse fails
+ * @returns {number[]|null}  T60 array [7] in seconds, or null if unparseable
  */
 export function parseRewRT60(rewRT60Response) {
   try {
@@ -101,7 +161,6 @@ export function parseRewRT60(rewRT60Response) {
     for (let i = 0; i < bin.length; i++) buf[i] = bin.charCodeAt(i);
     const view = new DataView(buf.buffer);
 
-    // Pairs of [freq_Hz, T60_s] as big-endian float32
     const points = [];
     for (let i = 0; i + 7 < buf.length; i += 8) {
       points.push({ freq: view.getFloat32(i, false), t60: view.getFloat32(i + 4, false) });