Files
micarrsim/packages/acoustics-core/propagation.js
T

504 lines
18 KiB
JavaScript

import { interpolateFrequencyResponse } from "./frequency-response.js";
import {
ceilSampleIndex,
ceilSamplePosition,
dbToLinear,
distance3,
TAU,
} from "./math.js";
import {
addFrequencyDomainDelay,
addFrequencyDomainPaths,
createDelayWorkspace,
nextPowerOfTwo,
prepareDelaySpectrum,
} from "./fractional-delay.js";
import { buildAcousticPaths } from "./environment.js";
import { renderSignal } from "./signal.js";
import { describeFrequencyAliasing } from "./sampling.js";
export function propagationDelaySeconds(sourcePosition, receiverPosition, speedOfSound) {
return distance3(sourcePosition, receiverPosition) / speedOfSound;
}
export function sphericalAttenuation(distance, referenceDistance = 1, minimumDistance = 0.01) {
return referenceDistance / Math.max(distance, minimumDistance);
}
function addToneEvent(
events,
sampleIndex,
sineDelta,
cosineDelta,
activePathDelta,
) {
const event = events.get(sampleIndex);
if (event) {
event.sine += sineDelta;
event.cosine += cosineDelta;
event.activePaths += activePathDelta;
} else {
events.set(sampleIndex, {
sine: sineDelta,
cosine: cosineDelta,
activePaths: activePathDelta,
});
}
}
function carrierAtSample(signedFrequencyHz, sampleRate, absoluteSampleIndex) {
const cycles = (signedFrequencyHz * absoluteSampleIndex) / sampleRate;
const fractionalCycles = cycles - Math.floor(cycles);
const phase = TAU * fractionalCycles;
return { cosine: Math.cos(phase), sine: Math.sin(phase) };
}
/**
* Apply propagation and microphone response to an analog tone before sampling.
*
* This order is essential above Nyquist: temporal samples fold to the ADC
* alias, while inter-microphone phase still comes from the original acoustic
* frequency and path delay.
*/
function addSampledAnalogTonePaths(options) {
const {
output,
entry,
scene,
microphone,
paths,
timeMapping = {},
} = options;
const signal = entry.signal;
const frequencyHz = signal.frequency;
const sampling = describeFrequencyAliasing(frequencyHz, scene.sampleRate);
const sceneSampleOffset = timeMapping.sceneSampleOffset ?? 0;
const signalStartTime = timeMapping.signalStartTime ?? entry.startTime ?? 0;
const signalDuration = timeMapping.signalDuration ?? entry.duration ?? scene.duration;
if (!Number.isSafeInteger(sceneSampleOffset) || sceneSampleOffset < 0) {
throw new RangeError("sceneSampleOffset must be a non-negative safe integer.");
}
if (!Number.isFinite(signalStartTime)) {
throw new TypeError("signalStartTime must be finite.");
}
if (!Number.isFinite(signalDuration) || signalDuration <= 0) {
throw new RangeError("signalDuration must be a positive finite number.");
}
const sourcePhase = ((signal.phaseDeg ?? 0) * Math.PI) / 180;
const response = interpolateFrequencyResponse(
microphone.frequencyResponse,
frequencyHz,
{},
);
const microphonePhase =
((microphone.phaseOffsetDeg ?? 0) * Math.PI) / 180 + response.phaseRadians;
const baseGain =
dbToLinear(entry.gainDb ?? 0) *
dbToLinear(microphone.sensitivityDb ?? 0) *
response.linearGain;
const events = new Map();
for (const path of paths) {
const firstSample = Math.max(
0,
ceilSampleIndex(signalStartTime + path.delaySeconds, scene.sampleRate) -
sceneSampleOffset,
);
const lastSample = Math.min(
output.length,
ceilSampleIndex(
signalStartTime + signalDuration + path.delaySeconds,
scene.sampleRate,
) - sceneSampleOffset,
);
if (firstSample >= lastSample) continue;
const amplitude = baseGain * path.gain;
const delayedCycles = (frequencyHz * (signalStartTime + path.delaySeconds)) % 1;
const phase = sourcePhase + microphonePhase - TAU * delayedCycles;
const sineCoefficient = amplitude * Math.cos(phase);
const cosineCoefficient = amplitude * Math.sin(phase);
addToneEvent(events, firstSample, sineCoefficient, cosineCoefficient, 1);
addToneEvent(events, lastSample, -sineCoefficient, -cosineCoefficient, -1);
}
if (events.size === 0) return;
const orderedEvents = [...events.entries()].sort((left, right) => left[0] - right[0]);
const phaseStep =
(TAU * sampling.signedAdcFrequencyHz) / scene.sampleRate;
const stepReal = Math.cos(phaseStep);
const stepImaginary = Math.sin(phaseStep);
const firstSample = orderedEvents[0][0];
const finalEventSample = orderedEvents.at(-1)[0];
const initialCarrier = carrierAtSample(
sampling.signedAdcFrequencyHz,
scene.sampleRate,
sceneSampleOffset + firstSample,
);
let carrierCosine = initialCarrier.cosine;
let carrierSine = initialCarrier.sine;
let sineCoefficient = 0;
let cosineCoefficient = 0;
let activePathCount = 0;
let eventIndex = 0;
const phaseReanchorInterval = 1024;
let sampleIndex = firstSample;
while (sampleIndex < output.length && sampleIndex < finalEventSample) {
while (
eventIndex < orderedEvents.length &&
orderedEvents[eventIndex][0] === sampleIndex
) {
sineCoefficient += orderedEvents[eventIndex][1].sine;
cosineCoefficient += orderedEvents[eventIndex][1].cosine;
activePathCount += orderedEvents[eventIndex][1].activePaths;
eventIndex += 1;
}
if (activePathCount === 0 && eventIndex < orderedEvents.length) {
sampleIndex = orderedEvents[eventIndex][0];
const carrier = carrierAtSample(
sampling.signedAdcFrequencyHz,
scene.sampleRate,
sceneSampleOffset + sampleIndex,
);
carrierCosine = carrier.cosine;
carrierSine = carrier.sine;
continue;
}
output[sampleIndex] +=
carrierSine * sineCoefficient + carrierCosine * cosineCoefficient;
const nextSample = sampleIndex + 1;
if (nextSample % phaseReanchorInterval === 0) {
const carrier = carrierAtSample(
sampling.signedAdcFrequencyHz,
scene.sampleRate,
sceneSampleOffset + nextSample,
);
carrierCosine = carrier.cosine;
carrierSine = carrier.sine;
} else {
const nextCosine =
carrierCosine * stepReal - carrierSine * stepImaginary;
carrierSine =
carrierSine * stepReal + carrierCosine * stepImaginary;
carrierCosine = nextCosine;
}
sampleIndex = nextSample;
}
}
export function synthesizeMicrophoneSignals(options) {
const {
scene,
arrayDefinition,
assetBuffers = new Map(),
referenceDistance = 1,
minimumDistance = 0.01,
maximumBufferBytes = 512 * 1024 * 1024,
signalOffsets,
} = options;
const sampleCount = ceilSampleIndex(scene.duration, scene.sampleRate);
const sceneSampleOffset = signalOffsets?.sceneSampleOffset ?? 0;
const cropStartSample = signalOffsets?.cropStartSample ?? 0;
const cropSampleCount = signalOffsets?.cropSampleCount ?? sampleCount;
if (!Number.isSafeInteger(sceneSampleOffset) || sceneSampleOffset < 0) {
throw new RangeError("sceneSampleOffset must be a non-negative safe integer.");
}
if (
!Number.isSafeInteger(cropStartSample) ||
cropStartSample < 0 ||
!Number.isSafeInteger(cropSampleCount) ||
cropSampleCount <= 0 ||
cropStartSample + cropSampleCount > sampleCount
) {
throw new RangeError("signal output crop must be a non-empty range inside the synthesis buffer.");
}
const cropEndSample = cropStartSample + cropSampleCount;
const needsOutputCrop = cropStartSample !== 0 || cropSampleCount !== sampleCount;
const sourcePathSets = scene.sources.map((definition) => ({
definition,
microphonePaths: arrayDefinition.microphones.map((microphone) =>
buildAcousticPaths({
scene,
sourcePosition: definition.position,
receiverPosition: microphone.position,
sourceId: definition.id,
receiverId: microphone.id,
referenceDistance,
minimumDistance,
}),
),
}));
const acousticPathCounts = {};
let acousticPathCount = 0;
let maximumDelaySamples = 0;
for (const { microphonePaths } of sourcePathSets) {
for (const paths of microphonePaths) {
for (const path of paths) {
acousticPathCount += 1;
acousticPathCounts[path.kind] = (acousticPathCounts[path.kind] ?? 0) + 1;
maximumDelaySamples = Math.max(
maximumDelaySamples,
path.delaySeconds * scene.sampleRate,
);
}
}
}
const sampledSourceCount = scene.sources
.filter((source) => source.signal.type !== "tone")
.length;
const hasSampledAmbient = (scene.ambient ?? [])
.some((ambient) => ambient.signal.type !== "tone");
const retainedBufferCount =
arrayDefinition.microphones.length +
sampledSourceCount +
(hasSampledAmbient ? 1 : 0);
const delayFftSize = nextPowerOfTwo(
sampleCount + ceilSamplePosition(maximumDelaySamples) + 2,
);
const delayWorkingBytes =
sampledSourceCount > 0 || hasSampledAmbient
? delayFftSize * Float64Array.BYTES_PER_ELEMENT * 4
: 0;
const outputCropBytes = needsOutputCrop
? arrayDefinition.microphones.length * cropSampleCount * Float32Array.BYTES_PER_ELEMENT
: 0;
const estimatedBufferBytes =
retainedBufferCount * sampleCount * Float32Array.BYTES_PER_ELEMENT +
delayWorkingBytes +
outputCropBytes;
if (estimatedBufferBytes > maximumBufferBytes) {
throw new Error(
`예상 signal buffer ${(estimatedBufferBytes / 1024 ** 2).toFixed(1)} MiB가 ` +
`현재 한도 ${(maximumBufferBytes / 1024 ** 2).toFixed(1)} MiB를 넘습니다. duration 또는 sample rate를 줄이세요.`,
);
}
const synthesisChannels = arrayDefinition.microphones.map(
() => new Float32Array(sampleCount),
);
const renderedSources = sourcePathSets.map(({ definition, microphonePaths }, sourceIndex) => ({
definition,
microphonePaths,
timeMapping: {
sceneSampleOffset,
...(signalOffsets?.sources?.[sourceIndex] ?? {}),
},
...(definition.signal.type === "tone"
? {}
: {
samples: renderSignal(definition, scene, assetBuffers, 0, {
sceneSampleOffset,
...(signalOffsets?.sources?.[sourceIndex] ?? {}),
}),
}),
}));
const outputStartSample = sceneSampleOffset + cropStartSample;
const outputEndSample = outputStartSample + cropSampleCount;
const sourceArrivals = renderedSources.map(({ definition, microphonePaths, timeMapping }) => {
const emissionStartTime = timeMapping.signalStartTime ?? definition.startTime ?? 0;
const emissionDuration = timeMapping.signalDuration ?? definition.duration ?? scene.duration;
let earliestArrivalSample = Infinity;
let latestArrivalEndSample = -Infinity;
let overlapsOutput = false;
const overlappingPathKinds = new Set();
for (const paths of microphonePaths) {
for (const path of paths) {
const arrivalStartSample = ceilSampleIndex(
emissionStartTime + path.delaySeconds,
scene.sampleRate,
);
const arrivalEndSample = ceilSampleIndex(
emissionStartTime + emissionDuration + path.delaySeconds,
scene.sampleRate,
);
earliestArrivalSample = Math.min(earliestArrivalSample, arrivalStartSample);
latestArrivalEndSample = Math.max(latestArrivalEndSample, arrivalEndSample);
if (arrivalStartSample < outputEndSample && arrivalEndSample > outputStartSample) {
overlapsOutput = true;
overlappingPathKinds.add(path.kind);
}
}
}
return {
sourceId: definition.id,
emissionStartSample: ceilSampleIndex(emissionStartTime, scene.sampleRate),
emissionEndSample: ceilSampleIndex(
emissionStartTime + emissionDuration,
scene.sampleRate,
),
earliestArrivalSample: Number.isFinite(earliestArrivalSample)
? earliestArrivalSample
: null,
latestArrivalEndSample: Number.isFinite(latestArrivalEndSample)
? latestArrivalEndSample
: null,
overlapsOutput,
overlappingPathKinds: [...overlappingPathKinds].sort(),
};
});
renderedSources.forEach(({ definition, microphonePaths, samples, timeMapping }) => {
let sourceMaximumDelaySamples = 0;
for (const paths of microphonePaths) {
for (const path of paths) {
sourceMaximumDelaySamples = Math.max(
sourceMaximumDelaySamples,
path.delaySeconds * scene.sampleRate,
);
}
}
const spectrum = samples
? prepareDelaySpectrum(samples, sourceMaximumDelaySamples)
: undefined;
const workspace = spectrum
? createDelayWorkspace(spectrum.fftSize)
: undefined;
arrayDefinition.microphones.forEach((microphone, microphoneIndex) => {
const paths = microphonePaths[microphoneIndex];
const sensitivity = dbToLinear(microphone.sensitivityDb ?? 0);
if (definition.signal.type === "tone") {
addSampledAnalogTonePaths({
output: synthesisChannels[microphoneIndex],
entry: definition,
scene,
microphone,
paths,
timeMapping,
});
} else {
addFrequencyDomainPaths({
output: synthesisChannels[microphoneIndex],
spectrum,
paths: paths.map((path) => ({
delaySamples: path.delaySeconds * scene.sampleRate,
gain: path.gain,
})),
gain: sensitivity,
phaseOffsetDeg: microphone.phaseOffsetDeg ?? 0,
frequencyResponse: microphone.frequencyResponse,
sampleRate: scene.sampleRate,
workspace,
});
}
});
});
if (Array.isArray(scene.ambient)) {
scene.ambient.forEach((ambient, ambientIndex) => {
synthesisChannels.forEach((channel, microphoneIndex) => {
const microphone = arrayDefinition.microphones[microphoneIndex];
const ambientTimeMapping = {
sceneSampleOffset,
...(signalOffsets?.ambient?.[ambientIndex] ?? {}),
};
if (ambient.signal.type === "tone") {
addSampledAnalogTonePaths({
output: channel,
entry: ambient,
scene,
microphone,
paths: [{ delaySeconds: 0, gain: 1 }],
timeMapping: ambientTimeMapping,
});
return;
}
const seedSalt = Math.imul(ambientIndex + 1, 0x45d9f3b) ^ Math.imul(microphoneIndex + 1, 0x27d4eb2d);
const ambientSamples = renderSignal(ambient, scene, assetBuffers, seedSalt, {
...ambientTimeMapping,
});
const spectrum = prepareDelaySpectrum(ambientSamples, 0);
addFrequencyDomainDelay({
output: channel,
spectrum,
delaySamples: 0,
gain: dbToLinear(microphone.sensitivityDb ?? 0),
phaseOffsetDeg: microphone.phaseOffsetDeg ?? 0,
frequencyResponse: microphone.frequencyResponse,
sampleRate: scene.sampleRate,
workspace: createDelayWorkspace(spectrum.fftSize),
});
});
});
}
const channels = needsOutputCrop
? synthesisChannels.map((channel) => channel.slice(cropStartSample, cropEndSample))
: synthesisChannels;
let preNormalizationPeak = 0;
for (const channel of channels) {
for (const value of channel) {
preNormalizationPeak = Math.max(preNormalizationPeak, Math.abs(value));
}
}
const appliedScale =
scene.mix?.normalization === "peak" && preNormalizationPeak > 0
? scene.mix.targetPeak / preNormalizationPeak
: 1;
let peak = 0;
let clippedSampleCount = 0;
for (const channel of channels) {
for (let sampleIndex = 0; sampleIndex < channel.length; sampleIndex += 1) {
channel[sampleIndex] *= appliedScale;
const magnitude = Math.abs(channel[sampleIndex]);
peak = Math.max(peak, magnitude);
if (magnitude > 1) clippedSampleCount += 1;
}
}
return {
channels,
sampleRate: scene.sampleRate,
peak,
preNormalizationPeak,
appliedScale,
clippedSampleCount,
estimatedBufferBytes,
acousticModel: scene.acoustics?.model ?? "free-field",
acousticPathCount,
acousticPathCounts,
windowing: {
model: signalOffsets?.windowingModel ?? "whole-scene",
renderStartSample: sceneSampleOffset,
outputStartSample,
outputSampleCount: cropSampleCount,
preRollSamples: cropStartSample,
activePathMaximumDelaySeconds: maximumDelaySamples / scene.sampleRate,
requestedMaximumPropagationDelaySeconds:
signalOffsets?.maximumPropagationDelaySeconds ?? maximumDelaySamples / scene.sampleRate,
geometryModel: Number.isFinite(signalOffsets?.sourcePositionSampleTime)
? "source-position-fixed-at-window-midpoint"
: "static-source-position",
sourcePositionSampleTime: signalOffsets?.sourcePositionSampleTime ?? null,
dopplerMode: "not-modeled",
delayedTailPolicy: needsOutputCrop
? "propagation-preroll-then-output-crop"
: "whole-render-buffer",
arrivalOverlapModel: "path-gated-half-open-sample-interval",
arrivingSourceIds: sourceArrivals
.filter(({ overlapsOutput }) => overlapsOutput)
.map(({ sourceId }) => sourceId),
sourceArrivals,
},
sampling: {
model: "ideal-real-sampler-no-antialias-filter",
sampleRate: scene.sampleRate,
nyquistHz: scene.sampleRate / 2,
toneObservations: [
...scene.sources
.filter((source) => source.signal.type === "tone")
.map((source) => ({
scope: "source",
id: source.id,
...describeFrequencyAliasing(source.signal.frequency, scene.sampleRate),
})),
...(scene.ambient ?? [])
.filter((ambient) => ambient.signal.type === "tone")
.map((ambient) => ({
scope: "ambient",
id: ambient.id,
...describeFrequencyAliasing(ambient.signal.frequency, scene.sampleRate),
})),
],
},
};
}