#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
使用adt_evaluation核心组件实现AllRAD双耳解码
"""
import os
import numpy as np
import soundfile as sf
from scipy import signal
import spaudiopy as spa
import spherical_grids as sg
import spherical_data as sd
import adt_basic_decoders as bd
def allrad_binaural_adt(amb_path, sofa_path, output_path, N_sph=None):
"""
使用ADT核心组件的AllRAD双耳解码
步骤：
1. 使用ADT的t-design5200生成虚拟扬声器
2. 使用ADT的inversion解码器解码到虚拟扬声器
3. 使用ADT的allrad_v2rp映射到HRTF方向
4. HRTF卷积生成双耳音频
"""
print("=== ADT组件 AllRAD双耳解码 ===")
# 1. 加载Ambisonics音频
amb_signal, fs_amb = sf.read(amb_path)
amb_signal = amb_signal.T
num_channels = amb_signal.shape[0]
if N_sph is None:
N_sph = int(np.sqrt(num_channels)) - 1
if (N_sph + 1) ** 2 != num_channels:
raise ValueError(f"通道数{num_channels}不符合{N_sph}阶Ambisonics要求")
print(f"输入: {N_sph}阶Ambisonics, {num_channels}通道, {fs_amb}Hz")
# 2. 加载SOFA文件
hrirs_obj = spa.io.load_sofa_hrirs(sofa_path)
print(f"HRTF: {len(hrirs_obj.azi)}个方向, {hrirs_obj.fs}Hz")
# 3. 使用ADT的t-design生成虚拟扬声器
print("生成ADT虚拟t-design扬声器...")
td = sg.t_design5200()
virtual_az = td.az
virtual_el = td.el
print(f"虚拟扬声器: {len(virtual_az)}个")
# 4. 第一步：使用ADT解码器解码到虚拟扬声器
print("第一步: ADT解码到虚拟扬声器...")
sh_l, sh_m = zip(*[(l, m) for l in range(N_sph + 1) for m in range(-l, l + 1)])
# 使用ADT的逆变换解码器
M_virtual = bd.inversion(sh_l, sh_m, virtual_az, virtual_el)
virtual_signals = M_virtual @ amb_signal
print(f"虚拟扬声器信号: {virtual_signals.shape}")
# 5. 第二步：使用ADT的VBAP映射到HRTF方向
print("第二步: ADT VBAP映射到HRTF方向...")
# 转换为ADT所需的格式
Su_hrtf = np.array(sd.sph2cart(hrirs_obj.azi, hrirs_obj.zen)) # HRTF方向
Vu = np.array(sd.sph2cart(virtual_az, virtual_el)) # 虚拟扬声器
# 使用ADT的allrad_v2rp函数进行映射
V2H, _, _ = bd.allrad_v2rp(Su_hrtf, Vu, vbap_norm=True)
print(f"ADT VBAP映射矩阵: {V2H.shape}")
hrtf_signals = V2H @ virtual_signals
print(f"HRTF方向信号: {hrtf_signals.shape}")
# 6. 第三步：HRTF卷积
print("第三步: HRTF卷积...")
binaural_left = np.sum(
signal.oaconvolve(hrtf_signals, hrirs_obj.left, axes=-1),
axis=0
)
binaural_right = np.sum(
signal.oaconvolve(hrtf_signals, hrirs_obj.right, axes=-1),
axis=0
)
binaural = np.array([binaural_left, binaural_right])
# 7. 保存结果
sf.write(output_path, binaural.T, fs_amb)
print(f"✓ ADT AllRAD完成: {output_path}")
return binaural
def allrad_binaural_adt_direct(amb_path, sofa_path, output_path, N_sph=None):
"""
直接使用ADT的AllRAD函数的简化版本
"""
print("=== ADT直接 AllRAD双耳解码 ===")
# 1. 加载数据
amb_signal, fs_amb = sf.read(amb_path)
amb_signal = amb_signal.T
num_channels = amb_signal.shape[0]
if N_sph is None:
N_sph = int(np.sqrt(num_channels)) - 1
if (N_sph + 1) ** 2 != num_channels:
raise ValueError(f"通道数{num_channels}不符合{N_sph}阶Ambisonics要求")
print(f"输入: {N_sph}阶Ambisonics, {num_channels}通道, {fs_amb}Hz")
# 2. 加载HRTF
hrirs_obj = spa.io.load_sofa_hrirs(sofa_path)
print(f"HRTF: {len(hrirs_obj.azi)}个方向, {hrirs_obj.fs}Hz")
# 3. 直接使用ADT的AllRAD函数
print("直接使用ADT AllRAD...")
sh_l, sh_m = zip(*[(l, m) for l in range(N_sph + 1) for m in range(-l, l + 1)])
M_allrad = bd.allrad(
degree=sh_l,
order=sh_m,
speakers_azimuth=hrirs_obj.azi,
speakers_elevation=hrirs_obj.zen
)
print(f"ADT AllRAD矩阵: {M_allrad.shape}")
# 4. 解码和卷积
hrtf_signals = M_allrad @ amb_signal
print(f"HRTF方向信号: {hrtf_signals.shape}")
binaural_left = np.sum(
signal.oaconvolve(hrtf_signals, hrirs_obj.left, axes=-1),
axis=0
)
binaural_right = np.sum(
signal.oaconvolve(hrtf_signals, hrirs_obj.right, axes=-1),
axis=0
)
binaural = np.array([binaural_left, binaural_right])
# 5. 保存
sf.write(output_path, binaural.T, fs_amb)
print(f"ADT直接AllRAD完成: {output_path}")
return binaural
def main():
"""演示两种ADT实现方法"""
amb_path = "D:/projects/A2B/overall_Amb.wav"
sofa_path = "D:/projects/A2B/3d3a_Subject1_BIRs.sofa"
output_dir = "D:/projects/A2B/"
try:
os.makedirs(output_dir, exist_ok=True)
# 方法1：ADT组件完整实现
print("=== 方法1: ADT组件完整实现 ===")
output1 = os.path.join(output_dir, "adt_allrad_full.wav")
allrad_binaural_adt(amb_path, sofa_path, output1, N_sph=5)
# 方法2：ADT直接实现
print("\n=== 方法2: ADT直接实现 ===")
output2 = os.path.join(output_dir, "adt_allrad_direct.wav")
allrad_binaural_adt_direct(amb_path, sofa_path, output2, N_sph=5)
print("\n✓ ADT AllRAD双耳解码完成！")
except Exception as e:
print(f"错误: {e}")
if __name__ == '__main__':
main()
这个算法每一步的原理是什么，就是对什么形式的参数做了什么操作有什么物理意义，越详细越好，还有声源在扬声器的凸包外是什么意思为什么会报错，第二步: ADT VBAP映射到HRTF方向...
错误: 原点 [0 0 0] 不在扬声器阵列的凸包内
最新发布
本文介绍如何通过修改PHP配置文件来开启或关闭错误报告，适用于前台和后台界面的排错过程。通过简单的步骤即可查看详细的错误信息，便于快速定位并解决问题。
扫一扫
3704
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
