对头部相关冲激响应进行建模的制作方法

1.本公开涉及渲染空间声音。


背景技术：

2.我们有捕获朝着我们传播的声波的两只耳朵。图1示出了从由球坐标系中的仰角角度和方位角角度对指定的到达方向(doa)朝着收听者传播的声波。在朝着我们的传播路径上，在到达我们的左耳鼓和右耳鼓之前，每个声波与我们的上躯干、头部、外耳和环境介质相互作用。这个相互作用导致到达左耳鼓和右耳鼓的波形的时间改变和频谱改变，其中的一些改变是doa相关的。我们的听觉系统已经学会解释这些改变，以推断声波本身的各种空间特性以及收听者发现他自身/她自身所处的声音环境。这个能力称为空间听觉，其涉及我们如何评估嵌入在双耳信号(即右耳道和左耳道中的声音信号)中的空间线索，以推断由声音事件(物理声源)引发的听觉事件的位置和由我们所处的物理环境(例如，小房间、贴砖的浴室、礼堂、洞穴)导致的声音特性。人类的这个能力(空间听觉)进而可以被开发用于通过在双耳信号中重新引入将导致对声音的空间感知的空间线索来创建空间声音场景。
3.主空间线索包括：1)角度相关的线索：双耳线索(即耳间电平差(ild)和耳间时差(itd))和单耳(或频谱)线索；2)距离相关的线索：强度以及直接-混响(d/r)能量比。图2示出了朝着收听者传播的声波的itd和频谱线索的示例。这两个图示出了在0度仰角和40度方位角处获得的hr滤波器对的幅度响应(所述数据来自cipic数据库：对象id 28。所述数据库是公开可用的，其可以通过urlwww.ece.ucdavis.edu/cipic/spatial-sound/hrtf-data/访问)。在图1和图2中，使用正方位角方向向右的惯例，并且这也是在本文的其余部分中使用的惯例。然而，一些hr滤波器集合的确使用另一惯例(其中的正方位角方向向左)。波形的短时间的doa相关的时间和频谱改变(1-5毫秒)的数学表示是所谓的头部相关(hr)滤波器。这些滤波器的频域(fd)表示是所谓的头部相关传递函数(hrtf)并且时域(td)表示是头部相关冲激响应(hrir)。已经逐渐建立了基于hr滤波器的双耳渲染方法，其中，通过直接利用希望位置的hr滤波器对对声源信号滤波而生成空间声音场景。这个方法对于许多新兴应用特别有吸引力，例如虚拟现实(vr)、增强现实(ar)、混合现实(mr)或扩展现实(xr)、以及移动通信系统，其中通常都使用耳机。
4.通常根据作为线性动态系统的冲激响应的测量结果来估计头部相关的(hr)滤波器，该线性动态系统将原始声音信号(输入信号)转换成可以在收听对象的耳道内测量到的左耳信号和右耳信号(输出信号)，其中该收听对象在距该收听对象(例如，人造头部、人体模型或人体)恒定半径的球表面上的预定义仰角角度和方位角角度集合处。已估计的hr滤波器通常设置为fir滤波器并且可以在该格式下直接使用。为了实现高效的双耳渲染，hrtf对可以转换成耳间传递函数(itf)或经修改的itf以防止陡峭的频谱峰。备选地，hrtf可以通过参数表示进行描述。这些参数化的hrtf易于与参数多通道声音编码器(例如mpeg环绕声和空间声音对象编码(saoc))集成。
5.渲染空间声音信号以提供针对空间中任意位置的声音的令人信服的空间感知需
要针对对应位置的hr滤波器对，因此需要针对在2d球体上的精细采样位置的hr滤波器集合。最小可听角度(maa)将我们的听觉系统的灵敏度表征为声音事件的角位移。关于方位角的位置，所报告的是，对于宽带噪声突发，maa在前向和后向上最小(约1度)并且对于侧向声源要大得多(约10度)。正中面中的maa随仰角增大。所报告的是，对于宽带噪声突发，仰角上的平均maa与4度一样小。当前，存在一些公开可用的在空间上密集采样的hr滤波器数据库，例如sadie数据库、cipic数据库。然而，特别是在仰角的样本上，它们都不能完全满足maa要求。即使人造头部neumann(纽曼)ku100和kemar人体模型的saide数据集包含多于8000个测量结果，但是其在-15度与15度之间的仰角上的采样分辨率是15度，而根据maa研究则需要4度。不可避免的是，需要hr滤波器的角度插值，从而可以在没有经测量的实际滤波器的位置渲染声源。图3示出了在2d球体上对网格进行采样的示例，其中的各个点指示对hr滤波器进行测量的位置。
6.已经针对hr滤波器的角度插值开发了多个不同的插值方案。通常，根据在球上处的测量结果估计m个hr滤波器对其中，r表示右耳，l表示左耳，表示仰角，表示方位角。任务是找到函数其中其中其在未经采样的角度处提供传递具有良好感知准确度的声音渲染的左滤波器和右滤波器。一旦获得就可以在由指定的任何任意位置生成左耳和右耳hr滤波器。要注意，为了简单而不造成混淆有时省略了上标l或r。
7.用于hrtf角度插值的两个主要的方法如下：
8.(1)本地邻居法：通常采用的方法是线性插值，其中，通过在其周围最近位置测量的hrtf的贡献进行加权而推断缺失的hrtf。在插值之前可以对hrtf进行预处理，例如首先将在两个或更多个最近位置测量的hrtf转换成最小相位，然后应用线性插值。
9.(2)变化的方法：更复杂的数据驱动的方法是将已测量的hrtf线性变换到由基函数集合定义的另一个空间，其中，基函数的一个集合覆盖仰角角度维度和方位角角度维度并且另一个集合覆盖频率维度。可以通过对已测量的hrtf的协方差矩阵进行特征分解来获得基函数[1,2]。在[3]中，已经将完备并且在2d球体上正交的球谐函数(sh)用于覆盖仰角角度维度和方位角角度维度，并且已经将复指数函数用于覆盖频率维度。基于sh的hrtf模型在模型的平均均方误差(mse)和感知响度稳定性方面产生了令人鼓舞的性能水平[4]。


技术实现要素：

[0010]
准确和高效地渲染空间位置的声源的能力是基于hr滤波器的空间声音渲染器的主要特征之一。渲染器中使用的hr滤波器集合的空间分辨率确定经渲染的声源的空间分辨率。使用在2d球体上进行粗采样的hr滤波器集合，vr/ar/mr/xr用户经常报告移动声音的空间不连续性。这些空间不连续性导致音频视频同步误差，其明显地降低了沉浸感。使用在球体上经精细采样的hr滤波器集合是一个解决方案。然而，根据在满足maa要求的精细网格上的输入-输出测量结果估计hr滤波器集合可能非常耗时并且对于对象和实验人员两者而言很枯燥。备选地，更高效的是在给定经稀疏采样的hr滤波器数据集时推断与缺失的hr滤波器有关的空间相关信息。
[0011]
最近邻hr滤波器插值法假设在每个采样位置hr滤波器仅最多影响某个有限距离的区域。然后将未采样位置的hr滤波器近似为在某个截止距离之内位置的hr滤波器的加权平均，或基于直线2d网格上的指定数量的最近点对未采样位置的hr滤波器进行近似，例如，其中，是在未采样位置处估计的hr滤波器向量，并且并且这个方法简单，并且计算复杂度低，其可以导致高效的实现。然而，插值准确度可能不足以产生令人信服的空间声音场景。这仅仅是由于样本点之间的情况的变化比滤波器的加权平均可以产生的更复杂的事实。
[0012]
所述变化的方法将hr滤波器表示为基函数集合的线性组合(即其中，ω
p
是第p基函数的系数)。无论基函数是什么，系数通常都是通过使已测量的点集合上的平方估计误差之和最小化而获得的最小二乘估计，即给定基函数集合，认为所述系数是在求解二次最小问题意义下的“最佳”拟合。在原理上，对基函数的选择没有限制。然而，实际上，实际选择能够在估计准确度方面有效地表示以及在基函数的数量和基函数的复杂度方面有效地表示hr滤波器集合的基函数集合。
[0013]
关于对hrtf幅度响应进行建模的早期工作使用主成分(pc)作为基函数，其中的pc是通过对根据10个收听者在265个源位置测量的hrtf幅度响应的协方差矩阵进行特征分解获得的。仅使用五个pc，所产生的模型达到了原始数据库中的方差的接近90％。这个模型是高效的。其很好地表示了原始数据集同时无需在缺失的位置对hrtf进行插值的机制。最近，提出了将主成分分析(pca)与最近邻方法组合的混合方法，其中通过偏导数来近似模型系数。然而，该混合方法仅实现了与基于最近邻的双线性插值类似的结果。
[0014]
已经将sh用于对hrtf集合的角度相关性进行建模。所产生的模型在模型的平均均方误差(mse)方面产生了令人鼓舞的性能水平。然而，与基于特征分解的模型中的基函数(其是固定的pc向量)不同，sh基函数复杂并且评估代价高。将p次q阶的sh函数写为是相关联的勒让德多项式，其实质上是p次三角多项式。对于整个模型，需要评估最高为p阶的(p 1)2个sh。
[0015]
为了实现高空间分辨率，sh表示的阶应该尽可能的高。已经在所感知的空间响度稳定性的背景下研究了sh阶对空间混叠的影响，所感知的空间响度稳定性被定义为在不同的头部定向上感知的经渲染的声音场景的响度有多稳定。主观结果显示需要将高阶(p》10)sh hrtf表示用于促进高质量动态虚拟声音场景。这导致了l(p 1)2＝15488个系数，其中，l＝128与频率仓(bin)的数量相对应。另一个研究还利用复指数对hrtf频率部分进行了建模，并且系数的总数量是l(p 1)2，其中，l是频率部分表示的截短数量。结果显示，为了在mse方面表示高达20khz的整个频率范围上的hrtf，sh的阶需要达到与p＝30一样大，并且经截短的频率部分是l＝40。系数的数量是38440。使用这样的高阶sh hrtf模型来评估hrtf基本上不可能在实时vr/ar/mr/xr系统中实现。
[0016]
本公开提供用于以下的过程：在空间中的任何任意位置生成对于实时vr/ar/mr/xr系统足够准确和高效的hr滤波器。在一个实施例中，采用变化的方法，其中，利用b样条基函数对hr滤波器集合的空间变化进行建模，以及根据时域fir滤波器或频域上的某种映射对滤波器进行参数化，其中，dft是这样的一个映射。所产生的模型在mse测量和感知评估方面是准确的。它在如下方面是高效的：基于该模型评估hr滤波器所需要的基函数的总数量和计算量远小于使用球谐函数或其他这样的复杂基函数的模型的量。
[0017]
因此，在一个方面，提供了一种用于声音信号滤波的方法。所述方法包括：生成用于由仰角角度和方位角角度指定的特定位置的滤波器对，所述滤波器对由右滤波器和左滤波器构成。所述方法还包括：使用所述右滤波器对声音信号滤波以及使用所述左滤波器对所述声音信号滤波。生成所述滤波器对包括：i)在所述仰角角度处至少获得第一仰角基函数值集合；ii)在所述方位角角度处至少获得第一方位角基函数值集合；iii)使用以下生成所述右滤波器：a)至少所述第一仰角基函数值集合，b)至少所述第一方位角基函数值集合，以及c)右滤波器模型参数；以及iv)使用以下生成所述左滤波器：a)至少所述第一仰角基函数值集合，b)至少所述第一方位角基函数值集合，以及c)左滤波器模型参数。
[0018]
在另一个方面，提供了一种用于声音信号滤波的滤波装置。所述滤波装置被适配为执行方法，所述方法包括：生成由仰角角度和方位角角度指定的特定位置的滤波器对，所述滤波器对由右滤波器和左滤波器构成。所述方法还包括：使用所述右滤波器对声音信号滤波以及使用所述左滤波器对所述声音信号滤波。生成所述滤波器对包括：i)在所述仰角角度处至少获得第一仰角基函数值集合；ii)在所述方位角角度处至少获得第一方位角基函数值集合；iii)使用以下生成所述右滤波器：a)至少所述第一仰角基函数值集合，b)至少所述第一方位角基函数值集合，以及c)右滤波器模型参数；以及iv)使用以下生成所述左滤波器：a)至少所述第一仰角基函数值集合，b)至少所述第一方位角基函数值集合，以及c)左滤波器模型参数。
[0019]
所提出的过程的主要优点包括：a)比基于双线性pc的解更准确，b)比基于sh的解更高效，c)建立模型无需对hr滤波器数据库密集采样，以及d)相对于原始hr滤波器数据库，该模型在存储器中占用更少的空间。以上优点使所提出的实施例对实时vr/ar/mr/xr系统具有吸引力。
附图说明
[0020]
被并入本文且形成说明书的部分的附图示出各种实施例。
[0021]
图1示出了从由球坐标系中的仰角角度和方位角角度对指定的到达方向(doa)朝着收听者传播的声波。
[0022]
图2示出了朝着收听者传播的声波的itd和频谱线索的示例。
[0023]
图3示出了在2d球体上对网格进行采样的示例。
[0024]
图4示出了根据实施例的hr滤波单元。
[0025]
图5是示出hr滤波器建模的一个实施例的流程图。
[0026]
图6是根据实施例的描述用于获得零时延hr滤波器和itd的预处理的过程的流程
图。
[0027]
图7a示出了具有0度仰角和从0度到360度方位角的水平平面上的右耳hrtf(实线)和左耳hrtf(虚线)的延迟估计。
[0028]
图7b示出了90度方位角处的对应的右耳hrtf(实线)和左耳hrtf(虚线)。
[0029]
图8描绘了根据实施例的建模过程的框图。
[0030]
图9示出了b样条基函数的示例。
[0031]
图10示出了周期性基函数的示例。
[0032]
图11示出了根据实施例的过程。
[0033]
图12示出了周期性b样条基函数的示例。
[0034]
图13a示出了b样条基函数的示例。
[0035]
图13b示出了标准b样条基函数的示例。
[0036]
图14a示出了b样条基函数的另一个示例。
[0037]
图14b示出了没有在节点(knot-point)0/180度处的平滑条件的标准b样条基函数。
[0038]
图15示出了根据一个实施例的hr滤波器数据集的模型表示。
[0039]
图16是根据一个实施例的用于生成零时延hr滤波器对和对应itd的系统的框图。
[0040]
图17示出了根据一个实施例的用于在给定hr滤波器模型表示时在位置生成零时延hr滤波器对的过程。
[0041]
图18示出了根据一个实施例的用于在给定itd模型表示时在位置生成itd的过程。
[0042]
图19是示出根据实施例的过程的流程图。
[0043]
图20是示出根据实施例的过程的流程图。
[0044]
图21是根据一个实施例的hr滤波装置2100的框图。
具体实施方式
[0045]
图4示出了根据实施例的hr滤波单元400。hr滤波单元400包括渲染单元402。单元400还包括：在由渲染单元402实时请求的任意仰角和方位角角度处分别生成hr滤波器以及itd的hr滤波器生成器404以及itd生成器406。这需要基于已经加载进单元400的hr滤波器模型高效地评估左右hr滤波器对。这还需要基于已经加载进该单元的itd模型高效地评估itd。这个hr滤波单元400因此将具有用于从这些模型的数据库410加载hr滤波器模型和itd模型的接口408。hr滤波器模型的数据库通过估计不同的hr滤波器数据库的hr滤波器模型而离线地生成。
[0046]
1、hr滤波器集合建模
[0047]
如上所述，hr滤波器是包括itd、ild和频谱线索的角度相关的空间线索的数学表示。itd被定义为声音信号到达两耳的时间上的差异，如图2所示。我们从hr滤波器中移除与频率不相关的时间延迟，并将其单独保留为每对hr滤波器的纯延迟。其余的零时延hr滤波器包含耳间相位差(ipd)、ild和频谱线索。滤波器和itd根据方位角和仰角分别建模。
[0048]
使用大多数都是公开可用的现有hr滤波器数据库来执行变化的方法。这些数据库
中的hr滤波器是基于在不同的空间采样网格上完成的声音测量结果估计的，并且通常将它们以不同的文件格式进行存储，其对于提供数据库的实验室自然是有利的。最近，针对具有一致定义的自描述数据开发了面向空间的声学格式(sofa)的格式，其统一了不同的hr滤波器数据库的表示。因此，在一个实施例中，采用了sofa格式，因此无需在建模之前用于交换数据格式所需的额外的努力。可以在www.sofaconventions.org/mediawiki/index.php处找到关于sofa格式的更多信息。
[0049]
图5描述了hr滤波器建模的一个实施例的流程图，其中，经由sofa api加载sofa格式的hr滤波器集合。在预处理单元中，如果未在原始数据库中提供与频率不相关的时延，则针对每个hr滤波器估计这些信息。然后将hr滤波器分成零时延hr滤波器和itd。最后，在对单元进行建模时，将零时延hr滤波器和itd分别建模为仰角角度和方位角角度的连续基函数的线性求和。
[0050]
在前三个小节中更详细地描述了预处理、hr滤波器模型估计和itd模型估计的步骤。在其之后给出了整个模型表示的描述。
[0051]
1.1预处理
[0052]
用于基于测量结果估计hr滤波器集合的基本过程包括以下步骤：
[0053]
(1)经由放置在指定的仰角方位角以及离对象的头部固定距离处的扬声器发射已知信号；
[0054]
(2)使用放置在对象的耳道中或其入口处的麦克风记录对象的左耳信号和右耳信号；
[0055]
(3)对已记录的原始数据进行后处理，其主要用于移除测量系统的响应；以及
[0056]
(4)将已知的扬声器信号用作输入信号以及将预处理的耳信号用作输出信号，基于预处理的数据来估计作为线性动态系统的冲激响应的hr滤波器。
[0057]
在冲激响应开始(onset)之前通常存在与频率不相关的延迟。一些数据库(例如，cipic数据库)提供开始信息。然而，大多数数据库不提供这样的信息。如上文提及的，hr滤波器集合可以被建模为类似于最小相位的系统和纯延迟线的组合。在这个情况下需要延迟估计。给定延迟信息，通过从右耳hr滤波器的延迟减去左耳hr滤波器的延迟简单地计算itd。其次，通过对hr滤波器加窗移除延迟并获得零时延hr滤波器。在图6中示出了描述用于获得零时延hr滤波器和itd的预处理过程的流程图。
[0058]
在hrir的时间结构中，很容易观察到在开始(onset)出现之后表现出幅度的突然的增大。基于这个时间特征，用于估计延迟的一个方法是使用遵从冲激响应(ir)的能量包络的开始(onset)检测函数。这种开始检测函数可以被构造为其中，{w(l):l＝1,
…
,l}是l个样本长的加窗函数并且r是两个窗之间的以样本为单位的时间步长。为了不导致模糊，为了简单在此处省略了角度参数和耳的表示。可以将覆盖hrir的总能量的90％的部分的长度选择为窗的长度l。当hrir中存在强冲击瞬态时，以上解决方案产生了满意的结果。然而，不总是这个情况，因此通过使用累积能量与总能量之比n＝1,
···
n优化所述解决方案，其中，n是hrir的长度。累积能量被定义为其中，w(l)是n个点的窗。总能量是
另一种优化采用比的导数，并且找到当导数超过某个阈值时的第一个样本的索引作为开始(onset)的索引。以样本为单位的时延τ
td
可以写为其中η是阈值。通常，相对于对侧hrtf，同侧hrtf的阈值较高。图7a和图7b示出使用普林斯顿hrtf数据集——对象id 27(数据库的url是www.princeton.edu/3d3a/hrtfmeasurements.html)估计的hrtf的延迟的示例。图7a中的曲线示出了具有0度仰角和从0度到360度方位角的水平平面上的右耳hrtf(实线)和左耳hrtf(虚线)的延迟估计。在数据提示中示出了hrtf在90度方位角处的延迟。在图7b中示出了在90度方位角处对应的右耳hrtf(实线)和左耳hrtf(虚线)。星号强调了已检测到的开始(onset)。
[0059]
当给定延迟估计时，可以通过对原始hr滤波器加窗而获得零时延hr滤波器。已知对hr滤波器的频谱内容的最明显的位置相关的影响可以追踪到外耳或耳廓，其持续约0.3毫秒。“肩部弹跳”影响稍后出现。位置相关的ir的总长度通常不会超过1毫秒。因此，1毫秒矩形窗对于保留主要的与频谱有关的线索而言足够长。如果没有增加另外的位置相关的信息，则可以不需要更长的窗。
[0060]
1.2 hr滤波器模型估计
[0061]
右耳和左耳的hr滤波器被分别建模。在下文以两个可能的扩展形式(仰角扩展形式和方位角扩展形式)给出了长度n的hr滤波器的一般时域截短(td)fir模型，其中针对仰角和方位角使用单独的基函数。
[0062][0063]
在仰角扩展形式下，存在针对仰角维度的单个基函数集合{θ
p
:p＝1,
…
,p}和针对方位角维度的p个基函数集合，每个集合针对一个仰角索引p，{φ
p,q
:q＝1,
…
,q
p
}。k《n是滤波器参数向量的n维向量空间的基向量的数量，并且ek是长度为n的标准正交基向量：
[0064][0065]
α＝{α
p,q,k
:p＝1,
…
,p；q＝1,
…
,q
p
；k＝1,
…
,k}是需要估计的模型参数集合。因为fir模型被截短，k《n，所以hr滤波器模型值等于0。
[0066]
方位角扩展形式是仰角扩展形式的镜像形式，具有对应的镜像术语。从现在开始
我们将示出仰角扩展形式的性质。这些性质对于方位角扩展形式在镜像意义上也成立并且本领域普通技术人员可以基于仰角扩展形式的这些性质推导出那些镜像性质。
[0067]
仰角扩展形式非常灵活之处在于其支持针对每个仰角索引p的方位角基函数的各个集合。不是始终需要这种全面的灵活性，但是使用多于一个方位角基函数集合的确是好主意。在分别在收听者正上方和正下方的 /-90度的仰角角度处，在不同的方位角角度处的hr滤波器都一样。这可以通过对仰角索引p使用等于1的单个方位角基函数来处理，该仰角索引p具有对仰角 /-90度有贡献的基函数。其他仰角索引可以共享单个不同的方位角基函数集合，其中基函数的数量q》1，或者共享经仔细选择的方位角基函数的几个集合，以捕获所建模的滤波器集合的仰角-方位角变化。
[0068]
在下文我们将得出一般仰角扩展形式的性质。然而，本领域技术人员将清楚，当方位角基函数的不同集合的数量小于p时如何修改这些性质。
[0069]
为了估计模型参数{α
p,q,k
}，需要两个东西。
[0070]
(1)需要指定最小化准则，其通常采用时域、频域或两者的组合上的建模误差的测量的形式，并且这个准则甚至可以包括：用于减少过拟合所建模的数据的趋势的正则化项。
[0071]
(2)用于估计使最小化准则最小化的参数的优化方法。
[0072]
图8描绘了给定与对应的仰角角度和方位角角度(即)相关联的零时延hr滤波器集合时的建模过程的框图。给定仰角和方位角的列表，分别构造仰角角度和方位角角度上的基函数。然后采用最小二乘方法来估计模型参数。
[0073]
在小节1.2.1中更详细地描述了这个模型估计过程。
[0074]
在这个阶段，模型规定是相当通用的，因为还没有规定基函数的两个集合和获得可以输出准确的建模并且高效地评估hr滤波器模型的关键在于基函数的这两个集合的选择。在使用不同类型的函数试验之后，我们选择将我们称为周期性b样条基函数的基函数用作方位角基函数以及将标准b样条函数用作仰角基函数。在小节1.2.2中更详细地说明已选择的这些基函数。
[0075]
1.2.1模型参数估计
[0076]
给定基函数集合和以及在m个不同的角度位置上采样的右耳或左耳的零时延hr滤波器集合时域上典型的最小化准则是m个hr滤波器构成的集合(右耳或左耳)上的建模误差的范数之和：
[0077][0078]
其中，
[0079]
并且
[0080][0081]
所估计的参数的数量是其应该远少于可用数据样本的数量m*n，以避免不确定系统。
[0082]
因为标准正交基向量ek是正交向量，所以参数αk＝{α
p,q,k
:p＝1,
…
,p；q＝1,
…
,q
p
}可以针对每个样本k独立求解。对于每个样本k，最小化准则变为：
[0083][0084]
其可以用矩阵形式表示为：
[0085][0086]
其中，
[0087][0088]
j(αk)是线性最小二乘准则。通过标准方程的求解αk＝(b
t
b)-1bthk
获得使j(αk)最小化的解。然而，直接使以上成本函数最小化导致对线性系统的准确解。该解对数据中的噪声敏感并且可以导致过拟合。然后应用tikhonov正则化，并且最小化准则变为：
[0089][0090]
其中i是尺寸为的单位矩阵并且0是具有个元素的0列向量。
[0091]
j(αk)也是线性最小二乘准则。类似地，通过对标准方程的求解获得使最小化的解，其中，可以将λ的值确定为使矩阵的条件数小于10或导致良好模型准确度的其他某个值。
[0092]
为了更好的数值准确度，借助于的奇异值分解(svd)对αk进行实际求解：
basics(b样条基础知识)”(ftp://ftp.cs.wisc.edu/approx/bsplbasic.pdf)中找到这个过程的细节。
[0103]
在图9中示出了使用节点序列θ＝(-90,-60,-30,0,30,60,90)和多重序列m＝(4,1,1,1,1,1,3)评估的针对3度的仰角角度的b样条基函数的示例。
[0104]
在度的方位角角度在空间上是与(对于任意整数值κ)度的方位角角度相同的点的意义上方位角角度是周期性的(例如，循环的)，并且为了在方位角维度上获得高效的建模，以同样方式使用周期性的基函数(即)很重要。在图10中示出了这样的周期性的基函数的示例，其中，使用实线绘制了函数在从0到360的角度范围的部分并且使用虚线绘制了函数在该范围之外的部分。
[0105]
我们已经设计了用于生成0到360度的方位角范围上的周期性b样条基函数集合的方法。在图11中示出了该方法，其包括以下步骤。
[0106]
(步骤1)指定从0到360度范围上的节点序列。将该节点序列的长度表示为l。
[0107]
(步骤2)使用小于0度的j个值和大于360度的j-1值以周期性方式扩展该节点序列。
[0108]
(步骤3)使用该扩展的节点序列和该节点的扩展的多重序列来生成扩展的b样条基函数集合，使用标准方法生成b样条函数集合。
[0109]
(步骤4)选择扩展的基函数中自索引2开始的l-1个连续的扩展的基函数，以及以周期性方式将它们映射到0到360度的方位角范围。
[0110]
这个方法提供了0到360度范围上的l-1个周期性的基函数集合。
[0111]
方位角角度上的每个基函数也是j-1次多项式函数，并且被写为：
[0112][0113]
图12中示出了使用长度为l＝11的节点序列φ＝(0,30,70,110,150,180,210,250,290,330,360)进行评估的针对3度的方位角角度的周期性b样条基函数的示例。
[0114]
1.3 itd模型估计
[0115]
itd模型的一般形式由下式给出：
[0116][0117]
和是分别针对仰角角度和方位角角度的b样条基函数。是模型参数集合。
[0118]
1.3.1模型参数估计
[0119]
通过使最小二乘准则最小化而获得模型参数{c
p
′
,q
′
}，
[0120][0121]
其中，
[0122][0123]
是针对的itd。是由原始数据库提供的或使用小节1.1中描述的方法所估计的与频率不相关的时延。
[0124]
应用tikhonov正则化以避免过拟合，最小化准则变为：
[0125][0126]
其中，是尺寸为的单位矩阵并且0是具有个元素的0列向量。
[0127]
可以将的值确定为使矩阵的条件数小于10或导致良好模型准确度的其他某个值。如小节1.2.1中所述，借助于时的svd，通过c＝v
′s′-1u′
t
τ获得模型参数，其是具有个元素的列向量。
[0128]
1.3.2.仰角和方位角基函数的规定
[0129]
当仰角从-90度向上移动到90度时，itd从零增大到在仰角0度处的最大值并且随后减小到零。基于此，很自然地使用在 /-90度仰角处为零的基函数。这个要求相当于在 /-90度仰角处的至少一个平滑条件。如在小节1.2.2中所说明的，函数在节点处的平滑由多重序列m控制。将每个基函数写为：
[0130][0131]
在图13a中示出了使用节点序列和多重序列m＝(3,1,1,1,2)评估的针对3度的仰角角度的b样条基函数的示例。
[0132]
考虑到由于测量设置中的不对称以及对象而导致itd在 /-90度仰角处可能不正
好为零，所以使用标准b样条基函数而无需在节点 /-90度处的平滑条件仍然是好的选择。在图13b中示出了使用节点序列θ＝(-90,-45,0,45,90)和多重序列m＝(4,1,1,1,3)评估的针对3度的仰角角度的标准b样条基函数的示例。
[0133]
当方位角沿着圆周移动时，itd的改变呈类似正弦的形状，其中，零itd出现在0/180/360度的方位角处并且最大itd出现在90/270度的方位角处。类似地，在0/180/360度的方位角处可以满足一个平滑条件。此外，可以认为在180度与360度之间的方位角处的itd是在0与180度之间的方位角处的itd的镜像。因此，我们将基函数的一个集合用于两个区间[0,180]和[180,360]的方位角角度。将每个基函数写为：
[0134][0135]
在图14a中示出了最初使用节点序列和多重序列m＝(3,1,
···
,1,2)评估的针对3度的方位角角度的b样条基函数的示例。
[0136]
考虑到itd在0/180/360度的方位角处可能不正好是零，可以使用标准b样条基函数而无需在节点0/180度处的平滑条件。在图14b中示出了这样的基函数的示例。
[0137]
1.4模型表示
[0138]
图15示出了hr滤波器数据集的模型表示。所述表示包括：均包括基函数和模型参数的一个零时延hr滤波器模型表示和一个itd模型表示。建模解决方案的建模准确度和计算效率的关键是仔细地构造用于对hr滤波器集合的角度变化进行建模的b样条基函数集合，其足够简单地给出了良好计算效率并且足够丰富地给出了良好建模准确度。
[0139]
对于零延迟hr滤波器模型，存在：p个仰角b样条基函数；均包含q
p
个函数的方位角b样条基函数的p个集合；以及均为乘k的矩阵的模型参数的两个集合。对于itd模型，存在：个仰角b样条基函数；均包含个函数的方位角b样条基函数的个集合；以及作为具有个元素的向量的模型参数的一个集合。
[0140]
b样条基函数的每个集合由节点序列和作为三维数组的多项式模型系数γ表示。第一维度与b样条的阶相对应，第二维度与节点区间的数量相对应，并且第三维度与基函数的数量相对应。
[0141]
p或远小于原始hr滤波器数据集中的仰角角度的数量。q或远小于数据集中的方位角角度的数量。k也小于原始滤波器的频率仓的长度或数量。因此，该模型表示在表示hr滤波器数据集方面是高效的。
[0142]
此外，因为角度基函数连续，所以模型表示可以用于在由仰角和方位角指定的任何任意的位置生成hr滤波器对。
[0143]
2、hr滤波器生成
[0144]
图16是用于在给定模型表示时生成零时延hr滤波器对(即右耳滤波器和左耳滤波器)和对应的itd的系统的框图。所述模型表示可以写到二进制文件或文本文件中。经由api加载文件以获取模型结构。下面描述如何使用该模型表示在指定的位置获得hr滤波器对和itd。
[0145]
2.1生成零时延hr滤波器
[0146]
图17示出了在给定hr滤波器模型表示时在位置生成零时延hr滤波器对的过程。如1.2.2节中所说明的，仰角b样条基函数集合的模型{θ
p
:p＝1,
···
,p}包括：节点序列θ＝(θ1,
…
,θu)，其指定子区间所述函数是所述子区间上的多项式；以及指示模型参数的3维数组述子区间上的多项式；以及指示模型参数的3维数组涉及评估p个仰角基函数在仰角角度处的值的步骤如下：
[0147]
(1)找到满足的索引u；以及
[0148]
(2)根据下式评估第p仰角b样条基函数在仰角角度处的值：
[0149][0150]
将类似的过程用于评估方位角b样条基函数集合在给定方位角角度处的值
[0151]
一旦获得基函数值的这些集合，就根据下式获得在位置处的右耳零时延hr滤波器：
[0152][0153]
基于此，根据下式获得对的评估也是清楚的：
[0154][0155]
根据下式获得在位置处的左耳零时延hr滤波器：
[0156][0157]
根据下式获得对的评估：
[0158][0159]
2.2生成itd
[0160]
图18示出了根据一个实施例的用于在给定itd模型表示时在位置处生成itd的过程。
[0161]
在小节2.1中描述的过程之后，分别评估个仰角基函数在仰角角度处的值和方位角b样条基函数集合在给定方位角角度处的值一旦对仰角基函数和方位角基函数的值进行了评估，就根据下式获得itd：
[0162][0163]
在小节5.1.1中所提及的，我们将hr滤波器集合建模为类似于最小相位的系统和纯延迟线的组合。右耳hr滤波器的延迟是：
[0164][0165]
左耳hr滤波器的延迟是：
[0166][0167]
要注意，和的计算应该与所使用的itd的定义和坐标系一致。
[0168]
图19是示出根据实施例的过程1900的流程图。过程1900可以从步骤s1902开始。
[0169]
步骤s1902包括：生成用于由仰角角度和方位角角度指定的特定位置的滤波器对，所述滤波器对由右滤波器和左滤波器构成。
[0170]
步骤s1904包括：使用右滤波器对声音信号滤波。
[0171]
步骤s1906包括：使用左滤波器对声音信号滤波。
[0172]
如图20所示，步骤s1902包括：i)在所述仰角角度处至少获得第一仰角基函数值集合(步骤s2002)；ii)在所述方位角角度处至少获得第一方位角基函数值集合(步骤s2004)；iii)使用以下生成所述右滤波器：a)至少所述第一仰角基函数值集合，b)至少所述第一方位角基函数值集合，以及c)右滤波器模型参数(步骤s2006)；以及iv)使用以下生成所述左滤波器：a)至少所述第一仰角基函数值集合，b)至少所述第一方位角基函数值集合，以及c)左滤波器模型参数(步骤s2008)。
[0173]
在一些实施例中，获得所述第一方位角基函数值集合包括：获得方位角基函数值的p个集合，其中，所述方位角基函数值的p个集合包括所述第一方位角基函数值集合以及以及其中，(p＝1至p、q＝1至q
p
、以及k＝1至k)是左模型参数集合，(p＝1至p、q＝1至q
p
、以及k＝1至k)是右模型参数集合，(p＝1至p)定义在仰角角度处的第一仰角基函数值集合，并且(p＝1至p、以及q＝
1至q
p
)定义在方位角角度处的方位角基函数值的p个集合；并且ek(k＝1至k)是长度为n的标准正交基向量集合。
[0174]
在一些实施例中，获得所述第一仰角基函数值集合包括：获得仰角基函数值的q个集合，其中，所述仰角基函数值的q个集合包括第一仰角基函数值集合以及其中，(p＝1至pq、q＝1至q、以及k＝1至k)是左模型参数集合，(p＝1至pq、q＝1至q、以及k＝1至k)是右模型参数集合，(q＝1至q、以及p＝1至pq)定义在仰角角度处的所述仰角基函数值的q个集合，并且(q＝1至q)定义方位角角度处的第一方位角基函数值集合；并且ek(k＝1至k)是长度为n的标准正交基向量集合。
[0175]
在一些实施例中，获得第一仰角基函数值集合包括：针对第一仰角基函数集合中包括的每个仰角基函数，在所述仰角角度处评估该仰角基函数，以产生与所述仰角角度和该仰角基函数相对应的仰角基函数值，并且获得所述第一方位角基函数值集合包括：针对第一方位角基函数集合中包括的每个方位角基函数，在所述方位角角度处评估该方位角基函数，以产生与所述方位角角度和该方位角基函数相对应的方位角基函数值。
[0176]
在一些实施例中，所述第一仰角基函数集合中包括的每个仰角基函数是b样条基函数，并且所述第一方位角基函数集合中包括的每个方位角基函数是周期性b样条基函数。
[0177]
在一些实施例中，所述过程还包括：获得至少表示所述第一仰角基函数集合的模型，其中，所述模型包括：序列(θ)，其中θ＝(θ1,
…
,θu)，其指定子区间所述仰角基函数是所述子区间上的多项式；以及模型参数的三维数组型参数的三维数组
[0178]
在一些实施例中，所述第一仰角基函数集合包括第p仰角基函数，在所述仰角角度处评估所述第一仰角基函数集合中包括的每个仰角基函数包括：在所述仰角角度处评估第p仰角基函数，并且在所述仰角角度评估第p仰角基函数包括以下步骤：找到满足的索引u；以及根据评估所述第p仰角基函数在仰角角度处的值。
[0179]
在一些实施例中，所述过程还包括：获得至少表示所述第一方位角基函数集合的模型，其中，所述模型包括：序列(φ1)，其中其指定子区间其指定子区间所述方位角基函数是所述子区间上的多项式；以及模型参数的三维数组
[0180]
在一些实施例中，所述第一方位角基函数集合包括第q方位角基函数，在方位角角度处评估所述第一方位角基函数集合中包括的每个方位角基函数包括：在方位角角度处评估所述第q方位角基函数，并且在方位角角度处评估所述第q方位角基函数包括以下
2141包括计算机可读介质(crm)2142，计算机可读介质(crm)2142存储计算机程序(cp)2143，计算机程序(cp)2143包括计算机可读指令(cri)2144。crm 2142可以是非暂时性计算机可读介质(例如，磁介质(例如硬盘)、光介质、存储器设备(例如随机存取存储器、闪存)等)。在一些实施例中，计算机程序2143的cri 2144被配置为使得：当由pc 2102执行时，cri导致装置2100执行本文描述的步骤(例如，参考流程图在本文描述的步骤)。在其他实施例中，装置2100可被配置为在不需要代码的情况下执行本文描述的步骤。也即是说，例如，pc 2102可以仅由一个或多个asic构成。因此，本文描述的实施例的特征可以在硬件和/或软件中实现。
[0187]
下面是本文描述的各种实施例的总结：
[0188]
a1、一种用于声音信号滤波的方法，所述方法包括：生成用于由仰角角度和方位角角度指定的特定位置的滤波器对，所述滤波器对由右滤波器和左滤波器构成；使用所述右滤波器对声音信号滤波；以及使用所述左滤波器对所述声音信号滤波，其中，生成所述滤波器对包括：i)在所述仰角角度处至少获得第一仰角基函数值集合；ii)在所述方位角角度处至少获得第一方位角基函数值集合；iii)使用以下生成所述右滤波器：a)至少所述第一仰角基函数值集合，b)至少所述第一方位角基函数值集合，以及c)右滤波器模型参数；以及iv)使用以下生成所述左滤波器：a)至少所述第一仰角基函数值集合，b)至少所述第一方位角基函数值集合，以及c)左滤波器模型参数。
[0189]
a2、根据权利要求a1所述的方法，其中，获得所述第一方位角基函数值集合包括：获得方位角基函数值的p个集合，其中所述方位角基函数值的p个集合包括所述第一方位角基函数值集合。
[0190]
a3、根据权利要求a1所述的方法，其中，生成所述右滤波器包括计算并且生成所述左滤波器包括计算其中，(p＝1至p、q＝1至q
p
、以及k＝1至k)是右模型参数集合，(p＝1至p、q＝1至q
p
、以及k＝1至k)是左模型参数集合，(p＝1至p)定义在仰角角度处的所述第一仰角基函数值集合，并且(p＝1至p、以及q＝1至q
p
)定义在所述方位角角度处的方位角基函数值的p个集合，并且ek(k＝1至k)是长度为n的标准正交基向量集合。
[0191]
a4、根据权利要求a1所述的方法，其中，获得所述第一仰角基函数值集合包括：获得仰角基函数值的q个集合，其中所述仰角基函数值的q个集合包括所述第一仰角基函数值集合。
[0192]
a5、根据权利要求a1所述的方法，其中，生成所述右滤波器包括计算并且生成所述左滤波器包括计算其中，(p＝1至pq、q＝1
至q、以及k＝1至k)是右模型参数集合，(p＝1至pq、q＝1至q、以及k＝1至k)是左模型参数集合，(q＝1至q、以及p＝1至pq)定义在仰角角度处的仰角基函数值的q个集合，并且(q＝1至q)定义在方位角角度处的所述第一方位角基函数值集合，并且ek(k＝1至k)是长度为n的标准正交基向量集合。
[0193]
a6、根据权利要求a1至a5中任一项所述的方法，其中，每个所述仰角基函数值取决于所述方位角角度，和/或每个所述方位角基函数值取决于所述仰角角度。
[0194]
a7、根据权利要求a1至a5中任一项所述的方法，其中，获得所述第一仰角基函数值集合包括：针对第一仰角基函数集合中包括的每个仰角基函数，在所述仰角角度处评估该仰角基函数，以产生与所述仰角角度和该仰角基函数相对应的仰角基函数值，并且获得所述第一方位角基函数值集合包括：针对第一方位角基函数集合中包括的每个方位角基函数，在所述方位角角度处评估该方位角基函数，以产生与所述方位角角度和该方位角基函数相对应的方位角基函数值。
[0195]
a8、根据权利要求a7所述的方法，其中，所述第一仰角基函数集合中包括的每个仰角基函数是b样条基函数，并且所述第一方位角基函数集合中包括的每个方位角基函数是周期性b样条基函数。
[0196]
a9、根据权利要求a7或a8所述的方法，还包括：获得至少表示所述第一仰角基函数集合的模型，其中，所述模型包括：序列(θ)，其中θ＝(θ1,
…
,θu)，其指定子区间所述仰角基函数是所述子区间上的多项式；以及模型参数的三维数组
[0197]
a10、根据权利要求a9所述的方法，其中，所述第一仰角基函数集合包括第p仰角基函数，在所述仰角角度处评估所述第一仰角基函数集合中包括的每个仰角基函数包括：在所述仰角角度处评估所述第p仰角基函数，并且在所述仰角角度处评估所述第p仰角基函数包括以下步骤：找到满足的索引u；以及根据评估所述第p仰角基函数在该仰角角度处的值。
[0198]
a11、根据权利要求a7或a8所述的方法，还包括：获得至少表示所述第一方位角基函数集合的模型，其中，所述模型包括：序列(φ1)，其中其指定子区间子区间所述方位角基函数是所述子区间上的多项式；以及模型参数的三维数组
[0199]
a12、根据权利要求a11所述的方法，其中，所述第一方位角基函数集合包括第q方位角基函数，在所述方位角角度处评估所述第一方位角基函数集合中包括的每个方位角基函数包括：在所述方位角角度处评估所述第q方位角基函数，并且在所述方位角角度处评估所述第q方位角基函数包括以下步骤：找到满足的索引l；以及根据评估所述第q方位角基函数在所述方位角角度处的值。
[0200]
a13、根据权利要求a7至a12中任一项所述的方法，其中获得所述第一方位角基函数值集合的步骤还包括生成所述第一方位角基函数集合。
[0201]
a14、根据权利要求a13所述的方法，其中，生成所述第一方位角基函数集合包括：生成0到360度方位角范围上的周期性b样条基函数集合。
[0202]
a15、根据权利要求a14所述的方法，其中，生成所述0到360度方位角范围上的周期性b样条基函数集合包括：指定0到360度范围上长度为l的节点序列；基于长度为l的节点序列生成扩展的节点序列，其中，生成所述扩展的节点序列包括：使用小于0度的j个值和大于360度的j-1个值以周期性方式扩展长度为l的节点序列；获得节点的扩展的多重序列；使用所述扩展的节点序列和所述扩展的多重序列生成扩展的b样条基函数集合；选择扩展的基函数中自索引2开始的l-1个连续的扩展的基函数；以及以周期性方式将已选择的扩展的基函数映射到0到360度的方位角范围。
[0203]
a16、根据权利要求a1至a15中任一项所述的方法，还包括：针对仰角-方位角角度确定耳间时差
[0204]
a17、根据权利要求a16所述的方法，还包括：基于确定右延迟以及基于确定左延迟
[0205]
a18、根据权利要求a17所述的方法，其中使用所述右滤波器对所述声音信号滤波包括：使用所述右滤波器和右延迟对所述声音信号滤波；并且使用所述左滤波器对所述声音信号滤波包括：使用所述左滤波器和所述左延迟对所述声音信号滤波。
[0206]
a19、根据权利要求a18所述的方法，其中，使用所述右滤波器和对所述声音信号滤波包括计算使用所述左滤波器和对所述声音信号滤波包括计算声音信号滤波包括计算其中，u(n)是声音信号。
[0207]
a20、根据权利要求a17至a19中任一项所述的方法，其中，
[0208]
并且
[0209]
a21、根据权利要求a7至a15中任一项所述的方法，其中，所述方位角基函数具有以360度为周期的周期性。
[0210]
虽然本文(包括存在的任意附录)描述了各种实施例，但应该理解，它们仅通过示例方式提出而非限制。因此，本公开的宽度和范围不应该由上述示例性实施例中的任何示例性实施例所限制。此外，按照上述要素的所有可能的变化的其任意组合被本公开包括，除非在本文另外指示或在上下文中另外明显矛盾。
[0211]
此外，虽然上文所述和在附图中示出的过程被示出为一系列步骤，但这样做仅为了说明。因此，要考虑的是，可以增加一些步骤、可以省略一些步骤、可以重新安排步骤的顺序、并且可以并行地执行一些步骤。
[0212]
缩写词：
[0213]
ar增强现实
[0214]
doa到达方向
[0215]
fir有限冲激响应
[0216]
hr头部相关的
[0217]
hrir头部相关的冲激响应
[0218]
hrtf头部相关的传递函数
[0219]
ild耳间强度差
[0220]
ipd耳间相位差
[0221]
itd耳间时间差
[0222]
itf耳间传递函数
[0223]
maa最小可听角度
[0224]
mpeg活动图像专家组
[0225]
mr混合现实
[0226]
mse均方误差
[0227]
pca主成分分析
[0228]
saoc空间声音对象编码
[0229]
sh球谐函数
[0230]
sofa面向空间的声学格式
[0231]
svd奇异值分解
[0232]
vr虚拟现实
[0233]
xr扩展现实
[0234]
参考文献：
[0235]
[1]dorisj.kistler，fredericl.wightman，“amodelofhead-relatedtransferfunctionsbasedonprincipalcomponentsanalysisandminimum-phasereconstruction(基于主成分分析和最小相位重构的头部相关传递函数的模型)”，journal of acoustical society of america(美国声学学会学报)，91(3):1637-1647，1992年3月。
[0236]
[2]f
á
bio p.freeland，luiz w.p.biscainho以及paulo s.r.diniz，"interpolation of head-related transfer functions(hrtfs):a multi-source approach(头部相关的传递函数(hrtfs)的插值：多源法)"，in 12th european signal processing conference(第12次欧洲信号处理会议)，第1761-1764页，维也纳，2004年9月。
[0237]
[3]mengqiu zhang，rodney a.kennedy以及thushara d.abhayapala，"empirical determination of frequency representation in spherical harmonics-based hrtf functional modeling(基于球谐函数的hrtf函数建模中的频率表示的经验确定)"，ieee/acm transactions on audio，speech，和language processing(ieee/acm声音、语音及语言处理学报)，第23(2)卷，第351-360页，2015年2月。
[0238]
[4]zamir ben-hur，david lou alon，boaz rafaely以及ravish mehra，"loudness stability of binaural sound with spherical harmonic representation of sparse head-related transfer functions(使用稀疏头部相关传递函数的球谐函数表示的双耳声音的响度稳定性)"，eurasip journal on audio，speech，and music processing(eurasip声音、语音及音乐处理期刊)2019，2019年5月。