83 lines
3.4 KiB
Python
83 lines
3.4 KiB
Python
"""Numerical parity check between the standalone C++ CLI and the Python path.
|
|
|
|
Build first with ``cmake -S . -B build && cmake --build build --config Release``.
|
|
The script intentionally uses only the repository's existing NumPy and
|
|
Pyroomacoustics source tree for the reference calculation.
|
|
"""
|
|
|
|
from __future__ import annotations
|
|
|
|
import json
|
|
import subprocess
|
|
import sys
|
|
import tempfile
|
|
from pathlib import Path
|
|
|
|
import numpy as np
|
|
|
|
ROOT = Path(__file__).resolve().parents[1]
|
|
sys.path.insert(0, str(ROOT))
|
|
from src.doa_estimator import DoaConfig, DoaEstimator # noqa: E402
|
|
|
|
|
|
def plane_wave(config: DoaConfig, angle: float) -> np.ndarray:
|
|
rng = np.random.default_rng(1234)
|
|
source = rng.standard_normal(config.window_samples)
|
|
frequencies = np.fft.rfftfreq(config.window_samples, 1.0 / config.sample_rate)
|
|
source_spectrum = np.fft.rfft(source)
|
|
azimuths = np.deg2rad(config.microphone_azimuths_deg)
|
|
positions = config.radius_m * np.vstack((np.cos(azimuths), np.sin(azimuths)))
|
|
direction = np.array([np.cos(np.deg2rad(angle)), np.sin(np.deg2rad(angle))])
|
|
delays = direction @ positions / config.sound_speed_m_s
|
|
return np.stack(
|
|
[np.fft.irfft(source_spectrum * np.exp(2j * np.pi * frequencies * delay),
|
|
n=config.window_samples) for delay in delays], axis=1
|
|
).astype(np.float64) * 0.05
|
|
|
|
|
|
def main() -> int:
|
|
config = DoaConfig()
|
|
audio = plane_wave(config, 73.0)
|
|
build = ROOT / "build" / "cpp" / "doa_cli"
|
|
if sys.platform == "win32":
|
|
build = build.with_suffix(".exe")
|
|
candidate = ROOT / "build" / "cpp" / "Release" / "doa_cli.exe"
|
|
if candidate.exists():
|
|
build = candidate
|
|
if not build.exists():
|
|
raise SystemExit(f"CLI not found: {build}; build it with CMake first")
|
|
|
|
with tempfile.TemporaryDirectory() as directory:
|
|
matrix = Path(directory) / "audio.txt"
|
|
with matrix.open("w", encoding="ascii") as output:
|
|
output.write(f"{audio.shape[0]} {audio.shape[1]}\n")
|
|
np.savetxt(output, audio, fmt="%.17g")
|
|
estimator = DoaEstimator(config=config, start_worker=False)
|
|
try:
|
|
for algorithm, command in (("SRP-PHAT", "srp"), ("NormMUSIC", "normmusic")):
|
|
reference = estimator.estimate(audio, algorithm=algorithm)
|
|
raw = subprocess.check_output(
|
|
[str(build), command, str(matrix), str(config.sample_rate), str(config.nfft)],
|
|
text=True,
|
|
)
|
|
result = json.loads(raw)
|
|
error = abs((result["angle_degrees"] - reference.angle_degrees + 180) % 360 - 180)
|
|
print(f"{algorithm}: Python={reference.angle_degrees:.3f} C++={result['angle_degrees']:.3f} error={error:.3f} deg")
|
|
if error > 1.0:
|
|
return 1
|
|
reference_mvdr = estimator._mvdr.process(audio, 73.0).audio
|
|
raw = subprocess.check_output(
|
|
[str(build), "mvdr", str(matrix), "73", str(config.sample_rate), str(config.nfft)],
|
|
text=True,
|
|
)
|
|
actual_mvdr = np.fromstring(raw, sep="\n")
|
|
error = np.max(np.abs(reference_mvdr.astype(float) - actual_mvdr))
|
|
print(f"MVDR: max absolute sample error={error:.3e}")
|
|
return int(error > 1e-5)
|
|
finally:
|
|
estimator.close()
|
|
|
|
|
|
if __name__ == "__main__":
|
|
raise SystemExit(main())
|