diff --git a/packages/acoustics-core/beam-pattern.js b/packages/acoustics-core/beam-pattern.js index 6bba96c..3fe8547 100644 --- a/packages/acoustics-core/beam-pattern.js +++ b/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) { - const sampleCount = responseLinear.length; +function sampledLocalMaxima(responsePower, anglesDeg, offsetsDeg, microphoneCount) { + 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; - if (maximum - minimum <= comparisonTolerance) return []; + const comparisonTolerance = + 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 current = responseLinear[index]; - const next = responseLinear[(index + 1) % sampleCount]; - if ( - current >= previous - comparisonTolerance && - 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, - }); + const previousIndex = (index - 1 + sampleCount) % sampleCount; + const difference = Math.abs(responsePower[index] - responsePower[previousIndex]); + if (difference > largestBoundaryDifference) { + largestBoundaryDifference = difference; + largestBoundaryIndex = index; } } - 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: { - status: metricResolutionSatisfied ? "sampled" : "resolution-unresolved", - peaks, - }, + sampledLocalMaxima: peaks, }, }; } diff --git a/tests/unit/beam-pattern.test.js b/tests/unit/beam-pattern.test.js index 31117ca..6a2fccd 100644 --- a/tests/unit/beam-pattern.test.js +++ b/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"); });