fix: stabilize localization comparison semantics

packages/acoustics-core/das.js

@@ -601,14 +601,23 @@ export function estimateDasWorkload(options) {
   };
 }
 
-export function finalizeDasBands(rawMaps, workload, experiment, rawScale = "power") {
+export function finalizeDasBands(
+  rawMaps,
+  workload,
+  experiment,
+  rawScale = "power",
+  options = {},
+) {
   if (rawScale !== "power" && rawScale !== "score") {
     throw new RangeError('rawScale must be "power" or "score".');
   }
+  const averageBandBins = options.averageBandBins !== false;
   return experiment.frequencyBands.map((band, bandIndex) => {
     const rawMap = rawMaps[bandIndex];
-    for (let pointIndex = 0; pointIndex < workload.pointCount; pointIndex += 1) {
-      rawMap[pointIndex] /= workload.bandBins[bandIndex].length;
+    if (averageBandBins) {
+      for (let pointIndex = 0; pointIndex < workload.pointCount; pointIndex += 1) {
+        rawMap[pointIndex] /= workload.bandBins[bandIndex].length;
+      }
     }
     const rawPeak = findPeak(rawMap);
     const displayMap = Float32Array.from(rawMap);

packages/acoustics-core/gcc-phat.js

@@ -14,7 +14,10 @@ import {
 
 export const GCC_PHAT_ANALYSIS_PLAN_KIND = "micarrsim.gcc-phat-analysis-plan/1";
 export const DEFAULT_GCC_PHAT_FLOOR_RATIO = 1e-6;
+export const DEFAULT_GCC_PHAT_POINT_BATCH_SIZE = 4096;
+export const DEFAULT_GCC_PHAT_CROSS_SPECTRUM_BATCH_SIZE = 16384;
 export const DEFAULT_MAXIMUM_GCC_PHAT_PLAN_BYTES = 512 * 1024 * 1024;
+export const DEFAULT_MAXIMUM_GCC_PHAT_WORKING_BYTES = 512 * 1024 * 1024;
 
 function createAbortError(signal) {
   if (signal?.reason instanceof Error && signal.reason.name === "AbortError") {
@@ -241,6 +244,60 @@ function requirePhatFloorRatio(value) {
   return ratio;
 }
 
+function requireMaximumWorkingBytes(value) {
+  const maximumWorkingBytes = value ?? DEFAULT_MAXIMUM_GCC_PHAT_WORKING_BYTES;
+  if (!Number.isSafeInteger(maximumWorkingBytes) || maximumWorkingBytes <= 0) {
+    throw new RangeError("maximumWorkingBytes must be a positive safe integer.");
+  }
+  return maximumWorkingBytes;
+}
+
+/**
+ * Reject an oversized no-plan analysis before createDistanceMatrix allocates
+ * its P x M table. This mirrors estimateWithGeometry's retained/spectrum/
+ * calibration/output accounting without requiring geometry to exist first.
+ */
+function preflightUnplannedWorkingSet(options, configuration) {
+  const microphoneCount = configuration.arrayDefinition.microphones.length;
+  const sampleCount = validateChannels(options.channels, microphoneCount);
+  const { bandBins, selectedBins } = createBandGeometry(configuration);
+  const frameCount = getFrameCount(
+    sampleCount,
+    configuration.fftSize,
+    configuration.hopSize,
+  );
+  const bandCount = bandBins.length;
+  const distanceBytes =
+    configuration.pointCount * microphoneCount * Float32Array.BYTES_PER_ELEMENT;
+  const stftWorkspaceBytes =
+    configuration.fftSize * Float64Array.BYTES_PER_ELEMENT * 3;
+  const correlationWorkspaceBytes =
+    configuration.fftSize * Float64Array.BYTES_PER_ELEMENT * 2;
+  const spectrumBytes =
+    selectedBins.length * frameCount * microphoneCount *
+    Float32Array.BYTES_PER_ELEMENT * 2;
+  const calibrationPhaseBytes =
+    selectedBins.length * microphoneCount * Float64Array.BYTES_PER_ELEMENT;
+  const outputBytes =
+    bandCount * configuration.pointCount * Float32Array.BYTES_PER_ELEMENT * 2;
+  const estimatedWorkingBytes =
+    distanceBytes + stftWorkspaceBytes + correlationWorkspaceBytes +
+    spectrumBytes + calibrationPhaseBytes + outputBytes;
+  if (!Number.isSafeInteger(estimatedWorkingBytes)) {
+    throw new RangeError("GCC-PHAT estimated working set exceeds the safe integer range.");
+  }
+  const maximumWorkingBytes = requireMaximumWorkingBytes(options.maximumWorkingBytes);
+  if (estimatedWorkingBytes > maximumWorkingBytes) {
+    throw new RangeError(
+      "GCC-PHAT estimated working set " +
+      (estimatedWorkingBytes / 1024 ** 2).toFixed(1) +
+      " MiB exceeds maximumWorkingBytes " +
+      (maximumWorkingBytes / 1024 ** 2).toFixed(1) + " MiB before geometry allocation.",
+    );
+  }
+  return { estimatedWorkingBytes, maximumWorkingBytes };
+}
+
 function estimateWithGeometry(options, configuration, geometry) {
   const sampleCount = validateChannels(options.channels, geometry.microphoneCount);
   const frameCount = getFrameCount(
@@ -332,6 +389,7 @@ export function estimateGccPhatWorkload(options) {
     assertCompatibleAnalysisPlan(options.analysisPlan, options);
     geometry = options.analysisPlan;
   } else {
+    preflightUnplannedWorkingSet(options, configuration);
     geometry = createGeometry(configuration);
   }
   return estimateWithGeometry(options, configuration, geometry);
@@ -467,7 +525,7 @@ function circularLinearSample(values, lagSamples) {
   return values[firstIndex] * (1 - fraction) + values[secondIndex] * fraction;
 }
 
-function fillBandPairCorrelation({
+async function fillBandPairCorrelation({
   correlationReal,
   correlationImaginary,
   bins,
@@ -476,11 +534,14 @@ function fillBandPairCorrelation({
   referenceIndex,
   phasesByBin,
   phatFloorRatio,
+  signal,
+  crossSpectrumBatchSize,
 }) {
   correlationReal.fill(0);
   correlationImaginary.fill(0);
   const microphoneCount = stft.microphoneCount;
   let maximumMagnitude = 0;
+  let operationsUntilYield = crossSpectrumBatchSize;
 
   for (const bin of bins) {
     const spectrum = stft.spectra.get(bin);
@@ -498,6 +559,12 @@ function fillBandPairCorrelation({
         maximumMagnitude,
         Math.hypot(crossReal, crossImaginary),
       );
+      operationsUntilYield -= 1;
+      if (operationsUntilYield === 0) {
+        await cooperativeYield();
+        throwIfAborted(signal);
+        operationsUntilYield = crossSpectrumBatchSize;
+      }
     }
   }
 
@@ -535,6 +602,12 @@ function fillBandPairCorrelation({
         normalizedReal * correctionReal - normalizedImaginary * correctionImaginary;
       accumulatedImaginary +=
         normalizedReal * correctionImaginary + normalizedImaginary * correctionReal;
+      operationsUntilYield -= 1;
+      if (operationsUntilYield === 0) {
+        await cooperativeYield();
+        throwIfAborted(signal);
+        operationsUntilYield = crossSpectrumBatchSize;
+      }
     }
 
     const averagedReal = accumulatedReal / stft.frameCount;
@@ -552,13 +625,38 @@ function fillBandPairCorrelation({
   fftInPlace(correlationReal, correlationImaginary, true);
 }
 
-function shiftAndSquareScores(rawMaps) {
-  for (const rawMap of rawMaps) {
-    let minimum = Infinity;
-    for (const value of rawMap) minimum = Math.min(minimum, value);
+async function normalizePositiveScores(
+  rawMaps,
+  bandBins,
+  fftSize,
+  pairCount,
+  pointBatchSize,
+  signal,
+) {
+  for (let bandIndex = 0; bandIndex < rawMaps.length; bandIndex += 1) {
+    const rawMap = rawMaps[bandIndex];
+    const coefficientCount = bandBins[bandIndex].reduce(
+      (count, bin) => count + (bin === 0 || bin === fftSize / 2 ? 1 : 2),
+      0,
+    );
+    const normalization = fftSize / (pairCount * coefficientCount);
     for (let pointIndex = 0; pointIndex < rawMap.length; pointIndex += 1) {
-      const shifted = rawMap[pointIndex] - minimum;
-      rawMap[pointIndex] = shifted * shifted;
+      // PHAT coefficients have magnitude <= 1. After Hermitian completion and
+      // inverse-FFT scaling, coefficientCount / fftSize is therefore a strict
+      // per-pair bound. Pair averaging makes the positive score stable across
+      // microphone count, band width, and scan extent without implying a
+      // probability. Negative correlation is evidence against the hypothesis.
+      rawMap[pointIndex] = Math.min(
+        1,
+        Math.max(0, rawMap[pointIndex] * normalization),
+      );
+      if (
+        (pointIndex + 1) % pointBatchSize === 0 &&
+        pointIndex + 1 < rawMap.length
+      ) {
+        await cooperativeYield();
+        throwIfAborted(signal);
+      }
     }
   }
 }
@@ -568,16 +666,16 @@ async function computeGccPhatMapInternal(options, configuration, geometry, plan)
     channels,
     signal,
     onProgress,
-    maximumWorkingBytes = 512 * 1024 * 1024,
+    maximumWorkingBytes = DEFAULT_MAXIMUM_GCC_PHAT_WORKING_BYTES,
     maximumFftButterflies = 1_000_000_000,
-    pointBatchSize = 256,
+    pointBatchSize = DEFAULT_GCC_PHAT_POINT_BATCH_SIZE,
+    crossSpectrumBatchSize = DEFAULT_GCC_PHAT_CROSS_SPECTRUM_BATCH_SIZE,
     stftTransformBatchSize = 128,
   } = options;
   requirePositiveInteger(pointBatchSize, "pointBatchSize");
+  requirePositiveInteger(crossSpectrumBatchSize, "crossSpectrumBatchSize");
   requirePositiveInteger(stftTransformBatchSize, "stftTransformBatchSize");
-  if (!Number.isSafeInteger(maximumWorkingBytes) || maximumWorkingBytes <= 0) {
-    throw new RangeError("maximumWorkingBytes must be a positive safe integer.");
-  }
+  requireMaximumWorkingBytes(maximumWorkingBytes);
   const phatFloorRatio = requirePhatFloorRatio(
     options.phatFloorRatio ??
     options.phatRelativeFloor ??
@@ -666,7 +764,7 @@ async function computeGccPhatMapInternal(options, configuration, geometry, plan)
     const rawMap = rawMaps[bandIndex];
     for (const microphoneIndex of geometry.pairMicrophoneIndexes) {
       throwIfAborted(signal);
-      fillBandPairCorrelation({
+      await fillBandPairCorrelation({
         correlationReal,
         correlationImaginary,
         bins,
@@ -675,6 +773,8 @@ async function computeGccPhatMapInternal(options, configuration, geometry, plan)
         referenceIndex: geometry.reference.index,
         phasesByBin,
         phatFloorRatio,
+        signal,
+        crossSpectrumBatchSize,
       });
       if (plan) {
         plan.statistics.inverseFftTransforms += 1;
@@ -715,12 +815,20 @@ async function computeGccPhatMapInternal(options, configuration, geometry, plan)
   }
 
   throwIfAborted(signal);
-  shiftAndSquareScores(rawMaps);
+  await normalizePositiveScores(
+    rawMaps,
+    geometry.bandBins,
+    configuration.fftSize,
+    geometry.pairCount,
+    pointBatchSize,
+    signal,
+  );
   const bands = finalizeDasBands(
     rawMaps,
     workload,
     configuration.experiment,
     "score",
+    { averageBandBins: false },
   );
   return {
     backend: "cpu",
@@ -753,6 +861,7 @@ export async function computeGccPhatMap(options) {
   const configuration = validateConfiguration(options);
   const plan = options.analysisPlan;
   if (!plan) {
+    preflightUnplannedWorkingSet(options, configuration);
     const geometry = createGeometry(configuration);
     return computeGccPhatMapInternal(options, configuration, geometry, undefined);
   }

packages/acoustics-core/localization.js

@@ -25,9 +25,9 @@ function requestedAlgorithm(options) {
  */
 export function estimateLocalizationWorkload(options) {
   const algorithm = requestedAlgorithm(options);
-  return isDasAlgorithm(algorithm.id)
-    ? estimateDasWorkload(options)
-    : estimateGccPhatWorkload(options);
+  if (isDasAlgorithm(algorithm.id)) return estimateDasWorkload(options);
+  if (algorithm.id === "gcc-phat") return estimateGccPhatWorkload(options);
+  throw new RangeError("No workload estimator is registered for " + algorithm.id + ".");
 }
 
 /**
@@ -35,9 +35,9 @@ export function estimateLocalizationWorkload(options) {
  */
 export function createLocalizationAnalysisPlan(options) {
   const algorithm = requestedAlgorithm(options);
-  return isDasAlgorithm(algorithm.id)
-    ? createDasAnalysisPlan(options)
-    : createGccPhatAnalysisPlan(options);
+  if (isDasAlgorithm(algorithm.id)) return createDasAnalysisPlan(options);
+  if (algorithm.id === "gcc-phat") return createGccPhatAnalysisPlan(options);
+  throw new RangeError("No analysis-plan factory is registered for " + algorithm.id + ".");
 }
 
 /**
@@ -47,7 +47,7 @@ export function createLocalizationAnalysisPlan(options) {
  */
 export async function computeLocalizationMap(options) {
   const algorithm = requestedAlgorithm(options);
-  return isDasAlgorithm(algorithm.id)
-    ? computeDasMap(options)
-    : computeGccPhatMap(options);
+  if (isDasAlgorithm(algorithm.id)) return computeDasMap(options);
+  if (algorithm.id === "gcc-phat") return computeGccPhatMap(options);
+  throw new RangeError("No CPU localization implementation is registered for " + algorithm.id + ".");
 }

packages/acoustics-core/realtime-localizer.js

@@ -3,6 +3,10 @@ import {
   computeLocalizationMap,
   createLocalizationAnalysisPlan,
 } from "./localization.js";
+import {
+  DEFAULT_GCC_PHAT_CROSS_SPECTRUM_BATCH_SIZE,
+  DEFAULT_GCC_PHAT_POINT_BATCH_SIZE,
+} from "./gcc-phat.js";
 import {
   isDasAlgorithm,
   requireLocalizationAlgorithmDefinition,
@@ -355,8 +359,13 @@ export function createLocalizationPlan(config) {
       ? { dasPlan: analysisPlan }
       : { gccPhatPlan: analysisPlan }),
     computeOptions: Object.freeze({
-      pointBatchSize: DEFAULT_REALTIME_POINT_BATCH_SIZE,
+      pointBatchSize: isDasAlgorithm(algorithm.id)
+        ? DEFAULT_REALTIME_POINT_BATCH_SIZE
+        : DEFAULT_GCC_PHAT_POINT_BATCH_SIZE,
       stftTransformBatchSize: DEFAULT_REALTIME_STFT_TRANSFORM_BATCH_SIZE,
+      ...(!isDasAlgorithm(algorithm.id)
+        ? { crossSpectrumBatchSize: DEFAULT_GCC_PHAT_CROSS_SPECTRUM_BATCH_SIZE }
+        : {}),
       ...(phatRelativeFloor !== undefined ? { phatRelativeFloor } : {}),
       ...(phatFloorRatio !== undefined ? { phatFloorRatio } : {}),
       ...(config.computeOptions ?? {}),