fix: harden beam pattern numerical diagnostics

packages/acoustics-core/beam-pattern.js

@@ -6,6 +6,8 @@ export const DEFAULT_RECEIVE_BEAM_PATTERN_DB_FLOOR = -60;
 export const MAXIMUM_RECEIVE_BEAM_PATTERN_EVALUATIONS = 16_000_000;
 export const MAXIMUM_RECEIVE_BEAM_PATTERN_SAMPLES = 14_400;
 export const MAXIMUM_RECEIVE_BEAM_PATTERN_MICROPHONES = 1_024;
+export const MAXIMUM_RECEIVE_BEAM_PATTERN_FREQUENCY_HZ = 1_000_000;
+export const MAXIMUM_RECEIVE_BEAM_PATTERN_BASELINE_COMPARISONS = 600_000;
 export const RECEIVE_BEAM_PATTERN_METRIC_PHASE_STEP_LIMIT_RAD = Math.PI / 8;
 
 const HALF_POWER = 0.5;
@@ -107,6 +109,7 @@ function halfPowerBeamwidth(
   stepDeg,
   steeringAngleDeg,
   evaluatePowerAtOffset,
+  analyticallyConstant,
 ) {
   const sampleCount = responsePower.length;
   let allAboveHalfPower = true;
@@ -115,11 +118,14 @@ function halfPowerBeamwidth(
   }
   if (allAboveHalfPower) {
     return {
-      status: "full-circle-above-half-power",
+      status: analyticallyConstant
+        ? "full-circle-above-half-power"
+        : "no-sampled-half-power-crossing",
       widthDeg: null,
       startAngleDeg: null,
       endAngleDeg: null,
-      fullCircleAboveHalfPower: true,
+      fullCircleAboveHalfPower: analyticallyConstant ? true : null,
+      bracketSource: null,
     };
   }
 
@@ -137,11 +143,14 @@ function halfPowerBeamwidth(
   }
   if (leftSteps + rightSteps + 1 >= sampleCount) {
     return {
-      status: "full-circle-above-half-power",
+      status: analyticallyConstant
+        ? "full-circle-above-half-power"
+        : "no-sampled-half-power-crossing",
       widthDeg: null,
       startAngleDeg: null,
       endAngleDeg: null,
-      fullCircleAboveHalfPower: true,
+      fullCircleAboveHalfPower: analyticallyConstant ? true : null,
+      bracketSource: null,
     };
   }
 
@@ -156,49 +165,112 @@ function halfPowerBeamwidth(
     evaluatePowerAtOffset,
   );
   return {
-    status: "bounded",
+    status: "sampled-bracketed-bisection",
     widthDeg: rightBoundaryOffset - leftBoundaryOffset,
     startAngleDeg: normalizeBeamPatternAngle(steeringAngleDeg + leftBoundaryOffset),
     endAngleDeg: normalizeBeamPatternAngle(steeringAngleDeg + rightBoundaryOffset),
     fullCircleAboveHalfPower: false,
+    bracketSource: "adjacent-angular-samples",
+    bisectionIterations: HPBW_BISECTION_ITERATIONS,
   };
 }
 
-function sampledLocalMaxima(responseLinear, anglesDeg, offsetsDeg) {
+function sampledLocalMaxima(responsePower, anglesDeg, offsetsDeg, microphoneCount) {
-  const sampleCount = responseLinear.length;
+  const sampleCount = responsePower.length;
   let minimum = Infinity;
   let maximum = -Infinity;
-  for (const value of responseLinear) {
+  for (const value of responsePower) {
     minimum = Math.min(minimum, value);
     maximum = Math.max(maximum, value);
   }
-  const comparisonTolerance = Math.max(1, maximum) * 64 * Number.EPSILON;
+  const comparisonTolerance =
-  if (maximum - minimum <= comparisonTolerance) return [];
+    Math.max(1, maximum) * 256 * Number.EPSILON * microphoneCount;
+  if (maximum - minimum <= comparisonTolerance) {
+    return {
+      status: "flat-on-grid",
+      comparisonTolerancePower: comparisonTolerance,
+      peaks: [],
+      globalMaximumRegions: [],
+      selected: null,
+    };
+  }
 
-  const peaks = [];
+  let largestBoundaryIndex = 0;
+  let largestBoundaryDifference = -Infinity;
   for (let index = 0; index < sampleCount; index += 1) {
-    const previous = responseLinear[(index - 1 + sampleCount) % sampleCount];
+    const previousIndex = (index - 1 + sampleCount) % sampleCount;
-    const current = responseLinear[index];
+    const difference = Math.abs(responsePower[index] - responsePower[previousIndex]);
-    const next = responseLinear[(index + 1) % sampleCount];
+    if (difference > largestBoundaryDifference) {
-    if (
+      largestBoundaryDifference = difference;
-      current >= previous - comparisonTolerance &&
+      largestBoundaryIndex = index;
-      current >= next - comparisonTolerance &&
-      (current > previous + comparisonTolerance || current > next + comparisonTolerance)
-    ) {
-      peaks.push({
-        sampleIndex: index,
-        angleDeg: anglesDeg[index],
-        offsetDeg: offsetsDeg[index],
-        responseLinear: current,
-        responseDb: current > 0 ? 20 * Math.log10(current) : -Infinity,
-      });
     }
   }
-  peaks.sort((left, right) =>
+
-    right.responseLinear - left.responseLinear ||
+  const runs = [];
+  let run = null;
+  for (let step = 0; step < sampleCount; step += 1) {
+    const index = (largestBoundaryIndex + step) % sampleCount;
+    const power = responsePower[index];
+    if (!run || Math.abs(power - run.anchorPower) > comparisonTolerance) {
+      run = {
+        anchorPower: power,
+        indices: [],
+        maximumPower: power,
+        representativeIndex: index,
+      };
+      runs.push(run);
+    }
+    run.indices.push(index);
+    if (power > run.maximumPower) {
+      run.maximumPower = power;
+      run.representativeIndex = index;
+    }
+  }
+
+  const peakRegions = [];
+  for (let runIndex = 0; runIndex < runs.length; runIndex += 1) {
+    const current = runs[runIndex];
+    const previous = runs[(runIndex - 1 + runs.length) % runs.length];
+    const next = runs[(runIndex + 1) % runs.length];
+    if (
+      current.maximumPower <= previous.maximumPower + comparisonTolerance ||
+      current.maximumPower <= next.maximumPower + comparisonTolerance
+    ) {
+      continue;
+    }
+    const representativeIndex = current.representativeIndex;
+    const firstIndex = current.indices[0];
+    const lastIndex = current.indices.at(-1);
+    peakRegions.push({
+      sampleIndices: [...current.indices],
+      startAngleDeg: anglesDeg[firstIndex],
+      endAngleDeg: anglesDeg[lastIndex],
+      representativeIndex,
+      angleDeg: anglesDeg[representativeIndex],
+      offsetDeg: offsetsDeg[representativeIndex],
+      responsePower: current.maximumPower,
+      responseLinear: Math.sqrt(current.maximumPower),
+      responseDb: current.maximumPower > 0
+        ? 10 * Math.log10(current.maximumPower)
+        : -Infinity,
+    });
+  }
+  peakRegions.sort((left, right) =>
+    right.responsePower - left.responsePower ||
     Math.abs(left.offsetDeg) - Math.abs(right.offsetDeg) ||
     left.angleDeg - right.angleDeg);
-  return peaks;
+  const globalMaximumRegions = peakRegions.filter(
+    (peak) => maximum - peak.responsePower <= comparisonTolerance,
+  );
+  return {
+    status: globalMaximumRegions.length === 1
+      ? "unique-global-maximum-on-grid"
+      : "tied-global-maxima-on-grid",
+    comparisonTolerancePower: comparisonTolerance,
+    peaks: peakRegions,
+    globalMaximumRegions,
+    selected: globalMaximumRegions.length === 1 ? globalMaximumRegions[0] : null,
+  };
 }
 
 function maximumProjectedBaseline(coordinates) {
@@ -240,6 +312,12 @@ export function computeReceiveBeamPattern(options) {
   const array = options?.array;
   const positions = requireMicrophones(array);
   const frequencyHz = requirePositive(options?.frequencyHz, "frequencyHz");
+  if (frequencyHz > MAXIMUM_RECEIVE_BEAM_PATTERN_FREQUENCY_HZ) {
+    throw new RangeError(
+      "frequencyHz must not exceed " +
+      MAXIMUM_RECEIVE_BEAM_PATTERN_FREQUENCY_HZ + " Hz.",
+    );
+  }
   const speedOfSound = requirePositive(options?.speedOfSound ?? 343, "speedOfSound");
   const cut = requireCut(options?.cut ?? "xz");
   const steeringAngleDeg = normalizeBeamPatternAngle(options?.steeringAngleDeg ?? 0);
@@ -268,7 +346,42 @@ export function computeReceiveBeamPattern(options) {
   if (!Number.isSafeInteger(maximumEvaluations) || maximumEvaluations < 1) {
     throw new RangeError("maximumEvaluations must be a positive safe integer.");
   }
+  const maximumBaselineComparisons = options?.maximumBaselineComparisons ??
+    MAXIMUM_RECEIVE_BEAM_PATTERN_BASELINE_COMPARISONS;
+  if (
+    !Number.isSafeInteger(maximumBaselineComparisons) ||
+    maximumBaselineComparisons < 0
+  ) {
+    throw new RangeError(
+      "maximumBaselineComparisons must be a non-negative safe integer.",
+    );
+  }
   const microphoneCount = positions.length;
+  const projectedBaselineComparisons =
+    microphoneCount * (microphoneCount - 1) / 2;
+  if (projectedBaselineComparisons > maximumBaselineComparisons) {
+    throw new RangeError(
+      "receive beam pattern requires " + projectedBaselineComparisons +
+      " projected-baseline comparisons; limit is " +
+      maximumBaselineComparisons + ".",
+    );
+  }
+  const preliminarySampleCount =
+    requestedSampleCount === "auto" ? 36 : requestedSampleCount;
+  const preliminaryMaximumEvaluationCount = microphoneCount * (
+    preliminarySampleCount + HPBW_BISECTION_ITERATIONS * 2
+  );
+  if (
+    !Number.isSafeInteger(preliminaryMaximumEvaluationCount) ||
+    preliminaryMaximumEvaluationCount > maximumEvaluations
+  ) {
+    throw new RangeError(
+      "receive beam pattern requires at most " +
+      preliminaryMaximumEvaluationCount +
+      " microphone-angle evaluations before geometry audit; limit is " +
+      maximumEvaluations + ".",
+    );
+  }
   const centroid = [0, 0, 0];
   for (const position of positions) {
     centroid[0] += position[0] / microphoneCount;
@@ -283,7 +396,12 @@ export function computeReceiveBeamPattern(options) {
     ], cut));
   const projectedApertureM = maximumProjectedBaseline(centeredCutCoordinates);
   const waveNumber = TAU * frequencyHz / speedOfSound;
-  if (!Number.isFinite(waveNumber) || !Number.isFinite(waveNumber * projectedApertureM)) {
+  const maximumPhaseMagnitudeRadBound =
+    2 * waveNumber * projectedApertureM;
+  if (
+    !Number.isFinite(waveNumber) ||
+    !Number.isFinite(maximumPhaseMagnitudeRadBound)
+  ) {
     throw new RangeError("frequency, speed of sound, and array coordinates overflow phase.");
   }
   const recommendedMetricSampleCount =
@@ -373,6 +491,9 @@ export function computeReceiveBeamPattern(options) {
     const real = summedReal[index] / microphoneCount;
     const imaginary = summedImaginary[index] / microphoneCount;
     const power = real * real + imaginary * imaginary;
+    if (!Number.isFinite(real) || !Number.isFinite(imaginary) || !Number.isFinite(power)) {
+      throw new RangeError("receive beam-pattern phasor accumulation became non-finite.");
+    }
     const magnitude = Math.sqrt(power);
     const decibels = power > 0 ? 10 * Math.log10(power) : -Infinity;
     responseReal[index] = real;
@@ -403,8 +524,12 @@ export function computeReceiveBeamPattern(options) {
       real += Math.cos(phase);
       imaginary += Math.sin(phase);
     }
-    return (real * real + imaginary * imaginary) /
+    const power = (real * real + imaginary * imaginary) /
       (microphoneCount * microphoneCount);
+    if (!Number.isFinite(power)) {
+      throw new RangeError("receive beam-pattern metric evaluation became non-finite.");
+    }
+    return power;
   };
   const halfPower = metricResolutionSatisfied
     ? halfPowerBeamwidth(
@@ -413,6 +538,7 @@ export function computeReceiveBeamPattern(options) {
       stepDeg,
       steeringAngleDeg,
       evaluatePowerAtOffset,
+      projectedApertureM === 0,
     )
     : {
       status: "resolution-unresolved",
@@ -420,10 +546,17 @@ export function computeReceiveBeamPattern(options) {
       startAngleDeg: null,
       endAngleDeg: null,
       fullCircleAboveHalfPower: null,
+      bracketSource: null,
     };
   const peaks = metricResolutionSatisfied
-    ? sampledLocalMaxima(responseLinear, anglesDeg, offsetsDeg)
+    ? sampledLocalMaxima(responsePower, anglesDeg, offsetsDeg, microphoneCount)
-    : [];
+    : {
+      status: "resolution-unresolved",
+      comparisonTolerancePower: null,
+      peaks: [],
+      globalMaximumRegions: [],
+      selected: null,
+    };
   const steeringLinear = responseLinear[steeringIndex];
   return {
     model: "far-field-narrowband-conventional-delay-and-sum-receive-array-factor",
@@ -448,7 +581,8 @@ export function computeReceiveBeamPattern(options) {
       metricDirectionEvaluationCount * microphoneCount,
     workload: {
       projectedBaselineComparisons:
-        microphoneCount * (microphoneCount - 1) / 2,
+        projectedBaselineComparisons,
+      maximumBaselineComparisons,
       sampleEvaluationCount,
       maximumMetricEvaluationCount,
       maximumEvaluationCount: evaluationCount,
@@ -457,6 +591,7 @@ export function computeReceiveBeamPattern(options) {
       mode: requestedSampleCount === "auto" ? "auto" : "fixed",
       criterion: "maximum-adjacent-pair-phase-step",
       projectedApertureM,
+      maximumPhaseMagnitudeRadBound,
       maximumAdjacentPhaseStepRadBound,
       metricPhaseStepLimitRad:
         RECEIVE_BEAM_PATTERN_METRIC_PHASE_STEP_LIMIT_RAD,
@@ -485,10 +620,7 @@ export function computeReceiveBeamPattern(options) {
       steeringLinear,
       steeringDb: steeringLinear > 0 ? 20 * Math.log10(steeringLinear) : -Infinity,
       halfPowerBeamwidth: halfPower,
-      sampledLocalMaxima: {
+      sampledLocalMaxima: peaks,
-        status: metricResolutionSatisfied ? "sampled" : "resolution-unresolved",
-        peaks,
-      },
     },
   };
 }

tests/unit/beam-pattern.test.js

@@ -95,13 +95,28 @@ test("half-wavelength two-microphone broadside cut has nulls, 60-degree HPBW, an
   assertNear(sampleAtAngle(pattern, 270), 0, 1e-12, "endfire null at -x");
   assertNear(sampleAtAngle(pattern, 0), 1, 1e-12, "steered broadside response");
   assertNear(sampleAtAngle(pattern, 180), 1, 1e-12, "planar front-back response");
-  assertEqual(pattern.metrics.halfPowerBeamwidth.status, "bounded", "HPBW status");
+  assertEqual(
+    pattern.metrics.halfPowerBeamwidth.status,
+    "sampled-bracketed-bisection",
+    "HPBW status",
+  );
   assertNear(pattern.metrics.halfPowerBeamwidth.widthDeg, 60, 1e-10, "power HPBW");
   const zeroDbAngles = pattern.metrics.sampledLocalMaxima.peaks
     .filter((peak) => Math.abs(peak.responseDb) < 1e-10)
     .map((peak) => peak.angleDeg)
     .sort((left, right) => left - right);
   assertEqual(zeroDbAngles.join(","), "0,180", "front-back maxima remain tied");
+  assertEqual(
+    pattern.metrics.sampledLocalMaxima.status,
+    "tied-global-maxima-on-grid",
+    "global tie is explicit",
+  );
+  assertEqual(pattern.metrics.sampledLocalMaxima.selected, null, "tie has no selected peak");
+  assertEqual(
+    pattern.metrics.sampledLocalMaxima.globalMaximumRegions.length,
+    2,
+    "both tied regions are preserved",
+  );
   assertNear(pattern.responseReal[Math.floor(pattern.sampleCount / 2)], 1, 1e-12, "raw real");
   assertNear(pattern.responseImaginary[Math.floor(pattern.sampleCount / 2)], 0, 1e-12, "raw imaginary");
   assertNear(pattern.responsePower[Math.floor(pattern.sampleCount / 2)], 1, 1e-12, "raw power");
@@ -258,6 +273,33 @@ test("coarse high-frequency cuts preserve raw samples but withhold unresolved me
   assertEqual(pattern.responsePower.length, 360, "raw sampled power remains available");
 });
 
+test("non-constant cuts without sampled half-power crossings do not claim 360-degree HPBW", () => {
+  const pattern = computeReceiveBeamPattern({
+    array: arrayFromPositions([[-0.0005, 0, 0], [0.0005, 0, 0]]),
+    frequencyHz: 1_000,
+    speedOfSound: 343,
+    cut: "xz",
+    steeringAngleDeg: 0,
+    sampleCount: 720,
+  });
+  assert(pattern.angularSampling.projectedApertureM > 0, "array is not constant on cut");
+  assert(
+    Math.min(...pattern.responsePower) > 0.5,
+    "sampled cut stays above half power",
+  );
+  assertEqual(
+    pattern.metrics.halfPowerBeamwidth.status,
+    "no-sampled-half-power-crossing",
+    "missing crossing remains unresolved",
+  );
+  assertEqual(pattern.metrics.halfPowerBeamwidth.widthDeg, null, "no fabricated width");
+  assertEqual(
+    pattern.metrics.halfPowerBeamwidth.fullCircleAboveHalfPower,
+    null,
+    "full-circle claim is withheld",
+  );
+});
+
 test("auto angular sampling increases resolution and reports when the safety cap is insufficient", () => {
   const moderate = computeReceiveBeamPattern({
     array: arrayFromPositions([[-0.225, 0, 0], [0.225, 0, 0]]),
@@ -324,4 +366,26 @@ test("beam-pattern preflight rejects invalid input and excessive work", () => {
     microphoneLimitError.message.includes("1024"),
     "microphone cap is explicit",
   );
+
+  let baselineError = null;
+  try {
+    computeReceiveBeamPattern({
+      array,
+      frequencyHz: 1_000,
+      maximumBaselineComparisons: 0,
+    });
+  } catch (error) {
+    baselineError = error;
+  }
+  assert(baselineError instanceof RangeError, "baseline pair cap is checked");
+  assert(baselineError.message.includes("1 projected-baseline"), "pair work is explicit");
+
+  let frequencyLimitError = null;
+  try {
+    computeReceiveBeamPattern({ array, frequencyHz: 1_000_001 });
+  } catch (error) {
+    frequencyLimitError = error;
+  }
+  assert(frequencyLimitError instanceof RangeError, "extreme analog frequency is bounded");
+  assert(frequencyLimitError.message.includes("1000000"), "frequency cap is explicit");
 });