micarrsim/packages/acoustics-core/cfar.js

const CFAR_FRAME_RESULT_KIND = "micarrsim.cfar-frame-result/1";

export const DEFAULT_OS_CFAR_MAXIMUM_BISECTION_ITERATIONS = 256;

function requirePositiveSafeInteger(value, label) {
  if (!Number.isSafeInteger(value) || value < 1) {
    throw new TypeError(`${label} must be a positive safe integer.`);
  }
  return value;
}

function requireNonNegativeSafeInteger(value, label) {
  if (!Number.isSafeInteger(value) || value < 0) {
    throw new TypeError(`${label} must be a non-negative safe integer.`);
  }
  return value;
}

function requirePositiveFinite(value, label) {
  if (!Number.isFinite(value) || value <= 0) {
    throw new TypeError(`${label} must be a positive finite number.`);
  }
  return value;
}

function requireProbability(value, label) {
  if (!Number.isFinite(value) || value <= 0 || value >= 1) {
    throw new RangeError(`${label} must be finite and strictly inside 0..1.`);
  }
  return value;
}

function requireTrainingCellCount(trainingCellCount) {
  return requirePositiveSafeInteger(trainingCellCount, "trainingCellCount");
}

function requireOsRank(rank, trainingCellCount) {
  if (
    !Number.isSafeInteger(rank) ||
    rank < 1 ||
    rank > trainingCellCount
  ) {
    throw new RangeError(
      `rank must be a 1-based ascending order-statistic index inside 1..${trainingCellCount}.`,
    );
  }
  return rank;
}

/**
 * CA-CFAR multiplier for an arithmetic-mean reference statistic.
 *
 * Under independent, identically distributed exponential CUT/reference
 * powers, Pfa = (1 + alpha / N)^(-N). The `expm1` form preserves accuracy
 * when the requested probability is close to one or N is large.
 */
export function computeCaCfarThresholdMultiplier({
  trainingCellCount,
  nominalFalseAlarmProbability,
}) {
  const count = requireTrainingCellCount(trainingCellCount);
  const probability = requireProbability(
    nominalFalseAlarmProbability,
    "nominalFalseAlarmProbability",
  );
  const multiplier = count * Math.expm1(-Math.log(probability) / count);
  if (!Number.isFinite(multiplier) || multiplier <= 0) {
    throw new RangeError(
      "CA-CFAR threshold multiplier exceeds the finite positive Float64 range.",
    );
  }
  return multiplier;
}

/**
 * Log false-alarm probability for OS-CFAR with an ascending 1-based rank.
 *
 * log(Pfa) = sum_(j=0)^(k-1) -log(1 + alpha / (N - j)).
 * The log product avoids factorial/gamma overflow and is also used by the
 * bracketed multiplier solver below.
 */
export function computeOsCfarLogFalseAlarmProbability({
  trainingCellCount,
  rank,
  thresholdMultiplier,
}) {
  const count = requireTrainingCellCount(trainingCellCount);
  const selectedRank = requireOsRank(rank, count);
  if (!Number.isFinite(thresholdMultiplier) || thresholdMultiplier < 0) {
    throw new RangeError(
      "thresholdMultiplier must be finite and non-negative.",
    );
  }

  let sum = 0;
  let compensation = 0;
  for (let offset = 0; offset < selectedRank; offset += 1) {
    const term = -Math.log1p(thresholdMultiplier / (count - offset));
    const next = sum + term;
    const correction = Math.abs(sum) >= Math.abs(term)
      ? (sum - next) + term
      : (term - next) + sum;
    sum = next;
    compensation += correction;
  }
  return sum + compensation;
}

function solveOsCfarThresholdMultiplier({
  trainingCellCount,
  rank,
  nominalFalseAlarmProbability,
  maximumIterations,
}) {
  const count = requireTrainingCellCount(trainingCellCount);
  const selectedRank = requireOsRank(rank, count);
  const probability = requireProbability(
    nominalFalseAlarmProbability,
    "nominalFalseAlarmProbability",
  );
  const iterationLimit = requirePositiveSafeInteger(
    maximumIterations,
    "maximumIterations",
  );
  const targetLogProbability = Math.log(probability);
  const evaluate = (thresholdMultiplier) =>
    computeOsCfarLogFalseAlarmProbability({
      trainingCellCount: count,
      rank: selectedRank,
      thresholdMultiplier,
    });

  let lowerMultiplier = 0;
  let upperMultiplier = 1;
  let upperLogProbability = evaluate(upperMultiplier);
  let bracketDoublings = 0;
  while (upperLogProbability > targetLogProbability) {
    if (upperMultiplier >= Number.MAX_VALUE / 2) {
      upperMultiplier = Number.MAX_VALUE;
      upperLogProbability = evaluate(upperMultiplier);
      if (upperLogProbability > targetLogProbability) {
        throw new RangeError(
          "OS-CFAR threshold multiplier exceeds the finite Float64 range.",
        );
      }
      break;
    }
    upperMultiplier *= 2;
    upperLogProbability = evaluate(upperMultiplier);
    bracketDoublings += 1;
  }

  const logarithmTolerance =
    64 * Number.EPSILON *
    Math.max(1, selectedRank, Math.abs(targetLogProbability));
  let iterations = 0;
  let multiplier = upperMultiplier;
  let logarithmResidual = upperLogProbability - targetLogProbability;
  let converged = Math.abs(logarithmResidual) <= logarithmTolerance;

  while (!converged && iterations < iterationLimit) {
    const midpoint =
      lowerMultiplier + (upperMultiplier - lowerMultiplier) / 2;
    if (midpoint === lowerMultiplier || midpoint === upperMultiplier) {
      multiplier = Math.abs(logarithmResidual) <= logarithmTolerance
        ? upperMultiplier
        : midpoint;
      const midpointLogProbability = evaluate(multiplier);
      logarithmResidual = midpointLogProbability - targetLogProbability;
      converged =
        Math.abs(logarithmResidual) <= logarithmTolerance ||
        midpoint === lowerMultiplier ||
        midpoint === upperMultiplier;
      break;
    }

    const midpointLogProbability = evaluate(midpoint);
    logarithmResidual = midpointLogProbability - targetLogProbability;
    multiplier = midpoint;
    iterations += 1;
    if (Math.abs(logarithmResidual) <= logarithmTolerance) {
      converged = true;
      break;
    }
    if (midpointLogProbability > targetLogProbability) {
      lowerMultiplier = midpoint;
    } else {
      upperMultiplier = midpoint;
      upperLogProbability = midpointLogProbability;
    }
  }

  if (!converged) {
    throw new RangeError(
      `OS-CFAR threshold multiplier did not converge within ${iterationLimit} bisection iterations.`,
    );
  }
  if (!Number.isFinite(multiplier) || multiplier <= 0) {
    throw new RangeError(
      "OS-CFAR threshold multiplier exceeds the finite positive Float64 range.",
    );
  }

  return Object.freeze({
    multiplier,
    iterations,
    bracketDoublings,
    targetLogProbability,
    logarithmResidual,
    logarithmTolerance,
    converged: true,
  });
}

/** Solve the standard exponential-background OS-CFAR multiplier by bisection. */
export function computeOsCfarThresholdMultiplier({
  trainingCellCount,
  rank,
  nominalFalseAlarmProbability,
  maximumIterations = DEFAULT_OS_CFAR_MAXIMUM_BISECTION_ITERATIONS,
}) {
  return solveOsCfarThresholdMultiplier({
    trainingCellCount,
    rank,
    nominalFalseAlarmProbability,
    maximumIterations,
  }).multiplier;
}

function requirePowerSpectrum(powerSpectrum, fftSize) {
  if (powerSpectrum == null || !Number.isSafeInteger(powerSpectrum.length)) {
    throw new TypeError("powerSpectrum must be a finite numeric array-like value.");
  }
  const expectedBinCount = fftSize / 2 + 1;
  if (powerSpectrum.length !== expectedBinCount) {
    throw new RangeError(
      `powerSpectrum must contain exactly ${expectedBinCount} one-sided bins for fftSize ${fftSize}.`,
    );
  }
  for (let bin = 0; bin < expectedBinCount; bin += 1) {
    const power = powerSpectrum[bin];
    if (typeof power !== "number" || !Number.isFinite(power)) {
      throw new TypeError(`powerSpectrum[${bin}] must be finite.`);
    }
    if (power < 0) {
      throw new RangeError(`powerSpectrum[${bin}] must be non-negative.`);
    }
  }
  return expectedBinCount;
}

function scaledArithmeticMean(values) {
  let scale = 0;
  for (const value of values) scale = Math.max(scale, value);
  if (scale === 0) return 0;

  let sum = 0;
  let compensation = 0;
  for (const value of values) {
    const normalized = value / scale;
    const next = sum + normalized;
    const correction = Math.abs(sum) >= Math.abs(normalized)
      ? (sum - next) + normalized
      : (normalized - next) + sum;
    sum = next;
    compensation += correction;
  }
  return ((sum + compensation) / values.length) * scale;
}

function freezeArray(values) {
  return Object.freeze(values);
}

function createFrequencyMetadata(bin, binWidthHz, nyquistHz) {
  return {
    adcBin: bin,
    adcCenterFrequencyHz: bin * binWidthHz,
    adcFftCellHz: Object.freeze({
      minimumHz: Math.max(0, (bin - 0.5) * binWidthHz),
      maximumHz: Math.min(nyquistHz, (bin + 0.5) * binWidthHz),
    }),
  };
}

function createUntestedBin({
  bin,
  cutPower,
  binWidthHz,
  nyquistHz,
  reason,
}) {
  return Object.freeze({
    ...createFrequencyMetadata(bin, binWidthHz, nyquistHz),
    status: "untested",
    reason,
    cutPower,
    referencePower: null,
    thresholdPower: null,
    thresholdRatio: null,
    thresholdRatioStatus: "unavailable-untested-bin",
  });
}

function resolveDetectorConfiguration(options, nyquistBin) {
  const detector = options.detector;
  if (detector !== "ca-cfar" && detector !== "os-cfar") {
    throw new RangeError('detector must be "ca-cfar" or "os-cfar".');
  }
  const trainingCellsPerSide = requirePositiveSafeInteger(
    options.trainingCellsPerSide,
    "trainingCellsPerSide",
  );
  const guardCellsPerSide = requireNonNegativeSafeInteger(
    options.guardCellsPerSide,
    "guardCellsPerSide",
  );
  if (
    trainingCellsPerSide >
    Number.MAX_SAFE_INTEGER - guardCellsPerSide
  ) {
    throw new RangeError("CFAR reference radius exceeds the safe integer range.");
  }
  const referenceRadius = trainingCellsPerSide + guardCellsPerSide;
  if (referenceRadius > (nyquistBin - 2) / 2) {
    throw new RangeError(
      "CFAR training and guard cells leave no interior bin with a full symmetric reference window.",
    );
  }
  const trainingCellCount = 2 * trainingCellsPerSide;
  if (!Number.isSafeInteger(trainingCellCount)) {
    throw new RangeError("trainingCellCount exceeds the safe integer range.");
  }
  const nominalFalseAlarmProbability = requireProbability(
    options.nominalFalseAlarmProbability,
    "nominalFalseAlarmProbability",
  );

  if (detector === "ca-cfar") {
    if (options.rank !== undefined) {
      throw new RangeError("rank is only valid for OS-CFAR.");
    }
    return {
      detector,
      trainingCellsPerSide,
      guardCellsPerSide,
      referenceRadius,
      trainingCellCount,
      nominalFalseAlarmProbability,
      rank: null,
      thresholdMultiplier: computeCaCfarThresholdMultiplier({
        trainingCellCount,
        nominalFalseAlarmProbability,
      }),
      thresholdMultiplierDerivation:
        "N*expm1(-log(nominalPfa)/N)-for-arithmetic-mean-reference",
      osSolver: null,
    };
  }

  const rank = requireOsRank(options.rank, trainingCellCount);
  const osSolver = solveOsCfarThresholdMultiplier({
    trainingCellCount,
    rank,
    nominalFalseAlarmProbability,
    maximumIterations:
      options.maximumOsBisectionIterations ??
      DEFAULT_OS_CFAR_MAXIMUM_BISECTION_ITERATIONS,
  });
  return {
    detector,
    trainingCellsPerSide,
    guardCellsPerSide,
    referenceRadius,
    trainingCellCount,
    nominalFalseAlarmProbability,
    rank,
    thresholdMultiplier: osSolver.multiplier,
    thresholdMultiplierDerivation:
      "log-product-bracketed-bisection-for-ascending-order-statistic",
    osSolver,
  };
}

/**
 * Detect narrowband activity in one complete, single-channel, one-sided power
 * spectrum. No temporal/microphone averaging, bin grouping, peak selection,
 * persistence, clustering, or detector fusion is performed.
 */
export function detectCfarFrame(options) {
  if (!options || typeof options !== "object") {
    throw new TypeError("CFAR options must be an object.");
  }
  const sampleRate = requirePositiveFinite(options.sampleRate, "sampleRate");
  const fftSize = requirePositiveSafeInteger(options.fftSize, "fftSize");
  if (fftSize < 4 || fftSize % 2 !== 0) {
    throw new RangeError("fftSize must be an even integer greater than or equal to four.");
  }
  if (options.window !== "hann") {
    throw new RangeError('window must be "hann".');
  }
  if (options.frameComplete !== true) {
    throw new RangeError("frameComplete must be true; partial frames are not tested.");
  }
  if (typeof options.channelId !== "string" || options.channelId.length === 0) {
    throw new TypeError("channelId must be a non-empty string.");
  }
  const frameIndex = requireNonNegativeSafeInteger(
    options.frameIndex,
    "frameIndex",
  );
  const frameStartSample = requireNonNegativeSafeInteger(
    options.frameStartSample,
    "frameStartSample",
  );
  const temporalAverageCount = options.temporalAverageCount ?? 1;
  const microphoneAverageCount = options.microphoneAverageCount ?? 1;
  if (temporalAverageCount !== 1) {
    throw new RangeError(
      "temporalAverageCount must equal one; the exponential CFAR multiplier is not valid after temporal averaging.",
    );
  }
  if (microphoneAverageCount !== 1) {
    throw new RangeError(
      "microphoneAverageCount must equal one; microphone averaging or fusion is not performed.",
    );
  }

  const binCount = requirePowerSpectrum(options.powerSpectrum, fftSize);
  const nyquistBin = fftSize / 2;
  const configuration = resolveDetectorConfiguration(options, nyquistBin);
  const {
    detector,
    trainingCellsPerSide,
    guardCellsPerSide,
    referenceRadius,
    trainingCellCount,
    nominalFalseAlarmProbability,
    rank,
    thresholdMultiplier,
    thresholdMultiplierDerivation,
    osSolver,
  } = configuration;
  const nyquistHz = sampleRate / 2;
  const binWidthHz = sampleRate / fftSize;
  const firstTestedBin = referenceRadius + 1;
  const lastTestedBin = nyquistBin - referenceRadius - 1;
  const referencePowers = new Float64Array(trainingCellCount);
  const bins = new Array(binCount);
  let testedBinCount = 0;
  let detectedBinCount = 0;

  for (let bin = 0; bin < binCount; bin += 1) {
    const cutPower = options.powerSpectrum[bin];
    let untestedReason = null;
    if (bin === 0) {
      untestedReason = "real-only-dc-bin";
    } else if (bin === nyquistBin) {
      untestedReason = "real-only-nyquist-bin";
    } else if (bin < firstTestedBin || bin > lastTestedBin) {
      untestedReason = "insufficient-full-symmetric-reference-window";
    }
    if (untestedReason !== null) {
      bins[bin] = createUntestedBin({
        bin,
        cutPower,
        binWidthHz,
        nyquistHz,
        reason: untestedReason,
      });
      continue;
    }

    let referenceIndex = 0;
    const leftStart = bin - guardCellsPerSide - trainingCellsPerSide;
    const leftEnd = bin - guardCellsPerSide;
    for (let referenceBin = leftStart; referenceBin < leftEnd; referenceBin += 1) {
      referencePowers[referenceIndex] = options.powerSpectrum[referenceBin];
      referenceIndex += 1;
    }
    const rightStart = bin + guardCellsPerSide + 1;
    const rightEnd = rightStart + trainingCellsPerSide;
    for (let referenceBin = rightStart; referenceBin < rightEnd; referenceBin += 1) {
      referencePowers[referenceIndex] = options.powerSpectrum[referenceBin];
      referenceIndex += 1;
    }
    if (referenceIndex !== trainingCellCount) {
      throw new Error("CFAR reference-window construction violated fixed N=2R.");
    }

    let referencePower;
    if (detector === "ca-cfar") {
      referencePower = scaledArithmeticMean(referencePowers);
    } else {
      referencePowers.sort();
      referencePower = referencePowers[rank - 1];
    }
    const thresholdPower = thresholdMultiplier * referencePower;
    if (!Number.isFinite(thresholdPower)) {
      throw new RangeError(
        `CFAR threshold exceeds the finite Float64 range at bin ${bin}.`,
      );
    }
    const detected = cutPower > thresholdPower;
    let thresholdRatio;
    let thresholdRatioStatus;
    if (thresholdPower === 0) {
      if (cutPower === 0) {
        thresholdRatio = null;
        thresholdRatioStatus = "indeterminate-zero-over-zero";
      } else {
        thresholdRatio = null;
        thresholdRatioStatus = "positive-infinity-zero-threshold";
      }
    } else {
      thresholdRatio = cutPower / thresholdPower;
      thresholdRatioStatus = "finite";
    }
    bins[bin] = Object.freeze({
      ...createFrequencyMetadata(bin, binWidthHz, nyquistHz),
      status: detected ? "detected" : "not-detected",
      reason: null,
      cutPower,
      referencePower,
      thresholdPower,
      thresholdRatio,
      thresholdRatioStatus,
    });
    testedBinCount += 1;
    if (detected) detectedBinCount += 1;
  }

  return {
    kind: CFAR_FRAME_RESULT_KIND,
    detector,
    frame: Object.freeze({
      channelId: options.channelId,
      frameIndex,
      frameStartSample,
      complete: true,
      temporalAverageCount: 1,
      microphoneAverageCount: 1,
    }),
    spectrum: Object.freeze({
      domain: "one-sided-real-adc-power-spectrum",
      statistic: "caller-supplied-non-negative-power-per-bin",
      sampleRate,
      fftSize,
      binCount,
      binWidthHz,
      nyquistHz,
      window: "hann",
      windowDefinition: "symmetric-hann-N-minus-one",
      analogFrequencyOrAliasOrderInferred: false,
    }),
    parameters: Object.freeze({
      trainingCellsPerSide,
      guardCellsPerSide,
      trainingCellCount,
      referenceStatistic: detector === "ca-cfar"
        ? "arithmetic-mean-power"
        : "ascending-order-statistic-power",
      rank,
      rankConvention: detector === "os-cfar" ? "1-based-ascending" : null,
      nominalFalseAlarmProbability,
      nominalPfaPerTestedCutUnderIidExponentialModel:
        nominalFalseAlarmProbability,
      thresholdMultiplier,
      thresholdMultiplierDerivation,
      thresholdComparison: "cut-power-strictly-greater-than-threshold-power",
      firstTestedBin: firstTestedBin <= lastTestedBin
        ? firstTestedBin
        : null,
      lastTestedBin: firstTestedBin <= lastTestedBin
        ? lastTestedBin
        : null,
      testedBinCount,
      untestedBinCount: binCount - testedBinCount,
      detectedBinCount,
      osSolver,
    }),
    calibration: Object.freeze({
      status: "nominal-model-only-not-empirically-calibrated",
      falseAlarmScope: "one-tested-bin-in-one-frame",
      model: "independent-identically-distributed-exponential-power-with-common-mean",
      calibratedFalseAlarmProbability: null,
      globalFrameFalseAlarmProbability: null,
      assumptions: freezeArray([
        "cut-and-training-powers-are-independent",
        "cut-and-training-powers-have-one-common-local-mean",
        "powers-are-exponential-square-law-complex-gaussian-statistics",
      ]),
      notCalibratedReasons: freezeArray([
        "hann-window-bin-independence-is-not-guaranteed",
        "colored-acoustic-noise-and-microphone-response-are-not-verified-homogeneous",
        "overlapping-sliding-reference-windows-share-cells",
      ]),
    }),
    processing: Object.freeze({
      completeFramesOnly: true,
      temporalAveraging: false,
      microphoneAveraging: false,
      peakSelection: false,
      topKSelection: false,
      clustering: false,
      contiguousBinMerging: false,
      temporalPersistence: false,
      detectorFusion: false,
      microphoneFusion: false,
      bandFusion: false,
    }),
    bins: Object.freeze(bins),
  };
}

export function detectCaCfarFrame(options) {
  return detectCfarFrame({ ...options, detector: "ca-cfar" });
}

export function detectOsCfarFrame(options) {
  return detectCfarFrame({ ...options, detector: "os-cfar" });
}