Files

623 lines
19 KiB
JavaScript

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" });
}