diff --git a/packages/acoustics-core/das.js b/packages/acoustics-core/das.js index 0a63760..7847f93 100644 --- a/packages/acoustics-core/das.js +++ b/packages/acoustics-core/das.js @@ -5,6 +5,11 @@ import { requireMaxCompensationDb, } from "./frequency-response.js"; import { cooperativeYield } from "./cooperative-yield.js"; +import { + DEFAULT_PHAT_RELATIVE_FLOOR, + isDasAlgorithm, + requireLocalizationAlgorithmDefinition, +} from "./localization-algorithms.js"; import { binsForBand, createStftWorkspace, extractStftBins, getFrameCount } from "./stft.js"; function throwIfAborted(signal) { @@ -48,20 +53,24 @@ function createBandIndexLookup(bandBins) { return lookup; } -function createBinCalibration( - arrayDefinition, - bin, - sampleRate, - fftSize, - compensateFrequencyResponse, - maximumCompensationDb, -) { +export function createBinCalibration(options) { + const { + arrayDefinition, + bin, + sampleRate, + fftSize, + compensateFrequencyResponse = true, + maxCompensationDb: maximumCompensationDb, + includeAmplitude = true, + } = options; const microphoneCount = arrayDefinition.microphones.length; const gains = new Float64Array(microphoneCount); const phases = new Float64Array(microphoneCount); const frequencyHz = (bin * sampleRate) / fftSize; const response = {}; - const maxCompensationDb = requireMaxCompensationDb(maximumCompensationDb); + const maxCompensationDb = includeAmplitude + ? requireMaxCompensationDb(maximumCompensationDb) + : 0; for (let microphoneIndex = 0; microphoneIndex < microphoneCount; microphoneIndex += 1) { const microphone = arrayDefinition.microphones[microphoneIndex]; let responseGainDb = 0; @@ -71,23 +80,72 @@ function createBinCalibration( responseGainDb = response.gainDb; responsePhaseDeg = response.phaseDeg; } - const responseCompensationDb = compensateFrequencyResponse + const responseCompensationDb = includeAmplitude && compensateFrequencyResponse ? inverseResponseGainDb(responseGainDb, maxCompensationDb) : 0; - gains[microphoneIndex] = dbToLinear( - -(microphone.sensitivityDb ?? 0) + responseCompensationDb, - ); + gains[microphoneIndex] = includeAmplitude + ? dbToLinear(-(microphone.sensitivityDb ?? 0) + responseCompensationDb) + : 1; phases[microphoneIndex] = ((microphone.phaseOffsetDeg ?? 0) + responsePhaseDeg) * Math.PI / 180; } return { gains, phases }; } +function requirePhatRelativeFloor(value) { + if (!Number.isFinite(value) || value < 0 || value > 1) { + throw new RangeError("phatRelativeFloor must lie inside 0..1."); + } + return value; +} + +export function phaseNormalizeSpectrum( + spectrum, + frameCount, + microphoneCount, + relativeFloor = DEFAULT_PHAT_RELATIVE_FLOOR, +) { + const expectedLength = frameCount * microphoneCount; + if ( + !spectrum || + spectrum.real?.length !== expectedLength || + spectrum.imaginary?.length !== expectedLength + ) { + throw new RangeError("spectrum length must match frameCount and microphoneCount."); + } + requirePhatRelativeFloor(relativeFloor); + const real = new Float32Array(expectedLength); + const imaginary = new Float32Array(expectedLength); + for (let frameIndex = 0; frameIndex < frameCount; frameIndex += 1) { + const frameOffset = frameIndex * microphoneCount; + let maximumMagnitude = 0; + for (let microphoneIndex = 0; microphoneIndex < microphoneCount; microphoneIndex += 1) { + const spectrumIndex = frameOffset + microphoneIndex; + maximumMagnitude = Math.max( + maximumMagnitude, + Math.hypot(spectrum.real[spectrumIndex], spectrum.imaginary[spectrumIndex]), + ); + } + const magnitudeFloor = maximumMagnitude * relativeFloor; + for (let microphoneIndex = 0; microphoneIndex < microphoneCount; microphoneIndex += 1) { + const spectrumIndex = frameOffset + microphoneIndex; + const inputReal = spectrum.real[spectrumIndex]; + const inputImaginary = spectrum.imaginary[spectrumIndex]; + const magnitude = Math.hypot(inputReal, inputImaginary); + if (!(magnitude > magnitudeFloor)) continue; + real[spectrumIndex] = inputReal / magnitude; + imaginary[spectrumIndex] = inputImaginary / magnitude; + } + } + return { real, imaginary }; +} + function fillSteeringVector({ angularWaveNumber, distances, distanceOffset, calibration, + applyCalibrationGains = true, outputReal, outputImaginary, outputOffset, @@ -96,10 +154,9 @@ function fillSteeringVector({ const phase = angularWaveNumber * distances[distanceOffset + microphoneIndex] - calibration.phases[microphoneIndex]; - outputReal[outputOffset + microphoneIndex] = - Math.cos(phase) * calibration.gains[microphoneIndex]; - outputImaginary[outputOffset + microphoneIndex] = - Math.sin(phase) * calibration.gains[microphoneIndex]; + const gain = applyCalibrationGains ? calibration.gains[microphoneIndex] : 1; + outputReal[outputOffset + microphoneIndex] = Math.cos(phase) * gain; + outputImaginary[outputOffset + microphoneIndex] = Math.sin(phase) * gain; } } @@ -113,6 +170,7 @@ async function createSteeringMatrix(options) { fftSize, speedOfSound, calibration, + applyCalibrationGains = true, signal, pointBatchSize, } = options; @@ -127,6 +185,7 @@ async function createSteeringMatrix(options) { distances, distanceOffset: offset, calibration, + applyCalibrationGains, outputReal: real, outputImaginary: imaginary, outputOffset: offset, @@ -158,21 +217,32 @@ export async function computeBinMap(options) { steeringImaginary: providedSteeringImaginary, signal, pointBatchSize = 64, + algorithm: algorithmId = "das", + phatRelativeFloor = DEFAULT_PHAT_RELATIVE_FLOOR, } = options; + const algorithm = requireLocalizationAlgorithmDefinition(algorithmId); + if (!isDasAlgorithm(algorithm.id)) { + throw new RangeError(`computeBinMap does not implement ${algorithm.id}.`); + } + const phaseOnly = algorithm.spectralWeighting === "phat"; const microphoneCount = arrayDefinition.microphones.length; const output = providedOutput ?? new Float32Array(pointCount); const steeringReal = providedSteeringReal ?? new Float64Array(microphoneCount); const steeringImaginary = providedSteeringImaginary ?? new Float64Array(microphoneCount); const angularWaveNumber = (TAU * bin * sampleRate) / (fftSize * speedOfSound); const microphoneScale = 1 / microphoneCount; - const calibration = providedCalibration ?? createBinCalibration( + const calibration = providedCalibration ?? createBinCalibration({ arrayDefinition, bin, sampleRate, fftSize, compensateFrequencyResponse, maxCompensationDb, - ); + includeAmplitude: !phaseOnly, + }); + const inputSpectrum = phaseOnly + ? phaseNormalizeSpectrum(spectrum, frameCount, microphoneCount, phatRelativeFloor) + : spectrum; const phaseRealValues = steeringMatrix?.real ?? steeringReal; const phaseImaginaryValues = steeringMatrix?.imaginary ?? steeringImaginary; if (output.length !== pointCount) throw new RangeError("output length must match pointCount."); @@ -208,6 +278,7 @@ export async function computeBinMap(options) { outputReal: steeringReal, outputImaginary: steeringImaginary, outputOffset: 0, + applyCalibrationGains: !phaseOnly, }); } @@ -218,8 +289,8 @@ export async function computeBinMap(options) { let sumImaginary = 0; for (let microphoneIndex = 0; microphoneIndex < microphoneCount; microphoneIndex += 1) { const spectrumIndex = spectrumOffset + microphoneIndex; - const inputReal = spectrum.real[spectrumIndex]; - const inputImaginary = spectrum.imaginary[spectrumIndex]; + const inputReal = inputSpectrum.real[spectrumIndex]; + const inputImaginary = inputSpectrum.imaginary[spectrumIndex]; const phaseReal = phaseRealValues[steeringOffset + microphoneIndex]; const phaseImaginary = phaseImaginaryValues[steeringOffset + microphoneIndex]; sumReal += inputReal * phaseReal - inputImaginary * phaseImaginary; @@ -265,6 +336,12 @@ export function createDasAnalysisPlan(options) { if (!Number.isSafeInteger(maximumPlanBytes) || maximumPlanBytes <= 0) { throw new RangeError("maximumPlanBytes must be a positive safe integer."); } + const algorithm = requireLocalizationAlgorithmDefinition( + experiment.beamforming?.algorithm ?? "das", + ); + if (!isDasAlgorithm(algorithm.id)) { + throw new RangeError(`createDasAnalysisPlan does not implement ${algorithm.id}.`); + } const { fftSize } = experiment.stft; const bandBins = experiment.frequencyBands.map((band) => binsForBand(band, sampleRate, fftSize)); @@ -319,6 +396,7 @@ export function createDasAnalysisPlan(options) { arrayDefinition, scene, experiment, + algorithm, microphoneCount, pointCount, bandBins, @@ -394,14 +472,16 @@ function planCalibration(plan, bin) { plan.statistics.calibrationCacheHits += 1; return calibration; } - calibration = createBinCalibration( - plan.arrayDefinition, + calibration = createBinCalibration({ + arrayDefinition: plan.arrayDefinition, bin, - plan.sampleRate, - plan.experiment.stft.fftSize, - plan.experiment.beamforming.compensateFrequencyResponse !== false, - plan.experiment.beamforming.maxCompensationDb, - ); + sampleRate: plan.sampleRate, + fftSize: plan.experiment.stft.fftSize, + compensateFrequencyResponse: + plan.experiment.beamforming.compensateFrequencyResponse !== false, + maxCompensationDb: plan.experiment.beamforming.maxCompensationDb, + includeAmplitude: plan.algorithm.spectralWeighting !== "phat", + }); plan.calibrations.set(bin, calibration); plan.statistics.calibrationBuilds += 1; plan.statistics.retainedBytes += calibration.gains.byteLength + calibration.phases.byteLength; @@ -431,6 +511,7 @@ async function planSteeringMatrix(plan, bin, calibration, signal, pointBatchSize fftSize: plan.experiment.stft.fftSize, speedOfSound: plan.scene.speedOfSound, calibration, + applyCalibrationGains: plan.algorithm.spectralWeighting !== "phat", signal, pointBatchSize, }); @@ -520,7 +601,10 @@ export function estimateDasWorkload(options) { }; } -export function finalizeDasBands(rawMaps, workload, experiment) { +export function finalizeDasBands(rawMaps, workload, experiment, rawScale = "power") { + if (rawScale !== "power" && rawScale !== "score") { + throw new RangeError('rawScale must be "power" or "score".'); + } return experiment.frequencyBands.map((band, bandIndex) => { const rawMap = rawMaps[bandIndex]; for (let pointIndex = 0; pointIndex < workload.pointCount; pointIndex += 1) { @@ -528,7 +612,7 @@ export function finalizeDasBands(rawMaps, workload, experiment) { } const rawPeak = findPeak(rawMap); const displayMap = Float32Array.from(rawMap); - let displayScale = "power"; + let displayScale = rawScale; if (experiment.beamforming.normalization === "peak" && rawPeak.value > 0) { displayScale = "normalized"; for (let pointIndex = 0; pointIndex < workload.pointCount; pointIndex += 1) { @@ -540,12 +624,17 @@ export function finalizeDasBands(rawMaps, workload, experiment) { color: band.color, bins: workload.bandBins[bandIndex], rawMap, + rawScale, displayMap, displayScale, peak: { index: rawPeak.index, position: scanIndexToPosition(rawPeak.index, experiment.scan), - rawPower: rawPeak.value, + rawValue: rawPeak.value, + rawScale, + ...(rawScale === "power" + ? { rawPower: rawPeak.value } + : { rawScore: rawPeak.value }), }, }; }); @@ -566,7 +655,14 @@ async function computeDasMapInternal(options, analysisPlan) { maximumSpectrumBatchBins = Number.MAX_SAFE_INTEGER, pointBatchSize = 64, stftTransformBatchSize = 32, + phatRelativeFloor = DEFAULT_PHAT_RELATIVE_FLOOR, } = options; + const algorithm = analysisPlan?.algorithm ?? requireLocalizationAlgorithmDefinition( + experiment.beamforming?.algorithm ?? "das", + ); + if (!isDasAlgorithm(algorithm.id)) { + throw new RangeError(`computeDasMap does not implement ${algorithm.id}.`); + } if (channels.length !== arrayDefinition.microphones.length) { throw new Error(`채널 수 ${channels.length}와 마이크 수 ${arrayDefinition.microphones.length}가 다릅니다.`); } @@ -685,6 +781,8 @@ async function computeDasMapInternal(options, analysisPlan) { steeringImaginary: analysisPlan?.steeringImaginary, signal, pointBatchSize, + algorithm: algorithm.id, + phatRelativeFloor, }); if (analysisPlan) { analysisPlan.statistics.beamformingComplexMacs += @@ -701,9 +799,10 @@ async function computeDasMapInternal(options, analysisPlan) { } } - const bands = finalizeDasBands(rawMaps, workload, experiment); + const bands = finalizeDasBands(rawMaps, workload, experiment, algorithm.rawScale); return { + algorithm: algorithm.id, backend: "cpu", frameCount: workload.frameCount, microphoneCount: workload.microphoneCount, diff --git a/packages/acoustics-core/evaluation.js b/packages/acoustics-core/evaluation.js index a7d651b..49045e4 100644 --- a/packages/acoustics-core/evaluation.js +++ b/packages/acoustics-core/evaluation.js @@ -50,8 +50,12 @@ function normalizeRawMap(rawMap, scan) { if (!Array.isArray(values) && !ArrayBuffer.isView(values)) { throw new TypeError("rawMap must be an Array, a typed array, or a map descriptor with values."); } - if (descriptor?.scale !== undefined && descriptor.scale !== "power") { - throw new RangeError("rawMap scale must be power."); + if ( + descriptor?.scale !== undefined && + descriptor.scale !== "power" && + descriptor.scale !== "score" + ) { + throw new RangeError("rawMap scale must be power or score."); } if (descriptor?.rows !== undefined && descriptor.rows !== scan.rows) { throw new RangeError("rawMap rows do not match scan.rows."); @@ -69,6 +73,13 @@ function normalizeRawMap(rawMap, scan) { return values; } +function rawMapScale(rawMap) { + if (rawMap && !Array.isArray(rawMap) && !ArrayBuffer.isView(rawMap)) { + return rawMap.scale ?? "power"; + } + return "power"; +} + function normalizeSources(sources) { if (!Array.isArray(sources)) throw new TypeError("sources must be an array."); const ids = new Set(); @@ -420,12 +431,14 @@ function collectHalfPowerIndices(values, scan, peakIndex, halfPowerRatio) { return { indices: queue.subarray(0, queueEnd), thresholdValue }; } -function measureBeamSpreadFromValues(values, scan, peakIndex, halfPowerRatio) { +function measureBeamSpreadFromValues(values, scan, peakIndex, halfPowerRatio, mapScale = "power") { const component = collectHalfPowerIndices(values, scan, peakIndex, halfPowerRatio); const peakPosition = scanIndexToPosition(peakIndex, scan); const base = { calibrated: false, - method: "connected-half-power-power-weighted-covariance", + method: mapScale === "score" + ? "connected-half-maximum-score-weighted-covariance" + : "connected-half-power-power-weighted-covariance", halfPowerRatio, thresholdValue: component.thresholdValue, pixelCount: component.indices.length, @@ -506,18 +519,28 @@ export function measureBeamSpread(rawMap, scan, peakIndex, halfPowerRatio = 0.5) throw new RangeError("peakIndex is outside rawMap."); } requireRatio("halfPowerRatio", halfPowerRatio, { positive: true }); - return measureBeamSpreadFromValues(values, scan, peakIndex, halfPowerRatio); + return measureBeamSpreadFromValues( + values, + scan, + peakIndex, + halfPowerRatio, + rawMapScale(rawMap), + ); } /** - * Evaluate one linear-power DAS band. Source z coordinates are projected onto - * scan.z, matching the image-plane ground-truth convention. + * Evaluate one nonnegative spatial power or score map. Source z coordinates + * are projected onto scan.z, matching the image-plane ground-truth convention. */ -export function evaluateDasBand(options) { +export function evaluateLocalizationBand(options) { if (!options || typeof options !== "object") throw new TypeError("evaluation options are required."); const { rawMap, scan } = options; validateScan(scan); const values = normalizeRawMap(rawMap, scan); + const mapScale = options.mapScale ?? rawMapScale(rawMap); + if (mapScale !== "power" && mapScale !== "score") { + throw new RangeError("mapScale must be power or score."); + } const sources = normalizeSources(options.sources ?? options.scene?.sources ?? []); const truths = sources.map((source) => ({ sourceId: source.id, @@ -558,7 +581,13 @@ export function evaluateDasBand(options) { const detections = selectDetections(localPeaks, thresholdValue, minimumPeakDistanceM, maximumDetections).map( (detection) => ({ ...detection, - beamSpread: measureBeamSpreadFromValues(values, scan, detection.index, halfPowerRatio), + beamSpread: measureBeamSpreadFromValues( + values, + scan, + detection.index, + halfPowerRatio, + mapScale, + ), }), ); const { sourceMatches, truePositiveCount } = matchTruths(truths, detections, matchToleranceM); @@ -604,6 +633,7 @@ export function evaluateDasBand(options) { mainLobeWidthM: mainLobe?.equivalentDiameterM ?? 0, }, diagnostics: { + mapScale, mainLobe, mainLobeWidthM: mainLobe?.equivalentDiameterM ?? 0, gridResolutionM: resolution, @@ -611,6 +641,11 @@ export function evaluateDasBand(options) { }; } +/** Backwards-compatible name for callers that evaluate conventional DAS maps. */ +export function evaluateDasBand(options) { + return evaluateLocalizationBand(options); +} + /** Select only metric names accepted by result.schema.json. */ export function selectResultMetrics(evaluationOrMetrics, metricNames = RESULT_METRIC_KEYS) { const metrics = evaluationOrMetrics?.metrics ?? evaluationOrMetrics; diff --git a/packages/acoustics-core/gcc-phat.js b/packages/acoustics-core/gcc-phat.js new file mode 100644 index 0000000..88b2996 --- /dev/null +++ b/packages/acoustics-core/gcc-phat.js @@ -0,0 +1,783 @@ +import { + createBinCalibration, + createDistanceMatrix, + finalizeDasBands, +} from "./das.js"; +import { cooperativeYield } from "./cooperative-yield.js"; +import { fftInPlace } from "./fft.js"; +import { + binsForBand, + createStftWorkspace, + extractStftBins, + getFrameCount, +} from "./stft.js"; + +export const GCC_PHAT_ANALYSIS_PLAN_KIND = "micarrsim.gcc-phat-analysis-plan/1"; +export const DEFAULT_GCC_PHAT_FLOOR_RATIO = 1e-6; +export const DEFAULT_MAXIMUM_GCC_PHAT_PLAN_BYTES = 512 * 1024 * 1024; + +function createAbortError(signal) { + if (signal?.reason instanceof Error && signal.reason.name === "AbortError") { + return signal.reason; + } + const error = new Error("GCC-PHAT localization was canceled."); + error.name = "AbortError"; + return error; +} + +function throwIfAborted(signal) { + if (signal?.aborted) throw createAbortError(signal); +} + +function requirePositiveFinite(value, label) { + if (!Number.isFinite(value) || value <= 0) { + throw new RangeError(label + " must be a positive finite number."); + } + return value; +} + +function requirePositiveInteger(value, label) { + if (!Number.isInteger(value) || value < 1) { + throw new RangeError(label + " must be a positive integer."); + } + return value; +} + +function validateArrayDefinition(arrayDefinition) { + if (!Array.isArray(arrayDefinition?.microphones) || arrayDefinition.microphones.length < 2) { + throw new RangeError("GCC-PHAT requires at least two microphones."); + } + arrayDefinition.microphones.forEach((microphone, index) => { + if ( + !Array.isArray(microphone?.position) || + microphone.position.length !== 3 || + !microphone.position.every(Number.isFinite) + ) { + throw new TypeError( + "arrayDefinition.microphones[" + index + "].position must be a finite [x, y, z] array.", + ); + } + }); + return arrayDefinition; +} + +function validateConfiguration(options) { + const sampleRate = requirePositiveFinite(options.sampleRate, "sampleRate"); + const arrayDefinition = validateArrayDefinition(options.arrayDefinition); + const speedOfSound = requirePositiveFinite(options.scene?.speedOfSound, "scene.speedOfSound"); + const experiment = options.experiment; + if (!experiment?.stft || !experiment?.scan || !Array.isArray(experiment.frequencyBands)) { + throw new TypeError("experiment must contain stft, scan, and frequencyBands."); + } + const fftSize = requirePositiveInteger(experiment.stft.fftSize, "experiment.stft.fftSize"); + if ((fftSize & (fftSize - 1)) !== 0) { + throw new RangeError("experiment.stft.fftSize must be a power of two."); + } + const hopSize = requirePositiveInteger(experiment.stft.hopSize, "experiment.stft.hopSize"); + if (hopSize > fftSize) { + throw new RangeError("experiment.stft.hopSize must not exceed fftSize."); + } + const pointCount = + requirePositiveInteger(experiment.scan.columns, "experiment.scan.columns") * + requirePositiveInteger(experiment.scan.rows, "experiment.scan.rows"); + if (!Number.isSafeInteger(pointCount)) { + throw new RangeError("scan point count exceeds the safe integer range."); + } + if (experiment.frequencyBands.length === 0) { + throw new RangeError("experiment.frequencyBands must not be empty."); + } + return { sampleRate, arrayDefinition, speedOfSound, experiment, fftSize, hopSize, pointCount }; +} + +function validateChannels(channels, microphoneCount) { + if (!Array.isArray(channels) || channels.length !== microphoneCount) { + throw new Error( + "channel count " + (channels?.length ?? 0) + + " does not match microphone count " + microphoneCount + ".", + ); + } + if (!channels.every((channel) => channel instanceof Float32Array)) { + throw new TypeError("every channel must be a Float32Array."); + } + const sampleCount = channels[0].length; + if (sampleCount === 0) throw new RangeError("channels must contain at least one sample."); + if (channels.some((channel) => channel.length !== sampleCount)) { + throw new Error("all channels must have the same sample count."); + } + return sampleCount; +} + +function referenceTieBreak(left, right) { + if (left.id < right.id) return -1; + if (left.id > right.id) return 1; + return left.index - right.index; +} + +/** + * Select a stable reference without depending on microphone input order. + * + * A peripheral microphone generally provides a longer TDOA baseline. Exact + * geometric ties are resolved by microphone id and finally by input index. + */ +export function selectGccPhatReference(arrayDefinition) { + validateArrayDefinition(arrayDefinition); + const centroid = [0, 0, 0]; + for (const microphone of arrayDefinition.microphones) { + for (let axis = 0; axis < 3; axis += 1) centroid[axis] += microphone.position[axis]; + } + for (let axis = 0; axis < 3; axis += 1) { + centroid[axis] /= arrayDefinition.microphones.length; + } + + let selected; + arrayDefinition.microphones.forEach((microphone, index) => { + const descriptor = { + index, + id: String(microphone.id ?? ""), + microphone, + distanceSquared: + (microphone.position[0] - centroid[0]) ** 2 + + (microphone.position[1] - centroid[1]) ** 2 + + (microphone.position[2] - centroid[2]) ** 2, + }; + if (!selected) { + selected = descriptor; + return; + } + if ( + descriptor.distanceSquared > selected.distanceSquared || + ( + descriptor.distanceSquared === selected.distanceSquared && + referenceTieBreak(descriptor, selected) < 0 + ) + ) { + selected = descriptor; + } + }); + + return Object.freeze({ + index: selected.index, + id: selected.id, + microphone: selected.microphone, + }); +} + +function createBandGeometry(configuration) { + const bandBins = configuration.experiment.frequencyBands.map((band) => + binsForBand(band, configuration.sampleRate, configuration.fftSize) + ); + const selectedBins = [...new Set(bandBins.flat())].sort((left, right) => left - right); + return { bandBins, selectedBins }; +} + +function maximumGeometricLagSamples({ + distances, + pointCount, + microphoneCount, + referenceIndex, + sampleRate, + speedOfSound, +}) { + let maximumLagSamples = 0; + for (let pointIndex = 0; pointIndex < pointCount; pointIndex += 1) { + const distanceOffset = pointIndex * microphoneCount; + const referenceDistance = distances[distanceOffset + referenceIndex]; + for (let microphoneIndex = 0; microphoneIndex < microphoneCount; microphoneIndex += 1) { + if (microphoneIndex === referenceIndex) continue; + const lagSamples = + Math.abs(distances[distanceOffset + microphoneIndex] - referenceDistance) * + sampleRate / speedOfSound; + maximumLagSamples = Math.max(maximumLagSamples, lagSamples); + } + } + return maximumLagSamples; +} + +function assertUnaliasedLag(maximumLagSamples, fftSize) { + const maximumUnaliasedLag = fftSize / 2; + if (maximumLagSamples > maximumUnaliasedLag) { + throw new RangeError( + "GCC-PHAT maximum absolute geometric lag " + maximumLagSamples.toFixed(3) + + " samples exceeds fftSize/2 (" + maximumUnaliasedLag + + "); increase fftSize or reduce the array aperture or scan range.", + ); + } +} + +function createGeometry(configuration, distances) { + const microphoneCount = configuration.arrayDefinition.microphones.length; + const reference = selectGccPhatReference(configuration.arrayDefinition); + const pairMicrophoneIndexes = []; + for (let microphoneIndex = 0; microphoneIndex < microphoneCount; microphoneIndex += 1) { + if (microphoneIndex !== reference.index) pairMicrophoneIndexes.push(microphoneIndex); + } + const resolvedDistances = + distances ?? createDistanceMatrix(configuration.arrayDefinition, configuration.experiment.scan); + const maximumLagSamples = maximumGeometricLagSamples({ + distances: resolvedDistances, + pointCount: configuration.pointCount, + microphoneCount, + referenceIndex: reference.index, + sampleRate: configuration.sampleRate, + speedOfSound: configuration.speedOfSound, + }); + assertUnaliasedLag(maximumLagSamples, configuration.fftSize); + return { + microphoneCount, + reference, + pairMicrophoneIndexes, + pairCount: pairMicrophoneIndexes.length, + distances: resolvedDistances, + maximumLagSamples, + ...createBandGeometry(configuration), + }; +} + +function requirePhatFloorRatio(value) { + const ratio = value ?? DEFAULT_GCC_PHAT_FLOOR_RATIO; + if (!Number.isFinite(ratio) || ratio < 0 || ratio > 1) { + throw new RangeError("phatFloorRatio must be a finite number from 0 through 1."); + } + return ratio; +} + +function estimateWithGeometry(options, configuration, geometry) { + const sampleCount = validateChannels(options.channels, geometry.microphoneCount); + const frameCount = getFrameCount( + sampleCount, + configuration.fftSize, + configuration.hopSize, + ); + const pairCount = geometry.pairCount; + const bandCount = geometry.bandBins.length; + const bandBinCount = geometry.bandBins.reduce((sum, bins) => sum + bins.length, 0); + const selectedBinCount = geometry.selectedBins.length; + const stftTransformCount = frameCount * geometry.microphoneCount; + const inverseFftCount = pairCount * bandCount; + const fftButterfliesPerTransform = + (configuration.fftSize / 2) * Math.log2(configuration.fftSize); + const stftButterflyCount = stftTransformCount * fftButterfliesPerTransform; + const inverseFftButterflyCount = inverseFftCount * fftButterfliesPerTransform; + const crossSpectrumSampleCount = pairCount * frameCount * bandBinCount; + // The relative-floor pass finds the maximum magnitude, then the PHAT pass + // recomputes and accumulates each complex cross-spectrum sample. + const crossSpectrumComplexOpCount = crossSpectrumSampleCount * 2; + const scanLookupCount = configuration.pointCount * pairCount * bandCount; + const dasEquivalentComplexMacCount = + configuration.pointCount * geometry.microphoneCount * frameCount * selectedBinCount; + const distanceBytes = + configuration.pointCount * geometry.microphoneCount * Float32Array.BYTES_PER_ELEMENT; + const stftWorkspaceBytes = + configuration.fftSize * Float64Array.BYTES_PER_ELEMENT * 3; + const correlationWorkspaceBytes = + configuration.fftSize * Float64Array.BYTES_PER_ELEMENT * 2; + const spectrumBytes = + selectedBinCount * frameCount * geometry.microphoneCount * + Float32Array.BYTES_PER_ELEMENT * 2; + const calibrationPhaseBytes = + selectedBinCount * geometry.microphoneCount * Float64Array.BYTES_PER_ELEMENT; + const outputBytes = + bandCount * configuration.pointCount * Float32Array.BYTES_PER_ELEMENT * 2; + const retainedPlanBytes = + options.analysisPlan?.retainedBytes ?? + distanceBytes + stftWorkspaceBytes + correlationWorkspaceBytes; + const estimatedWorkingBytes = + retainedPlanBytes + spectrumBytes + calibrationPhaseBytes + outputBytes; + + return { + bandBins: geometry.bandBins, + selectedBins: geometry.selectedBins, + selectedBinCount, + bandBinCount, + frameCount, + microphoneCount: geometry.microphoneCount, + pointCount: configuration.pointCount, + bandCount, + pairCount, + referenceMicrophoneId: geometry.reference.id, + referenceMicrophoneIndex: geometry.reference.index, + maximumLagSamples: geometry.maximumLagSamples, + stftTransformCount, + inverseFftCount, + stftButterflyCount, + inverseFftButterflyCount, + crossSpectrumSampleCount, + crossSpectrumComplexOpCount, + crossSpectrumCount: crossSpectrumComplexOpCount, + complexMacCount: crossSpectrumComplexOpCount, + scanLookupCount, + operationCount: crossSpectrumComplexOpCount + scanLookupCount, + dasEquivalentComplexMacCount, + distanceBytes, + stftWorkspaceBytes, + correlationWorkspaceBytes, + spectrumBytes, + calibrationPhaseBytes, + outputBytes, + retainedPlanBytes, + estimatedWorkingBytes, + }; +} + +/** + * Estimate CPU work for one reference-pair GCC-PHAT analysis. + */ +export function estimateGccPhatWorkload(options) { + if (!options || typeof options !== "object") { + throw new TypeError("options must be an object."); + } + const configuration = validateConfiguration(options); + let geometry; + if (options.analysisPlan) { + assertCompatibleAnalysisPlan(options.analysisPlan, options); + geometry = options.analysisPlan; + } else { + geometry = createGeometry(configuration); + } + return estimateWithGeometry(options, configuration, geometry); +} + +/** + * Retain geometry and FFT workspaces for repeated analyses with one immutable + * array/scan/band configuration. Plans are stateful and cannot run concurrently. + */ +export function createGccPhatAnalysisPlan(options) { + if (!options || typeof options !== "object") { + throw new TypeError("options must be an object."); + } + const configuration = validateConfiguration(options); + const maximumPlanBytes = + options.maximumPlanBytes ?? DEFAULT_MAXIMUM_GCC_PHAT_PLAN_BYTES; + if (!Number.isSafeInteger(maximumPlanBytes) || maximumPlanBytes <= 0) { + throw new RangeError("maximumPlanBytes must be a positive safe integer."); + } + + const microphoneCount = configuration.arrayDefinition.microphones.length; + const pointMicrophonePairs = configuration.pointCount * microphoneCount; + const distanceBytes = pointMicrophonePairs * Float32Array.BYTES_PER_ELEMENT; + const stftWorkspaceBytes = + configuration.fftSize * Float64Array.BYTES_PER_ELEMENT * 3; + const correlationWorkspaceBytes = + configuration.fftSize * Float64Array.BYTES_PER_ELEMENT * 2; + const retainedBytes = distanceBytes + stftWorkspaceBytes + correlationWorkspaceBytes; + if ( + !Number.isSafeInteger(pointMicrophonePairs) || + !Number.isSafeInteger(retainedBytes) || + retainedBytes > maximumPlanBytes + ) { + throw new RangeError( + "GCC-PHAT analysis plan requires " + retainedBytes.toLocaleString() + + " bytes, above maximumPlanBytes " + maximumPlanBytes.toLocaleString() + ".", + ); + } + + const distances = createDistanceMatrix( + configuration.arrayDefinition, + configuration.experiment.scan, + ); + const geometry = createGeometry(configuration, distances); + const plan = { + kind: GCC_PHAT_ANALYSIS_PLAN_KIND, + sampleRate: configuration.sampleRate, + arrayDefinition: configuration.arrayDefinition, + scene: options.scene, + experiment: configuration.experiment, + fftSize: configuration.fftSize, + hopSize: configuration.hopSize, + pointCount: configuration.pointCount, + maximumPlanBytes, + retainedBytes, + distanceBytes, + stftWorkspaceBytes, + correlationWorkspaceBytes, + stftWorkspace: createStftWorkspace(configuration.fftSize), + correlationReal: new Float64Array(configuration.fftSize), + correlationImaginary: new Float64Array(configuration.fftSize), + busy: false, + statistics: { + analysisAttempts: 0, + completedAnalyses: 0, + distanceMatrixBuilds: 1, + distanceEvaluations: pointMicrophonePairs, + avoidedDistanceEvaluations: 0, + stftWorkspaceAllocations: 1, + stftWorkspaceUses: 0, + avoidedStftWorkspaceAllocations: 0, + correlationWorkspaceAllocations: 1, + correlationWorkspaceUses: 0, + avoidedCorrelationWorkspaceAllocations: 0, + stftTransforms: 0, + inverseFftTransforms: 0, + crossSpectrumOperations: 0, + scanLookups: 0, + }, + ...geometry, + }; + return plan; +} + +function assertCompatibleAnalysisPlan(plan, options) { + if (plan?.kind !== GCC_PHAT_ANALYSIS_PLAN_KIND) { + throw new TypeError( + "analysisPlan must be created by createGccPhatAnalysisPlan().", + ); + } + if ( + plan.sampleRate !== options.sampleRate || + plan.arrayDefinition !== options.arrayDefinition || + plan.scene !== options.scene || + plan.experiment !== options.experiment + ) { + throw new Error( + "analysisPlan must be used with the exact configuration objects that created it.", + ); + } +} + +function analysisCacheSummary(plan) { + const statistics = plan.statistics; + return { + analysisAttempts: statistics.analysisAttempts, + completedAnalyses: statistics.completedAnalyses, + distanceMatrixBuilds: statistics.distanceMatrixBuilds, + distanceEvaluations: statistics.distanceEvaluations, + avoidedDistanceEvaluations: statistics.avoidedDistanceEvaluations, + stftWorkspaceAllocations: statistics.stftWorkspaceAllocations, + stftWorkspaceUses: statistics.stftWorkspaceUses, + avoidedStftWorkspaceAllocations: statistics.avoidedStftWorkspaceAllocations, + correlationWorkspaceAllocations: statistics.correlationWorkspaceAllocations, + correlationWorkspaceUses: statistics.correlationWorkspaceUses, + avoidedCorrelationWorkspaceAllocations: + statistics.avoidedCorrelationWorkspaceAllocations, + stftTransforms: statistics.stftTransforms, + inverseFftTransforms: statistics.inverseFftTransforms, + crossSpectrumOperations: statistics.crossSpectrumOperations, + scanLookups: statistics.scanLookups, + retainedBytes: plan.retainedBytes, + maximumPlanBytes: plan.maximumPlanBytes, + }; +} + +function circularLinearSample(values, lagSamples) { + const length = values.length; + const wrapped = ((lagSamples % length) + length) % length; + const firstIndex = Math.floor(wrapped); + const fraction = wrapped - firstIndex; + const secondIndex = (firstIndex + 1) % length; + return values[firstIndex] * (1 - fraction) + values[secondIndex] * fraction; +} + +function fillBandPairCorrelation({ + correlationReal, + correlationImaginary, + bins, + stft, + microphoneIndex, + referenceIndex, + phasesByBin, + phatFloorRatio, +}) { + correlationReal.fill(0); + correlationImaginary.fill(0); + const microphoneCount = stft.microphoneCount; + let maximumMagnitude = 0; + + for (const bin of bins) { + const spectrum = stft.spectra.get(bin); + for (let frameIndex = 0; frameIndex < stft.frameCount; frameIndex += 1) { + const offset = frameIndex * microphoneCount; + const microphoneReal = spectrum.real[offset + microphoneIndex]; + const microphoneImaginary = spectrum.imaginary[offset + microphoneIndex]; + const referenceReal = spectrum.real[offset + referenceIndex]; + const referenceImaginary = spectrum.imaginary[offset + referenceIndex]; + const crossReal = + microphoneReal * referenceReal + microphoneImaginary * referenceImaginary; + const crossImaginary = + microphoneImaginary * referenceReal - microphoneReal * referenceImaginary; + maximumMagnitude = Math.max( + maximumMagnitude, + Math.hypot(crossReal, crossImaginary), + ); + } + } + + const magnitudeFloor = maximumMagnitude * phatFloorRatio; + if (!(maximumMagnitude > 0)) { + fftInPlace(correlationReal, correlationImaginary, true); + return; + } + + for (const bin of bins) { + const spectrum = stft.spectra.get(bin); + const phases = phasesByBin.get(bin); + const correctionPhase = phases[referenceIndex] - phases[microphoneIndex]; + const correctionReal = Math.cos(correctionPhase); + const correctionImaginary = Math.sin(correctionPhase); + let accumulatedReal = 0; + let accumulatedImaginary = 0; + + for (let frameIndex = 0; frameIndex < stft.frameCount; frameIndex += 1) { + const offset = frameIndex * microphoneCount; + const microphoneReal = spectrum.real[offset + microphoneIndex]; + const microphoneImaginary = spectrum.imaginary[offset + microphoneIndex]; + const referenceReal = spectrum.real[offset + referenceIndex]; + const referenceImaginary = spectrum.imaginary[offset + referenceIndex]; + const crossReal = + microphoneReal * referenceReal + microphoneImaginary * referenceImaginary; + const crossImaginary = + microphoneImaginary * referenceReal - microphoneReal * referenceImaginary; + const magnitude = Math.hypot(crossReal, crossImaginary); + const denominator = Math.max(magnitude, magnitudeFloor); + if (!(denominator > 0)) continue; + const normalizedReal = crossReal / denominator; + const normalizedImaginary = crossImaginary / denominator; + accumulatedReal += + normalizedReal * correctionReal - normalizedImaginary * correctionImaginary; + accumulatedImaginary += + normalizedReal * correctionImaginary + normalizedImaginary * correctionReal; + } + + const averagedReal = accumulatedReal / stft.frameCount; + const averagedImaginary = accumulatedImaginary / stft.frameCount; + correlationReal[bin] = averagedReal; + correlationImaginary[bin] = + bin === 0 || bin === correlationReal.length / 2 ? 0 : averagedImaginary; + if (bin > 0 && bin < correlationReal.length / 2) { + const negativeBin = correlationReal.length - bin; + correlationReal[negativeBin] = averagedReal; + correlationImaginary[negativeBin] = -averagedImaginary; + } + } + + fftInPlace(correlationReal, correlationImaginary, true); +} + +function shiftAndSquareScores(rawMaps) { + for (const rawMap of rawMaps) { + let minimum = Infinity; + for (const value of rawMap) minimum = Math.min(minimum, value); + for (let pointIndex = 0; pointIndex < rawMap.length; pointIndex += 1) { + const shifted = rawMap[pointIndex] - minimum; + rawMap[pointIndex] = shifted * shifted; + } + } +} + +async function computeGccPhatMapInternal(options, configuration, geometry, plan) { + const { + channels, + signal, + onProgress, + maximumWorkingBytes = 512 * 1024 * 1024, + maximumFftButterflies = 1_000_000_000, + pointBatchSize = 256, + stftTransformBatchSize = 128, + } = options; + requirePositiveInteger(pointBatchSize, "pointBatchSize"); + requirePositiveInteger(stftTransformBatchSize, "stftTransformBatchSize"); + if (!Number.isSafeInteger(maximumWorkingBytes) || maximumWorkingBytes <= 0) { + throw new RangeError("maximumWorkingBytes must be a positive safe integer."); + } + const phatFloorRatio = requirePhatFloorRatio( + options.phatFloorRatio ?? + options.phatRelativeFloor ?? + configuration.experiment.beamforming?.phatFloorRatio, + ); + const workload = estimateWithGeometry(options, configuration, geometry); + if (workload.estimatedWorkingBytes > maximumWorkingBytes) { + throw new RangeError( + "GCC-PHAT estimated working set " + + (workload.estimatedWorkingBytes / 1024 ** 2).toFixed(1) + + " MiB exceeds maximumWorkingBytes " + + (maximumWorkingBytes / 1024 ** 2).toFixed(1) + " MiB.", + ); + } + const totalFftButterflies = + workload.stftButterflyCount + workload.inverseFftButterflyCount; + if (totalFftButterflies > maximumFftButterflies) { + throw new RangeError( + "GCC-PHAT estimated FFT butterfly count " + + Math.ceil(totalFftButterflies).toLocaleString() + + " exceeds maximumFftButterflies " + + maximumFftButterflies.toLocaleString() + ".", + ); + } + + throwIfAborted(signal); + if (plan) { + plan.statistics.stftWorkspaceUses += 1; + plan.statistics.avoidedStftWorkspaceAllocations = Math.max( + 0, + plan.statistics.stftWorkspaceUses - 1, + ); + plan.statistics.correlationWorkspaceUses += 1; + plan.statistics.avoidedCorrelationWorkspaceAllocations = Math.max( + 0, + plan.statistics.correlationWorkspaceUses - 1, + ); + } + const stft = await extractStftBins(channels, { + fftSize: configuration.fftSize, + hopSize: configuration.hopSize, + bins: geometry.selectedBins, + workspace: plan?.stftWorkspace, + transformBatchSize: stftTransformBatchSize, + signal, + onProgress: ({ completedTransforms, totalTransforms }) => + onProgress?.({ + stage: "stft", + algorithm: "gcc-phat", + completedTransforms, + totalTransforms, + }), + }); + if (plan) { + plan.statistics.stftTransforms += workload.stftTransformCount; + } + await cooperativeYield(); + throwIfAborted(signal); + + const phasesByBin = new Map(); + for (const bin of geometry.selectedBins) { + const calibration = createBinCalibration({ + arrayDefinition: configuration.arrayDefinition, + bin, + sampleRate: configuration.sampleRate, + fftSize: configuration.fftSize, + compensateFrequencyResponse: + configuration.experiment.beamforming?.compensateFrequencyResponse !== false, + maxCompensationDb: + configuration.experiment.beamforming?.maxCompensationDb, + includeAmplitude: false, + }); + phasesByBin.set(bin, calibration.phases); + } + + const correlationReal = + plan?.correlationReal ?? new Float64Array(configuration.fftSize); + const correlationImaginary = + plan?.correlationImaginary ?? new Float64Array(configuration.fftSize); + const rawMaps = geometry.bandBins.map(() => new Float32Array(configuration.pointCount)); + const totalPairs = geometry.pairCount * geometry.bandBins.length; + let completedPairs = 0; + + for (let bandIndex = 0; bandIndex < geometry.bandBins.length; bandIndex += 1) { + const bins = geometry.bandBins[bandIndex]; + const rawMap = rawMaps[bandIndex]; + for (const microphoneIndex of geometry.pairMicrophoneIndexes) { + throwIfAborted(signal); + fillBandPairCorrelation({ + correlationReal, + correlationImaginary, + bins, + stft, + microphoneIndex, + referenceIndex: geometry.reference.index, + phasesByBin, + phatFloorRatio, + }); + if (plan) { + plan.statistics.inverseFftTransforms += 1; + plan.statistics.crossSpectrumOperations += + bins.length * stft.frameCount * 2; + } + + for (let pointIndex = 0; pointIndex < configuration.pointCount; pointIndex += 1) { + const distanceOffset = pointIndex * geometry.microphoneCount; + const predictedLagSamples = + ( + geometry.distances[distanceOffset + microphoneIndex] - + geometry.distances[distanceOffset + geometry.reference.index] + ) * + configuration.sampleRate / configuration.speedOfSound; + rawMap[pointIndex] += circularLinearSample( + correlationReal, + predictedLagSamples, + ); + if ( + (pointIndex + 1) % pointBatchSize === 0 && + pointIndex + 1 < configuration.pointCount + ) { + await cooperativeYield(); + throwIfAborted(signal); + } + } + if (plan) plan.statistics.scanLookups += configuration.pointCount; + completedPairs += 1; + onProgress?.({ + stage: "localization", + algorithm: "gcc-phat", + completedPairs, + totalPairs, + }); + if (completedPairs < totalPairs) await cooperativeYield(); + } + } + + throwIfAborted(signal); + shiftAndSquareScores(rawMaps); + const bands = finalizeDasBands( + rawMaps, + workload, + configuration.experiment, + "score", + ); + return { + backend: "cpu", + algorithm: "gcc-phat", + frameCount: workload.frameCount, + microphoneCount: workload.microphoneCount, + pointCount: workload.pointCount, + referenceMicrophoneId: geometry.reference.id, + referenceMicrophoneIndex: geometry.reference.index, + pairCount: geometry.pairCount, + phatFloorRatio, + maximumLagSamples: geometry.maximumLagSamples, + workload, + bands, + ...(plan ? { analysisCache: analysisCacheSummary(plan) } : {}), + }; +} + +/** + * Build band-limited reference-pair GCC-PHAT maps on the CPU. + * + * For each pair and band, corrected PHAT cross spectra are averaged over STFT + * frames, transformed once to circular correlation, and sampled at the + * near-field geometric delay of each scan point. + */ +export async function computeGccPhatMap(options) { + if (!options || typeof options !== "object") { + throw new TypeError("options must be an object."); + } + const configuration = validateConfiguration(options); + const plan = options.analysisPlan; + if (!plan) { + const geometry = createGeometry(configuration); + return computeGccPhatMapInternal(options, configuration, geometry, undefined); + } + + assertCompatibleAnalysisPlan(plan, options); + if (plan.busy) { + throw new Error("analysisPlan is already in use by another GCC-PHAT analysis."); + } + plan.busy = true; + plan.statistics.analysisAttempts += 1; + if (plan.statistics.analysisAttempts > 1) { + plan.statistics.avoidedDistanceEvaluations += + plan.pointCount * plan.microphoneCount; + } + try { + const result = await computeGccPhatMapInternal( + options, + configuration, + plan, + plan, + ); + plan.statistics.completedAnalyses += 1; + result.analysisCache = analysisCacheSummary(plan); + return result; + } finally { + plan.busy = false; + } +} diff --git a/packages/acoustics-core/index.js b/packages/acoustics-core/index.js index ef034eb..f131829 100644 --- a/packages/acoustics-core/index.js +++ b/packages/acoustics-core/index.js @@ -4,6 +4,9 @@ export * from "./environment.js"; export * from "./fft.js"; export * from "./frequency-response.js"; export * from "./fractional-delay.js"; +export * from "./gcc-phat.js"; +export * from "./localization.js"; +export * from "./localization-algorithms.js"; export * from "./math.js"; export * from "./mixdown.js"; export * from "./noise.js"; diff --git a/packages/acoustics-core/localization-algorithms.js b/packages/acoustics-core/localization-algorithms.js new file mode 100644 index 0000000..21dcee0 --- /dev/null +++ b/packages/acoustics-core/localization-algorithms.js @@ -0,0 +1,52 @@ +export const DEFAULT_LOCALIZATION_ALGORITHM = "das"; +export const DEFAULT_PHAT_RELATIVE_FLOOR = 1e-6; + +function definition(id, label, rawScale, spectralWeighting, backends) { + return Object.freeze({ + id, + label, + rawScale, + spectralWeighting, + backends: Object.freeze([...backends]), + }); +} + +const ALGORITHM_DEFINITIONS = Object.freeze([ + definition("das", "DAS", "power", "none", ["cpu", "webgpu"]), + definition("das-phat", "DAS-PHAT", "score", "phat", ["cpu", "webgpu"]), + definition("gcc-phat", "GCC-PHAT", "score", "phat", ["cpu"]), +]); + +const ALGORITHM_BY_ID = new Map( + ALGORITHM_DEFINITIONS.map((algorithm) => [algorithm.id, algorithm]), +); + +export const LOCALIZATION_ALGORITHM_IDS = Object.freeze( + ALGORITHM_DEFINITIONS.map((algorithm) => algorithm.id), +); + +export function localizationAlgorithmDefinition(id = DEFAULT_LOCALIZATION_ALGORITHM) { + return ALGORITHM_BY_ID.get(id); +} + +export function requireLocalizationAlgorithmDefinition( + id = DEFAULT_LOCALIZATION_ALGORITHM, +) { + const algorithm = localizationAlgorithmDefinition(id); + if (!algorithm) { + throw new RangeError(`unsupported localization algorithm: ${id}`); + } + return algorithm; +} + +export function isDasAlgorithm(id = DEFAULT_LOCALIZATION_ALGORITHM) { + return id === "das" || id === "das-phat"; +} + +export function localizationAlgorithmLabel(id = DEFAULT_LOCALIZATION_ALGORITHM) { + return localizationAlgorithmDefinition(id)?.label ?? String(id); +} + +export function localizationAlgorithmSupportsBackend(id, backend) { + return localizationAlgorithmDefinition(id)?.backends.includes(backend) ?? false; +} diff --git a/packages/acoustics-core/localization.js b/packages/acoustics-core/localization.js new file mode 100644 index 0000000..939f47c --- /dev/null +++ b/packages/acoustics-core/localization.js @@ -0,0 +1,53 @@ +import { + computeDasMap, + createDasAnalysisPlan, + estimateDasWorkload, +} from "./das.js"; +import { + computeGccPhatMap, + createGccPhatAnalysisPlan, + estimateGccPhatWorkload, +} from "./gcc-phat.js"; +import { + isDasAlgorithm, + requireLocalizationAlgorithmDefinition, +} from "./localization-algorithms.js"; + +function requestedAlgorithm(options) { + return requireLocalizationAlgorithmDefinition( + options?.experiment?.beamforming?.algorithm, + ); +} + +/** + * Estimate implementation-level work for whichever localization algorithm the + * experiment requests. + */ +export function estimateLocalizationWorkload(options) { + const algorithm = requestedAlgorithm(options); + return isDasAlgorithm(algorithm.id) + ? estimateDasWorkload(options) + : estimateGccPhatWorkload(options); +} + +/** + * Compile immutable geometry/workspace state for repeated CPU analyses. + */ +export function createLocalizationAnalysisPlan(options) { + const algorithm = requestedAlgorithm(options); + return isDasAlgorithm(algorithm.id) + ? createDasAnalysisPlan(options) + : createGccPhatAnalysisPlan(options); +} + +/** + * Compute a CPU localization map using the algorithm selected by the + * experiment. The returned `algorithm` and each band's `rawScale` make the + * result semantics explicit. + */ +export async function computeLocalizationMap(options) { + const algorithm = requestedAlgorithm(options); + return isDasAlgorithm(algorithm.id) + ? computeDasMap(options) + : computeGccPhatMap(options); +} diff --git a/packages/acoustics-core/realtime-localizer.js b/packages/acoustics-core/realtime-localizer.js index cfe3a43..3661e00 100644 --- a/packages/acoustics-core/realtime-localizer.js +++ b/packages/acoustics-core/realtime-localizer.js @@ -1,5 +1,12 @@ -import { computeDasMap, createDasAnalysisPlan } from "./das.js"; import { findLocalPeaks } from "./evaluation.js"; +import { + computeLocalizationMap, + createLocalizationAnalysisPlan, +} from "./localization.js"; +import { + isDasAlgorithm, + requireLocalizationAlgorithmDefinition, +} from "./localization-algorithms.js"; import { DEFAULT_MAX_COMPENSATION_DB, MAX_FREQUENCY_RESPONSE_HZ, @@ -214,6 +221,10 @@ function throwIfAborted(signal) { } export function detectBandPeaks(rawMap, scan, options = {}) { + const rawScale = options.rawScale ?? rawMap?.scale ?? "power"; + if (rawScale !== "power" && rawScale !== "score") { + throw new RangeError('peak map rawScale must be "power" or "score".'); + } const peakDetection = copyPeakDetection(options, scan); const localPeaks = findLocalPeaks(rawMap, scan); const threshold = (localPeaks[0]?.value ?? 0) * peakDetection.thresholdRatio; @@ -235,7 +246,9 @@ export function detectBandPeaks(rawMap, scan, options = {}) { rank: detections.length + 1, index: peak.index, position: [...peak.position], - rawPower: peak.value, + rawValue: peak.value, + rawScale, + ...(rawScale === "power" ? { rawPower: peak.value } : {}), }); if (detections.length === peakDetection.maximumDetections) break; } @@ -244,7 +257,10 @@ export function detectBandPeaks(rawMap, scan, options = {}) { function locationsFromBands(bands, scan, peakDetection) { return bands.flatMap((band) => - detectBandPeaks(band.rawMap, scan, peakDetection).map((peak) => ({ + detectBandPeaks(band.rawMap, scan, { + ...peakDetection, + rawScale: band.rawScale, + }).map((peak) => ({ bandId: band.id, ...peak, })), @@ -287,6 +303,13 @@ export function createLocalizationPlan(config) { config.beamforming?.maxCompensationDb ?? DEFAULT_MAX_COMPENSATION_DB, ); + const algorithm = requireLocalizationAlgorithmDefinition( + config.algorithm ?? config.beamforming?.algorithm, + ); + const phatRelativeFloor = + config.phatRelativeFloor ?? config.beamforming?.phatRelativeFloor; + const phatFloorRatio = + config.phatFloorRatio ?? config.beamforming?.phatFloorRatio; const scene = Object.freeze({ speedOfSound }); const experiment = Object.freeze({ @@ -294,12 +317,14 @@ export function createLocalizationPlan(config) { scan, frequencyBands, beamforming: Object.freeze({ + algorithm: algorithm.id, + backend: "cpu", normalization, compensateFrequencyResponse, maxCompensationDb, }), }); - const dasPlan = createDasAnalysisPlan({ + const analysisPlan = createLocalizationAnalysisPlan({ sampleRate, arrayDefinition, scene, @@ -310,6 +335,7 @@ export function createLocalizationPlan(config) { return Object.freeze({ kind: REALTIME_LOCALIZATION_PLAN_KIND, + algorithm: algorithm.id, sampleRate, channelCount: arrayDefinition.microphones.length, arrayDefinition, @@ -322,10 +348,15 @@ export function createLocalizationPlan(config) { maxCompensationDb, scene, experiment, - dasPlan, + analysisPlan, + ...(isDasAlgorithm(algorithm.id) + ? { dasPlan: analysisPlan } + : { gccPhatPlan: analysisPlan }), computeOptions: Object.freeze({ pointBatchSize: DEFAULT_REALTIME_POINT_BATCH_SIZE, stftTransformBatchSize: DEFAULT_REALTIME_STFT_TRANSFORM_BATCH_SIZE, + ...(phatRelativeFloor !== undefined ? { phatRelativeFloor } : {}), + ...(phatFloorRatio !== undefined ? { phatFloorRatio } : {}), ...(config.computeOptions ?? {}), }), }); @@ -350,7 +381,7 @@ export async function localizeMultichannel(options) { if (sampleCount === 0) throw new RangeError("channels must contain at least one sample."); throwIfAborted(options.signal); const processingStartedAt = monotonicNow(); - const das = await computeDasMap({ + const localization = await computeLocalizationMap({ ...plan.computeOptions, ...(options.computeOptions ?? {}), channels: options.channels, @@ -358,14 +389,18 @@ export async function localizeMultichannel(options) { arrayDefinition: plan.arrayDefinition, scene: plan.scene, experiment: plan.experiment, - analysisPlan: plan.dasPlan, + analysisPlan: plan.analysisPlan, signal: options.signal, onProgress: options.onProgress, }); - const locations = locationsFromBands(das.bands, plan.scan, plan.peakDetection); + const locations = locationsFromBands( + localization.bands, + plan.scan, + plan.peakDetection, + ); const processingTimeMs = Math.max(0, monotonicNow() - processingStartedAt); return { - ...das, + ...localization, sampleRate: plan.sampleRate, sampleCount, duration: sampleCount / plan.sampleRate, @@ -488,6 +523,7 @@ export class StreamingDasLocalizer { getState() { return { sampleRate: this.plan.sampleRate, + algorithm: this.plan.algorithm, channelCount: this.plan.channelCount, analysisWindowSamples: this.analysisWindowSamples, analysisHopSamples: this.analysisHopSamples, @@ -675,6 +711,9 @@ export class StreamingDasLocalizer { } } +/** Algorithm-neutral alias; the older class name remains for compatibility. */ +export { StreamingDasLocalizer as StreamingLocalizer }; + export function createStreamingLocalizer(config) { return new StreamingDasLocalizer(config); } diff --git a/packages/contracts/validation.js b/packages/contracts/validation.js index 49a67e5..b5aa2c3 100644 --- a/packages/contracts/validation.js +++ b/packages/contracts/validation.js @@ -6,6 +6,7 @@ const MAX_SOURCES = 32; const MAX_GRID_SIDE = 512; const MAX_FREQUENCY_RESPONSE_POINTS = 64; const MAX_FREQUENCY_RESPONSE_HZ = 1_000_000; +const LOCALIZATION_ALGORITHMS = ["das", "das-phat", "gcc-phat"]; /** @typedef {{ path: string, message: string }} ValidationIssue */ @@ -436,10 +437,26 @@ function validateExperimentRelationships(value, context = {}) { if (value.stft.window !== "hann") addIssue(issues, "$.stft.window", "hann이어야 합니다."); } if (requireRecord(value.beamforming, "$.beamforming", issues)) { - if (value.beamforming.algorithm !== "das") addIssue(issues, "$.beamforming.algorithm", "das이어야 합니다."); + if (!LOCALIZATION_ALGORITHMS.includes(value.beamforming.algorithm)) { + addIssue( + issues, + "$.beamforming.algorithm", + "das, das-phat 또는 gcc-phat이어야 합니다.", + ); + } if (!["auto", "cpu", "webgpu"].includes(value.beamforming.backend)) { addIssue(issues, "$.beamforming.backend", "auto, cpu 또는 webgpu여야 합니다."); } + if ( + value.beamforming.algorithm === "gcc-phat" && + value.beamforming.backend === "webgpu" + ) { + addIssue( + issues, + "$.beamforming.backend", + "gcc-phat는 현재 CPU backend만 지원합니다. auto 또는 cpu를 선택하세요.", + ); + } requireFinite(value.beamforming.minimumDistance, "$.beamforming.minimumDistance", issues, { min: 0, max: 10, @@ -475,6 +492,16 @@ function validateExperimentRelationships(value, context = {}) { } } if (context.arrayDefinition && context.scene) { + if ( + value.beamforming?.algorithm === "gcc-phat" && + context.arrayDefinition.microphones.length < 2 + ) { + addIssue( + issues, + "$.beamforming.algorithm", + "gcc-phat에는 서로 다른 두 개 이상의 마이크가 필요합니다.", + ); + } const retainedBufferCount = context.arrayDefinition.microphones.length + context.scene.sources.length + 1; const sampleCount = Math.ceil(context.scene.duration * context.scene.sampleRate); let maximumDelaySamples = 0; @@ -519,6 +546,12 @@ function validateResultRelationships(value, context = {}) { if (value.experimentId !== context.experiment.id) addIssue(issues, "$.experimentId", "불러온 실험 id와 다릅니다."); if (value.arrayId !== context.experiment.arrayId) addIssue(issues, "$.arrayId", "실험의 배열 id와 다릅니다."); if (value.sceneId !== context.experiment.sceneId) addIssue(issues, "$.sceneId", "실험의 장면 id와 다릅니다."); + if ( + value.algorithm !== undefined && + value.algorithm !== context.experiment.beamforming?.algorithm + ) { + addIssue(issues, "$.algorithm", "실험의 위치 추정 알고리즘과 다릅니다."); + } } const sourceIds = new Set(context.scene?.sources?.map((source) => source.id) ?? []); const expectedBandIds = new Set(context.experiment?.frequencyBands?.map((band) => band.id) ?? []); @@ -541,8 +574,8 @@ function validateResultRelationships(value, context = {}) { addIssue(issues, `${path}.maps[${mapIndex}]`, "map shape이 실험 scan grid와 다릅니다."); } }); - if (isRecord(band.rawMap) && band.rawMap.scale !== "power") { - addIssue(issues, `${path}.rawMap.scale`, "rawMap은 power scale이어야 합니다."); + if (isRecord(band.rawMap) && !["power", "score"].includes(band.rawMap.scale)) { + addIssue(issues, `${path}.rawMap.scale`, "rawMap은 power 또는 score scale이어야 합니다."); } if (Array.isArray(band.peaks)) { band.peaks.forEach((peak, peakIndex) => { diff --git a/packages/webgpu-kernels/compute.js b/packages/webgpu-kernels/compute.js index a7f1080..11ee408 100644 --- a/packages/webgpu-kernels/compute.js +++ b/packages/webgpu-kernels/compute.js @@ -1,4 +1,9 @@ import { computeDasMap } from "../acoustics-core/das.js"; +import { computeGccPhatMap } from "../acoustics-core/gcc-phat.js"; +import { + isDasAlgorithm, + requireLocalizationAlgorithmDefinition, +} from "../acoustics-core/localization-algorithms.js"; import { computeDasWebGpu } from "./das.js"; import { isRecoverableWebGpuError } from "./device.js"; @@ -17,7 +22,7 @@ function throwIfAborted(signal) { * Execute the requested DAS backend. `auto` falls back only for typed WebGPU * operational errors; invalid inputs and AbortError remain visible to callers. */ -export async function computeDasWithBackend(options, dependencies = {}) { +async function computeDasBackend(options, dependencies = {}) { const computeGpu = dependencies.computeGpu ?? computeDasWebGpu; const computeCpu = dependencies.computeCpu ?? computeDasMap; const requestedBackend = options?.experiment?.beamforming?.backend ?? "auto"; @@ -56,3 +61,48 @@ export async function computeDasWithBackend(options, dependencies = {}) { } } } + +/** + * Execute a DAS-family algorithm with the requested CPU/WebGPU policy. + * Conventional DAS callers can keep using this backwards-compatible entry. + */ +export async function computeDasWithBackend(options, dependencies = {}) { + const algorithm = requireLocalizationAlgorithmDefinition( + options?.experiment?.beamforming?.algorithm, + ); + if (!isDasAlgorithm(algorithm.id)) { + throw new RangeError( + `computeDasWithBackend does not implement ${algorithm.id}.`, + ); + } + return computeDasBackend(options, dependencies); +} + +/** + * Execute any registered localization algorithm. GCC-PHAT currently uses its + * portable CPU implementation for both `auto` and `cpu`; an explicit WebGPU + * request is rejected instead of silently changing the requested backend. + */ +export async function computeLocalizationWithBackend(options, dependencies = {}) { + const algorithm = requireLocalizationAlgorithmDefinition( + options?.experiment?.beamforming?.algorithm, + ); + if (isDasAlgorithm(algorithm.id)) { + return computeDasBackend(options, dependencies); + } + + const requestedBackend = options?.experiment?.beamforming?.backend ?? "auto"; + if (requestedBackend === "webgpu") { + throw new TypeError( + `${algorithm.label} currently supports only the CPU backend; choose auto or cpu.`, + ); + } + if (requestedBackend !== "auto" && requestedBackend !== "cpu") { + throw new TypeError(`unsupported localization backend: ${requestedBackend}`); + } + const computeCpu = + dependencies.computeGccCpu ?? + dependencies.computeCpu ?? + computeGccPhatMap; + return computeCpu(options); +} diff --git a/packages/webgpu-kernels/das.js b/packages/webgpu-kernels/das.js index 520fa86..48b9110 100644 --- a/packages/webgpu-kernels/das.js +++ b/packages/webgpu-kernels/das.js @@ -9,6 +9,11 @@ import { requireMaxCompensationDb, } from "../acoustics-core/frequency-response.js"; import { dbToLinear } from "../acoustics-core/math.js"; +import { + DEFAULT_PHAT_RELATIVE_FLOOR, + isDasAlgorithm, + requireLocalizationAlgorithmDefinition, +} from "../acoustics-core/localization-algorithms.js"; import { extractStftBins } from "../acoustics-core/stft.js"; import { DAS_PARAMS_BYTE_LENGTH, @@ -72,11 +77,60 @@ function minimumBufferSize(byteLength) { return Math.max(4, Math.ceil(byteLength / 4) * 4); } +function requirePhatRelativeFloor(value) { + if (!Number.isFinite(value) || value < 0 || value > 1) { + throw new RangeError("phatRelativeFloor must lie inside 0..1."); + } + return value; +} + +function phaseNormalizePackedStft( + packed, + binCount, + frameCount, + microphoneCount, + relativeFloor, +) { + requirePhatRelativeFloor(relativeFloor); + const valuesPerBin = frameCount * microphoneCount; + for (let binIndex = 0; binIndex < binCount; binIndex += 1) { + const binOffset = binIndex * valuesPerBin * 2; + for (let frameIndex = 0; frameIndex < frameCount; frameIndex += 1) { + const frameOffset = binOffset + frameIndex * microphoneCount * 2; + let maximumMagnitude = 0; + for (let microphoneIndex = 0; microphoneIndex < microphoneCount; microphoneIndex += 1) { + const valueOffset = frameOffset + microphoneIndex * 2; + maximumMagnitude = Math.max( + maximumMagnitude, + Math.hypot(packed[valueOffset], packed[valueOffset + 1]), + ); + } + const magnitudeFloor = maximumMagnitude * relativeFloor; + for (let microphoneIndex = 0; microphoneIndex < microphoneCount; microphoneIndex += 1) { + const valueOffset = frameOffset + microphoneIndex * 2; + const magnitude = Math.hypot(packed[valueOffset], packed[valueOffset + 1]); + if (!(magnitude > magnitudeFloor)) { + packed[valueOffset] = 0; + packed[valueOffset + 1] = 0; + continue; + } + packed[valueOffset] /= magnitude; + packed[valueOffset + 1] /= magnitude; + } + } + } +} + export function packStftBatch(stft, bins, options = {}) { if (!stft || !Number.isInteger(stft.frameCount) || !Number.isInteger(stft.microphoneCount)) { throw new TypeError("유효한 STFT 결과가 필요합니다."); } if (!Array.isArray(bins) || bins.length === 0) throw new TypeError("하나 이상의 FFT bin이 필요합니다."); + const algorithm = requireLocalizationAlgorithmDefinition(options.algorithm ?? "das"); + if (!isDasAlgorithm(algorithm.id)) { + throw new RangeError(`packStftBatch does not implement ${algorithm.id}.`); + } + const phaseOnly = algorithm.spectralWeighting === "phat"; const valuesPerBin = stft.frameCount * stft.microphoneCount; const packed = new Float32Array(bins.length * valuesPerBin * 2); const microphones = options.microphones; @@ -105,9 +159,9 @@ export function packStftBatch(stft, bins, options = {}) { microphones.forEach((microphone, microphoneIndex) => { if (!Array.isArray(microphone.frequencyResponse)) return; interpolateFrequencyResponse(microphone.frequencyResponse, frequencyHz, response); - const inverseGain = dbToLinear( - inverseResponseGainDb(response.gainDb, maxCompensationDb), - ); + const inverseGain = phaseOnly + ? 1 + : dbToLinear(inverseResponseGainDb(response.gainDb, maxCompensationDb)); const responseOffset = (binIndex * stft.microphoneCount + microphoneIndex) * 2; inverseResponses[responseOffset] = inverseGain * Math.cos(response.phaseRadians); inverseResponses[responseOffset + 1] = -inverseGain * Math.sin(response.phaseRadians); @@ -138,6 +192,15 @@ export function packStftBatch(stft, bins, options = {}) { inputReal * inverseImaginary + inputImaginary * inverseReal; } }); + if (phaseOnly) { + phaseNormalizePackedStft( + packed, + bins.length, + stft.frameCount, + stft.microphoneCount, + options.phatRelativeFloor ?? DEFAULT_PHAT_RELATIVE_FLOOR, + ); + } return packed; } @@ -445,9 +508,16 @@ export async function computeDasWebGpu(options) { maximumGpuComplexMacs = 20_000_000_000, maximumFftButterflies = 1_000_000_000, stftTransformBatchSize = 32, + phatRelativeFloor = DEFAULT_PHAT_RELATIVE_FLOOR, gpu, powerPreference, } = options; + const algorithm = requireLocalizationAlgorithmDefinition( + experiment.beamforming?.algorithm ?? "das", + ); + if (!isDasAlgorithm(algorithm.id)) { + throw new RangeError(`computeDasWebGpu does not implement ${algorithm.id}.`); + } if (!Array.isArray(channels) || channels.length !== arrayDefinition?.microphones?.length) { throw new Error(`채널 수 ${channels?.length ?? 0}와 마이크 수 ${arrayDefinition?.microphones?.length ?? 0}가 다릅니다.`); } @@ -501,7 +571,9 @@ export async function computeDasWebGpu(options) { label: "DAS beamforming pipeline", }), })); - const calibration = packMicrophoneCalibration(arrayDefinition.microphones); + const calibration = packMicrophoneCalibration(arrayDefinition.microphones, { + includeAmplitude: algorithm.spectralWeighting !== "phat", + }); const calibrationBuffer = await withWebGpuErrorScopes(context.device, async () => createBuffer(context.device, { label: "DAS microphone calibration", @@ -543,12 +615,14 @@ export async function computeDasWebGpu(options) { }), }); const packedSpectrum = packStftBatch(stft, batchBins, { + algorithm: algorithm.id, microphones: arrayDefinition.microphones, sampleRate, fftSize: experiment.stft.fftSize, compensateFrequencyResponse: experiment.beamforming.compensateFrequencyResponse !== false, maxCompensationDb: experiment.beamforming.maxCompensationDb, + phatRelativeFloor, }); stftMs += now() - stftStartedAt; const binIndices = Uint32Array.from(batchBins); @@ -688,9 +762,10 @@ export async function computeDasWebGpu(options) { completedBins += batchBins.length; } - const bands = finalizeDasBands(rawMaps, workload, experiment); + const bands = finalizeDasBands(rawMaps, workload, experiment, algorithm.rawScale); const totalMs = now() - totalStartedAt; return { + algorithm: algorithm.id, backend: "webgpu", frameCount: workload.frameCount, microphoneCount: workload.microphoneCount, diff --git a/packages/webgpu-kernels/das.wgsl.js b/packages/webgpu-kernels/das.wgsl.js index aa6ce87..1d11419 100644 --- a/packages/webgpu-kernels/das.wgsl.js +++ b/packages/webgpu-kernels/das.wgsl.js @@ -12,7 +12,7 @@ * Flat buffer layouts: * spectrum: [bin][frame][microphone][real, imaginary] * distances: [point][microphone] - * calibration: [microphone][inverse sensitivity gain, phase offset radians] + * calibration: [microphone][inverse sensitivity gain or unity, phase offset radians] * steering: [bin][point][microphone][real, imaginary] * output: [bin][point] */ @@ -197,9 +197,13 @@ export function packDasCalibration(calibration) { return packed; } -/** Convert array-definition microphone metadata to the shader calibration layout. */ -export function packMicrophoneCalibration(microphones) { +/** + * Convert array-definition microphone metadata to the shader calibration layout. + * Phase-only algorithms set includeAmplitude=false while retaining phase offsets. + */ +export function packMicrophoneCalibration(microphones, options = {}) { if (!Array.isArray(microphones)) throw new TypeError("microphones must be an array."); + const includeAmplitude = options.includeAmplitude !== false; const packed = new Float32Array(microphones.length * 2); for (let index = 0; index < microphones.length; index += 1) { const microphone = microphones[index]; @@ -213,7 +217,7 @@ export function packMicrophoneCalibration(microphones) { ); packed[index * 2] = requireFiniteF32( `microphones[${index}] inverse sensitivity gain`, - 10 ** (-sensitivityDb / 20), + includeAmplitude ? 10 ** (-sensitivityDb / 20) : 1, ); packed[index * 2 + 1] = requireFiniteF32( `microphones[${index}] phase offset`, diff --git a/schemas/experiment.schema.json b/schemas/experiment.schema.json index 8959ec6..4e79a00 100644 --- a/schemas/experiment.schema.json +++ b/schemas/experiment.schema.json @@ -42,7 +42,8 @@ "type": "object", "additionalProperties": false, "required": ["algorithm", "backend", "minimumDistance", "normalization", "dynamicRangeDb"], "properties": { - "algorithm": { "const": "das" }, "backend": { "enum": ["auto", "cpu", "webgpu"] }, + "algorithm": { "enum": ["das", "das-phat", "gcc-phat"] }, + "backend": { "enum": ["auto", "cpu", "webgpu"] }, "minimumDistance": { "type": "number", "exclusiveMinimum": 0, "maximum": 10 }, "normalization": { "enum": ["peak", "none"] }, "dynamicRangeDb": { "type": "number", "minimum": 6, "maximum": 120 }, diff --git a/schemas/result.schema.json b/schemas/result.schema.json index 96eff58..32b318c 100644 --- a/schemas/result.schema.json +++ b/schemas/result.schema.json @@ -15,6 +15,7 @@ "arrayId": { "type": "string", "minLength": 1 }, "sceneId": { "type": "string", "minLength": 1 }, "createdAt": { "type": "string", "format": "date-time" }, + "algorithm": { "enum": ["das", "das-phat", "gcc-phat"] }, "backend": { "enum": ["cpu", "webgpu"] }, "inputHash": { "type": "string", "pattern": "^[a-f0-9]{64}$" }, "signalStats": { @@ -43,6 +44,7 @@ "properties": { "total": { "type": "number", "minimum": 0 }, "signalGeneration": { "type": "number", "minimum": 0 }, + "localization": { "type": "number", "minimum": 0 }, "beamforming": { "type": "number", "minimum": 0 } } }, @@ -64,7 +66,7 @@ "properties": { "rows": { "type": "integer", "minimum": 1, "maximum": 512 }, "columns": { "type": "integer", "minimum": 1, "maximum": 512 }, - "scale": { "enum": ["power", "normalized", "db"] }, + "scale": { "enum": ["power", "score", "normalized", "db"] }, "values": { "type": "array", "minItems": 1, "maxItems": 262144, "items": { "type": "number" } @@ -91,7 +93,12 @@ ], "properties": { "calibrated": { "const": false }, - "method": { "const": "connected-half-power-power-weighted-covariance" }, + "method": { + "enum": [ + "connected-half-power-power-weighted-covariance", + "connected-half-maximum-score-weighted-covariance" + ] + }, "halfPowerRatio": { "type": "number", "exclusiveMinimum": 0, "maximum": 1 }, "thresholdValue": { "type": "number", "minimum": 0 }, "pixelCount": { "type": "integer", "minimum": 0, "maximum": 262144 }, diff --git a/tests/test-runner.js b/tests/test-runner.js index 8083379..4c38dc2 100644 --- a/tests/test-runner.js +++ b/tests/test-runner.js @@ -11,6 +11,8 @@ 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/gcc-phat.test.js"); + await import("/tests/unit/localization.test.js"); await import("/tests/unit/realtime-localizer.test.js"); await import("/tests/unit/builder.test.js"); await import("/tests/unit/evaluation.test.js"); diff --git a/tests/unit/acoustics.test.js b/tests/unit/acoustics.test.js index 272a596..84b64f0 100644 --- a/tests/unit/acoustics.test.js +++ b/tests/unit/acoustics.test.js @@ -12,6 +12,13 @@ import { prepareDelaySpectrum, } from "/packages/acoustics-core/fractional-delay.js"; import { dbToLinear, distance3 } from "/packages/acoustics-core/math.js"; +import { + DEFAULT_LOCALIZATION_ALGORITHM, + LOCALIZATION_ALGORITHM_IDS, + isDasAlgorithm, + localizationAlgorithmDefinition, + localizationAlgorithmLabel, +} from "/packages/acoustics-core/localization-algorithms.js"; import { renderSceneMix } from "/packages/acoustics-core/mixdown.js"; import { whiteNoiseAt } from "/packages/acoustics-core/noise.js"; import { sphericalAttenuation, synthesizeMicrophoneSignals } from "/packages/acoustics-core/propagation.js"; @@ -29,6 +36,115 @@ function assertNear(actual, expected, tolerance, label) { assert(Math.abs(actual - expected) <= tolerance, `${label}: expected ${expected} ± ${tolerance}, received ${actual}`); } +test("localization algorithm registry exposes stable ids and capabilities", () => { + assertEqual(DEFAULT_LOCALIZATION_ALGORITHM, "das", "default algorithm"); + assertEqual( + LOCALIZATION_ALGORITHM_IDS.join(","), + "das,das-phat,gcc-phat", + "stable algorithm ids", + ); + assertEqual(localizationAlgorithmDefinition("das").rawScale, "power", "DAS map scale"); + assertEqual( + localizationAlgorithmDefinition("das-phat").spectralWeighting, + "phat", + "DAS-PHAT weighting", + ); + assertEqual( + localizationAlgorithmDefinition("gcc-phat").backends.join(","), + "cpu", + "GCC-PHAT backend capability", + ); + assert(isDasAlgorithm("das") && isDasAlgorithm("das-phat"), "DAS engine algorithms"); + assert(!isDasAlgorithm("gcc-phat"), "GCC-PHAT uses a separate engine"); + assertEqual(localizationAlgorithmLabel("das-phat"), "DAS-PHAT", "algorithm label"); + assertEqual(localizationAlgorithmDefinition("unknown"), undefined, "unknown lookup"); +}); + +test("CPU DAS-PHAT is invariant to spectrum and calibration amplitude gains", async () => { + const sampleRate = 8_000; + const fftSize = 256; + const bin = 32; + const frequency = (bin * sampleRate) / fftSize; + const speedOfSound = 343; + const target = [0.2, 0, 1]; + const positions = [[-0.15, 0, 0], [-0.02, 0, 0], [0.12, 0, 0]]; + const phaseOffsets = [30, -45, 75]; + const responsePhases = [10, 20, -35]; + const createDefinition = (withAmplitudeMetadata) => ({ + microphones: positions.map((position, index) => ({ + id: `mic-${index}`, + position, + sensitivityDb: withAmplitudeMetadata ? [-30, 24, 60][index] : 0, + phaseOffsetDeg: phaseOffsets[index], + frequencyResponse: [ + { + frequencyHz: 0, + gainDb: withAmplitudeMetadata ? [-60, 40, 100][index] : 0, + phaseDeg: responsePhases[index], + }, + { + frequencyHz: 4_000, + gainDb: withAmplitudeMetadata ? [-60, 40, 100][index] : 0, + phaseDeg: responsePhases[index], + }, + ], + })), + }); + const createSpectrum = (magnitudes) => { + const real = new Float32Array(positions.length); + const imaginary = new Float32Array(positions.length); + positions.forEach((position, index) => { + const propagationPhase = + (2 * Math.PI * frequency * distance3(target, position)) / speedOfSound; + const calibrationPhase = + ((phaseOffsets[index] + responsePhases[index]) * Math.PI) / 180; + const observedPhase = -propagationPhase + calibrationPhase; + real[index] = Math.cos(observedPhase) * magnitudes[index]; + imaginary[index] = Math.sin(observedPhase) * magnitudes[index]; + }); + return { real, imaginary }; + }; + const scan = { + z: 1, + xMin: -0.2, + xMax: 0.2, + yMin: 0, + yMax: 0, + columns: 2, + rows: 1, + }; + const baselineArray = createDefinition(false); + const gainedArray = createDefinition(true); + const common = { + bin, + frameCount: 1, + pointCount: 2, + sampleRate, + fftSize, + speedOfSound, + algorithm: "das-phat", + pointBatchSize: 64, + }; + const baseline = await computeBinMap({ + ...common, + spectrum: createSpectrum([1, 1, 1]), + arrayDefinition: baselineArray, + distances: createDistanceMatrix(baselineArray, scan), + }); + const gained = await computeBinMap({ + ...common, + spectrum: createSpectrum([0.125, 8, 32]), + arrayDefinition: gainedArray, + distances: createDistanceMatrix(gainedArray, scan), + }); + + baseline.forEach((value, index) => { + assertNear(gained[index], value, 2e-6, `gain-invariant map value ${index}`); + }); + assertNear(baseline[1], 1, 2e-6, "phase-compensated target score"); + assert(baseline[1] > baseline[0] * 1.1, "phase compensation preserves target peak"); +}); + function rms(values, start = 0, end = values.length) { let sum = 0; for (let index = start; index < end; index += 1) sum += values[index] ** 2; @@ -521,12 +637,31 @@ test("CPU DAS peak가 calibration 후 단일 tone의 실제 grid 위치를 찾 onProgress: (progress) => progressEvents.push(progress), }); const peakPosition = result.bands[0].peak.position; + assertEqual(result.algorithm, "das", "default result algorithm"); + assertEqual(result.bands[0].rawScale, "power", "DAS raw map scale"); assert(distance3(peakPosition, scene.sources[0].position) <= 0.071, `peak=${peakPosition.join(",")}`); assert(result.bands[0].rawMap !== result.bands[0].displayMap, "raw/display buffers are separate"); assertNear(Math.max(...result.bands[0].rawMap), result.bands[0].peak.rawPower, 1e-4, "raw map peak"); assertNear(Math.max(...result.bands[0].displayMap), 1, 1e-6, "normalized display peak"); assert(progressEvents.some(({ stage }) => stage === "beamforming"), "beamforming progress"); + const phaseOnly = await computeDasMap({ + channels: synthesis.channels, + sampleRate: synthesis.sampleRate, + arrayDefinition, + scene, + experiment: { + ...experiment, + beamforming: { ...experiment.beamforming, algorithm: "das-phat" }, + }, + }); + assertEqual(phaseOnly.algorithm, "das-phat", "phase-only result algorithm"); + assertEqual(phaseOnly.bands[0].rawScale, "score", "DAS-PHAT raw map scale"); + assert( + distance3(phaseOnly.bands[0].peak.position, scene.sources[0].position) <= 0.071, + `DAS-PHAT peak=${phaseOnly.bands[0].peak.position.join(",")}`, + ); + const uncompensated = await computeDasMap({ channels: synthesis.channels, sampleRate: synthesis.sampleRate, diff --git a/tests/unit/contracts.test.js b/tests/unit/contracts.test.js index 325016c..641d346 100644 --- a/tests/unit/contracts.test.js +++ b/tests/unit/contracts.test.js @@ -157,6 +157,51 @@ test("단일 tone 장면과 실험이 관계 검증을 통과한다", async () = assertEqual((await validateExperiment(experiment, { arrayDefinition, scene })).length, 0, "experiment issues"); }); +test("세 위치 추정 알고리즘과 backend 제약을 계약이 검증한다", async () => { + const [arrayDefinition, scene, baseline] = await Promise.all([ + loadJson("/assets/array-geometries/rectangular-16x12.array.json"), + loadJson("/assets/sample-scenes/single-tone.scene.json"), + loadJson("/assets/sample-scenes/single-tone.experiment.json"), + ]); + + for (const algorithm of ["das", "das-phat", "gcc-phat"]) { + const experiment = structuredClone(baseline); + experiment.beamforming.algorithm = algorithm; + experiment.beamforming.backend = algorithm === "gcc-phat" ? "cpu" : "auto"; + assertEqual( + (await validateExperiment(experiment, { arrayDefinition, scene })).length, + 0, + `${algorithm} experiment issues`, + ); + } + + const incompatible = structuredClone(baseline); + incompatible.beamforming.algorithm = "gcc-phat"; + incompatible.beamforming.backend = "webgpu"; + let issues = await validateExperiment(incompatible, { arrayDefinition, scene }); + assert( + issues.some(({ path }) => path === "$.beamforming.backend"), + "GCC-PHAT rejects the WebGPU backend", + ); + + incompatible.beamforming.backend = "cpu"; + issues = await validateExperiment(incompatible, { + arrayDefinition: { ...arrayDefinition, microphones: arrayDefinition.microphones.slice(0, 1) }, + scene, + }); + assert( + issues.some(({ path }) => path === "$.beamforming.algorithm"), + "GCC-PHAT requires two microphones", + ); + + incompatible.beamforming.algorithm = "music"; + issues = await validateExperiment(incompatible, { arrayDefinition, scene }); + assert( + issues.some(({ path }) => path === "$.beamforming.algorithm"), + "unsupported algorithms are rejected", + ); +}); + test("shoebox 음향 환경의 room과 occluder 경계를 검증한다", async () => { const scene = await loadJson("/assets/sample-scenes/single-tone.scene.json"); scene.acoustics = { @@ -357,6 +402,7 @@ test("복수 음원 대응 result map 계약과 shape 관계를 검증한다", a sceneId: experiment.sceneId, createdAt: "2026-08-28T00:00:00.000Z", backend: "cpu", + algorithm: "das", inputHash: "a".repeat(64), signalStats: { preNormalizationPeak: 0.5, appliedScale: 1, clippedSampleCount: 0 }, timingMs: { total: 10, signalGeneration: 3, beamforming: 7 }, @@ -369,10 +415,25 @@ test("복수 음원 대응 result map 계약과 shape 관계를 검증한다", a }], }; assertEqual((await validateResult(result, { experiment, scene })).length, 0, "valid result issues"); + + experiment.beamforming.algorithm = "das-phat"; + result.algorithm = "das-phat"; + result.bands[0].rawMap.scale = "score"; + assertEqual( + (await validateResult(result, { experiment, scene })).length, + 0, + "phase-weighted score result issues", + ); + + result.algorithm = "das"; + let issues = await validateResult(result, { experiment, scene }); + assert(issues.some(({ path }) => path === "$.algorithm"), "result algorithm must match its experiment"); + + experiment.beamforming.algorithm = "das"; result.bands[0].rawMap.values.pop(); result.bands[0].rawMap.scale = "normalized"; result.bands[0].peaks[0].sourceId = "missing-source"; - const issues = await validateResult(result, { experiment, scene }); + issues = await validateResult(result, { experiment, scene }); assert(issues.some(({ message }) => message.includes("rows × columns"))); assert(issues.some(({ path }) => path.endsWith("rawMap.scale"))); assert(issues.some(({ path }) => path.endsWith("sourceId"))); diff --git a/tests/unit/evaluation.test.js b/tests/unit/evaluation.test.js index eb65636..71444f6 100644 --- a/tests/unit/evaluation.test.js +++ b/tests/unit/evaluation.test.js @@ -1,6 +1,7 @@ import { RESULT_METRIC_KEYS, evaluateDasBand, + evaluateLocalizationBand, exportBandEvaluation, findLocalPeaks, measureBeamSpread, @@ -213,6 +214,44 @@ test("beam spread가 회전된 비등방성 half-power lobe의 covariance ellips ); }); +test("generic localization 평가는 score map의 척도와 beam spread 방식을 보존한다", () => { + const scan = { xMin: -1, xMax: 1, yMin: -1, yMax: 1, z: 2, columns: 3, rows: 3 }; + const rawMap = { + rows: 3, + columns: 3, + scale: "score", + values: Float32Array.of( + 0, 1, 0, + 1, 4, 1, + 0, 1, 0, + ), + }; + const evaluation = evaluateLocalizationBand({ + rawMap, + scan, + sources: [{ id: "center", position: [0, 0, 2] }], + matchToleranceM: 0, + truthRadiusM: 0, + minimumPeakDistanceM: 0, + }); + + assertEqual(evaluation.diagnostics.mapScale, "score", "score scale diagnostic"); + assertEqual( + evaluation.detections[0].beamSpread.method, + "connected-half-maximum-score-weighted-covariance", + "score beam-spread method", + ); + assertEqual(evaluation.sourceMatches[0].localizationErrorM, 0, "score peak localization"); + + const legacy = evaluateDasBand({ + rawMap: rawMap.values, + scan, + sources: [], + minimumPeakDistanceM: 0, + }); + assertEqual(legacy.diagnostics.mapScale, "power", "legacy DAS defaults to power"); +}); + test("beam spread가 1-row, 1-column, single-pixel lobe에서 유한하고 결정적이다", () => { const powers = Float32Array.of(5, 8, 10, 8, 5); const rowSpread = measureBeamSpread( @@ -337,6 +376,6 @@ test("evaluation이 non-finite, negative, shape 오류를 입력 경계에서 scan, sources: [], })) instanceof RangeError, - "non-power map rejected", + "display-only map scale rejected", ); }); diff --git a/tests/unit/gcc-phat.test.js b/tests/unit/gcc-phat.test.js new file mode 100644 index 0000000..ffba682 --- /dev/null +++ b/tests/unit/gcc-phat.test.js @@ -0,0 +1,314 @@ +import { + computeGccPhatMap, + createGccPhatAnalysisPlan, + estimateGccPhatWorkload, + selectGccPhatReference, +} from "/packages/acoustics-core/gcc-phat.js"; +import { + addFrequencyDomainDelay, + createDelayWorkspace, + prepareDelaySpectrum, +} from "/packages/acoustics-core/fractional-delay.js"; +import { distance3 } from "/packages/acoustics-core/math.js"; + +function assertNear(actual, expected, tolerance, label) { + assert( + Math.abs(actual - expected) <= tolerance, + label + ": expected " + expected + " ± " + tolerance + ", received " + actual, + ); +} + +function deterministicBroadbandSamples(length) { + let state = 0x6d2b79f5; + return Float32Array.from({ length }, (_, index) => { + state ^= state << 13; + state ^= state >>> 17; + state ^= state << 5; + const noise = ((state >>> 0) / 0xffff_ffff) * 2 - 1; + const taper = Math.min(1, index / 64, (length - 1 - index) / 64); + return noise * Math.max(0, taper); + }); +} + +function createBroadbandFixture() { + const sampleRate = 4096; + const speedOfSound = 343; + const target = [0.3, 0.2, 0.6]; + const arrayDefinition = { + id: "gcc-array", + microphones: [ + { + id: "a-ref", + position: [-0.45, -0.35, 0], + sensitivityDb: 9, + phaseOffsetDeg: 31, + frequencyResponse: [ + { frequencyHz: 0, gainDb: -8, phaseDeg: 19 }, + { frequencyHz: sampleRate / 2, gainDb: -8, phaseDeg: 19 }, + ], + }, + { + id: "b-ne", + position: [0.45, 0.35, 0], + sensitivityDb: -7, + phaseOffsetDeg: -47, + frequencyResponse: [ + { frequencyHz: 0, gainDb: 6, phaseDeg: -23 }, + { frequencyHz: sampleRate / 2, gainDb: 6, phaseDeg: -23 }, + ], + }, + { + id: "c-nw", + position: [-0.45, 0.35, 0], + sensitivityDb: 4, + phaseOffsetDeg: 73, + }, + { + id: "d-se", + position: [0.45, -0.35, 0], + sensitivityDb: -3, + phaseOffsetDeg: -11, + }, + ], + }; + const scene = { speedOfSound }; + const experiment = { + stft: { fftSize: 256, hopSize: 128, window: "hann" }, + scan: { + z: target[2], + xMin: -0.4, + xMax: 0.4, + yMin: -0.4, + yMax: 0.4, + columns: 9, + rows: 9, + }, + frequencyBands: [ + { id: "broadband", minHz: 128, maxHz: 1800, color: "#65f2bf" }, + ], + beamforming: { + algorithm: "gcc-phat", + backend: "cpu", + normalization: "peak", + compensateFrequencyResponse: true, + maxCompensationDb: 24, + }, + }; + + const emitted = deterministicBroadbandSamples(2048); + const distances = arrayDefinition.microphones.map((microphone) => + distance3(target, microphone.position) + ); + const maximumDelaySamples = Math.max(...distances) * sampleRate / speedOfSound; + const spectrum = prepareDelaySpectrum(emitted, maximumDelaySamples); + const workspace = createDelayWorkspace(spectrum.fftSize); + const channels = arrayDefinition.microphones.map((microphone, index) => { + const channel = new Float32Array(emitted.length); + addFrequencyDomainDelay({ + output: channel, + spectrum, + delaySamples: distances[index] * sampleRate / speedOfSound, + gain: 1 / distances[index], + phaseOffsetDeg: microphone.phaseOffsetDeg, + frequencyResponse: microphone.frequencyResponse, + sampleRate, + workspace, + }); + return channel; + }); + + return { + sampleRate, + speedOfSound, + target, + arrayDefinition, + scene, + experiment, + channels, + }; +} + +function scanIndexFor(scan, x, y) { + const column = Math.round((x - scan.xMin) / (scan.xMax - scan.xMin) * (scan.columns - 1)); + const row = Math.round((scan.yMax - y) / (scan.yMax - scan.yMin) * (scan.rows - 1)); + return row * scan.columns + column; +} + +function assertFloat32Exact(actual, expected, label) { + assertEqual(actual.length, expected.length, label + " length"); + for (let index = 0; index < actual.length; index += 1) { + assertEqual(actual[index], expected[index], label + "[" + index + "]"); + } +} + +test("GCC-PHAT reference selection uses the farthest centroid radius and id tie-break", () => { + const arrayDefinition = { + microphones: [ + { id: "z-left", position: [-1, 0, 0] }, + { id: "middle", position: [0, 0, 0] }, + { id: "a-right", position: [1, 0, 0] }, + ], + }; + const selected = selectGccPhatReference(arrayDefinition); + assertEqual(selected.id, "a-right", "lexicographically first id wins a radius tie"); + assertEqual(selected.index, 2, "selected microphone index"); + assertEqual(selected.microphone, arrayDefinition.microphones[2], "selected microphone identity"); +}); + +test("GCC-PHAT workload reports reference-pair scan reduction and rejects aliased lag", () => { + const fixture = createBroadbandFixture(); + const workload = estimateGccPhatWorkload(fixture); + assertEqual(workload.pairCount, fixture.arrayDefinition.microphones.length - 1, "reference pair count"); + assertEqual( + workload.scanLookupCount, + workload.pointCount * workload.pairCount * workload.bandCount, + "one correlation lookup per point, pair, and band", + ); + assert( + workload.scanLookupCount < workload.dasEquivalentComplexMacCount, + "reference-pair lookup count is below equivalent DAS complex MAC count", + ); + assert( + workload.maximumLagSamples > 0 && + workload.maximumLagSamples <= fixture.experiment.stft.fftSize / 2, + "maximum geometric lag is exposed and unaliased", + ); + + let memoryError; + try { + createGccPhatAnalysisPlan({ ...fixture, maximumPlanBytes: 1 }); + } catch (error) { + memoryError = error; + } + assert( + memoryError instanceof RangeError && memoryError.message.includes("maximumPlanBytes"), + "plan hard memory bound", + ); + + const aliasedArray = { + microphones: [ + { id: "far", position: [-1, 0, 0] }, + { id: "near", position: [1, 0, 0] }, + ], + }; + let lagError; + try { + createGccPhatAnalysisPlan({ + sampleRate: 48000, + arrayDefinition: aliasedArray, + scene: { speedOfSound: 343 }, + experiment: { + stft: { fftSize: 128, hopSize: 64, window: "hann" }, + scan: { + z: 0.1, + xMin: 0.9, + xMax: 0.9, + yMin: 0, + yMax: 0, + columns: 1, + rows: 1, + }, + frequencyBands: [{ id: "wide", minHz: 500, maxHz: 5000, color: "#65f2bf" }], + beamforming: { algorithm: "gcc-phat", backend: "cpu", normalization: "peak" }, + }, + }); + } catch (error) { + lagError = error; + } + assert( + lagError instanceof RangeError && + lagError.message.includes("maximum absolute geometric lag") && + lagError.message.includes("fftSize/2"), + "circular-lag alias guard", + ); +}); + +test("GCC-PHAT corrected lag sign finds a deterministic broadband source and emits score maps", async () => { + const fixture = createBroadbandFixture(); + const result = await computeGccPhatMap({ + ...fixture, + phatFloorRatio: 1e-5, + pointBatchSize: 1024, + stftTransformBatchSize: 1024, + }); + const band = result.bands[0]; + const peakPosition = band.peak.position; + assertEqual(result.backend, "cpu", "CPU backend"); + assertEqual(result.algorithm, "gcc-phat", "algorithm id"); + assertEqual(result.referenceMicrophoneId, "a-ref", "deterministic reference metadata"); + assertEqual(result.pairCount, 3, "pair metadata"); + assertEqual(result.phatFloorRatio, 1e-5, "PHAT floor metadata"); + assertNear(result.maximumLagSamples, result.workload.maximumLagSamples, 0, "lag metadata"); + assertEqual(band.rawScale, "score", "raw map scale"); + assert( + band.rawMap.every((value) => Number.isFinite(value) && value >= 0), + "shifted and squared score map is finite and nonnegative", + ); + assert( + distance3(peakPosition, fixture.target) <= 0.15, + "broadband GCC-PHAT peak is within one grid diagonal: " + peakPosition.join(","), + ); + const targetIndex = scanIndexFor( + fixture.experiment.scan, + fixture.target[0], + fixture.target[1], + ); + const wrongSignIndex = scanIndexFor( + fixture.experiment.scan, + -fixture.target[0], + fixture.target[1], + ); + assert( + band.rawMap[targetIndex] > band.rawMap[wrongSignIndex], + "X_m * conj(X_ref) correlation uses the physical arrival-lag sign", + ); + assertNear(Math.max(...band.displayMap), 1, 1e-6, "normalized score display"); +}); + +test("GCC-PHAT analysis plan reuses geometry and FFT workspaces and releases its busy lock", async () => { + const fixture = createBroadbandFixture(); + const plan = createGccPhatAnalysisPlan(fixture); + const retainedDistances = plan.distances; + const retainedStftWorkspace = plan.stftWorkspace; + const retainedCorrelationReal = plan.correlationReal; + const computeOptions = { + ...fixture, + analysisPlan: plan, + pointBatchSize: 1024, + stftTransformBatchSize: 1024, + }; + const first = await computeGccPhatMap(computeOptions); + const second = await computeGccPhatMap(computeOptions); + + assertEqual(plan.distances, retainedDistances, "distance matrix identity is retained"); + assertEqual(plan.stftWorkspace, retainedStftWorkspace, "STFT workspace identity is retained"); + assertEqual(plan.correlationReal, retainedCorrelationReal, "correlation workspace identity is retained"); + assertFloat32Exact(second.bands[0].rawMap, first.bands[0].rawMap, "reused plan raw map"); + assertEqual(second.analysisCache.analysisAttempts, 2, "analysis attempts"); + assertEqual(second.analysisCache.completedAnalyses, 2, "completed analyses"); + assertEqual(second.analysisCache.distanceMatrixBuilds, 1, "distance matrix built once"); + assertEqual( + second.analysisCache.avoidedDistanceEvaluations, + plan.pointCount * plan.microphoneCount, + "second geometry pass avoided", + ); + assertEqual(second.analysisCache.stftWorkspaceUses, 2, "STFT workspace reused"); + assertEqual(second.analysisCache.correlationWorkspaceUses, 2, "correlation workspace reused"); + assertEqual( + second.analysisCache.inverseFftTransforms, + plan.pairCount * plan.bandBins.length * 2, + "one inverse FFT per pair and band per analysis", + ); + + const controller = new AbortController(); + controller.abort(); + let abortError; + try { + await computeGccPhatMap({ ...computeOptions, signal: controller.signal }); + } catch (error) { + abortError = error; + } + assertEqual(abortError?.name, "AbortError", "AbortError propagated"); + assertEqual(plan.busy, false, "plan busy lock released after cancellation"); + assertEqual(plan.statistics.completedAnalyses, 2, "canceled analysis is not completed"); +}); diff --git a/tests/unit/localization.test.js b/tests/unit/localization.test.js new file mode 100644 index 0000000..660487a --- /dev/null +++ b/tests/unit/localization.test.js @@ -0,0 +1,198 @@ +import { + computeLocalizationMap, + createLocalizationAnalysisPlan, + estimateLocalizationWorkload, +} from "/packages/acoustics-core/localization.js"; +import { + DEFAULT_LOCALIZATION_ALGORITHM, + LOCALIZATION_ALGORITHM_IDS, + localizationAlgorithmDefinition, + localizationAlgorithmSupportsBackend, +} from "/packages/acoustics-core/localization-algorithms.js"; +import { + StreamingDasLocalizer, + StreamingLocalizer, + createLocalizationPlan, + detectBandPeaks, + localizeMultichannel, +} from "/packages/acoustics-core/realtime-localizer.js"; +import { + computeDasWithBackend, + computeLocalizationWithBackend, +} from "/packages/webgpu-kernels/compute.js"; + +function createFixture(algorithm = "das") { + const sampleRate = 4096; + const arrayDefinition = { + id: "algorithm-fixture", + microphones: [ + { id: "nw", position: [-0.1, 0.1, 0] }, + { id: "ne", position: [0.1, 0.1, 0] }, + { id: "sw", position: [-0.1, -0.1, 0] }, + { id: "se", position: [0.1, -0.1, 0] }, + ], + }; + let state = 0x9e3779b9; + const emitted = Float32Array.from({ length: 128 }, () => { + state ^= state << 13; + state ^= state >>> 17; + state ^= state << 5; + return ((state >>> 0) / 0xffff_ffff) * 2 - 1; + }); + const channels = arrayDefinition.microphones.map((_, microphoneIndex) => + Float32Array.from( + { length: emitted.length }, + (_, sampleIndex) => emitted[Math.max(0, sampleIndex - microphoneIndex)], + ) + ); + const scan = { + z: 1, + xMin: -0.2, + xMax: 0.2, + yMin: -0.2, + yMax: 0.2, + columns: 3, + rows: 3, + }; + const frequencyBands = [ + { id: "wide", minHz: 128, maxHz: 1800, color: "#65f2bf" }, + ]; + const experiment = { + stft: { fftSize: 64, hopSize: 32, window: "hann" }, + scan, + frequencyBands, + beamforming: { + algorithm, + backend: "cpu", + normalization: "peak", + compensateFrequencyResponse: true, + }, + }; + return { + channels, + sampleRate, + arrayDefinition, + scene: { speedOfSound: 343 }, + experiment, + }; +} + +test("localization algorithm registry exposes explicit scale and backend semantics", () => { + assertEqual(DEFAULT_LOCALIZATION_ALGORITHM, "das", "default algorithm"); + assertEqual(LOCALIZATION_ALGORITHM_IDS.join(","), "das,das-phat,gcc-phat", "registered ids"); + assertEqual(localizationAlgorithmDefinition("das").rawScale, "power", "DAS power scale"); + assertEqual(localizationAlgorithmDefinition("das-phat").rawScale, "score", "DAS-PHAT score scale"); + assertEqual(localizationAlgorithmDefinition("gcc-phat").rawScale, "score", "GCC-PHAT score scale"); + assert(localizationAlgorithmSupportsBackend("das-phat", "webgpu"), "DAS-PHAT WebGPU support"); + assert(!localizationAlgorithmSupportsBackend("gcc-phat", "webgpu"), "GCC-PHAT CPU-only support"); +}); + +test("generic CPU dispatcher runs all algorithms with explicit raw-map scales", async () => { + for (const algorithm of LOCALIZATION_ALGORITHM_IDS) { + const fixture = createFixture(algorithm); + const analysisPlan = createLocalizationAnalysisPlan(fixture); + const workload = estimateLocalizationWorkload({ ...fixture, analysisPlan }); + const result = await computeLocalizationMap({ + ...fixture, + analysisPlan, + pointBatchSize: 1024, + stftTransformBatchSize: 1024, + }); + const definition = localizationAlgorithmDefinition(algorithm); + + assertEqual(result.algorithm, algorithm, `${algorithm} result id`); + assertEqual(result.backend, "cpu", `${algorithm} backend`); + assertEqual(result.bands[0].rawScale, definition.rawScale, `${algorithm} raw scale`); + assertEqual(result.bands[0].rawMap.length, 9, `${algorithm} map shape`); + assertEqual(result.bands[0].peak.rawScale, definition.rawScale, `${algorithm} peak scale`); + assertEqual( + "rawPower" in result.bands[0].peak, + definition.rawScale === "power", + `${algorithm} peak power field`, + ); + assert( + result.bands[0].rawMap.every((value) => Number.isFinite(value) && value >= 0), + `${algorithm} finite nonnegative map`, + ); + assert(workload.pointCount === 9 && workload.microphoneCount === 4, `${algorithm} workload`); + } +}); + +test("backend dispatcher keeps GCC-PHAT on CPU and rejects an explicit GPU request", async () => { + const fixture = createFixture("gcc-phat"); + const calls = []; + const result = await computeLocalizationWithBackend(fixture, { + computeGpu: async () => { + calls.push("gpu"); + throw new Error("GCC-PHAT must not call WebGPU"); + }, + computeGccCpu: async () => { + calls.push("gcc-cpu"); + return { backend: "cpu", algorithm: "gcc-phat", bands: [] }; + }, + }); + assertEqual(calls.join(","), "gcc-cpu", "GCC auto routing"); + assertEqual(result.algorithm, "gcc-phat", "GCC result"); + + fixture.experiment.beamforming.backend = "webgpu"; + let backendError; + try { + await computeLocalizationWithBackend(fixture); + } catch (error) { + backendError = error; + } + assert(backendError instanceof TypeError, "explicit GCC WebGPU is rejected"); + + let legacyError; + try { + await computeDasWithBackend(fixture); + } catch (error) { + legacyError = error; + } + assert(legacyError instanceof RangeError, "legacy DAS entry rejects GCC"); +}); + +test("portable and streaming API accept an algorithm choice without Scene synthesis", async () => { + const fixture = createFixture("gcc-phat"); + const config = { + sampleRate: fixture.sampleRate, + array: fixture.arrayDefinition, + scan: fixture.experiment.scan, + bands: fixture.experiment.frequencyBands, + speedOfSound: fixture.scene.speedOfSound, + stft: fixture.experiment.stft, + algorithm: "gcc-phat", + phatFloorRatio: 1e-5, + normalization: "peak", + computeOptions: { pointBatchSize: 1024, stftTransformBatchSize: 1024 }, + }; + const plan = createLocalizationPlan(config); + const localized = await localizeMultichannel({ plan, channels: fixture.channels }); + assertEqual(plan.algorithm, "gcc-phat", "portable plan algorithm"); + assertEqual(plan.computeOptions.phatFloorRatio, 1e-5, "portable PHAT floor forwarding"); + assertEqual(plan.analysisPlan, plan.gccPhatPlan, "portable GCC plan alias"); + assertEqual(localized.algorithm, "gcc-phat", "portable result algorithm"); + assertEqual(localized.phatFloorRatio, 1e-5, "portable PHAT floor reaches GCC core"); + assertEqual(localized.bands[0].rawScale, "score", "portable score map"); + assertEqual(StreamingLocalizer, StreamingDasLocalizer, "streaming compatibility alias"); + + const stream = new StreamingLocalizer({ + ...config, + analysisWindowSamples: 128, + analysisHopSamples: 128, + maximumBufferedSamples: 256, + }); + const streamed = await stream.push(fixture.channels); + assertEqual(streamed.length, 1, "one GCC streaming result"); + assertEqual(streamed[0].algorithm, "gcc-phat", "streaming result algorithm"); + assertEqual(stream.getState().algorithm, "gcc-phat", "streaming state algorithm"); + assertEqual(stream.plan.gccPhatPlan.busy, false, "streaming GCC plan lock released"); + + const scorePeaks = detectBandPeaks( + Float32Array.of(0, 2, 0, 1), + { z: 1, xMin: 0, xMax: 1, yMin: 0, yMax: 1, columns: 2, rows: 2 }, + { rawScale: "score", minimumPeakDistanceM: 0 }, + ); + assertEqual(scorePeaks[0].rawScale, "score", "portable peak scale"); + assert(!("rawPower" in scorePeaks[0]), "score is not mislabeled as acoustic power"); +}); diff --git a/tests/unit/webgpu.test.js b/tests/unit/webgpu.test.js index fe88210..0091661 100644 --- a/tests/unit/webgpu.test.js +++ b/tests/unit/webgpu.test.js @@ -12,6 +12,48 @@ import { } from "/packages/webgpu-kernels/das.js"; import { computeDasWithBackend } from "/packages/webgpu-kernels/compute.js"; +test("WebGPU DAS-PHAT packing phase-normalizes above a relative floor", () => { + const stft = { + frameCount: 1, + microphoneCount: 3, + spectra: new Map([[ + 8, + { + real: Float32Array.of(3, 0, 1e-8), + imaginary: Float32Array.of(4, 2, 0), + }, + ]]), + }; + const microphones = [ + { + frequencyResponse: [ + { frequencyHz: 0, gainDb: 60, phaseDeg: 90 }, + { frequencyHz: 4_000, gainDb: 60, phaseDeg: 90 }, + ], + }, + { + frequencyResponse: [ + { frequencyHz: 0, gainDb: -60, phaseDeg: 0 }, + { frequencyHz: 4_000, gainDb: -60, phaseDeg: 0 }, + ], + }, + {}, + ]; + const packed = packStftBatch(stft, [8], { + algorithm: "das-phat", + microphones, + sampleRate: 8_000, + fftSize: 256, + }); + + assertNear(packed[0], 0.8, 1e-6, "response-phase corrected real"); + assertNear(packed[1], -0.6, 1e-6, "response-phase corrected imaginary"); + assertNear(packed[2], 0, 1e-7, "second microphone real"); + assertNear(packed[3], 1, 1e-7, "response gain is ignored"); + assertEqual(packed[4], 0, "below-floor real is suppressed"); + assertEqual(packed[5], 0, "below-floor imaginary is suppressed"); +}); + function assertNear(actual, expected, tolerance, label) { assert( Math.abs(actual - expected) <= tolerance, @@ -53,6 +95,17 @@ test("WebGPU DAS uniform과 microphone calibration을 WGSL layout으로 packing assertNear(calibration[3], Math.PI / 2, 1e-6, "+90 degree phase"); assertNear(calibration[4], 2, 1e-6, "-6 dB inverse gain"); assertNear(calibration[5], -Math.PI / 4, 1e-6, "-45 degree phase"); + + const phaseOnlyCalibration = packMicrophoneCalibration([ + { sensitivityDb: 0, phaseOffsetDeg: 0 }, + { sensitivityDb: 6.020599913, phaseOffsetDeg: 90 }, + { sensitivityDb: -6.020599913, phaseOffsetDeg: -45 }, + ], { includeAmplitude: false }); + assertNear(phaseOnlyCalibration[0], 1, 1e-7, "phase-only unity gain 0"); + assertNear(phaseOnlyCalibration[2], 1, 1e-7, "phase-only unity gain 1"); + assertNear(phaseOnlyCalibration[4], 1, 1e-7, "phase-only unity gain 2"); + assertNear(phaseOnlyCalibration[3], Math.PI / 2, 1e-6, "phase-only preserves +90 phase"); + assertNear(phaseOnlyCalibration[5], -Math.PI / 4, 1e-6, "phase-only preserves -45 phase"); }); test("STFT batch를 bin-frame-microphone-complex 순서로 packing한다", () => {