micarrsim/packages/acoustics-core/beam-pattern.js

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 RECEIVE_BEAM_PATTERN_METRIC_PHASE_STEP_LIMIT_RAD = Math.PI / 8;

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];
}

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,
) {
  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: "full-circle-above-half-power",
      widthDeg: null,
      startAngleDeg: null,
      endAngleDeg: null,
      fullCircleAboveHalfPower: true,
    };
  }

  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: "full-circle-above-half-power",
      widthDeg: null,
      startAngleDeg: null,
      endAngleDeg: null,
      fullCircleAboveHalfPower: true,
    };
  }

  const leftBoundaryOffset = bisectHalfPowerBoundary(
    -leftSteps * stepDeg,
    -(leftSteps + 1) * stepDeg,
    evaluatePowerAtOffset,
  );
  const rightBoundaryOffset = bisectHalfPowerBoundary(
    rightSteps * stepDeg,
    (rightSteps + 1) * stepDeg,
    evaluatePowerAtOffset,
  );
  return {
    status: "bounded",
    widthDeg: rightBoundaryOffset - leftBoundaryOffset,
    startAngleDeg: normalizeBeamPatternAngle(steeringAngleDeg + leftBoundaryOffset),
    endAngleDeg: normalizeBeamPatternAngle(steeringAngleDeg + rightBoundaryOffset),
    fullCircleAboveHalfPower: false,
  };
}

function sampledLocalMaxima(responseLinear, anglesDeg, offsetsDeg) {
  const sampleCount = responseLinear.length;
  let minimum = Infinity;
  let maximum = -Infinity;
  for (const value of responseLinear) {
    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 peaks = [];
  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,
      });
    }
  }
  peaks.sort((left, right) =>
    right.responseLinear - left.responseLinear ||
    Math.abs(left.offsetDeg) - Math.abs(right.offsetDeg) ||
    left.angleDeg - right.angleDeg);
  return peaks;
}

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 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");
  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 microphoneCount = positions.length;
  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;
  if (!Number.isFinite(waveNumber) || !Number.isFinite(waveNumber * projectedApertureM)) {
    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;
    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);
    }
    return (real * real + imaginary * imaginary) /
      (microphoneCount * microphoneCount);
  };
  const halfPower = metricResolutionSatisfied
    ? halfPowerBeamwidth(
      responsePower,
      steeringIndex,
      stepDeg,
      steeringAngleDeg,
      evaluatePowerAtOffset,
    )
    : {
      status: "resolution-unresolved",
      widthDeg: null,
      startAngleDeg: null,
      endAngleDeg: null,
      fullCircleAboveHalfPower: null,
    };
  const peaks = metricResolutionSatisfied
    ? sampledLocalMaxima(responseLinear, anglesDeg, offsetsDeg)
    : [];
  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:
        microphoneCount * (microphoneCount - 1) / 2,
      sampleEvaluationCount,
      maximumMetricEvaluationCount,
      maximumEvaluationCount: evaluationCount,
    },
    angularSampling: {
      mode: requestedSampleCount === "auto" ? "auto" : "fixed",
      criterion: "maximum-adjacent-pair-phase-step",
      projectedApertureM,
      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: {
        status: metricResolutionSatisfied ? "sampled" : "resolution-unresolved",
        peaks,
      },
    },
  };
}