基于时域谱元法的高效头相关传输函数仿真技术制造技术

基于时域谱元法的高效头相关传输函数仿真技术，属于三维空间声场仿真领域。建立三维声场模型，赋予材料参数，设置边界条件；模型网格剖分；时域谱元法求解三维声波方程；计算声压是否收敛，若收敛，则生成时域头相关脉冲响应曲线(HRIR)；HRIR变换至频域形成头相关传输函数(HRTF)。时域谱元法应用于三维空间声场计算HRIR，运用快速傅里叶变换(FFT)计算HRTF。利用时域谱元法的高精度特性提高计算速度和减少计算时间。其计算结果与商业软件COMSOL Multiphysics计算结果进行对比验证其数值精度和可靠性。其中，HRIR应用于虚拟声源重现，构建三维空间声场。建三维空间声场。建三维空间声场。

【技术实现步骤摘要】
基于时域谱元法的高效头相关传输函数仿真技术


[0001]本专利技术属于三维空间声场仿真领域，涉及时域谱元法仿真技术，头相关脉冲响应(Head
‑
Related Impulse Response,HRIR)计算和头相关传输函数(Head Related Transfer Function,HRTF)计算。尤其是涉及一种基于时域谱元法的高效头相关传输函数仿真技术。

技术介绍

[0002]随着虚拟现实(VR)、增强现实(AR)以及空间音频环绕等技术的应用与发展，声场感知计算在相关
发挥重要作用。三维空间声场仿真是虚拟音频重现、空间声场音频等诸多领域中应用广泛的声场模拟技术。其中，HRTF是虚拟声源重现与空间音频环绕得以实现基础。目前由于HRTF实验测算需要耗费相当的成本，仅有部分实验室进行过实验测算。大多数实验室采用数值方法对人头模型进行处理求解。Mokhtari团队(Mokhtari P,Takemoto H,Nishimura R,et al.Computer simulation of KEMAR's head
‑
related transfer functions:Verificationwith measurements and acoustic effects of modifying head shape and pinna concavity[M].Principles and Applications of Spatial Hearing.2011:205
>‑
215.)建立KEMAR头模网格模型，使用FDTD计算双耳HRTF，结果证明该算法更接近实测值；Kahana等人(Kahana Y,Nelson P A.Boundary element simulations of the transfer function of human heads and baffledpinnae using accurate geometric models[J].Journal of sound and vibration,2007,300(3
‑
5):552
‑
579.)使用边界元法降低网格单元数量，减少矩阵求解时间，以快速获得相对精确的HRTF曲线。
[0003]三维空间声场仿真的商业软件普遍采用边界元法求解。边界元法基于边界积分方程理论的数值方法，只计算边界处的场且自动满足无穷远处边界条件。但处理人头模型时耳朵部分需要生成大量的网格单元。尤其对于高自由度的模型，随着扫频点增多，所需计算时间也随之增加。
[0004]声源定位、空间音频环绕技术在虚拟现实等领域的实现，使得声场计算算法成为技术广泛应用的关键。其中，HRTF计算是多个应用领域的重点研究对象。

技术实现思路

[0005]本专利技术的目的是提供可提高三维空间声场中求解HRTF曲线速度的一种基于时域谱元法的高效头相关传输函数仿真技术，使用时域谱元法来模拟三维空间中的声压场分布，可替代HRTF曲线的实验的直接测算，同时提高HRTF曲线求解精度与计算速度；最后生成声压级分布图。
[0006]本专利技术使用时域谱元法应用于三维空间声场计算HRIR曲线，接着运用快速傅里叶变换(Fast Fourier Transform，FFT)计算HRTF曲线。利用时域谱元算法的高精度特性提高计算速度和减少计算时间。其计算结果与商业软件COMSOL Multiphysics计算结果进行对
比验证其数值精度和可靠性。
[0007]本专利技术包括以下步骤：
[0008]1)构建仿真的三维人头，用其与三维空间模型做差集，生成去空人头的三维声场模型(以下简称：三维声场模型)；
[0009]2)将三维声场模型进行六面体网格剖分；
[0010]3)读取六面体网格单元，对三维声场模型进行预处理，处理过程包括设定边界条件、介质材料参数等；同时，预计算网格单元的雅可比行列式；
[0011]4)设置时间步长、总时长、Ricker子波声源点以及接收点；
[0012]5)建立时域有源声场
‑
三维声波方程(以下简称：三维声波方程)，选取合适的基函数和物理单元映射；使用时域谱元法需要将三维声波方程离散为矩阵方程；运用高斯积分法求解矩阵方程；
[0013]6)时域谱元法计算过程中进行时间迭代以求解时域总声压；同时，生成头相关脉冲响应HRIR曲线；利用快速傅里叶变换FFT计算频域总声压；
[0014]7)时域入射声压通过快速傅里叶变换为频域入射声压，作为参考声压；频域中计算每个频点的总声压与入射声压之比、取对数，生成头相关传输函数HRTF曲线；
[0015]8)上述步骤6)中若时域总声压计算不收敛，时间迭代中止；需要优化网格单元，重新读取模型计算。
[0016]在步骤1)中，所述生成去空人头的三维声场模型，选定去空人头的三维声场模型的计算区域时，确定其计算区域尺寸。
[0017]在步骤2)中，所述六面体网格剖分，选定的区域包含的声场边界是已知的，且区域生成的是单元雅可比矩阵大于零的任意六面体网格。
[0018]在步骤3)中，所述介质材料参数，包括声场介质密度和声速。
[0019]在步骤4)中，所述Ricker子波声源点需要确定声源点参数，即空间位置、数量与幅度，所述接收点需确定接收点参数，即空间位置与数量；
[0020]Ricker子波声源点的声场表达式如下：
[0021][0022]其中，f表示频率，f
M
表示峰值频率，ρ0表示静态介质密度，ω表示角频率，r表示源点到接收点空间距离，Q表示声源幅值(一般幅值取1)。
[0023]在步骤5)中，建立时域有源声场
‑
三维声波方程的步骤为：
[0024]三维运动方程和连续性方程：
[0025][0026][0027]其中，ρ(x,y,z,t)表示介质密度，p(x,y,z,t)表示压强，v(x,y,z,t)表示声传播速度，t表示时间；
[0028]下列各式中简记：ρ(x,y,z,t)
→
ρ，v(x,y,z,t)
→
v，p(x,y,z,t)
→
p,f(x,y,z,t)
→
f；
[0029]小振幅的振动情况下线性化近似ρ≈ρ0+Δρ，前二式替代为线性声学方程：
[0030][0031][0032]其中，ρ0表示静态介质密度(常数)，Δρ表示密度的微小变化量；
[0033]上式对t求导：
[0034][0035]小振幅振动情况下声学物态方程：
[0036][0037]其中，下标s表示绝热过程，下标0表示小振幅声波扰动平衡态时的数值，c0表示声速；
[0038]联立前二式，建立理想流媒体介质中的三维声波方程：
[0039][0040]其中，
[00本文档来自技高网...

【技术保护点】

【技术特征摘要】
1.基于时域谱元法的高效头相关传输函数仿真技术，包括以下步骤：1)构建仿真的三维人头，用其与三维空间模型做差集，生成去空人头的三维声场模型；2)将步骤1)生成的三维声场模型进行六面体网格剖分；3)读取六面体网格单元，对三维声场模型进行预处理，处理过程包括设定边界条件、介质材料参数；预计算网格单元的雅可比行列式；4)设置时间步长、总时长、Ricker子波声源点以及接收点；5)建立时域有源声场
‑
三维声波方程，选取合适的基函数和物理单元映射；运用高斯积分法离散时域有源声场
‑
三维声波方程并组合为矩阵方程；6)时域谱元法计算过程中进行时间迭代以求解时域总声压；同时，生成头相关脉冲响应HRIR曲线；利用快速傅里叶变换计算频域总声压；7)时域入射声压通过快速傅里叶变换为频域入射声压，作为参考声压；频域中计算每个频点的总声压与入射声压之比、取对数，生成头相关传输函数HRTF曲线；8)上述步骤6)中若时域总声压计算不收敛，时间迭代中止；需要优化网格单元，重新读取模型计算。2.如权利要求1所述基于时域谱元法的高效头相关传输函数仿真技术，其特征在于在步骤1)中，所述生成去空人头的三维声场模型，选定去空人头的三维声场模型的计算区域时，确定其计算区域尺寸。3.如权利要求1所述基于时域谱元法的高效头相关传输函数仿真技术，其特征在于在步骤2)中，所述六面体网格剖分，选定的区域包含的声场边界是已知的，且区域生成的是单元雅可比矩阵大于零的任意六面体网格。4.如权利要求1所述基于时域谱元法的高效头相关传输函数仿真技术，其特征在于在步骤3)中，所述介质材料参数，包括声场介质密度和声速。5.如权利要求1所述基于时域谱元法的高效头相关传输函数仿真技术，其特征在于在步骤4)中，所述Ricker子波声源点需要确定声源点参数，即空间位置、数量与幅度，所述接收点需确定接收点参数，即空间位置与数量；Ricker子波声源点的声场表达式如下：其中，f表示频率，f
M
表示峰值频率，ρ0表示静态介质密度，ω表示角频率，r表示源点到接收点空间距离，Q表示声源幅值。6.如权利要求1所述基于时域谱元法的高效头相关传输函数仿真技术，其特征在于步骤5)中，建立时域有源声场
‑
三维声波方程的步骤为：三维运动方程和连续性方程：三维运动方程和连续性方程：
其中，ρ(x,y,z,t)表示介质密度，p(x,y,z,t)表示压强，v(x,y,z,t)表示声传播速度，t表示时间；下列各式中简记：ρ(x,y,z,t)
→
ρ，v(x,y,z,t)
→
v，p(x,y,z,t)
→
p,f(x,y,z,t)
→
f；小振幅的振动情况下线性化近似ρ≈ρ0+Δρ，前二式替代为线性声学方程：+Δρ，前二式替代为线性声学方程：其中，ρ0表示静态介质密度，Δρ表示密度的微小变化量；上式对t求导：小振幅振动情况下声学物态方程：其中，下标s表示绝热过程，下标0表示小振幅声波扰动平衡态时的数值，c0表示声速；联立前二式，建立理想流媒体介质中的三维声波方程：其中，加入声源项，建立时域有源声场
‑
三维声波方程：所述三维声波方程可选用的边界条件：硬声场边界p＝0软声场边界其中，f表示声源，q
d
表示单极子源；确定N阶高斯
‑
勒让德
‑
洛巴托基(GLL)函数作为基函数，其形式如下：其中，L
N
(ξ
r
)与L'
N
(ξ
r
)为N点勒让德多项式及其导数；且...

【专利技术属性】
技术研发人员：刘娜，程文烁，万领，庄明伟，李长火，张志，蔡国雄，
申请(专利权)人：厦门大学，
类型：发明
国别省市：