feat: add realistic portable localization core

packages/acoustics-core/fractional-delay.js

@@ -1,4 +1,5 @@
 import { fftInPlace } from "./fft.js";
+import { interpolateFrequencyResponse } from "./frequency-response.js";
 import { TAU } from "./math.js";
 
 export function nextPowerOfTwo(value) {
@@ -21,25 +22,96 @@ export function createDelayWorkspace(fftSize) {
   return { real: new Float64Array(fftSize), imaginary: new Float64Array(fftSize) };
 }
 
-export function addFrequencyDomainDelay(options) {
-  const { output, spectrum, delaySamples, gain = 1, phaseOffsetDeg = 0, workspace } = options;
-  if (!(delaySamples >= 0) || !Number.isFinite(delaySamples)) throw new Error(`유효하지 않은 delay: ${delaySamples}`);
+export function addFrequencyDomainPaths(options) {
+  const {
+    output,
+    spectrum,
+    paths,
+    gain = 1,
+    phaseOffsetDeg = 0,
+    frequencyResponse,
+    sampleRate,
+    workspace,
+  } = options;
+  if (!Array.isArray(paths) || paths.length === 0) {
+    throw new Error("frequency-domain path가 하나 이상 필요합니다.");
+  }
+  for (const path of paths) {
+    if (!(path.delaySamples >= 0) || !Number.isFinite(path.delaySamples)) {
+      throw new Error("유효하지 않은 delay: " + path.delaySamples);
+    }
+    if (!Number.isFinite(path.gain)) throw new Error("path gain은 유한한 숫자여야 합니다.");
+  }
   if (workspace.real.length !== spectrum.fftSize || workspace.imaginary.length !== spectrum.fftSize) {
     throw new Error("fractional-delay workspace 크기가 spectrum과 다릅니다.");
   }
+  const hasFrequencyResponse = Array.isArray(frequencyResponse) && frequencyResponse.length > 0;
+  if (hasFrequencyResponse && (!(sampleRate > 0) || !Number.isFinite(sampleRate))) {
+    throw new Error("frequencyResponse 적용에는 유효한 sampleRate가 필요합니다.");
+  }
 
   const size = spectrum.fftSize;
   const calibrationPhase = (phaseOffsetDeg * Math.PI) / 180;
+  const response = {};
+  workspace.real.fill(0);
+  workspace.imaginary.fill(0);
+  const halfSize = size / 2;
+  const phaseReanchorInterval = 1024;
+  for (const path of paths) {
+    const phaseStep = (-TAU * path.delaySamples) / size;
+    const stepReal = Math.cos(phaseStep);
+    const stepImaginary = Math.sin(phaseStep);
+    let phaseReal = 1;
+    let phaseImaginary = 0;
+    for (let bin = 0; bin <= halfSize; bin += 1) {
+      workspace.real[bin] += phaseReal * path.gain;
+      workspace.imaginary[bin] += phaseImaginary * path.gain;
+      if (bin > 0 && bin < halfSize) {
+        const negativeBin = size - bin;
+        workspace.real[negativeBin] += phaseReal * path.gain;
+        workspace.imaginary[negativeBin] -= phaseImaginary * path.gain;
+      }
+      const nextBin = bin + 1;
+      if (nextBin % phaseReanchorInterval === 0) {
+        phaseReal = Math.cos(phaseStep * nextBin);
+        phaseImaginary = Math.sin(phaseStep * nextBin);
+      } else {
+        const nextReal = phaseReal * stepReal - phaseImaginary * stepImaginary;
+        phaseImaginary = phaseReal * stepImaginary + phaseImaginary * stepReal;
+        phaseReal = nextReal;
+      }
+    }
+  }
+
   for (let bin = 0; bin < size; bin += 1) {
     const signedBin = bin <= size / 2 ? bin : bin - size;
     const phaseSign = signedBin === 0 ? 0 : signedBin > 0 ? 1 : -1;
-    const phase = (-TAU * signedBin * delaySamples) / size + phaseSign * calibrationPhase;
-    const phaseReal = Math.cos(phase);
-    const phaseImaginary = Math.sin(phase);
+    let responseGain = 1;
+    let responsePhase = 0;
+    if (hasFrequencyResponse) {
+      interpolateFrequencyResponse(
+        frequencyResponse,
+        (Math.abs(signedBin) * sampleRate) / size,
+        response,
+      );
+      responseGain = response.linearGain;
+      responsePhase = response.phaseRadians;
+    }
+    const pathReal = workspace.real[bin];
+    const pathImaginary = workspace.imaginary[bin];
+    const responseAndCalibrationPhase = phaseSign * (calibrationPhase + responsePhase);
+    const calibrationReal = Math.cos(responseAndCalibrationPhase) * responseGain * gain;
+    const calibrationImaginary = Math.sin(responseAndCalibrationPhase) * responseGain * gain;
+    const transferReal =
+      pathReal * calibrationReal - pathImaginary * calibrationImaginary;
+    const transferImaginary =
+      pathReal * calibrationImaginary + pathImaginary * calibrationReal;
     const inputReal = spectrum.real[bin];
     const inputImaginary = spectrum.imaginary[bin];
-    workspace.real[bin] = inputReal * phaseReal - inputImaginary * phaseImaginary;
-    workspace.imaginary[bin] = inputReal * phaseImaginary + inputImaginary * phaseReal;
+    workspace.real[bin] =
+      inputReal * transferReal - inputImaginary * transferImaginary;
+    workspace.imaginary[bin] =
+      inputReal * transferImaginary + inputImaginary * transferReal;
   }
 
   // A real sequence has a real-valued Nyquist coefficient. Projecting this single self-conjugate
@@ -48,6 +120,14 @@ export function addFrequencyDomainDelay(options) {
   fftInPlace(workspace.real, workspace.imaginary, true);
   const outputLength = Math.min(output.length, workspace.real.length);
   for (let sampleIndex = 0; sampleIndex < outputLength; sampleIndex += 1) {
-    output[sampleIndex] += workspace.real[sampleIndex] * gain;
+    output[sampleIndex] += workspace.real[sampleIndex];
   }
 }
+
+export function addFrequencyDomainDelay(options) {
+  const { delaySamples, gain = 1, ...shared } = options;
+  return addFrequencyDomainPaths({
+    ...shared,
+    paths: [{ delaySamples, gain }],
+  });
+}