一种沉浸式宽带3D声场重放方法与流程

本发明属于声学技术领域，尤其涉及一种沉浸式宽带3d声场重放方法。

背景技术：

空间声重放是三维音频技术中最为关键的一项内容，目的是尽可能准确地重现一个期望的声场景，给听音者以临场感和沉浸感。主要的实现方式是利用一组扬声器阵列在指定的空间内重构一个与期望声场一致的声学环境。现有的声场重放系统中，扬声器阵列的布局方式通常是利用环形或球形阵列，在实际摆放时需要利用球架等大型设备，给用户带来了操作的困难与体验的不便。近年来，随着虚拟现实和增强现实的发展，用户对于沉浸式体验的需求越来越迫切，迎来了沉浸式空间声场重放的巨大发展机遇。沉浸式重放根本目的是给参与者创造一种共享同一个环境氛围的感觉，将一个与用户存在地理区隔的声源利用扬声器阵列进行重放，宛如用户处于声源所在环境当中，从而使用户达到身临其境的体验。对个人用户而言，复杂的环形/球形扬声器阵列布局方式已远远不能满足现实需要。针对此问题，本发明利用基于平面阵列的扬声器布局方式设计了一种面对沉浸式宽带语音信号场景下的3d声场重放系统。

在沉浸式交互中，3d视频在视觉角度为用户提供了声源的位置信息，而声场重放在这个过程可以作为视觉信息的补充从听觉角度加深用户的临场感。而这种沉浸感依赖于音视频的同步播放，这就要求重放算法的实时性要高，延迟要低。目前，大多数精度较高的声场重放算法是以场景中听音区域内的声压值为目标声场信息，利用扬声器阵列来恢复声压信息。重建过程中通常利用柱/球面坐标的球谐展开形式进行计算，辅以高阶汉克尔、贝塞尔函数及缔合勒让德函数。因此，其计算复杂度较高，在实际场景中应用较为困难。为此，本设计采用声压匹配技术来恢复目标声场，同时引入l1与l2正则化来控制对原始声压信息的拟合程度，用以提升算法的鲁棒性。在求解正则化问题的过程中，引入交替方向乘子法降低计算复杂度，提升算法的运行速度，达到较高的实时性要求。

技术实现要素：

本发明针对于现有多声道音频重放系统中扬声器布局繁琐及算法复杂度较高的问题，提供一种沉浸式宽带3d声场重放方法，适用于沉浸式应用场景下的宽带语音信号3d声场重放，利用一组布置在墙面上的扬声器通过声压匹配的方式来重建目标声场。

本发明用于解决实际场景中对目标声源进行精确感知重放困难的问题，包括以下步骤：：

步骤1，以扬声器所在墙面的中心位置为原点，建立空间直角坐标系，其中，x轴平行于地面，y轴垂直于地面，z轴垂直于墙面。指定待重放的虚拟声源的坐标位置为s，设所使用的扬声器数最多为n个，在听音区域内选取m个位置进行听音(即，听音点个数为m)。

步骤2，构建虚拟声源辐射声场。基于声学传输函数构建全指向性的虚拟声源到听音点的声学传播特征，将该函数值作为声源在指定区域内辐射声场的声压值。本发明考虑应用对象为宽带语音信号，故将所需计算的频率上限值定为16khz，下限值定为100hz，按一倍频程将100hz～16khz划分为8个频带，取每个频带的中心频率，故选定8个单频率虚拟声源进行计算。

步骤3，构建扬声器辐射声场。设定n个扬声器的位置，分别计算各扬声器在所有听音点位置处产生的声压值，作为扬声器辐射声场的声压值集合。

步骤4，针对虚拟声源求取扬声器权值向量。通过l1范数正则化，将虚拟声源辐射声场和扬声器辐射声场加以约束，把扬声器权值(也叫扬声器激励信号或扬声器驱动信号，即，扬声器播放的信号)的求取问题转换为最小绝对收缩和选择算子问题，并利用交替方向乘子法进行求解，得到重放各虚拟单频声源(共8个)的一组扬声器权值向量。

步骤5，确定激活扬声器的数量。汇总这个8个虚拟单频声源对应的扬声器权值向量，保留其中非零项内容，作为重放系统中选定的扬声器(称作，激活的扬声器)；同时去除其中的零项，即未被选择的扬声器(称作，未激活的扬声器)。

步骤6，针对待重放声源求取激活扬声器的权值信号。对拟重放声源和所有激活的扬声器重新按照步骤2与步骤3计算待重放声源辐射声场和激活扬声器的辐射声场，并利用l2范数正则化，在均方最小准则下计算出各个扬声器的权值，在重放中作为扬声器的驱动信号。

作为优选，步骤2的实现方式为：利用格林函数来表征频率为fb的虚拟单频声源到第m个听音点的声学传输函数：

式中，kb＝2πfb/c为波数；c为空气中的声速，约为340m/s；b为频率索引；s为虚拟单频声源的位置(三维坐标)；ym为第m个听音点的位置，m＝1,2,…,m。本设计考虑宽带语音声场，选择1倍频程从100hz至16khz中8个频带的中心频率进行联合优化(每一频带的上限频率与下限频率的比值为2，起始频率为88hz)，即8个频带分别为：88hz-177hz、177hz-355hz、355hz-710hz、710hz-1420hz、1420hz-2840hz、2840hz-5680hz、5680hz-11360hz和11360hz-33720hz，对应的中心频率分别为125hz、250hz、500hz、1khz、2khz、4khz、8khz和16khz，故设定f1＝125hz,…,f8＝16khz，即，b＝1,2,…,8。

此外，本发明中的虚拟声源用于系统设计过程(即：从所有扬声器中选取激活扬声器的过程)，激活的扬声器用于声场重放，所重放的声源与虚拟声源位置相同。本发明中的虚拟声源是8个频率分别为f1,f2,…,f8的单频声源。考虑到虚拟声源具有全指向性，其辐射声场的声压值与格林函数值相等，可记为：

式中，^t代表转置运算；ym为第m个听音点的位置，m＝1,2,…,m；是频率为fb的虚拟单频声源在m个听音点处的辐射声场集合，是扬声器阵列待拟合的声场集合。

作为优选，步骤3是为了求取扬声器辐射声场的声压值集合。首先计算第n个扬声器到第m个听音点处的声学传输函数，如下：

式中，xn为扬声器阵列中第n个扬声器的坐标，n＝1,2,…,n。则扬声器到所有听音点处的声传输函数集合为：

假定待求取的扬声器权值向量(也叫，扬声器激励信号向量或扬声器驱动信号向量)为wb＝[wb(1),wb(2),l,wb(n)]^t，n是扬声器总数。此时，扬声器辐射声场的声压值集合可表示为：

作为优选，步骤4是为了求取重建各虚拟声源(单频声源)时所需扬声器的权值向量(wb)。若直接使虚拟声源辐射声场与扬声器辐射声场相等，通过矩阵求逆法求解出wb，会造成病态解和模型过拟合问题。考虑到l1正则化具有变量选择的功能，会使得模型更稀疏，即在选取少量扬声器的情况下同时保证拟合的精确度。故本设计采用l1正则化方法对虚拟声源辐射声场和扬声器辐射声场进行约束，将扬声器权值的求取问题转换为最小绝对收缩和选择算子问题，如下式：

式中，||||1代表l1范数；||||2代表l2范数；λ1＝2为惩罚因子；是使最小时，计算出的wb值。即，在频率fb下求得的扬声器权值向量。

本设计采用交替方向乘子法对上式进行求解，首先，将上式写为交替方向乘子法的形式：

minimizef(wb)+g(zb)

s.t.wb-zb＝0

式中，g(zb)＝λ1||zb||1；zb是用于约束扬声器权值向量wb的变量。求解扬声器权值wb的过程是一个迭代的过程，第j+1次迭代得到的扬声器权值变量由如下公式求得：

式中，j+1为迭代次数，j以0为起点取值；为第j次迭代的对偶变量，用于判断算法的停止条件，为第j次迭代的约束变量，迭代初值初始化为i为单位矩阵；ρ为增广拉格朗日参数，其影响算法的收敛速度，在[0.1,10]内选取。约束变量迭代更新如下式所示：

式中，sκ(α)为软阈值判决函数，定义如下：

对偶变量的迭代更新过程如下式：

经j次迭代后，主残差为对偶残差为算法迭代的停止准则为：

其中，ε^pri是主残差的上限值，ε^dual是对偶残差的上限值，由下式计算：

本设计将迭代残差设定为：绝对误差ε^abs＝10^-4，相对误差ε^rel＝10^-2。当满足停止准则后，输出fb频点对应的一组扬声器权值向量wb。共需要求出8组不同频率下的扬声器权值向量wb，b＝1,2,…,8。

作为优选，步骤5是为了汇总步骤4所求得的所有扬声器权值向量，并统计保留其中的非零项。在8个频点下所有的扬声器权值向量w可表达为：

将上述矩阵逐列相加，最终求得的汇总扬声器权值矩阵为w∑：

w∑＝[w∑(1),w∑(2),l,w∑(n)]^t

其中，n为扬声器总数。w∑中非零元素所对应的扬声器判别为激活扬声器，可在扬声器阵列中选择出对应的扬声器位置信息，并将激活扬声器总数记为na。由此，就完成了从n个扬声器中选取出na个激活扬声器的过程，用这些激活扬声器进行后期声场重放。

作为优选，步骤6要计算在步骤5中所选取出的na个激活扬声器的驱动信号，用以合成待重放声源在听音点处的声压值。设待重放声源包含的频域成分为：f1,f2,…,fh,…,fh，对应波数为：k1,k2,…,kh,…,kh，其中，kh＝2πfh/c，h＝1,2,…,h代表了待重放声源的实际频率对应的索引值，h为最高频率索引。待重放的声源频率成分要比虚拟声源(用于扬声器选取过程)丰富，即h>>8。

按照步骤2计算待重放声源辐射的目标声场声压

其中，kh＝2πfh/c，s为实际待重放单频的位置(与虚拟声源位置相同)；ym为第m个听音点的位置，m＝1,2,…,m。

按照步骤3计算激活扬声器重建声场的声压

式中，wh＝[wh(1),l,wh(na)]^t为激活扬声器的权值向量；gh为所有激活扬声器到所有听音点处的声传输函数集合：

其中，代表激活扬声器阵列中的第n个扬声器到第m个听音点的声学传输函数，n＝1,2,…,na，m＝1,2,…,m。

为了最小化重放误差，本设计利用l2正则化在均方最小准则约束下求取激活扬声器权值：

式中，是使最小时计算出的wh值。即，在频率fh下求得的激活扬声器权值向量。λ2＞0为正则化惩罚参数，其控制扬声器的总功率，需要根据实际选择，在满足最小化误差的同时用尽可能小的总功率实现重放。上式的求解过程为：

其中，^h代表共轭转置运算，wh即为所求激活的na个扬声器在频率fh下的驱动信号。

最后，将w1,w2,…,wh叠加，合成待重放声源在所有频率成分下的激活扬声器的驱动信号。

本发明的基于扬声器阵列的沉浸式宽带3d声场重放方法；研究目的是设计一套将场景a中的声源信号通过扬声器阵列在场景b中精确重放的虚拟声系统。本发明利用凸优化模型理论，首先，计算在指定空间位置放置的场景a的虚拟声源(与待重放声源位置相同，频率不同)到场景b中各个听音点的声学传输函数，将该函数值作为虚拟声源辐射声场的声压值；其次，设定场景b中某一墙面的扬声器阵列为规则的矩形等间距布局，并基于声波的波动特性，利用格林函数来建模所有扬声器到听音点的声学传播特征；再次，基于线性凸优化理论，将l1范数作为稀疏规则算子利用交替方向乘子法进行正则化运算，选取1倍频程内八个频带的中心频率计算扬声器权值，进行激活扬声器选取。最后，利用l2范数正则化，计算重放系统中激活扬声器的权值信号，使得在最小均方准则下，待重放声源的辐射声场和激活扬声器的辐射声场最接近。本发明思路清晰，算法复杂度低、实时性高，解决了实际三维声场重放所需扬声器数量较多的问题，最重要的是本方法能够通过平面阵列布局的扬声器进行重放，更易于实际摆放，既可应用于家庭影院和网络播放平台，也可用于办公环境和音视频会议。

附图说明

图1是扬声器阵列重放示例图。

图2是本设计方法整体示例图。

具体实施方法

设计的总体流程中，首先计算虚拟声源(与待重放声源位置相同的单频声源)到听音点的声压值，随后计算所有扬声器到听音点的声学传输函数，并与对应的扬声器权值相乘，得到扬声器辐射声场的声压表达式。通过l1范数正则化结合交替方向乘子法求取并汇总所有频点下的扬声器权值，保留其中非零项作为激活扬声器。最后，分别计算拟重放声源和激活扬声器到听音点的声压值，以均方最小为准则计算出各扬声器驱动信号，得到每个扬声器的播放信号，完成整体设计。实施时，本发明是算法镶嵌进软件中实现各个流程的自动运行，下面以具体实施步骤结合图1、图2对本发明做进一步说明，具体的工作流程如下：

步骤1：以扬声器所在墙面的中心位置为原点，建立空间直角坐标系。其中，x轴平行于地面，y轴垂直于地面，z轴垂直于墙面。指定待重放的虚拟声源的坐标位置为s。设定听音区域内听音点的个数和最多使用的扬声器个数分别为m和n，并计算其空间位置坐标；

步骤2：构建虚拟声源辐射声场。

利用格林函数来表征频率为fb的虚拟单频声源到第m个听音点的声学传输函数：

式中，kb＝2πfb/c为波数；c为空气中的声速，约为340m/s；b为频率索引；s为虚拟单频声源的位置(三维坐标)；ym为第m个听音点的位置，m＝1,2,…,m。本设计考虑宽带语音声场，选择1倍频程从100hz至16khz中8个频带的中心频率进行联合优化(每一频带的上限频率与下限频率的比值为2，起始频率为88hz)，即8个频带分别为：88hz-177hz、177hz-355hz、355hz-710hz、710hz-1420hz、1420hz-2840hz、2840hz-5680hz、5680hz-11360hz和11360hz-33720hz，对应的中心频率分别为125hz、250hz、500hz、1khz、2khz、4khz、8khz和16khz，故设定f1＝125hz,…,f8＝16khz，即，b＝1,2,…,8。

此外，本发明中的虚拟声源用于系统设计过程(即：从所有扬声器中选取激活扬声器的过程)，激活的扬声器用于声场重放，所重放的声源与虚拟声源位置相同。本发明中的虚拟声源是8个频率分别为f1,f2,…,f8的单频声源。考虑到虚拟声源具有全指向性，其辐射声场的声压值与格林函数值相等，可记为：

式中，^t代表转置运算；ym为第m个听音点的位置，m＝1,2,…,m；是频率为fb的虚拟单频声源在m个听音点处的辐射声场集合，是扬声器阵列待拟合的声场集合。

步骤3：构建扬声器辐射声场。

首先计算第n个扬声器到第m个听音点处的声学传输函数：

式中，xn为扬声器阵列中第n个扬声器的坐标，n＝1,2,…,n。则扬声器到所有听音点处的声传输函数集合为：

假定待求取的扬声器权值向量(也叫，扬声器激励信号向量或扬声器驱动信号向量)为wb＝[wb(1),wb(2),l,wb(n)]^t，n是扬声器总数。此时，扬声器辐射声场的声压值集合可表示为：

步骤4：针对虚拟声源求取扬声器权值向量(wb)。

本设计采用l1正则化方法对虚拟声源辐射声场和扬声器辐射声场进行约束，将扬声器权值的求取问题转换为最小绝对收缩和选择算子问题，如下式：

式中，||||1代表l1范数；||||2代表l2范数；λ1＝2为惩罚因子；是使最小时，计算出的wb值。即，在频率fb下求得的扬声器权值向量。

本设计采用交替方向乘子法对上式进行求解，首先，将上式写为交替方向乘子法的形式：

minimizef(wb)+g(zb)

s.t.wb-zb＝0

式中，g(zb)＝λ1||zb||1；zb是用于约束扬声器权值向量wb的变量。求解扬声器权值wb的过程是一个迭代的过程，第j+1次迭代得到的扬声器权值变量由如下公式求得：

式中，j+1为迭代次数，j以0为起点取值；为第j次迭代的对偶变量，用于判断算法的停止条件，为第j次迭代的约束变量，迭代初值初始化为i为单位矩阵；ρ为增广拉格朗日参数，其影响算法的收敛速度，在[0.1,10]内选取。约束变量迭代更新如下式所示：

式中，sκ(α)为软阈值判决函数，定义如下：

对偶变量的迭代更新过程如下式：

经j次迭代后，主残差为对偶残差为算法迭代的停止准则为：

其中，ε^pri是主残差的上限值，ε^dual是对偶残差的上限值，由下式计算：

本设计将迭代残差设定为：绝对误差ε^abs＝10^-4，相对误差ε^rel＝10^-2。当满足停止准则后，输出fb频点对应的一组扬声器权值向量wb。共需要求出8组不同频率下的扬声器权值向量wb，b＝1,2,…,8。

步骤5：汇总步骤4所求得的所有扬声器权值向量，统计保留其中的非零项，作为激活扬声器，同时确定其数量和位置信息。

在8个频点下所有的扬声器权值向量w可表达为：

将上述矩阵逐列相加，最终求得的汇总扬声器权值矩阵为w∑：

w∑＝[w∑(1),w∑(2),l,w∑(n)]^t

其中，n为扬声器总数。w∑中非零元素所对应的扬声器判别为激活扬声器，可在扬声器阵列中选择出对应的扬声器位置信息，并将激活扬声器总数记为na。由此，就完成了从n个扬声器中选取出na个激活扬声器的过程，用这些激活扬声器进行后期声场重放。

步骤6：针对待重放声源求取激活扬声器的权值信号。

设待重放声源包含的频域成分为：f1,f2,…,fh,…,fh，对应波数为：k1,k2,…,kh,…,kh，其中，kh＝2πfh/c，h＝1,2,…,h代表了待重放声源的实际频率对应的索引值，h为最高频率索引。待重放的声源频率成分要比虚拟声源(用于扬声器选取过程)丰富，即h>>8。

按照步骤2计算待重放声源辐射的目标声场声压

其中，kh＝2πfh/c，s为实际待重放单频的位置(与虚拟声源位置相同)；ym为第m个听音点的位置，m＝1,2,…,m。

按照步骤3计算激活扬声器重建声场的声压

式中，wh＝[wh(1),l,wh(na)]^t为激活扬声器的权值向量；gh为所有激活扬声器到所有听音点处的声传输函数集合：

其中，代表激活扬声器阵列中的第n个扬声器到第m个听音点的声学传输函数，n＝1,2,…,na，m＝1,2,…,m。

为了最小化重放误差，本设计利用l2正则化在均方最小准则约束下求取激活扬声器权值：

式中，w％h是使最小时计算出的wh值。即，在频率fh下求得的激活扬声器权值向量。λ2>0为正则化惩罚参数，其控制扬声器的总功率，需要根据实际选择，在满足最小化误差的同时用尽可能小的总功率实现重放。上式的求解过程为：

其中，^h代表共轭转置运算，wh即为所求激活的na个扬声器在频率fh下的驱动信号。

最后，将w1,w2,…,wh叠加，合成待重放声源在所有频率成分下的激活扬声器的驱动信号。

本文中所描述的具体实施例仅仅是对本发明精神作举例说明。本发明所属技术领域的技术人员可以对所描述的具体实施例做各种各样的修改或补充或采用类似的方式替代，但并不会偏离本发明的精神或者超越所附权利要求书所定义的范围。