From 3ebc0f84e0256942dbd6c70fe9487149b798f97a Mon Sep 17 00:00:00 2001 From: MasterBabol Date: Sat, 29 Aug 2026 13:16:08 +0900 Subject: [PATCH] feat: add standard full-band activity detectors --- docs/activity-detection.md | 106 ++++ packages/acoustics-core/activity-monitor.js | 345 +++++++++++ packages/acoustics-core/cfar.js | 622 ++++++++++++++++++++ packages/acoustics-core/index.js | 3 + packages/acoustics-core/spectral-flux.js | 419 +++++++++++++ tests/test-runner.js | 3 + tests/unit/activity-monitor.test.js | 179 ++++++ tests/unit/cfar.test.js | 465 +++++++++++++++ tests/unit/spectral-flux.test.js | 372 ++++++++++++ 9 files changed, 2514 insertions(+) create mode 100644 docs/activity-detection.md create mode 100644 packages/acoustics-core/activity-monitor.js create mode 100644 packages/acoustics-core/cfar.js create mode 100644 packages/acoustics-core/spectral-flux.js create mode 100644 tests/unit/activity-monitor.test.js create mode 100644 tests/unit/cfar.test.js create mode 100644 tests/unit/spectral-flux.test.js diff --git a/docs/activity-detection.md b/docs/activity-detection.md new file mode 100644 index 0000000..c38df03 --- /dev/null +++ b/docs/activity-detection.md @@ -0,0 +1,106 @@ +# 전 대역 활동 검출 모델 + +이 문서는 위치 추정 전에 ADC 대역에서 관측되는 신호 변화를 찾는 보조 분석의 +의미와 한계를 정의한다. 검출기는 Scene의 음원 주파수·위치·개수를 입력으로 받지 +않으며, 검출 결과를 위치 정답이나 분석 대역으로 자동 변환하지 않는다. + +## 관측 경계 + +- 입력은 사용자가 고른 **물리 마이크 채널 하나**의 시간 샘플이다. 여러 채널의 + waveform 평균, power 평균, 최대값, 투표는 배열 형상과 TDOA에 따라 다른 결과를 + 만들기 때문에 기본 알고리즘에 포함하지 않는다. +- 실수 ADC의 관측 범위는 \(0\ldots F_s/2\)이다. anti-alias filter가 없는 현재 + 시뮬레이션에서는 \(F_s/2\)를 넘는 아날로그 성분도 샘플링되지만, 검출기가 보는 + 것은 접힌 ADC bin이다. +- 검출된 ADC bin만으로 아날로그 alias order, spectral inversion 전의 원 주파수, + 음원 수 또는 위치를 식별할 수 없다. 배열 설계 비교에서는 같은 ADC bin에 + 대응하는 아날로그 주파수 가설을 별도 Experiment band로 명시해 비교한다. +- 정확한 DC와 Nyquist bin은 실수 성분 하나뿐이어서 내부 complex FFT bin과 + 통계 모델이 다르다. CFAR의 CUT와 reference cell에서 제외한다. + +## CA-CFAR와 OS-CFAR + +CFAR는 한 STFT frame 안의 좁은 spectral outlier를 찾는다. frame 또는 마이크 +사이의 평균을 내지 않고 각 frame을 독립적으로 검사한다. 좌우 training cell을 +각각 \(R\)개, guard cell을 각각 \(G\)개 쓰며 \(N=2R\)이다. + +CA-CFAR는 training power의 평균 \(Z\)를 사용한다. + +\[ +Y_b > \alpha Z,\qquad +\alpha=N\left(P_{fa}^{-1/N}-1\right) +\] + +OS-CFAR는 training power를 오름차순 정렬한 뒤 명시한 1-based rank \(k\)의 +\(Y_{(k)}\)를 사용한다. + +\[ +Y_b > \alpha Y_{(k)},\qquad +P_{fa}=\prod_{j=0}^{k-1}\frac{N-j}{N-j+\alpha} +\] + +OS-CFAR multiplier는 이 단조식을 수치적으로 풀어 구한다. 다른 \(N,k\)의 값을 +재사용하지 않는다. 두 방식은 독립 결과이며 합의, winner 선택 또는 위치 map +융합을 하지 않는다. + +여기서 \(P_{fa}\)는 CUT와 모든 training cell의 power가 독립·동일 평균의 지수 +분포라는 square-law complex-Gaussian 모델 아래에서만 성립하는 **tested CUT +하나당 명목 오경보율**이다. posterior probability, confidence, frame 전체의 +오경보율 또는 실측 보정값이 아니다. Hann window의 인접 bin 상관, 마이크 주파수 +응답, colored noise, 잔향, noise-floor 경계와 off-bin leakage는 이 전제를 +깨뜨린다. 결과에는 원 power, reference statistic, threshold, ratio와 +`model-derived-not-empirically-calibrated` 상태를 그대로 남긴다. + +full symmetric reference window를 만들 수 없는 edge bin은 reference 수를 +줄이거나 wrap하지 않고 `untested`로 남긴다. 검출된 모든 bin을 반환하며 top-K, +NMS, gap merge, persistence 또는 임의 대역 padding을 적용하지 않는다. + +CA-CFAR는 [Finn–Johnson의 adaptive threshold +detector](https://www.rsp-italy.it/Electronics/Magazines/RCA%20Review/_contents/RCA%20Review%201968-09.pdf)를 +기준선으로 삼고, OS-CFAR는 다중 target/clutter 환경을 다룬 +[Rohling의 order-statistic CFAR](https://doi.org/10.1109/TAES.1983.309350)를 +따른다. 비균질 배경에서의 한계는 +[Gandhi–Kassam](https://doi.org/10.1109/7.7185), window leakage의 해석은 +[Harris](https://doi.org/10.1109/PROC.1978.10837)를 참고한다. + +## Spectral flux onset + +갑자기 시작되는 broadband 또는 tonal 변화를 보는 별도 지표는 +[Dixon의 linear-magnitude spectral +flux](https://www.dafx.de/paper-archive/2006/papers/p_133.pdf)를 사용한다. +Hamming-window STFT의 연속 frame에 대해 + +\[ +d[n,k]=\max(0,|X[n,k]|-|X[n-1,k]|),\qquad +SF[n]=\sum_k d[n,k] +\] + +를 계산한다. per-bin 양의 차이와 합계 원값을 모두 보존하며 log compression, +전역 정규화, SuperFlux 일부 또는 마이크 융합을 섞지 않는다. + +실시간 peak 선택은 +[Böck·Krebs·Schedl의 causal peak +picker](https://www.cp.jku.at/research/papers/Boeck_etal_ISMIR_2012.pdf)를 +별도 단계로 쓸 수 있다. 현재 frame이 과거 \(w_1\) frame의 local maximum이고, +과거 \(w_3\) frame 평균보다 사용자가 지정한 \(\delta\)만큼 크며, 마지막 +검출에서 \(w_5\) frame보다 멀 때 onset으로 표시한다. \(\delta\)는 확률이 아닌 +운용 임계값이고 환경 데이터로 보정해야 한다. 논문에서 음악 데이터로 최적화한 +값을 보편 default나 물리 상수로 취급하지 않는다. + +Spectral flux는 **onset novelty**이지 지속 activity detector가 아니다. 정상 +상태의 tone은 첫 상승 뒤 flux가 거의 0이 되며, slow attack은 여러 frame에 +분산되어 threshold를 넘지 못할 수 있다. vibrato, tremolo와 확률적 noise +변화는 false trigger를 만들 수 있다. STFT frame을 모두 받은 뒤에야 판단할 수 +있으므로 실제 availability latency는 최소 \(N/F_s\), 갱신 간격은 \(H/F_s\)다. + +## 배열 설계에서의 사용 순서 + +1. 물리 마이크 한 개와 검출 알고리즘·파라미터를 명시한다. +2. raw ADC-bin 검출 또는 raw onset curve와 모든 설정을 비교한다. +3. 관심 ADC bin을 선택하되, 이것을 주파수 정답이나 bandwidth 추정으로 부르지 + 않는다. +4. baseband 또는 명시한 초과-Nyquist 아날로그 가설을 Experiment band로 + 추가한다. +5. 같은 Scene·Array·가설에서 DAS, PHAT, GCC-PHAT, MUSIC의 raw map과 실패 + 양상을 비교한다. 서로 다른 band의 map은 알려진 별도 wideband 방법을 + 구현하기 전까지 자동 결합하지 않는다. diff --git a/packages/acoustics-core/activity-monitor.js b/packages/acoustics-core/activity-monitor.js new file mode 100644 index 0000000..8d1a732 --- /dev/null +++ b/packages/acoustics-core/activity-monitor.js @@ -0,0 +1,345 @@ +import { detectCfarFrame } from "./cfar.js"; +import { cooperativeYield } from "./cooperative-yield.js"; +import { fftInPlace } from "./fft.js"; +import { getCompleteFrameCount, hannWindow } from "./stft.js"; + +export const DEFAULT_MAXIMUM_ACTIVITY_FRAME_COUNT = 4_096; +export const DEFAULT_MAXIMUM_ACTIVITY_BIN_DECISIONS = 2_000_000; +export const DEFAULT_MAXIMUM_ACTIVITY_FFT_BUTTERFLIES = 100_000_000; + +function requirePositiveSafeInteger(value, label) { + if (!Number.isSafeInteger(value) || value < 1) { + throw new RangeError(label + " must be a positive safe integer."); + } + return value; +} + +function requireNonNegativeSafeInteger(value, label) { + if (!Number.isSafeInteger(value) || value < 0) { + throw new RangeError(label + " must be a non-negative safe integer."); + } + return value; +} + +function requirePositiveFinite(value, label) { + if (!Number.isFinite(value) || value <= 0) { + throw new RangeError(label + " must be a positive finite number."); + } + return value; +} + +function requireSingleChannel(samples) { + if ( + (!Array.isArray(samples) && !ArrayBuffer.isView(samples)) || + samples instanceof DataView + ) { + throw new TypeError( + "samples must be one numeric channel; channel collections and DataView are not accepted.", + ); + } + for (let index = 0; index < samples.length; index += 1) { + if (typeof samples[index] !== "number" || !Number.isFinite(samples[index])) { + throw new TypeError("samples[" + index + "] must be finite."); + } + } + return samples; +} + +function checkedProduct(values, label) { + let product = 1; + for (const value of values) { + product *= value; + if (!Number.isSafeInteger(product)) { + throw new RangeError(label + " exceeds the safe integer range."); + } + } + return product; +} + +function throwIfAborted(signal) { + if (!signal?.aborted) return; + const error = new Error("Full-band activity analysis was canceled."); + error.name = "AbortError"; + throw error; +} + +export function estimateCfarActivityWorkload({ + sampleCount, + fftSize, + hopSize, +}) { + requireNonNegativeSafeInteger(sampleCount, "sampleCount"); + requirePositiveSafeInteger(fftSize, "fftSize"); + requirePositiveSafeInteger(hopSize, "hopSize"); + if (fftSize < 4 || (fftSize & (fftSize - 1)) !== 0) { + throw new RangeError("fftSize must be a power of two greater than or equal to four."); + } + if (hopSize > fftSize) { + throw new RangeError("hopSize must not exceed fftSize."); + } + const frameCount = getCompleteFrameCount(sampleCount, fftSize, hopSize); + const binCount = fftSize / 2 + 1; + const binDecisionCount = checkedProduct( + [frameCount, binCount], + "CFAR bin decision count", + ); + const fftButterfliesPerFrame = + (fftSize / 2) * Math.log2(fftSize); + const fftButterflyCount = checkedProduct( + [frameCount, fftButterfliesPerFrame], + "CFAR FFT butterfly count", + ); + const workingBytes = + (3 * fftSize + binCount) * Float64Array.BYTES_PER_ELEMENT; + return Object.freeze({ + frameCount, + binCount, + binDecisionCount, + fftButterfliesPerFrame, + fftButterflyCount, + workingBytes, + outputRecordModel: "one-bin-object-per-frame-without-truncation", + }); +} + +function requireLimit(value, fallback, label) { + return requirePositiveSafeInteger(value ?? fallback, label); +} + +function compactFrameResult(frameResult, fftSize, hopSize) { + return Object.freeze({ + frameIndex: frameResult.frame.frameIndex, + frameStartSample: frameResult.frame.frameStartSample, + frameCenterSample: + frameResult.frame.frameStartSample + (fftSize - 1) / 2, + frameAvailableAfterSample: + frameResult.frame.frameStartSample + fftSize, + updateHopSamples: hopSize, + testedBinCount: frameResult.parameters.testedBinCount, + untestedBinCount: frameResult.parameters.untestedBinCount, + detectedBinCount: frameResult.parameters.detectedBinCount, + bins: frameResult.bins, + }); +} + +function sharedDetectorModel(frameResult, fftSize, hopSize, frameCount) { + const { + testedBinCount, + untestedBinCount, + detectedBinCount, + ...sharedParameters + } = frameResult.parameters; + void testedBinCount; + void untestedBinCount; + void detectedBinCount; + return Object.freeze({ + detector: frameResult.detector, + spectrum: Object.freeze({ + ...frameResult.spectrum, + statistic: + "symmetric-hann-windowed-dft-squared-magnitude-per-one-sided-bin", + powerNormalization: "none", + }), + parameters: Object.freeze(sharedParameters), + calibration: frameResult.calibration, + processing: frameResult.processing, + timeAnalysis: Object.freeze({ + frameCount, + fftSize, + hopSize, + completeFramesOnly: true, + partialFramePadding: false, + temporalAveraging: false, + temporalDecisionFusion: false, + overlappingFrames: + hopSize < fftSize, + frameIndependenceClaimed: false, + }), + }); +} + +/** + * Run standard CA-CFAR or OS-CFAR independently on every complete STFT frame + * from one explicitly selected physical channel. + * + * Every tested and untested ADC bin is retained. This function does not merge + * neighboring bins, select peaks, infer analog alias order, combine frames, + * combine microphones, or feed a localization plan. + */ +export async function analyzeCfarActivity(options) { + if (!options || typeof options !== "object") { + throw new TypeError("activity options must be an object."); + } + const samples = requireSingleChannel(options.samples); + const sampleRate = requirePositiveFinite(options.sampleRate, "sampleRate"); + const fftSize = requirePositiveSafeInteger(options.fftSize, "fftSize"); + const hopSize = requirePositiveSafeInteger(options.hopSize, "hopSize"); + if (typeof options.channelId !== "string" || options.channelId.length === 0) { + throw new TypeError("channelId must be a non-empty string."); + } + if (options.window !== "hann") { + throw new RangeError('window must be "hann".'); + } + const transformBatchSize = requirePositiveSafeInteger( + options.transformBatchSize ?? 16, + "transformBatchSize", + ); + const maximumFrameCount = requireLimit( + options.maximumFrameCount, + DEFAULT_MAXIMUM_ACTIVITY_FRAME_COUNT, + "maximumFrameCount", + ); + const maximumBinDecisionCount = requireLimit( + options.maximumBinDecisionCount, + DEFAULT_MAXIMUM_ACTIVITY_BIN_DECISIONS, + "maximumBinDecisionCount", + ); + const maximumFftButterflies = requireLimit( + options.maximumFftButterflies, + DEFAULT_MAXIMUM_ACTIVITY_FFT_BUTTERFLIES, + "maximumFftButterflies", + ); + const workload = estimateCfarActivityWorkload({ + sampleCount: samples.length, + fftSize, + hopSize, + }); + if (workload.frameCount === 0) { + throw new RangeError( + "samples must contain at least one complete activity-analysis frame.", + ); + } + if (workload.frameCount > maximumFrameCount) { + throw new RangeError( + "activity frame count " + workload.frameCount + + " exceeds maximumFrameCount " + maximumFrameCount + ".", + ); + } + if (workload.binDecisionCount > maximumBinDecisionCount) { + throw new RangeError( + "activity bin decision count " + workload.binDecisionCount + + " exceeds maximumBinDecisionCount " + maximumBinDecisionCount + ".", + ); + } + if (workload.fftButterflyCount > maximumFftButterflies) { + throw new RangeError( + "activity FFT butterfly count " + workload.fftButterflyCount + + " exceeds maximumFftButterflies " + maximumFftButterflies + ".", + ); + } + + const window = hannWindow(fftSize); + const real = new Float64Array(fftSize); + const imaginary = new Float64Array(fftSize); + const powerSpectrum = new Float64Array(workload.binCount); + const frames = new Array(workload.frameCount); + let detectedBinCount = 0; + let detectorModel = null; + throwIfAborted(options.signal); + + for (let frameIndex = 0; frameIndex < workload.frameCount; frameIndex += 1) { + const frameStartSample = frameIndex * hopSize; + imaginary.fill(0); + for (let offset = 0; offset < fftSize; offset += 1) { + real[offset] = samples[frameStartSample + offset] * window[offset]; + } + fftInPlace(real, imaginary); + for (let bin = 0; bin < workload.binCount; bin += 1) { + const power = + real[bin] * real[bin] + imaginary[bin] * imaginary[bin]; + if (!Number.isFinite(power)) { + throw new RangeError( + "STFT power became non-finite at frame " + frameIndex + + ", bin " + bin + ".", + ); + } + powerSpectrum[bin] = power; + } + const frameResult = detectCfarFrame({ + detector: options.detector, + powerSpectrum, + sampleRate, + fftSize, + window: "hann", + frameComplete: true, + channelId: options.channelId, + frameIndex, + frameStartSample, + temporalAverageCount: 1, + microphoneAverageCount: 1, + trainingCellsPerSide: options.trainingCellsPerSide, + guardCellsPerSide: options.guardCellsPerSide, + nominalFalseAlarmProbability: + options.nominalFalseAlarmProbability, + ...(options.detector === "os-cfar" + ? { rank: options.rank } + : {}), + }); + detectorModel ??= sharedDetectorModel( + frameResult, + fftSize, + hopSize, + workload.frameCount, + ); + frames[frameIndex] = compactFrameResult( + frameResult, + fftSize, + hopSize, + ); + detectedBinCount += frameResult.parameters.detectedBinCount; + + if ( + (frameIndex + 1) % transformBatchSize === 0 && + frameIndex + 1 < workload.frameCount + ) { + options.onProgress?.({ + stage: "activity-cfar", + completedFrames: frameIndex + 1, + totalFrames: workload.frameCount, + }); + await cooperativeYield(); + throwIfAborted(options.signal); + } + } + options.onProgress?.({ + stage: "activity-cfar", + completedFrames: workload.frameCount, + totalFrames: workload.frameCount, + }); + throwIfAborted(options.signal); + + return { + kind: "micarrsim.cfar-activity-analysis/1", + channel: Object.freeze({ + mode: "one-explicit-physical-channel", + channelId: options.channelId, + microphoneFusion: "none", + }), + input: Object.freeze({ + sampleRate, + sampleCount: samples.length, + adcObservableMinimumHz: 0, + adcObservableMaximumHz: sampleRate / 2, + analogFrequencyOrAliasOrderInferred: false, + sceneTruthUsed: false, + }), + detectorModel, + workload: Object.freeze({ + ...workload, + limits: Object.freeze({ + maximumFrameCount, + maximumBinDecisionCount, + maximumFftButterflies, + }), + }), + summary: Object.freeze({ + frameCount: workload.frameCount, + binRecordCount: workload.binDecisionCount, + detectedBinCount, + selectionApplied: false, + clusteringApplied: false, + localizationApplied: false, + }), + frames: Object.freeze(frames), + }; +} diff --git a/packages/acoustics-core/cfar.js b/packages/acoustics-core/cfar.js new file mode 100644 index 0000000..3522333 --- /dev/null +++ b/packages/acoustics-core/cfar.js @@ -0,0 +1,622 @@ +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" }); +} diff --git a/packages/acoustics-core/index.js b/packages/acoustics-core/index.js index 126dddf..0151532 100644 --- a/packages/acoustics-core/index.js +++ b/packages/acoustics-core/index.js @@ -1,3 +1,5 @@ +export * from "./activity-monitor.js"; +export * from "./cfar.js"; export * from "./das.js"; export * from "./evaluation.js"; export * from "./environment.js"; @@ -16,6 +18,7 @@ export * from "./covariance.js"; export * from "./propagation.js"; export * from "./realtime-localizer.js"; export * from "./signal.js"; +export * from "./spectral-flux.js"; export * from "./stft.js"; export * from "./trajectory.js"; export * from "./wav.js"; diff --git a/packages/acoustics-core/spectral-flux.js b/packages/acoustics-core/spectral-flux.js new file mode 100644 index 0000000..23e4bba --- /dev/null +++ b/packages/acoustics-core/spectral-flux.js @@ -0,0 +1,419 @@ +import { fftInPlace } from "./fft.js"; + +const DIXON_2006_URL = + "https://www.dafx.de/paper-archive/2006/papers/p_133.pdf"; +const BOECK_KREBS_SCHEDL_2012_URL = + "https://www.cp.jku.at/research/papers/Boeck_etal_ISMIR_2012.pdf"; +const DEFAULT_MAXIMUM_OUTPUT_BYTES = 64 * 1024 * 1024; + +function requireConfigurationObject(value, label) { + if (value === null || typeof value !== "object" || Array.isArray(value)) { + throw new TypeError(`${label} must be an object.`); + } + return value; +} + +function requireOwn(configuration, key) { + if (!Object.prototype.hasOwnProperty.call(configuration, key)) { + throw new TypeError(`${key} must be provided explicitly.`); + } + return configuration[key]; +} + +function requirePositiveFinite(value, label) { + if (!Number.isFinite(value) || value <= 0) { + throw new RangeError(`${label} must be finite and greater than zero.`); + } + return value; +} + +function requireNonNegativeFinite(value, label) { + if (!Number.isFinite(value) || value < 0) { + throw new RangeError(`${label} must be finite and non-negative.`); + } + return value; +} + +function requireNonNegativeSafeInteger(value, label) { + if (!Number.isSafeInteger(value) || value < 0) { + throw new RangeError(`${label} must be a non-negative safe integer.`); + } + return value; +} + +function requirePowerOfTwo(value, label) { + if ( + !Number.isSafeInteger(value) || + value < 2 || + value > 0x20000000 || + !Number.isInteger(Math.log2(value)) + ) { + throw new RangeError( + `${label} must be a power-of-two integer from 2 through 536870912.`, + ); + } + return value; +} + +function requireSingleChannelSamples(samples) { + const isTypedArray = ArrayBuffer.isView(samples) && !(samples instanceof DataView); + if (!Array.isArray(samples) && !isTypedArray) { + throw new TypeError("samples must be one numeric channel, not a channel collection."); + } + for (let sampleIndex = 0; sampleIndex < samples.length; sampleIndex += 1) { + if (typeof samples[sampleIndex] !== "number" || !Number.isFinite(samples[sampleIndex])) { + throw new TypeError(`samples[${sampleIndex}] must be finite.`); + } + } + return samples; +} + +function symmetricHammingWindow(size) { + const window = new Float64Array(size); + for (let index = 0; index < size; index += 1) { + window[index] = 0.54 - 0.46 * Math.cos((2 * Math.PI * index) / (size - 1)); + } + return window; +} + +function checkedElementCount(frameCount, binCount) { + const elementCount = frameCount * binCount; + if (!Number.isSafeInteger(elementCount)) { + throw new RangeError("spectral-flux output element count exceeds safe integer range."); + } + return elementCount; +} + +function checkedByteCount(elementCount, bytesPerElement, label) { + const byteCount = elementCount * bytesPerElement; + if (!Number.isSafeInteger(byteCount)) { + throw new RangeError(`${label} byte count exceeds safe integer range.`); + } + return byteCount; +} + +/** + * Compute Simon Dixon's 2006 magnitude spectral-flux onset detection function. + * + * This function accepts exactly one numeric channel. It uses a symmetric + * Hamming window, complete non-padded frames, linear DFT magnitudes, positive + * (half-wave-rectified) frame differences, and an L1 sum over the full DFT. + * No per-frame or global normalization, microphone fusion, filterbank, + * maximum filter, location inference, or band inference is performed. + */ +export function computeDixonSpectralFlux(options) { + const configuration = requireConfigurationObject(options, "options"); + const samples = requireSingleChannelSamples(requireOwn(configuration, "samples")); + const sampleRate = requirePositiveFinite( + requireOwn(configuration, "sampleRate"), + "sampleRate", + ); + const fftSize = requirePowerOfTwo( + requireOwn(configuration, "fftSize"), + "fftSize", + ); + const hopSize = requireNonNegativeSafeInteger( + requireOwn(configuration, "hopSize"), + "hopSize", + ); + if (hopSize < 1 || hopSize > fftSize) { + throw new RangeError("hopSize must be between 1 and fftSize, inclusive."); + } + const channelId = requireOwn(configuration, "channelId"); + if (typeof channelId !== "string" || channelId.length === 0) { + throw new TypeError("channelId must be a non-empty string."); + } + const maximumOutputBytes = configuration.maximumOutputBytes ?? DEFAULT_MAXIMUM_OUTPUT_BYTES; + if (!Number.isSafeInteger(maximumOutputBytes) || maximumOutputBytes <= 0) { + throw new RangeError("maximumOutputBytes must be a positive safe integer."); + } + + const frameCount = samples.length < fftSize + ? 0 + : 1 + Math.floor((samples.length - fftSize) / hopSize); + const binCount = fftSize; + const positiveDifferenceElements = checkedElementCount(frameCount, binCount); + const positiveDifferenceBytes = checkedByteCount( + positiveDifferenceElements, + Float64Array.BYTES_PER_ELEMENT, + "positiveDifferences", + ); + const rawOdfBytes = checkedByteCount( + frameCount, + Float64Array.BYTES_PER_ELEMENT, + "rawOdf", + ); + const frameTimingBytes = checkedByteCount( + frameCount * 3, + Float64Array.BYTES_PER_ELEMENT, + "frame timing", + ); + const outputBytes = positiveDifferenceBytes + rawOdfBytes + frameTimingBytes; + if (!Number.isSafeInteger(outputBytes)) { + throw new RangeError("spectral-flux output byte count exceeds safe integer range."); + } + if (outputBytes > maximumOutputBytes) { + throw new RangeError( + `spectral flux requires ${outputBytes.toLocaleString()} output bytes; ` + + `maximumOutputBytes is ${maximumOutputBytes.toLocaleString()}.`, + ); + } + const workspaceElements = frameCount === 0 ? 0 : fftSize * 4; + const workspaceBytes = checkedByteCount( + workspaceElements, + Float64Array.BYTES_PER_ELEMENT, + "spectral-flux workspace", + ); + const estimatedPeakTypedArrayBytes = outputBytes + workspaceBytes; + if (!Number.isSafeInteger(estimatedPeakTypedArrayBytes)) { + throw new RangeError("spectral-flux peak typed-array byte count exceeds safe integer range."); + } + + const positiveDifferences = new Float64Array(positiveDifferenceElements); + const rawOdf = new Float64Array(frameCount); + const frameStartSamples = new Float64Array(frameCount); + const frameCenterSamples = new Float64Array(frameCount); + const frameAvailableAfterSamples = new Float64Array(frameCount); + const workspaceLength = frameCount === 0 ? 0 : fftSize; + const window = symmetricHammingWindow(workspaceLength); + const real = new Float64Array(workspaceLength); + const imaginary = new Float64Array(workspaceLength); + const previousMagnitudes = new Float64Array(workspaceLength); + + for (let frameIndex = 0; frameIndex < frameCount; frameIndex += 1) { + const frameStart = frameIndex * hopSize; + frameStartSamples[frameIndex] = frameStart; + frameCenterSamples[frameIndex] = frameStart + (fftSize - 1) / 2; + frameAvailableAfterSamples[frameIndex] = frameStart + fftSize; + imaginary.fill(0); + for (let offset = 0; offset < fftSize; offset += 1) { + real[offset] = samples[frameStart + offset] * window[offset]; + } + fftInPlace(real, imaginary); + + let frameFlux = 0; + const differenceOffset = frameIndex * binCount; + for (let binIndex = 0; binIndex < binCount; binIndex += 1) { + const magnitude = Math.hypot(real[binIndex], imaginary[binIndex]); + if (!Number.isFinite(magnitude)) { + throw new RangeError( + `STFT magnitude became non-finite at frame ${frameIndex}, bin ${binIndex}.`, + ); + } + const difference = frameIndex === 0 + ? 0 + : Math.max(0, magnitude - previousMagnitudes[binIndex]); + positiveDifferences[differenceOffset + binIndex] = difference; + frameFlux += difference; + previousMagnitudes[binIndex] = magnitude; + } + if (!Number.isFinite(frameFlux)) { + throw new RangeError(`spectral flux became non-finite at frame ${frameIndex}.`); + } + rawOdf[frameIndex] = frameFlux; + } + + return { + rawOdf, + positiveDifferences, + positiveDifferenceLayout: { + order: "frame-major", + frameStride: binCount, + binCount, + firstFrameHasNoPredecessor: true, + }, + frameCount, + binCount, + sampleRate, + sampleCount: samples.length, + fftSize, + hopSize, + frameStartSamples, + frameCenterSamples, + frameAvailableAfterSamples, + workload: { + completeFrameCount: frameCount, + transformCount: frameCount, + binCount, + comparedBinCount: Math.max(0, frameCount - 1) * binCount, + positiveDifferenceElements, + positiveDifferenceBytes, + rawOdfBytes, + frameTimingBytes, + outputBytes, + maximumOutputBytes, + workspaceBytes, + estimatedPeakTypedArrayBytes, + }, + latency: { + completeFramesOnly: true, + analysisWindowSamples: fftSize, + analysisWindowSeconds: fftSize / sampleRate, + updateHopSamples: hopSize, + updateHopSeconds: hopSize / sampleRate, + firstSpectrumAvailableAfterSamples: frameCount > 0 ? fftSize : null, + firstDifferenceAvailableAfterSamples: frameCount > 1 ? fftSize + hopSize : null, + timestampConvention: + "frameStartSamples/frameCenterSamples are descriptive; frameAvailableAfterSamples is causal availability", + }, + provenance: { + algorithm: "dixon-2006-linear-magnitude-half-wave-l1-spectral-flux", + reference: { + title: "Onset Detection Revisited", + author: "Simon Dixon", + year: 2006, + url: DIXON_2006_URL, + }, + input: { + channelMode: "one-explicit-channel", + channelId, + microphoneFusion: "none", + }, + stft: { + window: "symmetric-hamming", + windowFormula: "0.54 - 0.46*cos(2*pi*n/(N-1))", + framePolicy: "complete-only-no-padding", + spectrumBins: "full-dft-bin-order-including-the-real-signal-negative-frequency-mirror", + }, + magnitudeScale: "linear", + difference: "positive-only-half-wave-rectified", + frequencyAggregation: "unweighted-l1-sum-over-full-dft", + spectrumNormalization: "none", + globalNormalization: "none", + }, + measurementSemantics: { + kind: "raw-onset-novelty", + sustainedActivityMeasure: false, + initialFrameValue: "zero-because-no-preceding-spectrum-exists", + }, + limitations: { + steadySound: + "After the initial spectral rise, an exactly steady spectrum produces zero flux; finite precision and noise can leave small positive values.", + excludedInferences: [ + "event probability", + "sustained activity duration", + "source count", + "source location", + "automatically selected frequency band", + ], + }, + }; +} + +/** + * Apply the causal online peak-picking conditions evaluated by Böck, Krebs, + * and Schedl (2012). All operating parameters are required from the caller. + * The local-maximum and adaptive-mean windows contain only the current and + * past frames; no future frames or global normalization are used. + */ +export function pickCausalSpectralFluxPeaks(rawOdfInput, options) { + const rawOdfSource = requireSingleChannelSamples(rawOdfInput); + const configuration = requireConfigurationObject(options, "options"); + const delta = requireNonNegativeFinite(requireOwn(configuration, "delta"), "delta"); + const preMaxFrames = requireNonNegativeSafeInteger( + requireOwn(configuration, "preMaxFrames"), + "preMaxFrames", + ); + const preMeanFrames = requireNonNegativeSafeInteger( + requireOwn(configuration, "preMeanFrames"), + "preMeanFrames", + ); + const refractoryFrames = requireNonNegativeSafeInteger( + requireOwn(configuration, "refractoryFrames"), + "refractoryFrames", + ); + + const rawOdf = Float64Array.from(rawOdfSource); + for (let frameIndex = 0; frameIndex < rawOdf.length; frameIndex += 1) { + if (rawOdf[frameIndex] < 0) { + throw new RangeError(`rawOdfInput[${frameIndex}] must be non-negative.`); + } + } + const adaptiveThresholds = new Float64Array(rawOdf.length); + const peakMask = new Uint8Array(rawOdf.length); + const peakFrameIndices = []; + const peaks = []; + let lastPeakFrame = null; + let runningMeanWindowSum = 0; + + for (let frameIndex = 0; frameIndex < rawOdf.length; frameIndex += 1) { + const current = rawOdf[frameIndex]; + runningMeanWindowSum += current; + const expiredMeanIndex = frameIndex - preMeanFrames - 1; + if (expiredMeanIndex >= 0) runningMeanWindowSum -= rawOdf[expiredMeanIndex]; + const meanWindowStart = Math.max(0, frameIndex - preMeanFrames); + const meanWindowLength = frameIndex - meanWindowStart + 1; + const localMean = runningMeanWindowSum / meanWindowLength; + const adaptiveThreshold = localMean + delta; + if (!Number.isFinite(adaptiveThreshold)) { + throw new RangeError(`adaptive threshold became non-finite at frame ${frameIndex}.`); + } + adaptiveThresholds[frameIndex] = adaptiveThreshold; + + const maxWindowStart = Math.max(0, frameIndex - preMaxFrames); + let isPastLocalMaximum = true; + for ( + let comparisonIndex = maxWindowStart; + comparisonIndex < frameIndex; + comparisonIndex += 1 + ) { + if (rawOdf[comparisonIndex] > current) { + isPastLocalMaximum = false; + break; + } + } + const clearsAdaptiveThreshold = current >= adaptiveThreshold; + const clearsRefractory = + lastPeakFrame === null || frameIndex - lastPeakFrame > refractoryFrames; + if (isPastLocalMaximum && clearsAdaptiveThreshold && clearsRefractory) { + peakMask[frameIndex] = 1; + peakFrameIndices.push(frameIndex); + peaks.push({ + frameIndex, + rawOdf: current, + adaptiveThreshold, + currentAndPastLocalMean: localMean, + }); + lastPeakFrame = frameIndex; + } + } + + return { + rawOdf, + adaptiveThresholds, + peakMask, + peakFrameIndices: Uint32Array.from(peakFrameIndices), + peaks, + configuration: { + delta, + preMaxFrames, + preMeanFrames, + refractoryFrames, + postMaxFrames: 0, + postMeanFrames: 0, + }, + latency: { + causal: true, + futureContextFrames: 0, + addedLookaheadFrames: 0, + }, + provenance: { + algorithm: "boeck-krebs-schedl-2012-causal-online-peak-picker", + reference: { + title: "Evaluating the Online Capabilities of Onset Detection Methods", + authors: ["Sebastian Böck", "Florian Krebs", "Markus Schedl"], + year: 2012, + url: BOECK_KREBS_SCHEDL_2012_URL, + }, + localMaximum: "current-equals-maximum-of-current-and-available-past-window", + adaptiveThreshold: "mean-of-current-and-available-past-window-plus-delta", + refractoryCondition: "currentFrame-minus-lastPeakFrame-strictly-greater-than-refractoryFrames", + edgePolicy: "truncate-windows-to-available-past", + globalNormalization: "none", + }, + measurementSemantics: { + kind: "thresholded-onset-events-from-raw-novelty", + deltaMeaning: "caller-selected-detection-operating-offset", + }, + }; +} diff --git a/tests/test-runner.js b/tests/test-runner.js index 3787330..8cde63a 100644 --- a/tests/test-runner.js +++ b/tests/test-runner.js @@ -11,6 +11,9 @@ globalThis.assertEqual = (actual, expected, message = "values differ") => { try { await import("/tests/unit/contracts.test.js"); await import("/tests/unit/acoustics.test.js"); + await import("/tests/unit/cfar.test.js"); + await import("/tests/unit/activity-monitor.test.js"); + await import("/tests/unit/spectral-flux.test.js"); await import("/tests/unit/gcc-phat.test.js"); await import("/tests/unit/music-numerics.test.js"); await import("/tests/unit/localization.test.js"); diff --git a/tests/unit/activity-monitor.test.js b/tests/unit/activity-monitor.test.js new file mode 100644 index 0000000..d377647 --- /dev/null +++ b/tests/unit/activity-monitor.test.js @@ -0,0 +1,179 @@ +import { + analyzeCfarActivity, + estimateCfarActivityWorkload, +} from "/packages/acoustics-core/activity-monitor.js"; + +function sinusoid(sampleCount, sampleRate, frequencyHz, amplitude = 1) { + return Float64Array.from( + { length: sampleCount }, + (_, index) => + amplitude * Math.sin((2 * Math.PI * frequencyHz * index) / sampleRate), + ); +} + +function detectorOptions(overrides = {}) { + return { + samples: sinusoid(128, 6_400, 1_200), + sampleRate: 6_400, + fftSize: 64, + hopSize: 64, + window: "hann", + channelId: "mic-reference", + detector: "os-cfar", + trainingCellsPerSide: 4, + guardCellsPerSide: 1, + rank: 6, + nominalFalseAlarmProbability: 0.01, + ...overrides, + }; +} + +test("CFAR activity workload counts complete frames and every one-sided bin", () => { + const workload = estimateCfarActivityWorkload({ + sampleCount: 160, + fftSize: 64, + hopSize: 32, + }); + assertEqual(workload.frameCount, 4, "complete frame count"); + assertEqual(workload.binCount, 33, "one-sided bin count"); + assertEqual(workload.binDecisionCount, 132, "all frame/bin records"); + assertEqual(workload.fftButterflyCount, 768, "radix-2 butterfly count"); +}); + +test("CFAR activity monitor preserves raw bins without time, microphone, or band fusion", async () => { + const result = await analyzeCfarActivity(detectorOptions()); + assertEqual(result.kind, "micarrsim.cfar-activity-analysis/1", "result kind"); + assertEqual(result.channel.channelId, "mic-reference", "explicit channel"); + assertEqual(result.channel.microphoneFusion, "none", "no microphone fusion"); + assertEqual(result.frames.length, 2, "every complete frame retained"); + assertEqual(result.frames[0].bins.length, 33, "every one-sided bin retained"); + assertEqual(result.frames[0].bins[0].status, "untested", "DC remains untested"); + assertEqual( + result.frames[0].bins[32].status, + "untested", + "Nyquist remains untested", + ); + assertEqual(result.summary.selectionApplied, false, "no selection"); + assertEqual(result.summary.clusteringApplied, false, "no clustering"); + assertEqual(result.summary.localizationApplied, false, "no localization"); + assertEqual( + result.detectorModel.processing.bandFusion, + false, + "no band fusion", + ); + assert( + result.frames.some((frame) => frame.bins[12].status === "detected"), + "the exact 1.2 kHz ADC bin is reported as detected", + ); +}); + +test("above-Nyquist analog tone is reported only at its observable ADC alias", async () => { + const sampleRate = 6_400; + const analogFrequencyHz = 7_600; + const result = await analyzeCfarActivity(detectorOptions({ + samples: sinusoid(128, sampleRate, analogFrequencyHz), + })); + assert( + result.frames.some((frame) => frame.bins[12].status === "detected"), + "7.6 kHz samples expose the 1.2 kHz ADC bin", + ); + assertEqual( + result.input.analogFrequencyOrAliasOrderInferred, + false, + "monitor does not infer analog alias order", + ); + assertEqual(result.input.sceneTruthUsed, false, "Scene truth is not an input"); + assertEqual( + Object.prototype.hasOwnProperty.call(result.input, "analogFrequencyHz"), + false, + "the result does not invent an analog carrier field", + ); + assertEqual( + Object.prototype.hasOwnProperty.call(result.frames[0].bins[12], "analogFrequencyHz"), + false, + "the detected ADC bin contains no hidden analog carrier", + ); +}); + +test("CA and OS activity analyses stay separate algorithms", async () => { + const ca = await analyzeCfarActivity(detectorOptions({ + detector: "ca-cfar", + rank: undefined, + })); + const os = await analyzeCfarActivity(detectorOptions()); + assertEqual(ca.detectorModel.detector, "ca-cfar", "CA identity"); + assertEqual(os.detectorModel.detector, "os-cfar", "OS identity"); + assertEqual( + ca.detectorModel.processing.detectorFusion, + false, + "CA has no detector fusion", + ); + assertEqual( + os.detectorModel.processing.detectorFusion, + false, + "OS has no detector fusion", + ); +}); + +test("CFAR activity monitor is deterministic and preserves raw powers", async () => { + const first = await analyzeCfarActivity(detectorOptions()); + const second = await analyzeCfarActivity(detectorOptions()); + for (let frameIndex = 0; frameIndex < first.frames.length; frameIndex += 1) { + for (let bin = 0; bin < first.frames[frameIndex].bins.length; bin += 1) { + const left = first.frames[frameIndex].bins[bin]; + const right = second.frames[frameIndex].bins[bin]; + assertEqual(left.status, right.status, "deterministic status"); + assertEqual(left.cutPower, right.cutPower, "deterministic raw CUT power"); + assertEqual( + left.referencePower, + right.referencePower, + "deterministic reference power", + ); + assertEqual( + left.thresholdPower, + right.thresholdPower, + "deterministic threshold", + ); + } + } +}); + +test("CFAR activity monitor rejects implicit channel fusion and hard-cap overflow", async () => { + let channelError = null; + try { + await analyzeCfarActivity(detectorOptions({ + samples: [ + sinusoid(64, 6_400, 1_200), + sinusoid(64, 6_400, 1_200), + ], + })); + } catch (error) { + channelError = error; + } + assert(channelError instanceof Error, "nested channel input is rejected"); + + let capError = null; + try { + await analyzeCfarActivity(detectorOptions({ + maximumBinDecisionCount: 65, + })); + } catch (error) { + capError = error; + } + assert( + capError?.message.includes("maximumBinDecisionCount"), + "bin decision cap rejects work instead of truncating results", + ); +}); + +test("CFAR activity monitor honors an already-aborted signal", async () => { + const controller = new AbortController(); + controller.abort(); + let abortError = null; + try { + await analyzeCfarActivity(detectorOptions({ signal: controller.signal })); + } catch (error) { + abortError = error; + } + assertEqual(abortError?.name, "AbortError", "typed abort"); +}); diff --git a/tests/unit/cfar.test.js b/tests/unit/cfar.test.js new file mode 100644 index 0000000..9071d0c --- /dev/null +++ b/tests/unit/cfar.test.js @@ -0,0 +1,465 @@ +import { + computeCaCfarThresholdMultiplier, + computeOsCfarLogFalseAlarmProbability, + computeOsCfarThresholdMultiplier, + detectCaCfarFrame, + detectCfarFrame, + detectOsCfarFrame, +} from "/packages/acoustics-core/cfar.js"; +import { fftInPlace } from "/packages/acoustics-core/fft.js"; +import { hannWindow } from "/packages/acoustics-core/stft.js"; + +function assertNear(actual, expected, tolerance, label) { + assert( + Math.abs(actual - expected) <= tolerance, + `${label}: expected ${expected} ± ${tolerance}, received ${actual}`, + ); +} + +function assertRelativeNear(actual, expected, relativeTolerance, label) { + const scale = Math.max(Number.MIN_VALUE, Math.abs(expected)); + assertNear(actual, expected, relativeTolerance * scale, label); +} + +function captureError(callback) { + try { + callback(); + } catch (error) { + return error; + } + return null; +} + +function assertError(callback, ErrorType, messageFragment, label) { + const error = captureError(callback); + assert(error instanceof ErrorType, `${label}: expected ${ErrorType.name}`); + assert( + error.message.includes(messageFragment), + `${label}: expected message containing ${messageFragment}, received ${error.message}`, + ); +} + +function createPowerFixture(fftSize, fill = 1) { + return new Float64Array(fftSize / 2 + 1).fill(fill); +} + +function createBaseOptions(overrides = {}) { + const fftSize = overrides.fftSize ?? 64; + return { + powerSpectrum: overrides.powerSpectrum ?? createPowerFixture(fftSize), + detector: "ca-cfar", + sampleRate: 8192, + fftSize, + window: "hann", + frameComplete: true, + channelId: "mic-0", + frameIndex: 0, + frameStartSample: 0, + temporalAverageCount: 1, + microphoneAverageCount: 1, + trainingCellsPerSide: 4, + guardCellsPerSide: 2, + nominalFalseAlarmProbability: 1e-3, + ...overrides, + }; +} + +function createHannTonePower({ + sampleRate, + fftSize, + frequencyHz, + amplitude = 1, +}) { + const window = hannWindow(fftSize); + const real = new Float64Array(fftSize); + const imaginary = new Float64Array(fftSize); + for (let sample = 0; sample < fftSize; sample += 1) { + real[sample] = + amplitude * + Math.sin((2 * Math.PI * frequencyHz * sample) / sampleRate) * + window[sample]; + } + fftInPlace(real, imaginary); + return Float64Array.from( + { length: fftSize / 2 + 1 }, + (_, bin) => real[bin] * real[bin] + imaginary[bin] * imaginary[bin], + ); +} + +test("CA-CFAR mean-reference multiplier recovers the requested nominal per-CUT Pfa", () => { + const trainingCellCount = 16; + const nominalFalseAlarmProbability = 1e-3; + const multiplier = computeCaCfarThresholdMultiplier({ + trainingCellCount, + nominalFalseAlarmProbability, + }); + const recoveredProbability = + (1 + multiplier / trainingCellCount) ** -trainingCellCount; + assertRelativeNear( + recoveredProbability, + nominalFalseAlarmProbability, + 2e-14, + "CA nominal false-alarm probability", + ); + + const nearOneProbability = 1 - 1e-12; + const nearOneMultiplier = computeCaCfarThresholdMultiplier({ + trainingCellCount: 1024, + nominalFalseAlarmProbability: nearOneProbability, + }); + assert( + nearOneMultiplier > 0 && Number.isFinite(nearOneMultiplier), + "expm1 form preserves a finite positive multiplier near Pfa=1", + ); +}); + +test("OS-CFAR log-product bisection recovers Pfa and the analytic rank-one case", () => { + const trainingCellCount = 16; + const rank = 12; + const nominalFalseAlarmProbability = 1e-4; + const multiplier = computeOsCfarThresholdMultiplier({ + trainingCellCount, + rank, + nominalFalseAlarmProbability, + }); + const recoveredLogProbability = computeOsCfarLogFalseAlarmProbability({ + trainingCellCount, + rank, + thresholdMultiplier: multiplier, + }); + assertNear( + recoveredLogProbability, + Math.log(nominalFalseAlarmProbability), + 2e-12, + "OS log false-alarm probability", + ); + + const rankOneProbability = 0.25; + const rankOneMultiplier = computeOsCfarThresholdMultiplier({ + trainingCellCount, + rank: 1, + nominalFalseAlarmProbability: rankOneProbability, + }); + assertNear( + rankOneMultiplier, + trainingCellCount * (1 / rankOneProbability - 1), + 1e-11, + "rank-one OS multiplier", + ); + assertError( + () => computeOsCfarThresholdMultiplier({ + trainingCellCount, + rank, + nominalFalseAlarmProbability, + maximumIterations: 1, + }), + RangeError, + "did not converge", + "explicit OS bisection iteration cap", + ); +}); + +test("CA-CFAR preserves every Hann-frame bin and detects an exact-bin tone", () => { + const sampleRate = 8192; + const fftSize = 256; + const toneBin = 43; + const powerSpectrum = createHannTonePower({ + sampleRate, + fftSize, + frequencyHz: toneBin * sampleRate / fftSize, + }); + const result = detectCaCfarFrame(createBaseOptions({ + powerSpectrum, + sampleRate, + fftSize, + trainingCellsPerSide: 12, + guardCellsPerSide: 3, + })); + const tone = result.bins[toneBin]; + + assertEqual(result.bins.length, fftSize / 2 + 1, "all one-sided bins retained"); + assertEqual(tone.status, "detected", "exact-bin Hann tone detected"); + assertEqual(tone.cutPower, powerSpectrum[toneBin], "raw CUT power retained"); + assert(tone.referencePower >= 0, "raw reference statistic retained"); + assert(tone.thresholdPower >= 0, "raw threshold retained"); + assert(tone.thresholdRatio > 1, "raw threshold ratio exceeds one"); + assertEqual( + result.parameters.referenceStatistic, + "arithmetic-mean-power", + "CA mean normalization is explicit", + ); + assertEqual( + result.calibration.calibratedFalseAlarmProbability, + null, + "nominal Pfa is not mislabeled as empirical calibration", + ); + assertEqual(result.processing.clustering, false, "no bin clustering"); + assertEqual(result.processing.temporalPersistence, false, "no persistence"); + assertEqual(result.processing.detectorFusion, false, "no CA/OS fusion"); +}); + +test("full symmetric references exclude DC, Nyquist, and every insufficient edge bin", () => { + const fftSize = 32; + const result = detectCaCfarFrame(createBaseOptions({ + fftSize, + powerSpectrum: createPowerFixture(fftSize), + trainingCellsPerSide: 2, + guardCellsPerSide: 1, + nominalFalseAlarmProbability: 0.1, + })); + + assertEqual(result.bins[0].status, "untested", "DC status"); + assertEqual(result.bins[0].reason, "real-only-dc-bin", "DC reason"); + assertEqual(result.bins[16].status, "untested", "Nyquist status"); + assertEqual( + result.bins[16].reason, + "real-only-nyquist-bin", + "Nyquist reason", + ); + for (const bin of [1, 2, 3, 13, 14, 15]) { + assertEqual(result.bins[bin].status, "untested", `edge bin ${bin}`); + assertEqual( + result.bins[bin].reason, + "insufficient-full-symmetric-reference-window", + `edge reason ${bin}`, + ); + assertEqual(result.bins[bin].cutPower, 1, `edge CUT power ${bin}`); + } + for (let bin = 4; bin <= 12; bin += 1) { + assertEqual(result.bins[bin].status, "not-detected", `tested bin ${bin}`); + assertEqual(result.bins[bin].referencePower, 1, `fixed-N mean ${bin}`); + } + assertEqual(result.parameters.trainingCellCount, 4, "N=2R"); + assertEqual(result.parameters.firstTestedBin, 4, "first fully supported bin"); + assertEqual(result.parameters.lastTestedBin, 12, "last fully supported bin"); + assertEqual(result.parameters.testedBinCount, 9, "tested bin count"); +}); + +test("CA uses a strict threshold comparison and does not promote a tie", () => { + const fftSize = 8; + const nominalFalseAlarmProbability = (1 + 1 / 2) ** -2; + const result = detectCaCfarFrame(createBaseOptions({ + fftSize, + powerSpectrum: createPowerFixture(fftSize), + trainingCellsPerSide: 1, + guardCellsPerSide: 0, + nominalFalseAlarmProbability, + })); + assertNear(result.parameters.thresholdMultiplier, 1, 1e-15, "unit alpha"); + assertEqual(result.bins[2].cutPower, 1, "tie CUT"); + assertNear(result.bins[2].thresholdPower, 1, 1e-15, "tie threshold"); + assertEqual(result.bins[2].status, "not-detected", "strict greater-than rule"); +}); + +test("zero thresholds preserve infinite-ratio meaning without non-JSON Infinity", () => { + const fftSize = 16; + const powerSpectrum = createPowerFixture(fftSize, 0); + powerSpectrum[4] = 1; + const result = detectCaCfarFrame(createBaseOptions({ + fftSize, + powerSpectrum, + trainingCellsPerSide: 1, + guardCellsPerSide: 0, + nominalFalseAlarmProbability: 0.1, + })); + const detected = result.bins[4]; + assertEqual(detected.status, "detected", "positive CUT exceeds zero threshold"); + assertEqual(detected.thresholdPower, 0, "zero threshold retained"); + assertEqual(detected.thresholdRatio, null, "non-finite ratio is JSON-safe null"); + assertEqual( + detected.thresholdRatioStatus, + "positive-infinity-zero-threshold", + "positive infinite limit remains explicit", + ); + assert( + JSON.stringify(result).includes('"positive-infinity-zero-threshold"'), + "JSON serialization preserves the ratio meaning", + ); +}); + +test("OS-CFAR remains a separate detector and resists one high training contaminant", () => { + const fftSize = 128; + const cutBin = 32; + const powerSpectrum = createPowerFixture(fftSize); + powerSpectrum[cutBin] = 100; + powerSpectrum[23] = 1000; + const common = createBaseOptions({ + fftSize, + powerSpectrum, + trainingCellsPerSide: 8, + guardCellsPerSide: 1, + nominalFalseAlarmProbability: 0.1, + }); + const ca = detectCaCfarFrame(common); + const os = detectOsCfarFrame({ ...common, rank: 12 }); + + assertEqual(ca.bins[cutBin].status, "not-detected", "CA masking exposed"); + assertEqual(os.bins[cutBin].status, "detected", "OS order statistic result"); + assertEqual(os.parameters.rank, 12, "explicit OS rank"); + assertEqual(os.parameters.rankConvention, "1-based-ascending", "rank convention"); + assertEqual(os.parameters.referencePower, undefined, "no invented fused reference"); + assertEqual(os.processing.detectorFusion, false, "detectors remain separate"); +}); + +test("a sampled above-Nyquist tone returns only its observed ADC bin", () => { + const sampleRate = 8192; + const fftSize = 256; + const analogFrequencyHz = 6144; + const expectedAdcBin = 64; + const expectedAdcFrequencyHz = 2048; + const powerSpectrum = createHannTonePower({ + sampleRate, + fftSize, + frequencyHz: analogFrequencyHz, + }); + const result = detectOsCfarFrame(createBaseOptions({ + powerSpectrum, + sampleRate, + fftSize, + trainingCellsPerSide: 12, + guardCellsPerSide: 3, + rank: 18, + })); + const observed = result.bins[expectedAdcBin]; + + assertEqual(observed.status, "detected", "sampled alias detected"); + assertNear( + observed.adcCenterFrequencyHz, + expectedAdcFrequencyHz, + 0, + "observed ADC center frequency", + ); + assertEqual( + result.spectrum.analogFrequencyOrAliasOrderInferred, + false, + "no analog alias-order inference", + ); + assertEqual( + Object.prototype.hasOwnProperty.call(observed, "analogFrequencyHz"), + false, + "bin result contains no invented analog frequency", + ); + assertEqual(result.bands, undefined, "no bins-to-band fusion"); +}); + +test("CFAR rejects nonfinite power, invalid configuration, partial frames, and averaging", () => { + const nonfinite = createPowerFixture(64); + nonfinite[10] = Number.NaN; + assertError( + () => detectCaCfarFrame(createBaseOptions({ powerSpectrum: nonfinite })), + TypeError, + "powerSpectrum[10] must be finite", + "NaN power", + ); + + const negative = createPowerFixture(64); + negative[10] = -1; + assertError( + () => detectCaCfarFrame(createBaseOptions({ powerSpectrum: negative })), + RangeError, + "must be non-negative", + "negative power", + ); + assertError( + () => detectCaCfarFrame(createBaseOptions({ + powerSpectrum: new Float64Array(32), + })), + RangeError, + "exactly 33", + "wrong one-sided shape", + ); + for (const probability of [0, 1, Number.NaN]) { + assertError( + () => computeCaCfarThresholdMultiplier({ + trainingCellCount: 8, + nominalFalseAlarmProbability: probability, + }), + RangeError, + "strictly inside 0..1", + `invalid Pfa ${probability}`, + ); + } + assertError( + () => detectOsCfarFrame(createBaseOptions({ rank: 0 })), + RangeError, + "1-based ascending", + "zero OS rank", + ); + assertError( + () => detectOsCfarFrame(createBaseOptions({ rank: 9 })), + RangeError, + "inside 1..8", + "OS rank above N", + ); + assertError( + () => detectCaCfarFrame(createBaseOptions({ trainingCellsPerSide: 0 })), + TypeError, + "trainingCellsPerSide must be a positive safe integer", + "zero training cells", + ); + assertError( + () => detectCaCfarFrame(createBaseOptions({ guardCellsPerSide: -1 })), + TypeError, + "guardCellsPerSide must be a non-negative safe integer", + "negative guard cells", + ); + assertError( + () => detectCaCfarFrame(createBaseOptions({ + trainingCellsPerSide: 16, + guardCellsPerSide: 0, + })), + RangeError, + "leave no interior bin", + "reference window larger than the usable spectrum", + ); + assertError( + () => detectCfarFrame(createBaseOptions({ window: "rectangular" })), + RangeError, + 'window must be "hann"', + "non-Hann window", + ); + assertError( + () => detectCaCfarFrame(createBaseOptions({ frameComplete: false })), + RangeError, + "partial frames are not tested", + "partial frame", + ); + assertError( + () => detectCaCfarFrame(createBaseOptions({ temporalAverageCount: 2 })), + RangeError, + "temporalAverageCount must equal one", + "temporal averaging", + ); + assertError( + () => detectCaCfarFrame(createBaseOptions({ microphoneAverageCount: 2 })), + RangeError, + "microphoneAverageCount must equal one", + "microphone averaging", + ); +}); + +test("OS-CFAR is deterministic and sorting never mutates the caller spectrum", () => { + const fftSize = 128; + const powerSpectrum = Float64Array.from( + { length: fftSize / 2 + 1 }, + (_, bin) => 1 + ((bin * 17) % 23) / 11, + ); + powerSpectrum[28] = 40; + powerSpectrum[41] = 75; + const before = Float64Array.from(powerSpectrum); + const options = createBaseOptions({ + fftSize, + powerSpectrum, + trainingCellsPerSide: 8, + guardCellsPerSide: 2, + rank: 12, + }); + const first = detectOsCfarFrame(options); + const second = detectOsCfarFrame(options); + + assertEqual(JSON.stringify(second), JSON.stringify(first), "repeatable result"); + assertEqual(powerSpectrum.length, before.length, "input length retained"); + for (let bin = 0; bin < before.length; bin += 1) { + assertEqual(powerSpectrum[bin], before[bin], `input power ${bin}`); + } +}); diff --git a/tests/unit/spectral-flux.test.js b/tests/unit/spectral-flux.test.js new file mode 100644 index 0000000..0fade05 --- /dev/null +++ b/tests/unit/spectral-flux.test.js @@ -0,0 +1,372 @@ +import { + computeDixonSpectralFlux, + pickCausalSpectralFluxPeaks, +} from "/packages/acoustics-core/spectral-flux.js"; + +function assertNear(actual, expected, tolerance, label) { + assert( + Math.abs(actual - expected) <= tolerance, + `${label}: expected ${expected} +/- ${tolerance}, received ${actual}`, + ); +} + +function assertArrayNear(actual, expected, tolerance, label) { + assertEqual(actual.length, expected.length, `${label} length`); + for (let index = 0; index < expected.length; index += 1) { + assertNear(actual[index], expected[index], tolerance, `${label}[${index}]`); + } +} + +function assertArrayExact(actual, expected, label) { + assertEqual(actual.length, expected.length, `${label} length`); + for (let index = 0; index < expected.length; index += 1) { + assertEqual(actual[index], expected[index], `${label}[${index}]`); + } +} + +function assertThrows(callback, ErrorType, expectedText, label) { + let caught = null; + try { + callback(); + } catch (error) { + caught = error; + } + assert(caught instanceof ErrorType, `${label}: expected ${ErrorType.name}`); + assert( + caught.message.includes(expectedText), + `${label}: expected message containing "${expectedText}", received "${caught.message}"`, + ); +} + +function fluxConfiguration(samples, overrides = {}) { + return { + samples, + sampleRate: 4, + fftSize: 4, + hopSize: 4, + channelId: "adc-channel-0", + ...overrides, + }; +} + +test("Dixon flux는 Hamming linear magnitude의 양의 차이를 full DFT에서 그대로 합산한다", () => { + const result = computeDixonSpectralFlux(fluxConfiguration( + Float64Array.of(0, 0, 0, 0, 1, 1, 1, 1), + { maximumOutputBytes: 128 }, + )); + const sideMagnitude = Math.hypot(-0.69, -0.69); + const expectedSecondFrame = [1.7, sideMagnitude, 0, sideMagnitude]; + + assertEqual(result.frameCount, 2, "frame count"); + assertEqual(result.binCount, 4, "full DFT bin count"); + assert(result.rawOdf instanceof Float64Array, "raw ODF type"); + assert(result.positiveDifferences instanceof Float64Array, "raw per-bin type"); + assertArrayExact(result.positiveDifferences.subarray(0, 4), [0, 0, 0, 0], "first frame"); + assertArrayNear( + result.positiveDifferences.subarray(4, 8), + expectedSecondFrame, + 1e-12, + "second-frame positive differences", + ); + assertNear( + result.rawOdf[1], + expectedSecondFrame.reduce((sum, value) => sum + value, 0), + 1e-12, + "raw L1 ODF", + ); + assertArrayExact(result.frameStartSamples, [0, 4], "frame starts"); + assertArrayExact(result.frameCenterSamples, [1.5, 5.5], "frame centers"); + assertArrayExact(result.frameAvailableAfterSamples, [4, 8], "causal availability"); + assertEqual(result.latency.firstDifferenceAvailableAfterSamples, 8, "first difference latency"); + assertEqual(result.provenance.input.channelId, "adc-channel-0", "explicit channel provenance"); + assertEqual(result.provenance.input.microphoneFusion, "none", "no microphone fusion"); + assertEqual(result.provenance.globalNormalization, "none", "no global normalization"); + assertEqual(result.provenance.magnitudeScale, "linear", "linear magnitude"); + assertEqual( + result.provenance.stft.spectrumBins, + "full-dft-bin-order-including-the-real-signal-negative-frequency-mirror", + "full DFT aggregation disclosure", + ); + assertEqual(result.workload.positiveDifferenceBytes, 64, "raw bin bytes"); + assertEqual(result.workload.outputBytes, 128, "total returned typed-array bytes"); + assertEqual(result.workload.maximumOutputBytes, 128, "caller output limit"); +}); + +test("Dixon flux는 silence와 incomplete tail을 패딩 없이 처리한다", () => { + const samples = new Float64Array(9); + samples[8] = 1_000_000; + const result = computeDixonSpectralFlux(fluxConfiguration(samples)); + + assertEqual(result.frameCount, 2, "only two complete frames"); + assertArrayExact(result.rawOdf, [0, 0], "silent raw ODF"); + assertArrayExact( + result.positiveDifferences, + new Float64Array(8), + "silent per-bin differences", + ); + assertEqual(result.latency.completeFramesOnly, true, "complete frame policy"); + assertEqual( + result.provenance.stft.framePolicy, + "complete-only-no-padding", + "no padding provenance", + ); + + const tooShort = computeDixonSpectralFlux(fluxConfiguration( + Float64Array.of(1, 2, 3), + { maximumOutputBytes: 1 }, + )); + assertEqual(tooShort.frameCount, 0, "short input has no frame"); + assertEqual(tooShort.workload.outputBytes, 0, "short input allocates no output data"); + assertEqual(tooShort.workload.workspaceBytes, 0, "short input allocates no FFT workspace"); +}); + +test("정상상태 tone은 최초 spectral rise 뒤 raw flux가 0에 가까워지는 한계를 보존한다", () => { + const tone = Float64Array.from( + { length: 8 }, + (_, index) => Math.cos((2 * Math.PI * index) / 8), + ); + const samples = new Float64Array(24); + samples.set(tone, 8); + samples.set(tone, 16); + const result = computeDixonSpectralFlux({ + samples, + sampleRate: 8, + fftSize: 8, + hopSize: 8, + channelId: "tone-channel", + }); + + assertEqual(result.frameCount, 3, "three complete frames"); + assertEqual(result.rawOdf[0], 0, "no predecessor for first frame"); + assert(result.rawOdf[1] > 0, "tone rise is nonzero"); + assertNear(result.rawOdf[2], 0, 1e-12, "unchanged tone flux"); + assertArrayNear( + result.positiveDifferences.subarray(16, 24), + new Float64Array(8), + 1e-12, + "unchanged tone per-bin differences", + ); + assertEqual( + result.measurementSemantics.sustainedActivityMeasure, + false, + "not a sustained activity measure", + ); + assert( + result.limitations.steadySound.includes("steady spectrum produces zero flux"), + "steady-sound limitation is structural metadata", + ); +}); + +test("Boeck causal peak picker는 과거 max/mean, delta와 strict refractory만 사용한다", () => { + const source = Float64Array.of(0, 1, 0, 10, 0); + const configuration = { + delta: 0.4, + preMaxFrames: 1, + preMeanFrames: 2, + refractoryFrames: 1, + }; + const result = pickCausalSpectralFluxPeaks(source, configuration); + + assertArrayExact(result.rawOdf, [0, 1, 0, 10, 0], "preserved raw ODF"); + assertArrayNear( + result.adaptiveThresholds, + [0.4, 0.9, 1 / 3 + 0.4, 11 / 3 + 0.4, 10 / 3 + 0.4], + 1e-12, + "causal thresholds", + ); + assertArrayExact(result.peakFrameIndices, [1, 3], "causal peaks"); + assertArrayExact(result.peakMask, [0, 1, 0, 1, 0], "peak mask"); + assertEqual(result.configuration.postMaxFrames, 0, "no future maximum frames"); + assertEqual(result.configuration.postMeanFrames, 0, "no future mean frames"); + assertEqual(result.latency.futureContextFrames, 0, "no future context"); + assertEqual(result.provenance.globalNormalization, "none", "no peak-picker normalization"); + assertEqual( + result.measurementSemantics.kind, + "thresholded-onset-events-from-raw-novelty", + "event semantics", + ); + + const prefix = pickCausalSpectralFluxPeaks(Float64Array.of(0, 1), configuration); + assertArrayExact(prefix.peakFrameIndices, [1], "prefix peak"); + assertEqual(prefix.adaptiveThresholds[1], result.adaptiveThresholds[1], "future invariance"); + + const refractory = pickCausalSpectralFluxPeaks( + Float64Array.of(0, 1, 0, 2), + { + delta: 0.4, + preMaxFrames: 1, + preMeanFrames: 1, + refractoryFrames: 2, + }, + ); + assertArrayExact(refractory.peakFrameIndices, [1], "strict refractory boundary"); + + source[1] = 99; + assertEqual(result.rawOdf[1], 1, "returned raw ODF is a preserved copy"); +}); + +test("spectral flux와 causal picker는 non-finite 및 음수 novelty를 거부한다", () => { + assertThrows( + () => computeDixonSpectralFlux(fluxConfiguration(Float64Array.of(0, Number.NaN, 0, 0))), + TypeError, + "samples[1] must be finite", + "NaN sample", + ); + assertThrows( + () => computeDixonSpectralFlux(fluxConfiguration( + Float64Array.of(Number.MAX_VALUE, Number.MAX_VALUE, Number.MAX_VALUE, Number.MAX_VALUE), + )), + RangeError, + "STFT magnitude became non-finite", + "overflowing transform", + ); + assertThrows( + () => pickCausalSpectralFluxPeaks( + Float64Array.of(0, Number.POSITIVE_INFINITY), + { delta: 1, preMaxFrames: 1, preMeanFrames: 1, refractoryFrames: 1 }, + ), + TypeError, + "samples[1] must be finite", + "infinite ODF", + ); + assertThrows( + () => pickCausalSpectralFluxPeaks( + Float64Array.of(0, -1), + { delta: 1, preMaxFrames: 1, preMeanFrames: 1, refractoryFrames: 1 }, + ), + RangeError, + "rawOdfInput[1] must be non-negative", + "negative ODF", + ); + assertThrows( + () => pickCausalSpectralFluxPeaks( + Float64Array.of(0, 1), + { delta: Number.NaN, preMaxFrames: 1, preMeanFrames: 1, refractoryFrames: 1 }, + ), + RangeError, + "delta must be finite", + "non-finite delta", + ); +}); + +test("spectral flux 설정은 single channel, FFT, provenance와 출력 상한을 명시적으로 검증한다", () => { + assertThrows( + () => computeDixonSpectralFlux({ + samples: Float64Array.of(0, 0, 0, 0), + sampleRate: 4, + fftSize: 4, + hopSize: 4, + }), + TypeError, + "channelId must be provided explicitly", + "missing channel id", + ); + assertThrows( + () => computeDixonSpectralFlux(fluxConfiguration(Float64Array.of(0, 0, 0, 0), { + channelId: "", + })), + TypeError, + "channelId must be a non-empty string", + "empty channel id", + ); + assertThrows( + () => computeDixonSpectralFlux(fluxConfiguration([[0], [0], [0], [0]])), + TypeError, + "samples[0] must be finite", + "channel collection", + ); + assertThrows( + () => computeDixonSpectralFlux(fluxConfiguration(new Float64Array(8), { fftSize: 3 })), + RangeError, + "power-of-two integer", + "non-radix FFT", + ); + assertThrows( + () => computeDixonSpectralFlux(fluxConfiguration(new Float64Array(8), { + fftSize: 2 ** 32, + })), + RangeError, + "through 536870912", + "FFT exceeds radix implementation limit", + ); + assertThrows( + () => computeDixonSpectralFlux(fluxConfiguration(new Float64Array(8), { hopSize: 0 })), + RangeError, + "hopSize must be between 1 and fftSize", + "zero hop", + ); + assertThrows( + () => computeDixonSpectralFlux(fluxConfiguration(new Float64Array(8), { + maximumOutputBytes: 127, + })), + RangeError, + "maximumOutputBytes is 127", + "pre-allocation output cap", + ); + assertThrows( + () => computeDixonSpectralFlux(fluxConfiguration(new Float64Array(4), { + maximumOutputBytes: 1.5, + })), + RangeError, + "maximumOutputBytes must be a positive safe integer", + "invalid output cap", + ); + assertThrows( + () => pickCausalSpectralFluxPeaks(Float64Array.of(0, 1), { + delta: 0.1, + preMaxFrames: 1, + preMeanFrames: 1, + }), + TypeError, + "refractoryFrames must be provided explicitly", + "missing peak-picker parameter", + ); + assertThrows( + () => pickCausalSpectralFluxPeaks(Float64Array.of(0, 1), { + delta: 0.1, + preMaxFrames: -1, + preMeanFrames: 1, + refractoryFrames: 1, + }), + RangeError, + "preMaxFrames must be a non-negative safe integer", + "negative peak-picker window", + ); +}); + +test("spectral flux와 causal peak picker는 같은 입력에서 결정적이다", () => { + const samples = Float64Array.of( + 0, 0.25, -0.5, 0.75, + 0.5, -0.25, 0.125, -0.625, + -0.5, 0.25, 0.75, -0.125, + ); + const untouched = Float64Array.from(samples); + const options = { + samples, + sampleRate: 8, + fftSize: 4, + hopSize: 2, + channelId: "determinism-channel", + }; + const first = computeDixonSpectralFlux(options); + const second = computeDixonSpectralFlux(options); + + assertArrayExact(samples, untouched, "input samples remain unchanged"); + assertArrayExact(first.rawOdf, second.rawOdf, "deterministic raw ODF"); + assertArrayExact( + first.positiveDifferences, + second.positiveDifferences, + "deterministic per-bin differences", + ); + assertEqual(JSON.stringify(first), JSON.stringify(second), "deterministic metadata"); + + const pickerOptions = { + delta: 0.01, + preMaxFrames: 2, + preMeanFrames: 3, + refractoryFrames: 1, + }; + const firstPeaks = pickCausalSpectralFluxPeaks(first.rawOdf, pickerOptions); + const secondPeaks = pickCausalSpectralFluxPeaks(first.rawOdf, pickerOptions); + assertEqual(JSON.stringify(firstPeaks), JSON.stringify(secondPeaks), "deterministic peaks"); +});