本发明属声场主动控制
技术领域：
，具体涉及一种基于广义极小极大凹罚函数的声场重构方法。
背景技术：
：目前，声场重构方法主要可分为三类：(1)重构声场的球谐函数，称为高保真立体声重现；(2)波场合成，基于惠更斯原理和基尔霍夫-赫姆霍兹积分公式；(3)声压匹配(pressurematching，pm)，利用传声器阵列在空间中对期望声场采样以设计扬声器幅度，使重构声场与期望声场之间的误差平方最小。其中，pm方法可以看作以最小二乘(leastsquare，ls)为损失函数的线性回归问题。然而，传统的pm方法会激励所有的重构扬声器，这在高频或欠采样情况下，往往会导致重构区域内的声场发生空间混叠。尽管传声器位置处的重构误差很小，但是整个区域内的重构效果很差。为此，人们借鉴稀疏信号处理的方法，通过构造稀疏线性回归问题，采用最小绝对收缩和选择算子(least-absoluteshrinkageandselectionoperator，lasso)来重构声场。对比结果表明，基于lasso的声场重构方法比传统的pm方法效果更好，特别在高频或欠采样情况下，同时还减少了重构扬声器的个数。由于lasso方法本身的原因，求解得到的扬声器幅度偏小，使得重构声场达不到期望声场的声级，从而导致误差。相反，若在重构误差后加上非凸罚函数进行正则化，可以更精确地得到解的幅度。但是，此时目标函数往往也非凸，导致出现局部极值，难以求解。技术实现要素：为了克服现有技术的不足，本发明提供一种基于广义极小极大凹罚函数的声场重构方法。利用selesnick提出的广义极小极大凹(generalizedminimax-concave，gmc)罚函数，完成二维自由空间某区域内的声场重构，改善lasso对重构扬声器幅度的过低估计，更精确地在区域内重构期望声场。一种基于广义极小极大凹罚函数的声场重构方法，其特征在于步骤如下：步骤一：在自由空间中任意高度的水平面上选择一块区域，作为待重构的二维区域a，扬声器阵列均匀布置在区域a的边界上，传声器阵列位于区域a内，分别得到每个扬声器和传声器的位置矢量；所述的扬声器阵列包含l个扬声器，所述的传声器阵列包含m个传声器；扬声器的位置矢量表示为rl，l为扬声器序号，l＝1,2,…,l；传声器的位置矢量表示为rm，m为传声器序号，m＝1,2,…,m；步骤二：测量期望声场在传声器位置处的声压，得到期望声场声压向量d(ω)＝[d1(ω),d2(ω),…,dm(ω)]，其中，d表示期望声场声压向量，di(ω)表示第i个传声器位置处的期望声场声压，i＝1,2,…m；ω表示圆频率；步骤三：利用计算得到传声器到扬声器之间的声传递函数矩阵g，其中，c是声速，j为虚数单位，||·||2表示范数；步骤四：按照构造gmc目标函数jgmc，其中，s是重构扬声器幅度向量，||·||1表示范数，λ是稀疏度调节参数，λ＞0；是广义huber函数，其中，为中间变量，表示l维复数空间，是参数矩阵，0.5≤γ≤0.8；步骤五：利用邻近点算法求解最小化问题得到扬声器幅度向量s；步骤六：利用p＝gs计算得到传声器处的重构声场声压向量p。所述的利用邻近点算法求解最小化问题的具体过程如下：步骤1：参数初始化，令s＝0，v＝0，η＝max{1,γ/(1-γ)}||ghg||2，0＜μ＜2/η；步骤2：分别按照w＝s-μgh{g[s+γ(v-s)]-d}和u＝v-μγghg(v-s)计算中间参量w和u；步骤3：分别按照s＝soft(w,μλ)和v＝soft(u,μλ)更新s和v；其中，soft(·,·)表示软阈值函数；步骤4：若达到最大迭代次数k，则s即为最终的扬声器幅度向量；否则，重复步骤2-3；k的取值范围为5000～10000；所述的稀疏度调节参数λ还可以采用交叉验证方法确定，具体过程如下：步骤1：将向量d任意划分为ng个非重叠的子集，即一般ng取值为2、3或4；步骤2：在区间[0,||ghd||∞]上均匀选择n个值作为λ的可能值，即其中||·||∞表示取向量元素的最大绝对值，n的取值范围为1000～3000；步骤3：按照计算得到第g个子集的误差其中，s(-g)为λ＝λn、d＝d(-g)、g＝g(-g)时利用邻近点算法求解最小化问题得到的解；d(g)表示d的第g个子集，g(g)为g中与第g个子集d(g)对应的列组成的矩阵，d(-g)表示d删除d(g)后得到的向量，g(-g)表示g删除g(g)后得到的矩阵；步骤4：按照计算得到λn对应的误差en；步骤5：按照λ＝λn′确定最终参数λ，其中本发明的有益效果是：由于采用gmc罚函数对重构误差进行正则化，可以改善一般稀疏正则化对扬声器幅度的过低估计，更精确地在区域内重构期望声场。因此，通过激励少量的扬声器，能够在传声器处获得比lasso方法更好的重构精度。同时，由于只激励少量的扬声器，能够避免产生空间混叠现象，特别是对于高频期望声场和当期望声场欠采样时。附图说明图1是本发明的一种基于广义极小极大凹罚函数的声场重构方法流程图图2是期望声场示意图图3是利用传统pm方法重构得到的声场示意图图4是利用lasso方法重构得到的声场示意图图5是利用本发明方法重构得到的声场示意图具体实施方式下面结合附图和实施例对本发明进一步说明，本发明包括但不仅限于下述实施例。本发明的一种基于广义极小极大凹罚函数的声场重构方法，其基本流程如图1所示，具体过程如下：1、在自由空间中，取z＝0平面上的一块圆形区域r＝1.5m作为待重构区域a，其圆心作为坐标系原点，64个扬声器均匀布置在区域a的边界上，即组成半径r＝1.5m的圆形阵列；12个传声器位于区域a内，组成半径r＝0.25m的圆形阵列。如图2所示，其中，方框表示扬声器，圆点表示传声器。从而可以分别确定得到每个扬声器的位置矢量rl和每个传声器的位置矢量rm，其中，l为扬声器序号，l＝1,2,…,64，m为传声器序号，m＝1,2,…,12。2、假设在频率f＝800hz下，由两个分别位于(4,0,0)、(0,4,0)的点声源发声产生期望声场，并测量期望声场在传声器阵列位置处的声压，得到期望声场声压向量d(ω)＝[d1(ω),d2(ω),…,d12(ω)]，其中，表示期望声场声压向量，di(ω)表示第i个传声器位置处的期望声场声压，i＝1,2,…,12；ω＝2πf为圆频率。3、利用下式计算传声器到扬声器之间的声传递函数矩阵其中，g(m,l)表示矩阵g中的第m行l列的值，c是声速，j为虚数单位，||·||2表示范数。4、构造gmc目标函数，如下式所示：其中，是重构扬声器幅度向量，||·||1表示范数。λ＞0是稀疏度调节参数，它直接控制解的稀疏度，λ值越大，解越稀疏，但当λ大于||ghd||∞时，计算结果s＝0，没有意义，因此，可以在区间[0,||ghd||∞]上确定一个合适的λ值。本发明采用交叉验证法确定参数λ，即将向量d分为ng＝2个非重叠的子集，在区间[0,||ghd||∞]上均匀选择n＝101个可能值按照计算得到第g个子集的误差其中，s(-g)为λ＝λn、d＝d(-g)、g＝g(-g)时利用邻近点算法求解最小化问题得到的解；d(g)表示d的第g个子集，g(g)为g中与第g个子集d(g)对应的列组成的矩阵，d(-g)表示d删除d(g)后得到的向量，g(-g)表示g删除g(g)后得到的矩阵；然后，按照计算得到λn对应的误差en；最后，按照λ＝λn′确定最终参数λ，其中本实施中利用交叉验证确定参数λ为3.5974×10-4。hb(s)是广义huber函数，定义为：式中，是参数矩阵。selesnick证明，若：则jgmc(s)为凸函数。因此，可以简单地令：其中，参数γ控制gmc罚函数的非凸性：若γ＝0，则b＝o，gmc目标函数退化为lasso目标函数若γ＝1，则gmc罚函数最大程度地非凸。实际中，一般取0.5≤γ≤0.8，本实施例中γ＝0.8。5、利用邻近点算法求解最小化问题计算得到扬声器幅度向量s。具体过程为：(1)参数初始化，令s＝0，v＝0，η＝max{1,γ/(1-γ)}||ghg||2，0＜μ＜2/η；(2)分别按照w＝s-μgh{g[s+γ(v-s)]-d}和u＝v-μγghg(v-s)计算中间参量w和u；(3)分别按照s＝soft(w,μλ)和v＝soft(u,μλ)更新s和v；其中，soft(·,·)表示软阈值函数；(4)若达到最大迭代次数k，则s即为最终的扬声器幅度向量；否则，重复步骤(2)-(3)；本实施例中k取值为10000。6、利用p＝gs计算得到传声器处的重构声场声压向量p。为了对比本发明方法效果，分别利用传统pm方法和lasso方法进行声场重构，其中，图3、图4分别为其重构后的声场示意图，图5为利用本发明方法重构后声场的示意图。图中，方框表示激励的重构扬声器，圆点表示传声器，圆圈表示扬声器的幅度。可以看出，传统pm方法激励了所有的重构扬声器，造成了空间混叠。表1给出了三种方法重构误差和激励扬声器个数的对比数据。可以看出，本发明方法的重构误差更小，激励的扬声器个数更少。表1本发明方法传统pm方法lasso方法重构误差e0.00180.00310.0046激励扬声器个数146415当前第1页12