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