CN103453980B - 一种基于压缩感知的声场参数获取方法 - Google Patents

Abstract

本发明涉及一种基于压缩感知的声场参数获取方法，属于数字信号处理技术领域。本发明包括球型传声器阵列模块，常数观测矩阵生成模块，观测信号向量生成模块，正交基构造模块，随机观测矩阵生成模块，球谐波基系数重构模块，目标区域声压分布重构模块。在所述的确定球型传声器阵列设计模块中，考虑到可实现性及阵列小型化，人为确定球型传声器阵列半径，在确定下的球型半径下，球谐波基系数具有稀疏性，从而依据压缩感知理论，分别在正交基构造模块和随机观测矩阵生成模块中，构造出正交基和随机观测矩阵，同时输入球谐波基系数重构模块中，重构出球谐波基系数，最后将这些系数输入目标区域声压分布重构模块中重构出目标区域声压分布。

Description

一种基于压缩感知的声场参数获取方法

技术领域

本发明涉及一种基于压缩感知的声场参数获取方法，属于数字信号处理技术领域。

背景技术

三维音频重现是产生能给听众带来真实、自然感受的声场，在通信、数字娱乐、人机交互等众多领域有着广泛的应用。此处说一下声场重现的含义与应用领域。Ambisonics方法是在物理意义上能重现三维音频声场的技术之一，其最初由Gerzon提出并用于二维声场复制。Ambisonics方法的基本思路是将声场表示成圆谐波的线性组合。近年来，随着技术的发展，多采用球谐波分解来表示声场，用于高阶声场的重现，如三维声场重现^[2]，即高阶Ambisonics方法。该方法将声场分布在无穷阶球谐波域中进行分解展开，这样只要采集到各球谐波域中的分量，理论上就可重现出声场。该方法不受声场重现系统的限制，可任意放置扬声器位置，但在声场数据采集时，其性能与麦克风数目以及麦克风指向性密切相关。

目前已有的采集方案中，采集N阶声场至少需要(N+1)²个传声器，随着声场阶数N增大，传声器数目也增大，大量传声器置于声场中易造成声场畸变。本发明提出使用压缩感知方法来解决该问题。

设在自由空间中，有一单极子声源位于S(1米,0米,0米)处，单极子声源作球的伸缩运动产生球面波，单极子声源引起的声压为

$P=\mathrm{Amp}\frac{\exp (-k|r-r_S\left|\right)}{|r-r_S|},---\left(1\right)$

式中，Amp为声源处声压幅度值，一般视为单位幅度值1；|r-r_s|为空间某点到声源S的距离；k＝Ω/c＝2πf/c为波数，c为声波在媒介中传播速度，Ω为角频率(单位弧度)，f为频率(单位Hz)。采集声场分布的区域为一半径为r，球心为空间坐标系原点的封闭球面A如说明书附图图1所示。

可用Fourier-Bessel级数将N阶声场声压分解展开为

$p(r,\mathrm{\Ω})=\underset{n=0}{\overset{N}{\mathrm{\Σ}}}\underset{m=-n}{\overset{n}{\mathrm{\Σ}}}{A}_{\mathrm{nm}}\left(f\right)j_n\left(\frac{\mathrm{\Ω}}{c}r\right)Y_{\mathrm{nm}}\left(\hat{r}\right)=\underset{n=0}{\overset{N}{\mathrm{\Σ}}}\underset{m=-n}{\overset{n}{\mathrm{\Σ}}}{B}_{\mathrm{nm}}(f,r)Y_{\mathrm{nm}}\left(\hat{r}\right),---\left(2\right)$

式中，只与方位角仰角θ有关，B_nm(f,r)为一组球谐波系数，与空间位置无关，j_n(·)为第一类n阶球Bessel函数，Y_nm(·)为球谐波基

式中，P_n|m|(·)为连带Legendre函数，n为球谐波展开阶数，它决定重现时需要的声道数，m＝-n,…,n为每一阶对应的球谐波基，球谐波基具有正交性

${\∫}_A{Y}_{\mathrm{bd}}\left(\hat{r}\right)Y_{\mathrm{nm}}^{*}\left(\hat{r}\right)\mathrm{dA}={\mathrm{\δ}}_{\mathrm{nb}}{\mathrm{\δ}}_{\mathrm{md}},---\left(4\right)$

式中δ_nb为δ函数。

通过式(2)将N阶声场分解到(N+1)²个球谐波基上，对应有(N+1)²个球谐波基系数B_nm。此时要获取球面A的完整声场声压，需先获取声场声压分布在各球谐波基上的系数，该系数可根据球谐波傅里叶变换得到

实际采集时在封闭球面A上只能布局离散的L个采集点来获取B_nm，采集N阶声场声压的传声器布局如说明书附图2所示。同时，每个采集位点上需要引入加权因子α_i才能保证式(4)在离散采样时仍然成立：

$\underset{i}{\overset{L}{\mathrm{\Σ}}}{\mathrm{\α}}_i{Y}_{\mathrm{bd}}\left({\hat{r}}_i\right)Y_{\mathrm{nm}}^{*}\left({\hat{r}}_i\right)={\mathrm{\δ}}_{\mathrm{bn}}{\mathrm{\δ}}_{\mathrm{dm}},,---\left(6\right)$

在系数B_nm已知下，可根据式(2)重构球面A声压分布，使用扬声器球型阵列重放时，先从球面A的声压分布中获取声源空间位置信息，然后使用这些空间信息控制扬声器重现三维音频声场。

R.Duraiswami,N.A.Gumerov,A.O'Donovan在专利《Audiocamerausingmicrophonearraysforrealtimecaptureofaudioimagesandmethodforjointlyprocessingtheaudioimageswithvideoimages》(US,381/92,US20120288114A1.2012)中直接根据球谐波变换式(5)以及球谐波基的正交性式(6)，在封闭球面A上等角布局离散的L个传声器，离散采集N阶声场，各球谐波基系数B_nm能根据加权求和法获取

$B_{\mathrm{nm}}=\underset{i}{\overset{L}{\mathrm{\Σ}}}{\mathrm{\α}}_{i}p(r_i,\mathrm{\Ω})Y_{\mathrm{nm}}^{*}\left({\hat{r}}_i\right),---\left(7\right)$

该技术方案能使重构声场达到较高精度，若用一个旋转的麦克风进行空间采样时，旋转步长与角度差异一致时，采集信息因旋转造成的误差会很小，但无法实现实时采集。若在球面A上等角布局传声器，可以实现实时采集，但在球面上布局密集的传声器会影响声场分布，甚至明显改变声场的分布。换句话说，若测量手段较大地影响测量对象，则无论测量手段多么精确，测量结果也不可能精确，这也使测量过程本身失去了意义。

J.Ahrens,S.Spors在文献《Amodalanalysisofspatialdiscretizationofsphericalloudspeakerdistributionsusedforsoundfieldsynthesis》(IEEETransactiononAudio,Speech,andLanguageProcessing.2012,20(9):2564-2574)中应用最小二乘法技术，提出了一种减少传声器数目的方法。在该方法法中，若L＝(N+1)²，那么(N+1)²个谐波系数可对应有(N+1)²个式(7)，便构成线性方程组，使用最小二乘法求解，此方法精度较高，但是随着阶数地增加，需要大量传声器，在实际采样点布局时，对于N阶声场一般使用1.5(N+1)²个传声器，通过增加采集点提高采样结果精度。这些采样点根据Euclid范数准则确定，约束条件为

$E({\hat{r}}_1,...,{\hat{r}}_K)=\underset{i=1}{\overset{K}{\mathrm{\Σ}}}\underset{j=i+1}{\overset{K}{\mathrm{\Σ}}}\frac{1}{\left|\right|{\hat{r}}_i-{\hat{r}}_j\left|\right|},$

使上式最小值所对应的采样点为最佳采样点。加权值可用J.FliegeandU.Maier在文献《Atwo-stageapproachforcomputingcubatureformulaeforthesphere》(Dortmund:Dortmund,1996)中使用的快速算法计算得到。此技术在Euclid范数准则下使用尽可能少的传声器采集声场信息，但是随着声场阶数N增大，需要的传声器数目越大，虽然采集到的信号精度能达到较高程度，但如此多的传声器仍使原始声场发生畸变。

发明内容

为了克服上述已有技术的不足，本发明的目的在于提供一种基于压缩感知的声场参数获取方法。

本发明采取的技术方案如下：

本发明包括包括球型传声器阵列模块，常数观测矩阵生成模块，观测信号向量生成模块，正交基构造模块，随机观测矩阵生成模块，球谐波基系数重构模块，目标区域声压分布重构模块。

在所述的确定球型传声器阵列设计模块中，考虑到物理可实现性以及阵列小型化便于操作，人为确定球型传声器阵列半径，观测信号向量生成模块根据球型传声器阵列模块输出的声压向量和常数观测矩阵生成模块生成的常数观测矩阵计算出观测信号向量，在确定下的球型半径下，球谐波基系数具有稀疏性，从而依据压缩感知理论，分别在正交基构造模块和随机观测矩阵生成模块中，构造出正交基和随机观测矩阵，同时输入球谐波基系数重构模块中，可重构出球谐波基系数，最后将这些系数输入目标区域声压分布重构模块中可以重构出目标区域声压分布。

本发明原理及有益效果：在三维音频中，高阶Ambisonics技术需要大量传声器来感知声场分布。声场分布的感知是重建声场的前提，在感知声场的过程中，必须引入传声器阵列。通常阵列中传声器分布越密集，感知到的声场分布越精确，但是过多的传声器会影响或改变声场分布。为将引入传声器阵列对声场分布的影响降到最低，希望传声器数目尽可能少。为此，本发明提出了应用压缩感知技术来减少传声器数目。

附图说明

图1本发明涉及的Ambisonics三维空间示例。

图2球型传声器在球形空间A上布局示意图。

图3基于压缩感知的声场参数获取方法功能框图。

图4A_nm(f)误差最小时贝塞耳函数的参数arg(N)的上限。

图5球形空间A上M个采集点分布图。

图6在f＝3kHz、r＝4cm、N＝6时，选择不同M×N的高斯随机观测矩阵结果相对误差。

图7在f＝3kHz、r＝4cm、N＝6、M＝32时，球面上理想声压分布图。

图8在f＝3kHz、r＝4cm、N＝6、M＝32时，球面上重构声压分布图。

图9在f＝630Hz、r＝4cm、N＝6、M＝32时，球面上重构声压分布图。

图10在f＝2kHz、r＝4cm、N＝6、M＝32时，球面上重构声压分布图。

图11在f＝4kHz、r＝4cm、N＝6、M＝32时，球面上重构声压分布图。

图12在f＝5kHz、r＝4cm、N＝6、M＝32时，球面上重构声压分布图。

图13球面上重构声压相对误差曲线。

具体实施方式

下述非限制性实施例可以使本领域的普通技术人员更全面地理解本发明，但不以任何方式限制本发明。

灰度图能够更好的说明本发明的技术效果，特提供灰度图来说明本发明的技术效果，同时为了更好的让审查员审查本发明的技术效果，提供灰度图以供参考。灰度图为：图1(图1带颜色)，图2，图5，图7至图12。

实施例1

在声场重建方案中，过多传声器会干扰声场分布，导致采集信号不准确，从而影响重建效果，压缩感知能将高维信号降维进行采集，这样会减少采集设备数量。本发明用压缩感知方法，减少传声器数目，以减少传声器对感知声场的干扰，在保证精度的情况下重构出三维音频声场。本发明首先依据A_nm(f)计算误差最小原则设计球型传声器阵列半径r；然后根据球面积分误差最小原则，在半径为r的球面上布局传声器，获取声压信号p，得到观测信号x；同时使用球谐波基构造出正交基Ψ，并采用高斯随机矩阵设计与之无关的观测矩阵Φ；其次，使用原对偶内点法重构球谐波基系数γ；最后利用重构的球谐波基系数γ来重构出观测面的声压分布。本发明的系统框图如图3所示：本发明的系统包括球型传声器阵列模块，常数观测矩阵生成模块，观测信号向量生成模块，正交基构造模块，随机观测矩阵生成模块，球谐波基系数重构模块，目标区域声压分布重构模块；将球型传声器阵列模块安放在球形采集面处，采集球形采集面上的声压信号，该信号与常数观测矩阵生成模块生成的常数观测矩阵在观测信号向量生成模块中相乘，得到观测信号向量，正交基构造模块输出的正交矩阵输入到随机观测矩阵生成模块，用于选择适当的随机观测矩阵，随机观测矩阵生成模块生成的的随机观测矩阵、观测信号向量生成模块生成的观测信号向量在球谐波基系数重构模块中构造出球谐波基系数，根据球谐波基系数重构模块构造的球谐波基系数在目标区域声压分布重构模块中计算出球形采集面处任意位置处的声压，达到重构出目标区域声压分布的目的。

下面对本发明涉及到的模块做说明：

球型传声器阵列模块

球型传声器阵列模块处理方法如下：

在声场采集区域A中，把M个传声器等立体角地均匀布置在给定半径r的球形采集面就构成了所述的球型传声器阵列模块；其中，M、r由采集声场的精度确定；传声器的位置由以球形采集面的球心为原点，极坐标系下不同的方位角仰角θ描述。

常数观测矩阵生成模块

常数观测矩阵生成模块设计如下：

根据式(2)，球谐波系数B_nm中含有第一类n阶球贝塞耳函数j_n(kr)，这里j_n(kr)与频率和球半径有关。当声场截断为N阶声场时，根据采集到的N阶声场声压信息，计算B_nm中只与频率有关的部分A_nm(f)。为使计算得到的A_nm(f)误差最小，应该保证kr≤arg(N)，其中arg(N)为一常数，选择K(K>N)阶j_K(kr)比N阶j_N(kr)衰减达到10dB以上时对应参数kr＝arg(N)，具体值可参见图4。

从图4中可见，当阶数N一定时，期望能准确获取的频率f越大，球半径r就需要越小。考虑到物理可实现以及阵列小型化便宜操作，本发明采集区域A半径r设置为r＝4厘米，此时当N越大，准确获取的频率f越大。

根据球形传声器阵列中每个传声器的方位角仰角θ；计算出的M×M维常数观测矩阵C为，

$C=\left[\begin{array}{ccc}{Y}_{00}\left({\hat{r}}_1\right)& \mathrm{\Λ}& Y_{\mathrm{NN}}\left({\hat{r}}_1\right)\\ M& O& M\\ Y_{00}\left({\hat{r}}_M\right)& \mathrm{\Λ}& Y_{\mathrm{NN}}\left({\hat{r}}_M\right)\end{array}\right].$

观测信号向量生成模块

观测信号向量生成模块处理方法如下，

在r＝4厘米的球面上等立体角地均匀布置M个传声器，M<(N+1)²，用于采集这M个位置处由不同θ所区分的声压p_m，m＝1,2,…,M，构成声压信号向量p＝[p₁,p₂,…,p_M]，然后转换为观测信号向量x，

x＝Cp，(8)

其中，C为常数观测矩阵生成模块设计出的M×M维的常数观测矩阵C。

正交基构造模块

正交基构造模块处理方法如下：

要用压缩感知理论，首先要构造出一个正交基，即在球面A上选择L＝(N+1)²点处采集声压，一次构成L维的声压信号p，其每一维声压即每一个点处的用式(2)表示；以矩阵形式表示p，则有

$p=\left[\begin{array}{c}\underset{n=0}{\overset{N}{\mathrm{\&Sigma;}}}\underset{m=-n}{\overset{n}{\mathrm{\&Sigma;}}}{B}_{\mathrm{nm}}{Y}_{\mathrm{nm}}\left({\hat{r}}_1\right)\\ M\\ \underset{n=0}{\overset{N}{\mathrm{\&Sigma;}}}\underset{m=-n}{\overset{n}{\mathrm{\&Sigma;}}}{B}_{\mathrm{nm}}{Y}_{\mathrm{nm}}\left({\hat{r}}_L\right)\end{array}\right]=\left[\begin{array}{ccc}{Y}_{00}\left({\hat{r}}_1\right)& \mathrm{\&Lambda;}& Y_{\mathrm{NN}}\left({\hat{r}}_1\right)\\ M& O& M\\ Y_{00}\left({\hat{r}}_L\right)& \mathrm{\&Lambda;}& Y_{\mathrm{NN}}\left({\hat{r}}_L\right)\end{array}\right]\left[\begin{array}{c}{B}_{00}\\ M\\ B_{\mathrm{NN}}\end{array}\right]=\mathrm{\&Psi;\&gamma;},---\left(9\right)$

其中， $\mathrm{\&Psi;}=\left[\begin{array}{ccc}{Y}_{00}\left({\hat{r}}_1\right)& \mathrm{\&Lambda;}& Y_{\mathrm{NN}}\left({\hat{r}}_1\right)\\ M& O& M\\ Y_{00}\left({\hat{r}}_L\right)& \mathrm{\&Lambda;}& Y_{\mathrm{NN}}\left({\hat{r}}_L\right)\end{array}\right]$ 为N阶球谐波基构成的L×L维矩阵， $\mathrm{\&gamma;}=\left[\begin{array}{c}{B}_{00}\\ M\\ B_{\mathrm{NN}}\end{array}\right]$ 为各球谐波基系数构成的L维矢量；在引入加权因子的条件下，有

$\underset{i}{\overset{L}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_i{Y}_{\mathrm{bd}}\left({\hat{r}}_i\right)Y_{\mathrm{nm}}^{*}\left({\hat{r}}_i\right)\&ap;{\mathrm{\&delta;}}_{\mathrm{bn}}{\mathrm{\&delta;}}_{\mathrm{dm}}$

写成矩阵形式

$\left[\begin{array}{ccc}{\mathrm{\&alpha;}}_1{Y}_{00}\left({\hat{r}}_1\right)& \mathrm{\&Lambda;}& {\mathrm{\&alpha;}}_L{Y}_{00}\left({\hat{r}}_L\right)\\ M& O& M\\ {\mathrm{\&alpha;}}_1{Y}_{\mathrm{NN}}\left({\hat{r}}_1\right)& \mathrm{\&Lambda;}& {\mathrm{\&alpha;}}_L{Y}_{\mathrm{NN}}\left({\hat{r}}_L\right)\end{array}\right]\left[\begin{array}{ccc}{Y}_{00}\left({\hat{r}}_1\right)& \mathrm{\&Lambda;}& Y_{\mathrm{NN}}\left({\hat{r}}_1\right)\\ M& O& M\\ Y_{00}\left({\hat{r}}_L\right)& \mathrm{\&Lambda;}& Y_{\mathrm{NN}}\left({\hat{r}}_L\right)\end{array}\right]{\widetilde{\mathrm{\&Psi;}}}^T\mathrm{\&Psi;}\&ap;E,---\left(10\right)$

其中，Ψ为近似的正交基组成的矩阵，E为单位对角阵。

随机观测矩阵生成模块

随机观测矩阵生成模块处理方法如下：

定义声压信号向量p在球谐波正交基Ψ上的投影为γ，使用随机数发生器，生成M×L维高斯随机矩阵Φ，高斯随机矩阵Φ均值为0，方差为Φ的选择应使其与Ψ之间具有较高不相关性，Φ的任意列向量与Ψ的任意列向量之间相关系数越小越好，小于0.1，用Φ将γ降维为M维观测信号x

x＝Φγ。(11)

球谐波基系数重构模块

球谐波基系数重构模块的处理方法如下：

根据压缩感知理论，在l₁-范数约束条件下，x能重构γ，即式(12)表示的优化问题有解，

Minimize[||γ||_1范数]s.t.Φγ＝x(12)

将式(12)表示的优化问题转化为线性规划问题：

Minimize[sum(u)]s.t.-u≤γ≤u，Φγ＝x(13)

其中，s.t.表示是subjectto的缩写；u＝0.95·abs(γ)+0.1·max[abs(γ)]；Φ为M×L维高斯随机矩阵，均值为0，方差为Minimize[h(·)]表示使函数h(·)最小化，max()表示取向量元素中的最大值，abs()表示向量元素求绝对值，sum()表示向量元素求和，式(13)能用原对偶内点算法求解，具体步骤如下：

(1)根据线性规划，设定线性函数为 $g\left(\mathrm{\&gamma;}\right)=\left[\begin{array}{c}\mathrm{\&gamma;}-u\\ -\mathrm{\&gamma;}-u\end{array}\right],$ 这里γ为原始变量向量，初始值为Φ^Tx，同时设置迭代误差阈值为ε＝10^-3，迭代最大次数为50；设定第k次迭代时的线性函数 $g\left({\mathrm{\&gamma;}}^{\left(k\right)}\right)=\left[\begin{array}{c}{\mathrm{\&gamma;}}^{\left(k\right)}-u^{\left(k\right)}\\ -{\mathrm{\&gamma;}}^{\left(k\right)}-u^{\left(k\right)}\end{array}\right],$ k表示第k次迭代；

(2)计算补偿间隙其中 ${\mathrm{\&lambda;}}^{\left(k\right)}=\left[\begin{array}{c}1/({\mathrm{\&gamma;}}^{\left(k\right)}-u^{\left(k\right)})\\ 1/(-{\mathrm{\&gamma;}}^{\left(k\right)}-u^{\left(k\right)})\end{array}\right]=\left[\begin{array}{c}{{\mathrm{\&lambda;}}_1}^{\left(k\right)}\\ {{\mathrm{\&lambda;}}_2}^{\left(k\right)}\end{array}\right]$ 为一对偶变量向量，同时计算出中心参数其中，μ＝10为障碍常数，m＝2L；

(3)依据Karush-Kuhn-Tucker条件，计算第k次迭代的原始残差向量r_pri ^(k)、对偶残差向量r_dual ^(k)以及中心残差向量r_cent ^(k)

r_pri ^(k)＝Φγ^(k)-x^(k)(14)

${r_{\mathrm{dual}}}^{\left(k\right)}={\&dtri;g}_0\left({\mathrm{\&gamma;}}^{\left(k\right)}\right)+\mathrm{Dg}\left({\mathrm{\&gamma;}}^{\left(k\right)}\right){\mathrm{\&lambda;}}^{\left(k\right)}+{\mathrm{\&Phi;}}^T{v}^{\left(k\right)}---\left(15\right)$

r_cent ^(k)＝-diag(λ^(k))g(γ^(k))-(1/t^(k))E(16)

其中，v^(k)＝Φ^T[λ₂ ^(k)-λ₁ ^(k)]，为一对偶变数向量；g₀(γ^(k))＝sum(u^(k))，g₁(γ^(k))＝γ^(k)-u^(k)，g₂(γ^(k))＝-γ^(k)-u^(k)；为对γ^(k)与u^(k)求偏导，即 $\&dtri;g_0\left({\mathrm{\&gamma;}}^{\left(k\right)}\right)=\left[\begin{array}{c}0\\ 1\end{array}\right];$ $\mathrm{Dg}\left({\mathrm{\&gamma;}}^{\left(k\right)}\right)=\left[\begin{array}{c}\&dtri;g_1{\left({\mathrm{\&gamma;}}^{\left(k\right)}\right)}^T\\ \&dtri;g_2{\left({\mathrm{\&gamma;}}^{\left(k\right)}\right)}^T\end{array}\right]=\left[\begin{array}{cc}1& -1\\ -1& -1\end{array}\right];\mathrm{diag}\left({\mathrm{\&lambda;}}^{\left(k\right)}\right)=\left[\begin{array}{cc}\mathrm{diag}\left({{\mathrm{\&lambda;}}_1}^{\left(k\right)}\right)& 0\\ 0& \mathrm{diag}\left({{\mathrm{\&lambda;}}_2}^{\left(k\right)}\right)\end{array}\right];$ E为单位阵；所以r_dual ^(k)与r_cent ^(k)具体计算式为，

${r_{\mathrm{dual}}}^{\left(k\right)}=\left[\begin{array}{c}0\\ 1\end{array}\right]+\left[\begin{array}{cc}{{\mathrm{\&lambda;}}_1}^{\left(k\right)}& -{{\mathrm{\&lambda;}}_2}^{\left(k\right)}\\ -{{\mathrm{\&lambda;}}_1}^{\left(k\right)}& -{{\mathrm{\&lambda;}}_2}^{\left(k\right)}\end{array}\right]+\left[\begin{array}{c}{\mathrm{\&Phi;}}^T[{{\mathrm{\&lambda;}}_1}^{\left(k\right)}-{{\mathrm{\&lambda;}}_2}^{\left(k\right)}]\\ 1\end{array}\right]---\left(17\right)$

${r_{\mathrm{cent}}}^{\left(k\right)}=-\left[\begin{array}{cc}{{\mathrm{\&lambda;}}_1}^{\left(k\right)}{g}_1\left({\mathrm{\&gamma;}}^{\left(k\right)}\right)& 0\\ 0& {{\mathrm{\&lambda;}}_2}^{\left(k\right)}{g}_2\left({\mathrm{\&gamma;}}^{\left(k\right)}\right)\end{array}\right]-(1/t^{\left(k\right)})E---\left(18\right)$

其中，‘0’为元素皆为0的向量，‘1’为元素皆为1的向量，‘-1’为元素皆为-1的向量；

原对偶搜索方向需要使原始残差向量r_pri ^(k)、对偶残差向量r_dual ^(k)以及中心残差向量r_cent ^(k)接近0，根据式(14)～(16)所示的泰勒展开，得到

$\left[\begin{array}{ccc}{\&dtri;}^2{g}_0\left({\mathrm{\&gamma;}}^{\left(k\right)}\right)+{\mathrm{\&Sigma;}}_{i=1}^2{{\mathrm{\&lambda;}}_i}^{\left(k\right)}{\&dtri;}^2{g_i}^{\left(k\right)}\left({\mathrm{\&gamma;}}^{\left(k\right)}\right)& \mathrm{Dg}{\left({\mathrm{\&gamma;}}^{\left(k\right)}\right)}^T& {\mathrm{\&Phi;}}^T\\ -\mathrm{diag}\left({\mathrm{\&lambda;}}^{\left(k\right)}\right)\mathrm{Dg}\left({\mathrm{\&gamma;}}^{\left(k\right)}\right)& -\mathrm{diag}\left(g\left({\mathrm{\&gamma;}}^{\left(k\right)}\right)\right)& 0\\ \mathrm{\&Phi;}& 0& 0\end{array}\right]\left[\begin{array}{c}\mathrm{\&Delta;}{\mathrm{\&gamma;}}^{\left(k\right)}\\ \mathrm{\&Delta;}{\mathrm{\&lambda;}}^{\left(k\right)}\\ \mathrm{\&Delta;}{v}^{\left(k\right)}\end{array}\right]=-\left[\begin{array}{c}{r_{\mathrm{dual}}}^{\left(k\right)}\\ {r_{\mathrm{cent}}}^{\left(k\right)}\\ {r_{\mathrm{pri}}}^{\left(k\right)}\end{array}\right]---\left(19\right)$

即，

$\left[\begin{array}{ccc}0& \left[\begin{array}{cc}1& -1\\ -1& -1\end{array}\right]& {\mathrm{\&Phi;}}^T\\ \left[\begin{array}{cc}-{{\mathrm{\&lambda;}}_1}^{\left(k\right)}& {{\mathrm{\&lambda;}}_1}^{\left(k\right)}\\ {{\mathrm{\&lambda;}}_2}^{\left(k\right)}& {{\mathrm{\&lambda;}}_2}^{\left(k\right)}\end{array}\right]& \left[\begin{array}{cc}-g_1\left({\mathrm{\&gamma;}}^{\left(k\right)}\right)& 0\\ 0& -g_2\left({\mathrm{\&gamma;}}^{\left(k\right)}\right)\end{array}\right]& 0\\ \mathrm{\&Phi;}& 0& 0\end{array}\right]\left[\begin{array}{c}\mathrm{\&Delta;}{\mathrm{\&gamma;}}^{\left(k\right)}\\ \mathrm{\&Delta;}{\mathrm{\&lambda;}}^{\left(k\right)}\\ \mathrm{\&Delta;}{v}^{\left(k\right)}\end{array}\right]=-\left[\begin{array}{c}{r_{\mathrm{dual}}}^{\left(k\right)}\\ {r_{\mathrm{cent}}}^{\left(k\right)}\\ {r_{\mathrm{pir}}}^{\left(k\right)}\end{array}\right]---\left(20\right)$

可计算出搜索方向(Δγ^(k),Δλ^(k),Δv^(k))，同时，设定线性搜索的步长s^(k)＝0.99min{1,min[-λ^(k)/Δλ^(k)|Δλ_i ^(k)<0]}；具体步骤如下，

(A)开始线性搜索，将γ^(k)+s^(k)Δγ^(k)、λ^(k)+s^(k)Δλ^(k)、v^(k)+s^(k)Δv^(k)替代γ^(k)、λ^(k)、v^(k)后计算新的g(γ^(k))、r_pri ^(k)、r_dual ^(k)、r_cent ^(k)，为了有所区别，分别用g'(λ^(k))、r'_pri ^(k)、r'_dual ^(k)、r'_cent ^(k)表示，若α的取值范围是0<α<1，本发明建议α＝0.01，则结束线性搜索，进入第(B)步，否则继续第(A)步，同时令βs^(k)取代s^(k)，β＝0.5；

(B)在线性搜索后，第k+1次迭代的初始值γ^(k+1)、λ^(k+1)、v^(k+1)分别为γ^(k)+s^(k)Δγ^(k)、λ^(k)+s^(k)Δλ^(k)、v^(k)+s^(k)Δv^(k)，重新计算补偿间隙若或迭代次数超过50次，则迭代截止，否则从第(3)步继续迭代。

目标区域声压分布重构模块

目标区域声压分布重构模块处理方法如下：

根据重构出的γ以及式(2)，可获得整个球面A上的声压分布，这样，可以使用M个传声器获取到的声场信息来重构N阶的声场信息，M<(N+1)²，从而可减少传感器的数目，更进一步，利用重构的球面A上的声压分布信息以及记录的音频信号，在重放空间可高度近似重现真实的声场。

本发明技术方案带来的有益效果：

以N＝6为例，根据图4，此时能较准确地采集到频率上限约为5967Hz，采用压缩感知方法进行三维音频声场采集。当N＝6阶时，将有49个不同的球谐波基，在声压采集时，在球面上只需布局M(M<49)个传声器，具体如图5所示，这样可以重构出49个球谐波基系数，从而获得整个球面上任意一点的声压分布。

首先观察f＝3kHz情况，通过实验，当观测矩阵中M的选择不同，重构声压与理想声压之间的相对误差关系如图5所示。

在图6中，随着M逐渐变大，声场的相对误差逐渐减小，当M＝32后，相对误差控制在1％之内，并趋于稳定。依据式(8)，M越小，实际使用传声器数目越少，故选择M为32。当观测信号维数M确定之后，理想声压与重构球面声压如图7和图8所示。

由图7和图8可见，在f＝3kHz，r＝4cm，N＝6，M＝32时，使用压缩感知方法可较好地重构球面上声压。若以相同的实验条件，在f＝630Hz、2kHz、4kHz、5kHz条件下分析，重构声压如图7所示。

由图9～图12可见，在相同条件下，越靠近上限频率，压缩感知方法重构采集区域A上的声压分布的误差也逐渐变大，但总体上，在上限频率范围内都能较好地重建。

分别使用本发明提出的压缩感知方法与传统的加权求和法、最小二乘法作比较，取各频率下重构声压的相对误差最大值作分析，相对误差曲线如图13所示。在图13中，方块点曲线为r＝4厘米，N＝6，M＝32时，用压缩感知方案重构球面A上声场声压与理想声压之间的各频点上最大相对误差；圆形点曲线则反映使用加权求和法时重构相同声场声压与理想声压之间的各频率上的最大相对误差；三角点曲线分反映最小二乘法，后两种方法都使用49个传声器采集信号。从图8中可看到，越靠近上限频率，误差越大，这因为频率越大，高阶项项对当前阶的球谐波系数影响大，增加了重构误差，而压缩感知方案可以在一定频率范围内，使重构声场的相对误差基本控制在10％以内，尤其在3kHz范围内能控制在1％内。本发明方案与最小二乘法精度相近，但是比最小二乘法所使用传声器少，少于后者的1/3。

从实验结果看出，使用32个传声器就能重构的49个球谐波基系数，同时在重构声压的相对误差为10％的范围内重构整个球面A的声压，传声器的数目比传统的减少了1/3。使用布局在半径为4厘米的球传声器阵列，压缩感知方法能将一定频带下的高阶声场声压分布低相对误差地重构，尤其对低频成分最为明显。

本发明涉及到的缩略语和关键术语定义如下：

Ambisonics：高保真立体声响复制。

Fourier-Bessel级数：傅里叶-贝塞尔级数。

Legendre函数：勒让德函数。

δ函数：狄拉克函数。

Bessel函数：贝塞耳函数。

Legendre函数：勒让德函数。

Claims (4)

1.一种基于压缩感知的声场参数获取方法，其特征在于：包括球型传声器阵列模块，常数观测矩阵生成模块，观测信号向量生成模块，正交基构造模块，随机观测矩阵生成模块，球谐波基系数重构模块，目标区域声压分布重构模块；所述的球型传声器阵列模块的输出端、常数观测矩阵生成模块的输出端连接到观测信号向量生成模块的输入端，所述的正交基构造模块的输出端连接到随机观测矩阵生成模块的输入端，所述的观测信号向量生成模块的输出端、正交基构造模块的输出端、随机观测矩阵生成模块的输出端连接到球谐波基系数重构模块的输入端，所述的球谐波基系数重构模块的输出端连接到目标区域声压分布重构模块的输入端；将所述的球型传声器阵列模块安放在球形采集面处，所述的球型传声器阵列模块采集球形采集面上的声压信号与所述的常数观测矩阵生成模块生成的常数观测矩阵在观测信号向量生成模块中相乘，得到观测信号向量，所述的正交基构造模块输出的正交矩阵输入到随机观测矩阵生成模块，用于选择适当的随机观测矩阵，所述的随机观测矩阵生成模块生成的随机观测矩阵、所述的观测信号向量生成模块生成的观测信号向量在球谐波基系数重构模块中构造出球谐波基系数，根据所述的球谐波基系数重构模块构造的球谐波基系数在所述的目标区域声压分布重构模块中计算出球形采集面处任意位置处的声压，达到重构出目标区域声压分布的目的；

所述球型传声器阵列模块处理方法如下：

在声场采集区域中，把M个传声器等立体角地均匀布置在给定球半径r的球形采集面就构成了所述的球型传声器阵列模块；其中，M、r由采集声场的精度确定；传声器的位置由以球形采集面的球心为原点，极坐标系下不同的方位角仰角θ描述；

常数观测矩阵生成模块设计如下：

对一封闭球体空间内任意一点r处的声压能用Fourier-Bessel级数展开为

其中，r处的极坐标各个分量定义为：球半径r，方位角仰角θ；n为球谐波展开阶数，它决定重现声场时所需的声道数目，m＝-n,…,n，为n阶球谐波展开阶数对应的球谐波基，只与方位角仰角θ有关，A_nm(f)为与频率f有关的展开系数，j_n()为第一类n阶贝塞耳函数，c为声速，Y_nm()为球谐波，表示为

其中，P_n|m|(x)为连带勒让德函数，

B_nm(f,r)＝A_nm(f)j_n(2πfr/c)为一组球谐波系数，只与r有关，与无关，球谐波基具有正交性

其中，δ_nb为狄拉克函数，Λ为单位球面；由此可得单位圆上可积分的任意函数在球谐波域上的Fourier变换存在，其变换对为

式中，γ_nm为球谐波域的Fourier变换结果，也称为球谐波域的谱；

对式(1)中球谐波展开阶数n进行截断，即n的上限为N，

式中，角频率Ω＝2πf，球谐波系数B_nm(f,r)中含有第一类n阶贝塞耳函数j_n(kr)，k＝2πf/c，这里j_n(kr)与频率f和球半径r有关；当声场截断为N阶声场时，根据采集到的N阶声场声压信息，计算B_nm(f,r)中只与频率f有关的部分A_nm(f)；为使计算得到的A_nm(f)误差最小，应该保证kr≤arg(N)，这里arg(N)为一常数，其计算方法为：当K>N时，K阶j_K(kr)比N阶j_N(kr)衰减达到10dB以上时对应参数kr设为arg(N)，N＝3～8，其中arg(N)的具体值为{arg(3)＝2.63,arg(4)＝3.21,arg(5)＝3.78,arg(6)＝4.36,arg(7)＝4.94,arg(8)＝5.51}；当阶数N一定时，期望能准确获取的频率f越大，所需球半径r就越小；根据球型传声器阵列中每i个传声器的方位角仰角θ；计算出的M×M维常数观测矩阵C为，

所述正交基构造模块处理方法如下：

要用压缩感知理论，首先要构造出一个正交基，即在球面上选择L＝(N+1)²点处采集声压，一次构成L维的声压信号向量p，其每一维即每一个点处的声压用式(2)表示；以矩阵形式表示p，则有

其中，为N阶球谐波基构成的L×L维矩阵，为各球谐波基系数构成的L维矢量；在引入加权因子的条件下，有

式中，α_i为加权因子,

写成矩阵形式

其中，Ψ为近似的球谐波基矩阵，E为单位对角阵；

所述随机观测矩阵生成模块处理方法如下：

定义声压信号向量p在球谐波基矩阵Ψ上的投影为γ，使用随机数发生器，生成M×L维高斯随机矩阵Φ，高斯随机矩阵Φ均值为0，方差为Φ的选择应使其与Ψ之间具有较高不相关性，Φ的任意列向量与Ψ的任意列向量之间相关系数越小越好，小于0.1；用Φ将γ降维为M维观测信号x，即

x＝Φγ。

2.根据权利要求1所述的一种基于压缩感知的声场参数获取方法，其特征在于：观测信号向量生成模块处理方法如下：

考虑到物理可实现性及阵列小型化便于操作，采集区域半径设定为r＝4厘米；同时，在r＝4厘米的球面上等立体角地均匀布置M个传声器，M<(N+1)²，用于采集这M个位置处由不同及θ所区分的声压p_m，m＝1,2,…,M，构成声压信号向量p＝[p₁,p₂,…,p_M]，然后转换为观测信号向量x，

x＝Cp，

其中，C为常数观测矩阵生成模块设计出的M×M维的常数观测矩阵C。

3.根据权利要求1所述的一种基于压缩感知的声场参数获取方法，其特征在于：球谐波基系数重构模块处理方法如下：

根据压缩感知理论，在l₁-范数约束条件下，x能重构γ，即式(4)表示的优化问题有解，

Minimize[||γ||_1范数]s.t.Φγ＝x，(4)

将式(4)表示的优化问题转化为线性规划问题：

Minimize[sum(u)]s.t.-u≤γ≤u，Φγ＝x，(5)

其中，s.t.表示是subjectto的缩写；u＝0.95·abs(γ)+0.1·max[abs(γ)]；Φ为M×L维高斯随机矩阵，均值为0，方差为Minimize[h(·)]表示使函数h(·)最小化，max(·)表示取向量元素中的最大值，abs(·)表示向量元素求绝对值，sum(·)表示向量元素求和，式(5)能用原对偶内点算法求解，具体步骤如下：

(1)根据线性规划，设定线性函数为这里γ为原始变量向量，初始值为Φ^Tx，同时设置迭代误差阈值为ε＝10^-3，迭代最大次数为50；设定第k次迭代时的线性函数k表示第k次迭代；

(2)计算补偿间隙其中为一对偶变量向量，同时计算出中心参数其中，μ＝10为障碍常数，m＝2L；

(3)依据Karush-Kuhn-Tucker条件，计算第k次迭代的原始残差向量r_pri ^(k)、对偶残差向量r_dual ^(k)以及中心残差向量r_cent ^(k)

r_pri ^(k)＝Φγ^(k)-x^(k)，(6)

r_cent ^(k)＝-diag(λ^(k))g(γ^(k))-(1/t^(k))E，(8)

其中，v^(k)＝Φ^T[λ₂ ^(k)-λ₁ ^(k)]，为一对偶变数向量；g₀(γ^(k))＝sum(u^(k))，g₁(γ^(k))＝γ^(k)-u^(k)，g₂(γ^(k))＝-γ^(k)-u^(k)；为对γ^(k)与u^(k)求偏导，即 E为单位阵；

所以r_dual ^(k)与r_cent ^(k)具体计算式为，

其中，‘0’为元素皆为0的向量，‘1’为元素皆为1的向量，‘-1’为元素皆为-1的向量；

原对偶搜索方向需要使原始残差向量r_pri ^(k)、对偶残差向量r_dual ^(k)以及中心残差向量r_cent ^(k)接近0，根据式(6)～(8)所示的泰勒展开，得到

即，

可计算出搜索方向(Δγ^(k),Δλ^(k),Δv^(k))，同时，设定线性搜索的步长s^(k)＝0.99min{1,min[-λ^(k)/Δλ^(k)|Δλ^(k)<0]}；具体步骤如下：

(A).开始线性搜索，将γ^(k)+s^(k)Δγ^(k)、λ^(k)+s^(k)Δλ^(k)、v^(k)+s^(k)Δv^(k)替代γ^(k)、λ^(k)、v^(k)后计算新的g(γ^(k))、r_pri ^(k)、r_dual ^(k)、r_cent ^(k)，为了有所区别，分别用g'(λ^(k))、r'_pri ^(k)、r'_dual ^(k)、r'_cent ^(k)表示，若α的取值范围是0<α<1，则结束线性搜索，进入第(B)步，否则继续第(A)步，同时令βs^(k)取代s^(k)，β＝0.5；

(B).在线性搜索后，第k+1次迭代的初始值γ^(k+1)、λ^(k+1)、v^(k+1)分别为γ^(k)+s^(k)Δγ^(k)、λ^(k)+s^(k)Δλ^(k)、v^(k)+s^(k)Δv^(k)，重新计算补偿间隙若或迭代次数超过50次，则迭代截止，否则从第(3)步继续迭代。

4.根据权利要求3所述的一种基于压缩感知的声场参数获取方法，其特征在于：目标区域声压分布重构模块处理方法如下：

根据重构出的γ以及式(2)，可获得整个球面上的声压分布，这样，使用M个传声器获取到的声场信息来重构N阶的声场信息，M<(N+1)²，从而可减少传感器的数目，更进一步，利用重构的球面上的声压分布信息以及记录的音频信号，在重放空间可高度近似重现真实的声场。