Files

1033 lines
37 KiB
JavaScript

import { TAU } from "./math.js";
export const RECEIVE_BEAM_PATTERN_CUTS = Object.freeze(["xy", "xz", "yz"]);
export const DEFAULT_RECEIVE_BEAM_PATTERN_SAMPLES = 720;
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;
export const DEFAULT_RECEIVE_BEAM_PATTERN_SPHERE_AZIMUTH_SAMPLES = 360;
export const DEFAULT_RECEIVE_BEAM_PATTERN_SPHERE_ELEVATION_SAMPLES = 181;
export const MAXIMUM_RECEIVE_BEAM_PATTERN_SPHERE_AZIMUTH_SAMPLES = 1_440;
export const MAXIMUM_RECEIVE_BEAM_PATTERN_SPHERE_ELEVATION_SAMPLES = 721;
export const MAXIMUM_RECEIVE_BEAM_PATTERN_SPHERE_SAMPLES = 262_144;
const HALF_POWER = 0.5;
const RADIANS_PER_DEGREE = Math.PI / 180;
const HPBW_BISECTION_ITERATIONS = 48;
function requireFinite(value, label) {
const number = Number(value);
if (!Number.isFinite(number)) throw new TypeError(label + " must be finite.");
return number;
}
function requirePositive(value, label) {
const number = requireFinite(value, label);
if (!(number > 0)) throw new RangeError(label + " must be greater than zero.");
return number;
}
function requireCut(cut) {
if (!RECEIVE_BEAM_PATTERN_CUTS.includes(cut)) {
throw new RangeError(
"cut must be one of " + RECEIVE_BEAM_PATTERN_CUTS.join(", ") + ".",
);
}
return cut;
}
function requireMicrophones(array) {
const microphones = array?.microphones;
if (!Array.isArray(microphones) || microphones.length === 0) {
throw new TypeError("array.microphones must contain at least one microphone.");
}
if (microphones.length > MAXIMUM_RECEIVE_BEAM_PATTERN_MICROPHONES) {
throw new RangeError(
"array.microphones exceeds the receive beam-pattern limit of " +
MAXIMUM_RECEIVE_BEAM_PATTERN_MICROPHONES + ".",
);
}
return microphones.map((microphone, microphoneIndex) => {
const position = microphone?.position;
if (!Array.isArray(position) || position.length !== 3) {
throw new TypeError(
"array.microphones[" + microphoneIndex + "].position must contain x, y, z.",
);
}
return position.map((coordinate, axis) =>
requireFinite(
coordinate,
"array.microphones[" + microphoneIndex + "].position[" + axis + "]",
));
});
}
export function normalizeBeamPatternAngle(angleDeg) {
const angle = requireFinite(angleDeg, "angleDeg");
const normalized = ((angle % 360) + 360) % 360;
return Object.is(normalized, -0) ? 0 : normalized;
}
/**
* Convert a polar angle on an orthogonal great-circle cut to an array-to-source
* unit vector. Zero degrees points along the first named axis and +90 degrees
* points along the second: XY=(+x,+y), XZ=(+z,+x), YZ=(+z,+y).
*/
export function receiveBeamPatternDirection(cut, angleDeg) {
const plane = requireCut(cut);
const angleRadians = requireFinite(angleDeg, "angleDeg") * RADIANS_PER_DEGREE;
const cosine = Math.cos(angleRadians);
const sine = Math.sin(angleRadians);
if (plane === "xy") return [cosine, sine, 0];
if (plane === "xz") return [sine, 0, cosine];
return [0, sine, cosine];
}
/**
* Convert standard spherical azimuth/elevation to an array-to-source unit
* vector. Azimuth is measured in the XY plane from +X toward +Y. Elevation is
* measured from the XY plane toward +Z, with +90 degrees at the north pole.
*/
export function receiveBeamPatternSphereDirection(azimuthDeg, elevationDeg) {
const azimuthRadians = normalizeBeamPatternAngle(azimuthDeg) *
RADIANS_PER_DEGREE;
const elevation = requireFinite(elevationDeg, "elevationDeg");
if (elevation < -90 || elevation > 90) {
throw new RangeError("elevationDeg must lie inside -90..90 degrees.");
}
if (elevation === 90) return [0, 0, 1];
if (elevation === -90) return [0, 0, -1];
const elevationRadians = elevation * RADIANS_PER_DEGREE;
const cosineElevation = Math.cos(elevationRadians);
return [
cosineElevation * Math.cos(azimuthRadians),
cosineElevation * Math.sin(azimuthRadians),
Math.sin(elevationRadians),
];
}
function cutCoordinates(position, cut) {
if (cut === "xy") return [position[0], position[1]];
if (cut === "xz") return [position[2], position[0]];
return [position[2], position[1]];
}
function bisectHalfPowerBoundary(
insideOffsetDeg,
outsideOffsetDeg,
evaluatePowerAtOffset,
) {
let inside = insideOffsetDeg;
let outside = outsideOffsetDeg;
for (let iteration = 0; iteration < HPBW_BISECTION_ITERATIONS; iteration += 1) {
const middle = (inside + outside) / 2;
if (evaluatePowerAtOffset(middle) >= HALF_POWER) inside = middle;
else outside = middle;
}
return (inside + outside) / 2;
}
function halfPowerBeamwidth(
responsePower,
steeringIndex,
stepDeg,
steeringAngleDeg,
evaluatePowerAtOffset,
analyticallyConstant,
) {
const sampleCount = responsePower.length;
let allAboveHalfPower = true;
for (let index = 0; index < sampleCount; index += 1) {
if (responsePower[index] < HALF_POWER) allAboveHalfPower = false;
}
if (allAboveHalfPower) {
return {
status: analyticallyConstant
? "full-circle-above-half-power"
: "no-sampled-half-power-crossing",
widthDeg: null,
startAngleDeg: null,
endAngleDeg: null,
fullCircleAboveHalfPower: analyticallyConstant ? true : null,
bracketSource: null,
};
}
let leftSteps = 0;
while (leftSteps < sampleCount - 1) {
const nextIndex = (steeringIndex - leftSteps - 1 + sampleCount) % sampleCount;
if (responsePower[nextIndex] < HALF_POWER) break;
leftSteps += 1;
}
let rightSteps = 0;
while (rightSteps < sampleCount - leftSteps - 1) {
const nextIndex = (steeringIndex + rightSteps + 1) % sampleCount;
if (responsePower[nextIndex] < HALF_POWER) break;
rightSteps += 1;
}
if (leftSteps + rightSteps + 1 >= sampleCount) {
return {
status: analyticallyConstant
? "full-circle-above-half-power"
: "no-sampled-half-power-crossing",
widthDeg: null,
startAngleDeg: null,
endAngleDeg: null,
fullCircleAboveHalfPower: analyticallyConstant ? true : null,
bracketSource: null,
};
}
const leftBoundaryOffset = bisectHalfPowerBoundary(
-leftSteps * stepDeg,
-(leftSteps + 1) * stepDeg,
evaluatePowerAtOffset,
);
const rightBoundaryOffset = bisectHalfPowerBoundary(
rightSteps * stepDeg,
(rightSteps + 1) * stepDeg,
evaluatePowerAtOffset,
);
return {
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(responsePower, anglesDeg, offsetsDeg, microphoneCount) {
const sampleCount = responsePower.length;
let minimum = Infinity;
let maximum = -Infinity;
for (const value of responsePower) {
minimum = Math.min(minimum, value);
maximum = Math.max(maximum, value);
}
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,
};
}
let largestBoundaryIndex = 0;
let largestBoundaryDifference = -Infinity;
for (let index = 0; index < sampleCount; index += 1) {
const previousIndex = (index - 1 + sampleCount) % sampleCount;
const difference = Math.abs(responsePower[index] - responsePower[previousIndex]);
if (difference > largestBoundaryDifference) {
largestBoundaryDifference = difference;
largestBoundaryIndex = index;
}
}
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);
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) {
let maximum = 0;
for (let left = 0; left < coordinates.length; left += 1) {
for (let right = left + 1; right < coordinates.length; right += 1) {
maximum = Math.max(
maximum,
Math.hypot(
coordinates[left][0] - coordinates[right][0],
coordinates[left][1] - coordinates[right][1],
),
);
}
}
return maximum;
}
function maximumSpatialBaseline(positions) {
let maximum = 0;
for (let left = 0; left < positions.length; left += 1) {
for (let right = left + 1; right < positions.length; right += 1) {
maximum = Math.max(
maximum,
Math.hypot(
positions[left][0] - positions[right][0],
positions[left][1] - positions[right][1],
positions[left][2] - positions[right][2],
),
);
}
}
return maximum;
}
function requireSphereSampleCount(value, label, minimum, maximum) {
if (!Number.isInteger(value) || value < minimum || value > maximum) {
throw new RangeError(
label + " must be an integer inside " + minimum + ".." + maximum + ".",
);
}
return value;
}
function recommendedSampleCount(waveNumber, projectedApertureM) {
if (projectedApertureM === 0) return 36;
const ratio = RECEIVE_BEAM_PATTERN_METRIC_PHASE_STEP_LIMIT_RAD /
(2 * waveNumber * projectedApertureM);
if (ratio >= 1) return 36;
const required = Math.ceil(Math.PI / Math.asin(ratio));
if (!Number.isSafeInteger(required)) return Number.MAX_SAFE_INTEGER;
return Math.ceil(required / 4) * 4;
}
/**
* Compute the ideal, calibrated, omnidirectional receive array factor for a
* narrowband far-field plane wave and conventional delay-and-sum steering.
*
* This function intentionally has no sample-rate input. frequencyHz is the
* analog acoustic hypothesis used for k=2πf/c, including frequencies above a
* later ADC Nyquist limit. The result is a great-circle slice, not a full 3-D
* directivity surface and not a probability distribution.
*/
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);
const requestedSampleCount =
options?.sampleCount ?? DEFAULT_RECEIVE_BEAM_PATTERN_SAMPLES;
if (
requestedSampleCount !== "auto" &&
(!Number.isInteger(requestedSampleCount) ||
requestedSampleCount < 36 ||
requestedSampleCount > MAXIMUM_RECEIVE_BEAM_PATTERN_SAMPLES)
) {
throw new RangeError(
"sampleCount must be 'auto' or an integer inside 36.." +
MAXIMUM_RECEIVE_BEAM_PATTERN_SAMPLES + ".",
);
}
const dbFloor = requireFinite(
options?.dbFloor ?? DEFAULT_RECEIVE_BEAM_PATTERN_DB_FLOOR,
"dbFloor",
);
if (dbFloor < -300 || dbFloor > -1) {
throw new RangeError("dbFloor must lie inside -300..-1 dB.");
}
const maximumEvaluations = options?.maximumEvaluations ??
MAXIMUM_RECEIVE_BEAM_PATTERN_EVALUATIONS;
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;
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;
centroid[1] += position[1] / microphoneCount;
centroid[2] += position[2] / microphoneCount;
}
const centeredCutCoordinates = positions.map((position) =>
cutCoordinates([
position[0] - centroid[0],
position[1] - centroid[1],
position[2] - centroid[2],
], cut));
const projectedApertureM = maximumProjectedBaseline(centeredCutCoordinates);
const waveNumber = TAU * frequencyHz / speedOfSound;
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 =
recommendedSampleCount(waveNumber, projectedApertureM);
const maximumMetricDirectionEvaluations = projectedApertureM > 0
? HPBW_BISECTION_ITERATIONS * 2
: 0;
const maximumSamplesByWork = Math.floor(
maximumEvaluations / microphoneCount - maximumMetricDirectionEvaluations,
);
const allowedSampleCount = Math.min(
MAXIMUM_RECEIVE_BEAM_PATTERN_SAMPLES,
maximumSamplesByWork,
);
if (allowedSampleCount < 36) {
throw new RangeError(
"maximumEvaluations cannot cover the minimum 36-angle receive pattern.",
);
}
const sampleCount = requestedSampleCount === "auto"
? Math.min(
allowedSampleCount,
Math.max(DEFAULT_RECEIVE_BEAM_PATTERN_SAMPLES, recommendedMetricSampleCount),
)
: requestedSampleCount;
const maximumMetricEvaluationCount =
microphoneCount * maximumMetricDirectionEvaluations;
const sampleEvaluationCount = microphoneCount * sampleCount;
const evaluationCount = sampleEvaluationCount + maximumMetricEvaluationCount;
if (!Number.isSafeInteger(evaluationCount) || evaluationCount > maximumEvaluations) {
throw new RangeError(
"receive beam pattern requires at most " + evaluationCount +
" microphone-angle evaluations; limit is " + maximumEvaluations + ".",
);
}
const stepDeg = 360 / sampleCount;
const stepRadians = TAU / sampleCount;
const maximumAdjacentPhaseStepRadBound =
2 * waveNumber * projectedApertureM * Math.sin(stepRadians / 2);
const metricResolutionSatisfied =
maximumAdjacentPhaseStepRadBound <=
RECEIVE_BEAM_PATTERN_METRIC_PHASE_STEP_LIMIT_RAD *
(1 + 64 * Number.EPSILON);
const steeringIndex = Math.floor(sampleCount / 2);
const anglesDeg = new Float64Array(sampleCount);
const offsetsDeg = new Float64Array(sampleCount);
const directionCosines = new Float64Array(sampleCount);
const directionSines = new Float64Array(sampleCount);
for (let index = 0; index < sampleCount; index += 1) {
const offsetDeg = (index - steeringIndex) * stepDeg;
const angleDeg = normalizeBeamPatternAngle(steeringAngleDeg + offsetDeg);
const radians = angleDeg * RADIANS_PER_DEGREE;
offsetsDeg[index] = offsetDeg;
anglesDeg[index] = angleDeg;
directionCosines[index] = Math.cos(radians);
directionSines[index] = Math.sin(radians);
}
const steeringRadians = steeringAngleDeg * RADIANS_PER_DEGREE;
const steeringCosine = Math.cos(steeringRadians);
const steeringSine = Math.sin(steeringRadians);
const summedReal = new Float64Array(sampleCount);
const summedImaginary = new Float64Array(sampleCount);
for (const [firstCoordinate, secondCoordinate] of centeredCutCoordinates) {
const steeringProjection =
firstCoordinate * steeringCosine + secondCoordinate * steeringSine;
for (let index = 0; index < sampleCount; index += 1) {
const arrivalProjection =
firstCoordinate * directionCosines[index] +
secondCoordinate * directionSines[index];
const phase = waveNumber * (arrivalProjection - steeringProjection);
summedReal[index] += Math.cos(phase);
summedImaginary[index] += Math.sin(phase);
}
}
const responseReal = new Float64Array(sampleCount);
const responseImaginary = new Float64Array(sampleCount);
const responsePower = new Float64Array(sampleCount);
const responseLinear = new Float64Array(sampleCount);
const responseDb = new Float64Array(sampleCount);
const displayDb = new Float64Array(sampleCount);
let maximumLinear = 0;
let maximumPower = 0;
for (let index = 0; index < sampleCount; index += 1) {
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;
responseImaginary[index] = imaginary;
responsePower[index] = power;
responseLinear[index] = magnitude;
responseDb[index] = decibels;
displayDb[index] = Math.max(dbFloor, decibels);
maximumLinear = Math.max(maximumLinear, magnitude);
maximumPower = Math.max(maximumPower, power);
}
let metricDirectionEvaluationCount = 0;
const evaluatePowerAtOffset = (offsetDeg) => {
metricDirectionEvaluationCount += 1;
const angleRadians =
(steeringAngleDeg + offsetDeg) * RADIANS_PER_DEGREE;
const directionCosine = Math.cos(angleRadians);
const directionSine = Math.sin(angleRadians);
let real = 0;
let imaginary = 0;
for (const [firstCoordinate, secondCoordinate] of centeredCutCoordinates) {
const steeringProjection =
firstCoordinate * steeringCosine + secondCoordinate * steeringSine;
const arrivalProjection =
firstCoordinate * directionCosine + secondCoordinate * directionSine;
const phase = waveNumber * (arrivalProjection - steeringProjection);
real += Math.cos(phase);
imaginary += Math.sin(phase);
}
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(
responsePower,
steeringIndex,
stepDeg,
steeringAngleDeg,
evaluatePowerAtOffset,
projectedApertureM === 0,
)
: {
status: "resolution-unresolved",
widthDeg: null,
startAngleDeg: null,
endAngleDeg: null,
fullCircleAboveHalfPower: null,
bracketSource: null,
};
const peaks = metricResolutionSatisfied
? 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",
assumptions: Object.freeze({
wavefront: "far-field-plane-wave",
bandwidth: "narrowband",
elements: "ideal-calibrated-omnidirectional-point-microphones",
environment: "free-field",
coordinates: "meters",
}),
cut,
frequencyHz,
speedOfSound,
wavelengthM: speedOfSound / frequencyHz,
waveNumberRadPerM: waveNumber,
steeringAngleDeg,
steeringDirection: receiveBeamPatternDirection(cut, steeringAngleDeg),
sampleCount,
angularResolutionDeg: stepDeg,
microphoneCount,
evaluationCount: sampleEvaluationCount +
metricDirectionEvaluationCount * microphoneCount,
workload: {
projectedBaselineComparisons:
projectedBaselineComparisons,
maximumBaselineComparisons,
sampleEvaluationCount,
maximumMetricEvaluationCount,
maximumEvaluationCount: evaluationCount,
},
angularSampling: {
mode: requestedSampleCount === "auto" ? "auto" : "fixed",
criterion: "maximum-adjacent-pair-phase-step",
projectedApertureM,
maximumPhaseMagnitudeRadBound,
maximumAdjacentPhaseStepRadBound,
metricPhaseStepLimitRad:
RECEIVE_BEAM_PATTERN_METRIC_PHASE_STEP_LIMIT_RAD,
metricResolutionSatisfied,
recommendedSampleCount: recommendedMetricSampleCount,
autoCapped:
requestedSampleCount === "auto" &&
sampleCount < Math.max(
DEFAULT_RECEIVE_BEAM_PATTERN_SAMPLES,
recommendedMetricSampleCount,
),
},
dbFloor,
anglesDeg,
offsetsDeg,
responseReal,
responseImaginary,
responsePower,
responseLinear,
responseDb,
displayDb,
metrics: {
maximumLinear,
maximumPower,
maximumDb: maximumPower > 0 ? 10 * Math.log10(maximumPower) : -Infinity,
steeringLinear,
steeringDb: steeringLinear > 0 ? 20 * Math.log10(steeringLinear) : -Infinity,
halfPowerBeamwidth: halfPower,
sampledLocalMaxima: peaks,
},
};
}
/**
* Sample the same ideal conventional receive array factor over a full sphere.
*
* The returned rectangular grid is equirectangular: azimuth samples span
* [0, 360) without a duplicate endpoint, while elevation samples include both
* poles. Pole rows therefore contain coincident directions for every azimuth.
* All response arrays are raw elevation-major samples; no interpolation,
* smoothing, frequency folding, or sampled-maximum normalization is applied.
*/
export function computeReceiveBeamPatternSphere(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 steeringAzimuthDeg = normalizeBeamPatternAngle(
options?.steeringAzimuthDeg ?? 0,
);
const steeringElevationDeg = requireFinite(
options?.steeringElevationDeg ?? 90,
"steeringElevationDeg",
);
if (steeringElevationDeg < -90 || steeringElevationDeg > 90) {
throw new RangeError(
"steeringElevationDeg must lie inside -90..90 degrees.",
);
}
const azimuthSampleCount = requireSphereSampleCount(
options?.azimuthSampleCount ??
DEFAULT_RECEIVE_BEAM_PATTERN_SPHERE_AZIMUTH_SAMPLES,
"azimuthSampleCount",
4,
MAXIMUM_RECEIVE_BEAM_PATTERN_SPHERE_AZIMUTH_SAMPLES,
);
const elevationSampleCount = requireSphereSampleCount(
options?.elevationSampleCount ??
DEFAULT_RECEIVE_BEAM_PATTERN_SPHERE_ELEVATION_SAMPLES,
"elevationSampleCount",
3,
MAXIMUM_RECEIVE_BEAM_PATTERN_SPHERE_ELEVATION_SAMPLES,
);
const sampleCount = azimuthSampleCount * elevationSampleCount;
if (
!Number.isSafeInteger(sampleCount) ||
sampleCount > MAXIMUM_RECEIVE_BEAM_PATTERN_SPHERE_SAMPLES
) {
throw new RangeError(
"sphere sample product exceeds the receive beam-pattern limit of " +
MAXIMUM_RECEIVE_BEAM_PATTERN_SPHERE_SAMPLES + ".",
);
}
const dbFloor = requireFinite(
options?.dbFloor ?? DEFAULT_RECEIVE_BEAM_PATTERN_DB_FLOOR,
"dbFloor",
);
if (dbFloor < -300 || dbFloor > -1) {
throw new RangeError("dbFloor must lie inside -300..-1 dB.");
}
const maximumEvaluations = options?.maximumEvaluations ??
MAXIMUM_RECEIVE_BEAM_PATTERN_EVALUATIONS;
if (!Number.isSafeInteger(maximumEvaluations) || maximumEvaluations < 1) {
throw new RangeError("maximumEvaluations must be a positive safe integer.");
}
const microphoneCount = positions.length;
const sampleEvaluationCount = microphoneCount * sampleCount;
if (
!Number.isSafeInteger(sampleEvaluationCount) ||
sampleEvaluationCount > maximumEvaluations
) {
throw new RangeError(
"receive beam-pattern sphere requires " + sampleEvaluationCount +
" microphone-direction evaluations; limit is " +
maximumEvaluations + ".",
);
}
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 spatialBaselineComparisons =
microphoneCount * (microphoneCount - 1) / 2;
if (spatialBaselineComparisons > maximumBaselineComparisons) {
throw new RangeError(
"receive beam-pattern sphere requires " + spatialBaselineComparisons +
" spatial-baseline comparisons; limit is " +
maximumBaselineComparisons + ".",
);
}
const centroid = [0, 0, 0];
for (const position of positions) {
centroid[0] += position[0] / microphoneCount;
centroid[1] += position[1] / microphoneCount;
centroid[2] += position[2] / microphoneCount;
}
const centeredPositions = positions.map((position) => [
position[0] - centroid[0],
position[1] - centroid[1],
position[2] - centroid[2],
]);
const spatialApertureM = maximumSpatialBaseline(centeredPositions);
const waveNumber = TAU * frequencyHz / speedOfSound;
const maximumPhaseMagnitudeRadBound = 2 * waveNumber * spatialApertureM;
if (
!Number.isFinite(waveNumber) ||
!Number.isFinite(maximumPhaseMagnitudeRadBound)
) {
throw new RangeError(
"frequency, speed of sound, and array coordinates overflow phase.",
);
}
const azimuthStepDeg = 360 / azimuthSampleCount;
const elevationStepDeg = 180 / (elevationSampleCount - 1);
const azimuthStepRad = TAU / azimuthSampleCount;
const elevationStepRad = Math.PI / (elevationSampleCount - 1);
const azimuthAdjacentDirectionChordBound =
2 * Math.sin(azimuthStepRad / 2);
const elevationAdjacentDirectionChordBound =
2 * Math.sin(elevationStepRad / 2);
const maximumCellDiagonalDirectionChordBound = 2 * Math.hypot(
Math.sin(elevationStepRad / 2),
Math.cos(elevationStepRad / 2) * Math.sin(azimuthStepRad / 2),
);
const maximumAzimuthAdjacentPhaseStepRadBound =
waveNumber * spatialApertureM * azimuthAdjacentDirectionChordBound;
const maximumElevationAdjacentPhaseStepRadBound =
waveNumber * spatialApertureM * elevationAdjacentDirectionChordBound;
const maximumCellDiagonalPhaseStepRadBound =
waveNumber * spatialApertureM *
maximumCellDiagonalDirectionChordBound;
const resolutionCriterionSatisfied =
maximumCellDiagonalPhaseStepRadBound <=
RECEIVE_BEAM_PATTERN_METRIC_PHASE_STEP_LIMIT_RAD *
(1 + 64 * Number.EPSILON);
const azimuthAxisDeg = new Float64Array(azimuthSampleCount);
const elevationAxisDeg = new Float64Array(elevationSampleCount);
const azimuthCosines = new Float64Array(azimuthSampleCount);
const azimuthSines = new Float64Array(azimuthSampleCount);
const elevationCosines = new Float64Array(elevationSampleCount);
const elevationSines = new Float64Array(elevationSampleCount);
for (let azimuthIndex = 0; azimuthIndex < azimuthSampleCount; azimuthIndex += 1) {
const azimuthDeg = azimuthIndex * azimuthStepDeg;
const azimuthRadians = azimuthIndex * azimuthStepRad;
azimuthAxisDeg[azimuthIndex] = azimuthDeg;
azimuthCosines[azimuthIndex] = Math.cos(azimuthRadians);
azimuthSines[azimuthIndex] = Math.sin(azimuthRadians);
}
for (
let elevationIndex = 0;
elevationIndex < elevationSampleCount;
elevationIndex += 1
) {
const elevationDeg = elevationIndex === elevationSampleCount - 1
? 90
: -90 + elevationIndex * elevationStepDeg;
const elevationRadians = elevationDeg * RADIANS_PER_DEGREE;
elevationAxisDeg[elevationIndex] = elevationDeg;
if (elevationIndex === 0 || elevationIndex === elevationSampleCount - 1) {
elevationCosines[elevationIndex] = 0;
elevationSines[elevationIndex] = elevationIndex === 0 ? -1 : 1;
} else {
elevationCosines[elevationIndex] = Math.cos(elevationRadians);
elevationSines[elevationIndex] = Math.sin(elevationRadians);
}
}
const steeringDirection = receiveBeamPatternSphereDirection(
steeringAzimuthDeg,
steeringElevationDeg,
);
const steeringProjections = new Float64Array(microphoneCount);
for (let microphoneIndex = 0; microphoneIndex < microphoneCount; microphoneIndex += 1) {
const position = centeredPositions[microphoneIndex];
steeringProjections[microphoneIndex] =
position[0] * steeringDirection[0] +
position[1] * steeringDirection[1] +
position[2] * steeringDirection[2];
}
const responseReal = new Float64Array(sampleCount);
const responseImaginary = new Float64Array(sampleCount);
const responsePower = new Float64Array(sampleCount);
const responseLinear = new Float64Array(sampleCount);
const responseDb = new Float64Array(sampleCount);
const displayDb = new Float64Array(sampleCount);
let sampledMaximumLinear = 0;
let sampledMaximumPower = 0;
for (
let elevationIndex = 0;
elevationIndex < elevationSampleCount;
elevationIndex += 1
) {
const cosineElevation = elevationCosines[elevationIndex];
const directionZ = elevationSines[elevationIndex];
for (
let azimuthIndex = 0;
azimuthIndex < azimuthSampleCount;
azimuthIndex += 1
) {
const directionX = cosineElevation * azimuthCosines[azimuthIndex];
const directionY = cosineElevation * azimuthSines[azimuthIndex];
let summedReal = 0;
let summedImaginary = 0;
for (
let microphoneIndex = 0;
microphoneIndex < microphoneCount;
microphoneIndex += 1
) {
const position = centeredPositions[microphoneIndex];
const arrivalProjection =
position[0] * directionX +
position[1] * directionY +
position[2] * directionZ;
const phase = waveNumber *
(arrivalProjection - steeringProjections[microphoneIndex]);
summedReal += Math.cos(phase);
summedImaginary += Math.sin(phase);
}
const flatIndex =
elevationIndex * azimuthSampleCount + azimuthIndex;
const real = summedReal / microphoneCount;
const imaginary = summedImaginary / microphoneCount;
const power = real * real + imaginary * imaginary;
if (
!Number.isFinite(real) ||
!Number.isFinite(imaginary) ||
!Number.isFinite(power)
) {
throw new RangeError(
"receive beam-pattern sphere phasor accumulation became non-finite.",
);
}
const magnitude = Math.sqrt(power);
const decibels = power > 0 ? 10 * Math.log10(power) : -Infinity;
responseReal[flatIndex] = real;
responseImaginary[flatIndex] = imaginary;
responsePower[flatIndex] = power;
responseLinear[flatIndex] = magnitude;
responseDb[flatIndex] = decibels;
displayDb[flatIndex] = Math.max(dbFloor, decibels);
sampledMaximumLinear = Math.max(sampledMaximumLinear, magnitude);
sampledMaximumPower = Math.max(sampledMaximumPower, power);
}
}
const returnedFloat64Count =
sampleCount * 6 + azimuthSampleCount + elevationSampleCount;
return {
model: "far-field-narrowband-conventional-delay-and-sum-receive-array-factor",
surface: "full-sphere-equirectangular-samples",
assumptions: Object.freeze({
wavefront: "far-field-plane-wave",
bandwidth: "narrowband",
elements: "ideal-calibrated-omnidirectional-point-microphones",
environment: "free-field",
coordinates: "meters",
}),
coordinateSystem: Object.freeze({
handedness: "right-handed",
azimuth: "XY-plane-from-+X-toward-+Y",
elevation: "from-XY-plane-toward-+Z",
directionEquation:
"[cos(elevation)*cos(azimuth), cos(elevation)*sin(azimuth), sin(elevation)]",
}),
frequencyHz,
speedOfSound,
wavelengthM: speedOfSound / frequencyHz,
waveNumberRadPerM: waveNumber,
steeringAzimuthDeg,
steeringElevationDeg,
steeringDirection,
microphoneCount,
azimuthSampleCount,
elevationSampleCount,
sampleCount,
gridShape: [elevationSampleCount, azimuthSampleCount],
flatIndexing: "elevation-major-row-major",
evaluationCount: sampleEvaluationCount,
dbFloor,
grid: {
topology: "equirectangular-azimuth-elevation",
shape: [elevationSampleCount, azimuthSampleCount],
flatIndex: "elevationIndex * azimuthSampleCount + azimuthIndex",
azimuthEndpointDuplicated: false,
elevationEndpointsIncluded: true,
poleAzimuthSamplesCoincident: true,
uniqueDirectionCount:
azimuthSampleCount * (elevationSampleCount - 2) + 2,
interpolation: "none",
smoothing: "none",
sampledMaximumNormalization: false,
},
workload: {
spatialBaselineComparisons,
maximumBaselineComparisons,
sampleEvaluationCount,
maximumEvaluationCount: sampleEvaluationCount,
maximumEvaluations,
returnedFloat64Count,
returnedFloat64Bytes:
returnedFloat64Count * Float64Array.BYTES_PER_ELEMENT,
maximumSurfaceSamples: MAXIMUM_RECEIVE_BEAM_PATTERN_SPHERE_SAMPLES,
},
angularSampling: {
mode: "fixed-equirectangular-grid",
criterion: "maximum-cell-diagonal-pair-phase-step",
azimuthStepDeg,
elevationStepDeg,
spatialApertureM,
maximumPhaseMagnitudeRadBound,
azimuthAdjacentDirectionChordBound,
elevationAdjacentDirectionChordBound,
maximumCellDiagonalDirectionChordBound,
maximumAzimuthAdjacentPhaseStepRadBound,
maximumElevationAdjacentPhaseStepRadBound,
maximumCellDiagonalPhaseStepRadBound,
phaseStepLimitRad:
RECEIVE_BEAM_PATTERN_METRIC_PHASE_STEP_LIMIT_RAD,
resolutionCriterionSatisfied,
interpretation:
"sampling-bound-only-not-proof-of-complete-lobe-detection",
},
azimuthAxisDeg,
elevationAxisDeg,
responseReal,
responseImaginary,
responsePower,
responseLinear,
responseDb,
displayDb,
sampledSummary: {
maximumLinearOnGrid: sampledMaximumLinear,
maximumPowerOnGrid: sampledMaximumPower,
maximumDbOnGrid: sampledMaximumPower > 0
? 10 * Math.log10(sampledMaximumPower)
: -Infinity,
},
};
}