CN102147458B - 一种针对宽带声源的波达方向估计方法及其装置 - Google Patents

Abstract

本发明涉及一种针对宽带声源的波达方向估计方法及装置，从传声器阵列采集的声音数据中，选取一段声音数据X(n)＝[x₁(n)，...，x_P(n)]^T，n∈Z^*，选择的声音数据X(n)均匀分成L帧，对每一帧声音数据X_l(n)做N点快速傅里叶变换，得到声音数据的频域表示X_l(k)＝[x_l1(k)，...，x_lP(k)]^T，k＝1，...，N/2；然后计算传声器阵列中的两个传声器间相干函数的模平方η(k)；从η(k)中随机挑选出q个相干函数的模平方值，仅使用这q个相干函数的模平方对应的频带进行波达方向估计。本发明在保证精度基本不降低的情况下极大降低算法的运算量。

Description

一种针对宽带声源的波达方向估计方法及其装置

技术领域

本发明涉及目标测向技术领域，特别涉及一种针对宽带声源的波达方向估计方法及其装置。

背景技术

目标测向在声呐、雷达、导航及无线传感器网络等领域都有重要的应用。而传感器阵列是目标测向的主要方式。

在目标测向技术中，一个重要技术问题是如何对宽带声源进行定位。对于这方面已经有很多的算法。宽带MUSIC算法是先将宽带信号通过傅里叶变换到频域上，然后对每个频带用窄带的MUSIC法处理，最后将各个频带的结果综合得到声源的方向估计。即非相干子空间方法(ISM，Incoherent Sub-space Method)。这类方法运算量小，但精度差，并且只能处理非相干信号；最大似然(ML，Maximum Likelihood)估计算法为代表的拟合类算法，以其良好的渐进性，灵活的适用性，可低的计算复杂度，得到了广泛关注。针对宽带信号，Yao Kung教授等人提出了针对宽带声源的AML(Approximated Maximum Likelihood)波达方向估计算法。这类算法精度高，能处理相干信号，但运算量很大。同时，对于相干信号源还可以用聚焦矩阵的方法处理。首先需要利用低复杂度算法对声源方向进行一个预估，通过这个预估的角度构造聚焦矩阵，通过聚焦矩阵把不同频带的信号聚焦到一个中心频率上，然后把它当成一个窄带信号来处理，这种方法的缺点是需要知道声源的预估角度，并且预估角度的选择对于算法的性能影响很大。

发明内容

本发明的目的在于，在保证精度基本不降低的情况下极大降低算法的运算量。

为达到上述目的，提出一种针对宽带声源的波达方向估计方法，该方法具体步骤包括：

步骤1)：从传声器阵列采集的声音数据中，选取一段声音数据X(n)＝[x₁(n)，...，x_P(n)]^T，其中，P表示传声器阵列中传声器的个数，P≥2；将该段声音数据X(n)均匀分成L帧，每帧的数据X_l(n)的长度为N，X(n)的数据长度为N×L，l＝1，...，L；n∈Z^*；

对每一帧声音数据X_l(n)做N点快速傅里叶变换，得到声音数据的频域表示X_l(k)＝[x_l1(k)，...，x_lP(k)]^T，其中k＝1，2，...，N/2；

步骤2)：根据所述的步骤1)得到的X_l(k)按照下式(1)计算传声器阵列中的两个传声器间相干函数的模平方η(k)；

$\mathrm{\η}\left(k\right)={\left|{\mathrm{\γ}}_{i,j}\left(k\right)\right|}^2=\frac{{\left|{\mathrm{\Φ}}_{i,j}\left(k\right)\right|}^2}{{\mathrm{\Φ}}_{i,i}\left(k\right){\mathrm{\Φ}}_{j,j}\left(k\right)}---\left(1\right)$

其中，k＝1，2，...，N/2；i，j∈[1，2，...，P]；Φ_i，j(k)为第i个传声器和第j个传声器的互功率谱，Φ_i，i(k)为第i个传声器的自功率谱，

Φ_j，j(k)为第j个传声器的自功率谱，

步骤3)：根据所述的步骤2)得到的η(k)随机选出q个相干函数模的平方；

步骤4)：根据所述的步骤3)选出的q个相干函数模的平方对应的声音数据X_l(k)进行波达方向估计，其中，k∈[k₁，...，k_q]，k_q表示频带序号，且[k₁，...，k_q]∈[1，...，N/2]。

所述的步骤1)还包括：保存所述的传声器阵列到当前为止采集的连续L帧声音数据；并保存来自所述传声器阵列最新的1帧声音数据后，将保存的连续L帧声音数据中最早的1帧声音数据删除。

所述的步骤2)中，两个传声器为传声器阵列中空间距离最大的两个传声器。

所述的步骤2)中，两个传声器之间连续(L-1)帧的数据作为数据缓存，根据当前获得第L帧的N个数据和缓存的N×(L-1)个数据来求取相干函数的模平方。

所述的步骤4)中，根据所述的步骤3)选出的q个相干函数模的平方对应的第L帧声音数据X_L(k)进行波达方向估计，其中，k∈[k₁，...，k_q]，k_q表示频带序号，且[k₁，...，k_q]∈[1，...，N/2]。

所述的步骤3)中，比较所述的步骤2)中得到的η(k)大小，从最大的开始依次挑选前q个η(k)值，即η(k₁)≥η(k₂)≥...≥η(k_q)≥η(k_m)，其中，k_q表示频带序号，[k₁，...，k_q]∈[1，2，....，N/2]，

为实现上述目的，提出一种针对宽带声源的波达方向估计装置，其特征在于，该装置包括：傅里叶变换模块、相干函数的模平方模块和波达方向估计模块；

所述的傅里叶变换模块，包括：数据选择单元、数据分帧单元和傅里叶变换单元；用于所述的数据选择单元从传声器阵列采集的声音数据中，选取一段声音数据X(n)＝[x₁(n)，...，x_P(n)]^T，所述的数据分帧单元将所述的数据选择单元选择的声音数据X(n)均匀分成L帧，所述的傅里叶变换单元对每一帧声音数据X_l(n)做N点快速傅里叶变换，得到声音数据的频域表示X_l(k)＝[x_l1(k)，...，x_lP(k)]^T；并将所述频域的声音数据输出至所述的相干函数的模平方模块；其中，P表示传声器阵列中传声器的个数，P≥2；每帧的数据X_l(n)的长度为N，l＝1，...，L，X(n)的数据长度为N×L；n∈Z^*；k＝1，2，...，N/2；

所述的相干函数的模平方模块，用于根据所述的傅里叶变换模块得到的X_l(k)按照式(2)计算传声器阵列中的两个传声器间相干函数的模平方η(k)；

$\mathrm{\η}\left(k\right)={\left|{\mathrm{\γ}}_{i,j}\left(k\right)\right|}^2=\frac{{\left|{\mathrm{\Φ}}_{i,j}\left(k\right)\right|}^2}{{\mathrm{\Φ}}_{i,i}\left(k\right){\mathrm{\Φ}}_{j,j}\left(k\right)}---\left(2\right)$

其中，k＝1，2，...，N/2；i，j∈[1，2，...，P]；Φ_i，j(k)为第i个传声器和第j个传声器的互功率谱，

Φ_i，i(k)为第i个传声器的自功率谱，

Φ_j，j(k)为第j个传声器的自功率谱，

所述的波达方向估计模块，用于根据从所述的相干函数的模平方模块得到的相干函数的模平方中随机选出q个相干函数模的平方对应的声音数据X_l(k)进行波达方向估计，其中，k∈[k₁，...，k_q]，k_q表示频带序号，且[k₁，...，k_q]∈[1，...，N/2]。

所述的相干函数的模平方模块中，两个传声器为传声器阵列中空间距离最大的两个传声器。

所述的装置还包括声音数据缓存模块；所述的声音数据缓存模块，用于保存传声器阵列到当前为止采集的连续L帧声音数据；并保存接收到来自传声器阵列最新的1帧声音数据后，将保存的连续L帧声音数据中最早的1帧声音数据删除；并将缓存的声音数据输出至所述的傅里叶变换模块。

所述的相干函数的模平方模块将选取的两个传声器之间连续(L-1)帧的数据作为数据缓存，根据当前获得第L帧的N个数据和缓存的N×(L-1)个数据来求取相干函数的模平方。

所述的波达方向估计模块根据所述的比较模块得到的q个相干函数模的平方对应的第L帧声音数据X_L(k)进行波达方向估计，其中，k∈[k₁，...，k_q]，k_q表示频带序号，且[k₁，...，k_q]∈[1，...，N/2]。

该装置还包括一比较模块，用于根据比较所述的相干函数的模平方模块得到的η(k)大小，从最大的开始依次挑选前q个η(k)值，即η(k₁)≥η(k₂)≥...≥η(k_q)≥η(k_m)，其中，k_q表示频带序号，[k₁，...，k_q]∈[1，2，....，N/2]，

；将挑选出的相干函数的模平方输出至所述的波达方向估计模块。

本发明的优点在于，本发明可以大大降低运算量。同时这种简化算法可适用于背景技术部分所讨论的任何算法。以AML算法为例，由于AML算法的运算量最大，通过验证可以得出在定向效果相当的情况下，大大降低了运算量。

附图说明

图1为传声器阵列的空间结构示意图；

图2为本发明提出的一种针对宽带声源的波达方向估计方法流程图；

图3为仿真试验中白噪声情况下本发明实施例方案与全频带方案进行波达方向估计的结果比较图；

图4为仿真试验中有色噪声情况下本发明实施例方案与全频带方案进行波达方向估计的结果比较图；

图5为野外试验中本发明实施例方案与全频带方案进行波达方向估计的结果比较图。

具体实施方式

下面结合附图和实施例对本发明进一步说明。

设有R个全向(omni-directional)传声器，第p个传声器接收信号为：

$x_p\left(n\right)=\underset{m=1}{\overset{M}{\mathrm{\Σ}}}{s}_0^{\left(m\right)}(n-t_p^{\left(m\right)})+w_p\left(n\right)---\left(1\right)$

其中：n＝0，...，L-1，L为一帧信号采样点个数；p＝1，...，R，R为接收阵列中传声器的个数；m＝1，...，M，M为声源数量(M＜R)。为源信号，为用采样点数表示的时间延迟。w_p(n)为阵列接收信号中噪声分量。其中，时间延迟

根据阵型及目标声源位置的几何关系获得。对于一个远场条件下的R元圆阵，即R个传声器等距排列在一个圆周上，第p个传声器的用采样点数表示的时间延迟为：

$t_p^{\left(m\right)}=\frac{f_{s}r}{c}\cos (\frac{2\mathrm{\π}(p-1)}{R}-{\mathrm{\θ}}_m)---\left(2\right)$

其中，f_s为信号采样频率，r为圆阵半径，c为各向同性介质中的波速，θ_m为第m个源的远场方位角。

对接收信号做N点的FFT变换，得频率域信号模型：

X(k)＝D(Θ，k)S₀(k)+η(k)    (3)

式中，k＝0，...，N/2；频域阵列输出信号：X(k)＝[X₁(k)，...，X_R(k)]^T；M个声源的频谱：

方向矩阵(Steering Matrix)为：D(Θ，k)＝[d⁽¹⁾(k)，...，d^(M)(k)]；方向矢量(Steering Vector)为：

η(k)为噪声频谱。

噪声谱向量η(k)是方差为Lσ²零均值复数高斯白噪声。注意到即使w(n)为任意独立同分布的随机变量，根据中心极限定理，都可推知η(k)服从高斯分布。

定义    S(k)＝D(k)S₀(k)     (4)

得到NR/2维空间-时频向量：

X＝G(θ)+ξ₁                (5)

其中，G(θ)＝[S(0)^T，...，S(N/2)^T]^T，并且有：

$R_{\mathrm{\ξ}}=E\left[{\mathrm{\ξ\ξ}}^H\right]=L\left[\begin{array}{cccc}{\mathrm{\σ}}^2{I}_R& & & 0\\ & {\mathrm{\σ}}^2{I}_R& & \\ & & L& \\ 0& & & {\mathrm{\σ}}^2{I}_R\end{array}\right]---\left(6\right)$

其中，上角标^H代表复矩阵的共轭转置。

假设初始未知参数空间为：

$\mathrm{\Θ}={[{\mathrm{\θ}}_s^T,S_0^{{\left(1\right)}^T},...,S_0^{{\left(M\right)}^T}]}^T---\left(7\right)$

其中，源方位角

源频谱

所以有最大似然函数：

$f\left(\mathrm{\Θ}\right)=\frac{1}{{\left(2\mathrm{\π}\right)}^{\mathrm{NR}/2/2}{\left|R_{\mathrm{\ξ}}\right|}^{1/2}}\exp \{-\frac{1}{2}{[X-G\left(\mathrm{\Θ}\right)]}^H{R_{\mathrm{\ξ}}}^{-1}[X-G\left(\mathrm{\Θ}\right)\left]\right\}---\left(8\right)$

对数化该似然函数并略去常数项得：

$L\left(\mathrm{\Θ}\right)=-\frac{{\left|\right|X-G\left(\mathrm{\Θ}\right)\left|\right|}^2}{R_{\mathrm{\ξ}}}---\left(9\right)$

可以得到最优化准则：

$\underset{\mathrm{\Θ}}{\max }L\left(\mathrm{\Θ}\right)=\underset{\mathrm{\Θ}}{\min }\underset{k=0}{\overset{N/2}{\mathrm{\Σ}}}\frac{{\left|\right|X\left(k\right)-D\left(k\right)S_0\left(k\right)\left|\right|}^2}{R_{\mathrm{\ξ}}\left(k\right)}$

$=\underset{\mathrm{\Θ}}{\min }\underset{k=0}{\overset{N/2}{\mathrm{\Σ}}}\frac{{\left|\right|X\left(k\right)-D\left(k\right)S_0\left(k\right)\left|\right|}^2}{{\mathrm{\σ}}^2}---\left(10\right)$

这等价于对所有k个函数，找到又由于σ²为常数，f(k)可简写为：

f(k)＝||X(k)-D(k)S₀(k)||²            (11)

使f(k)取的极小值的源信号矢量S₀(k)一定满足因此为了对任何源定位有最小的残差，源矢量的估计给定：

其中，

是方向矩阵D(k)的伪逆矩阵(pseudo-inverse)。定义正交规划那么AML算法的波达方向估计可以通过下式得到。

${\mathrm{\θ}}_s=\underset{{\mathrm{\θ}}_s}{\max }\underset{k=0}{\overset{N/2}{\mathrm{\Σ}}}{\left|\right|P(k,{\mathrm{\θ}}_s)X\left(k\right)\left|\right|}^2---\left(13\right)$

声音信号虽然是宽带的，但信号一般都不可能是白的，它一般都存在一些谐频结构，即信号主要集中在某些频段上；或由于噪声是有色的，某些频带的信号完全被噪声湮没掉了，因此在计算时，这些频带上的信息是冗余的，甚至会影响最后的结果。所以采用一定的方法来选取某些频带是必要的，本发明方案中，将频带上的相干函数作为选择频带的标准。

相干函数(coherence)用于描述两组数据之间的相关性。本发明方案中，相干函数的定义为：

${\mathrm{\γ}}_{\mathrm{xy}}=\frac{G_{\mathrm{xy}}\left(f\right)}{\sqrt{G_{\mathrm{xx}}\left(f\right)G_{\mathrm{yy}}\left(f\right)}}---\left(14\right)$

其中，x和y代表两组不同的数据，G_xy(f)为x和y的互谱，G_xx(f)和G_yy(f)分别为x和y的自谱。G_xy(f)＝E{X(f)Y^*(f)}，G_xx(f)＝E{X(f)X^*(f)}，G_yy(f)＝E{Y(f)Y^*(f)}。

幅度平方相干函数(MSC)为相干函数γ_xy的模的平方：

${\left|{\mathrm{\γ}}_{\mathrm{xy}}\left(f\right)\right|}^2=\frac{{\left|G_{\mathrm{xy}}\left(f\right)\right|}^2}{G_{\mathrm{xx}}\left(f\right)G_{\mathrm{yy}}\left(f\right)}---\left(15\right)$

由此，可见相干函数指的都是频域上的运算。

相干函数的取值为0到1之间，即：

0≤|γ_xy|²≤1                (16)

当取值为0时表示频段x和y上的频谱完全不相关，取值为1表示频段x和y上的频谱完全相关。

在特定的情况下有：信噪比

由此可见MSC可以用来衡量信号某个频带上的信噪比的大小。当|γ_xy(k)|²较大时，表明第k个频带上有用的信号较多；反之，则说明第k频带上噪声占主要部分。

传声器阵列的空间结构示意图，如图1所示。使用图1所示的十字阵型的传声器阵列进行风噪信号采集，圆圈101表示12枚B&K 4189型标准传声器，并配有直径为10cm的球形海绵风罩。102表示用于固定这些传声器的支架。十字阵上三个4元阵列的孔径分别是0.2米，0.6米，1.2米。采用SONY SIR-1000型16通道同步数字录音机，采样频率为48kHz。数据分析时，通过低通滤波、归一化、降采样到1kHz进行分析。

使用实测数据，统计阵列风噪声的空间相关性。采用下式计算通道间平均相关系数：

$\stackrel{\‾}{\mathrm{\γ}}\left(k\right)=\frac{2}{R(R-1)}\underset{i=1}{\overset{R}{\mathrm{\Σ}}}\underset{j=i+1}{\overset{R}{\mathrm{\Σ}}}\left|{\mathrm{\γ}}_{\mathrm{ij}}\left(k\right)\right|---\left(17\right)$

其中，k＝0，...，N/2；R为通道数，γ_ij(k)为第i通道与第j通道的互相关系数，

为平均相关系数。

表1  风噪声通道间平均相关系数

表1说明，与低频风噪声的空间相关性不同，中高频处的风噪声通道间相关系数小，空间相关性微弱，阵列风噪声基本可以认为通道间相互独立。

本发明方案是在采集了声音数据之后，根据相干系数对声音数据进行筛选，筛选出携带有用信号较多的频段的声音数据，而放弃其他频段的声音数据，这样一来，就减少了波达方向估计所需的数据量，从而减小后续的波达方向估计的运算量。

本发明提出的一种针对宽带声源的波达方向估计方法流程图；如图2所示，包括如下步骤：

步骤201：从由P个传声器组成的传声器阵列采集的声音数据中，选择一段声音数据x(n)＝[x₁(n)，...，x_P(n)]^T，其中，n∈Z^*。

步骤202：将所选择的声音数据x(n)均匀分成L段x_l(n)，l＝1，...，L，每段的数据x_l(n)的长度为N，x(n)的数据数据长度为N×L。

步骤203：对每一段声音数据x_l(n)做N点快速傅里叶变换(FFT)，得到声音数据的频域表示为X_l(k)＝[X_l1(k)，...，X_lP(k)]^T，其中，k＝1，2，...，N/2。

步骤204：计算传声器阵列中的两个传声器间的相干函数模的平方(MSC)

其中，i，j∈[1，2，...，P]，Φ_i，j(k)为第i个传声器和第j个传声器的互功率谱，Φ_i，i(k)和Φ_j，j(k)分别为第i个传声器和第j个传声器的自功率谱：

${\mathrm{\Φ}}_{i,i}\left(k\right)=\underset{l=1}{\overset{L}{\mathrm{\Σ}}}{\left|X_{\mathrm{li}}\left(k\right)\right|}^2;$ ${\mathrm{\Φ}}_{j,j}\left(k\right)=\underset{l=1}{\overset{L}{\mathrm{\Σ}}}{\left|X_{\mathrm{lj}}\left(k\right)\right|}^2;$ ${\mathrm{\Φ}}_{i,j}\left(k\right)=\underset{l=1}{\overset{L}{\mathrm{\Σ}}}{X}_{\mathrm{li}}\left(k\right)X_{\mathrm{lj}}^{*}\left(k\right).$

较佳地，所述两个传声器可以为传声器阵列中空间距离最大的两个传声器。

步骤205：从相干函数模的平方η(k)中挑选出较大的q个相干函数模的平方值，即η(k₁)≥η(k₂)≥...≥η(k_q)≥η(k_m)，其中，k_q为频带序号，[k₁，...，k_q]∈[1，2，....，N/2]，

步骤206：使用挑选出的q个相干函数模的平方对应的声音数据X_l(k)进行波达方向估计。其中，k∈[k₁，...，k_q]；以AML算法为例，则角度的计算公式由(18)可知，

${\mathrm{\θ}}_s=\underset{{\mathrm{\θ}}_s}{\max }\underset{k=k_1}{\overset{k_q}{\mathrm{\Σ}}}{\left|\right|P(k,{\mathrm{\θ}}_s)X_l\left(k\right)\left|\right|}^2---\left(18\right)$

由于相干函数的计算需要较多数据，即N×L要较大，因此可能影响方法的实时性。在相干函数的计算中，保存了之前连续(L-1)帧的数据作为数据缓存，根据当前获得的N个数据和缓存的N×(L-1)个数据求取MSC。并且每次获得新一帧数据之后，都丢弃最前一帧的数据，更新整个数据缓存。如***中的全部数据为x_l(n)，l＝1，...，L，其中x_L(n)为当前获得的一帧数据，而x_l(n)，l＝1，...，L-1为以前的数据。使用x_l(n)，l＝1，...，L算MSC，而仅使用x_L(n)进行角度估计。

本发明还提出一种针对宽带声源的波达方向估计的装置，包括：

傅立叶变换模块，用于对P个传声器组成的传声器阵列采集的声音数据进行快速傅立叶变换，得到声音数据的频域表示，并将所述频域声音数据输出至相干函数的模平方模块；

相干函数的模平方模块，用于计算传声器阵列中两个传声器间的声音数据的相干函数的模平方，并将计算得到的相干函数的模平方输出至比较模块；

比较模块，用于从相干函数的模平方中挑选出较大的q个值；

波达方向估计模块，用于根据比较模块所挑选出的相干函数的模平方对应的声音数据的频带进行波达方向估计。

较佳地，所述相干函数的模平方模块用于计算传声器阵列中距离最远的两个传声器间声音数据的相干函数的模平方，并将计算得到的相干函数的模平方输出至比较模块。

较佳地，所述傅立叶变换模块包括：

数据选择单元，用于从由P个传声器组成的传声器阵列采集的声音数据中选择长度为N×L的一段声音数据作为待处理声音数据x(n)＝[x₁(n)，...，x_P(n)]^T，其中，n∈Z^*；

数据分帧单元，用于将所述数据选择单元所选择的待处理声音数据x(n)分成L段，其中任一段的声音数据为x_l(n)，l＝1，...，L，并将分段后的声音数据输出；以及

傅立叶变换单元，用于对所述数据分帧单元输出的每一段声音数据x_l(n)做N点快速傅里叶变换，得到声音数据的频域表示为X_l(k)＝[X_l1(k)，...，X_lP(k)]^T，其中，k＝1，2，...，N/2。

较佳地，相干函数的模平方模块计算传声器阵列第i个传声器和第j个传声器之间的相干函数的模平方为：

其中，i，j∈[1，2，...，P]，Φ_i，j(k)为第i个传声器和第j个传声器的互功率谱，Φ_i，i(k)和Φ_j，j(k)分别为第i个传声器和第j个传声器的自功率谱。

较佳地，所述比较模块从N/2个相干函数的模平方η(k)中挑选出较大的q个相干函数的模平方值，即η(k₁)≥η(k₂)≥...≥η(k_q)，k_q表示频带序号，k＝1，2，...，N/2，[k₁，...，k_q]∈[1，2，....，N/2]。

较佳地，所述DOA估计模块根据比较模块所挑选出的相干函数的模平方对应的声音数据的频带进行波达方向估计。

较佳地，该装置进一步包括：

声音数据缓存模块，用于保存所述传声器阵列到当前为止采集的连续L帧声音数据；并当接收到来自所述传声器阵列最新的1帧声音数据后，将保存的连续L帧声音数据中最早的1帧声音数据删除。

所述傅立叶变换模块用于对声音数据缓存模块中存储的声音数据进行行快速傅立叶变换得到声音数据的频域表示，并将所述频域声音数据输出至相干函数的模平方模块。

以下通过仿真实验以及实地实验的结果对本发明方案的技术效果进行说明。在仿真试验中，传声器组阵列为圆阵，阵元数为4，圆阵的半径为0.2m，声源为白噪声，个数为1，入射角为70°。使用全频带的AML算法和频带选取的AML算法结果对比分别如图3和图4所示，其中，图3为白噪声的情况，图4为有色噪声的情况。其中有色噪声是由白噪声通过一个50阶的FIR滤波器所产生的。

从图3可以看出，当背景为白噪声的条件时，低信噪比下全频带的AML算法要优于频带选取的AML算法；高信噪比下这两种算法的偏差相近都趋向于零。从图4可以看出，当背景为有色噪声时，低信噪比使用全频带的AML算法已经不是最优的，从图中可以看出在不同信噪比下应选取不同的频带数，这是由于背景为有色噪声，在某些频带上信噪比会很大，而在某些频带上信噪比会很小，因此选取过多的频带会使这些信噪比低的频带也参与运算，从而使算法的偏差增大；高信噪比下，两种算法的偏差也都趋近于0。

在野外实际实验中，阵列为圆阵，阵元数为4，圆阵半径为0.2m，声源为坦克绕阵列做圆周运动。实验结果如图5所示，图标“*”是全频带搜索的结果，图标“o”是只取相干系数最大的20个频带计算所得的结果，从图中可以看出，两者定位精度基本一样，但使用相干系数选择频带后，运算量是原来的1/15。

通过以上的实施方式的描述，本领域的技术人员可以清楚地了解到本发明可借助软件加必需的硬件平台的方式来实现，当然也可以全部通过硬件来实施，但很多情况下前者是更佳的实施方式。基于这样的理解，本发明的技术方案对背景技术做出贡献的全部或者部分可以以软件产品的形式体现出来，该计算机软件产品可以存储在存储介质中，如ROM/RAM、磁碟、光盘等，包括若干指令用以使得一台计算机设备执行本发明各个实施例或者实施例的某些部分所述的方法。计算机设备如：个人计算机，服务器，或者网络设备等。

最后所应说明的是，以上实施例仅用以说明本发明的技术方案而非限制。尽管参照实施例对本发明进行了详细说明，本领域的普通技术人员应当理解，对本发明的技术方案进行修改或者等同替换，都不脱离本发明技术方案的精神和范围，其均应涵盖在本发明的权利要求范围当中。

Claims (12)

1.一种针对宽带声源的波达方向估计方法，该方法具体步骤包括：

步骤1）：从传声器阵列采集的声音数据中，选取一段声音数据X(n)=[x₁(n),...,x_P(n)]^T，其中，P表示传声器阵列中传声器的个数，P≥2；将该段声音数据X(n)均匀分成L帧，每帧的数据X_l(n)的长度为N，X(n)的数据长度为N×L，l=1,...,L；n∈Z^*；

对每一帧声音数据X_l(n)做N点快速傅里叶变换，得到声音数据的频域表示X_l(k)=[x_l1(k),...,x_lP(k)]^T，其中k=1,2,...,N/2；

步骤2）：根据所述的步骤1）得到的X_l(k)按照下式（1）计算传声器阵列中的两个传声器间相干函数的模平方η(k)；

$\mathrm{\η}\left(k\right)={\left|{\mathrm{\γ}}_{i,j\left(k\right)}\right|}^2=\frac{{\left|{\mathrm{\Φ}}_{i,j}\left(k\right)\right|}^2}{{\mathrm{\Φ}}_{i,i}\left(k\right){\mathrm{\Φ}}_{j,j}\left(k\right)}---\left(1\right)$

其中，k=1,2,...,N/2；i,j∈[1,2,...,P]；Φ_i,j(k)为第i个传声器和第j个传声器的互功率谱，

Φ_i,i(k)为第i个传声器的自功率谱， ${\mathrm{\Φ}}_{i,i}\left(k\right)=\underset{l=1}{\overset{L}{\mathrm{\Σ}}}{|x_{\mathrm{li}}\left(k\right)|}^2;$ Φ_j,j(k)为第j个传声器的自功率谱， ${\mathrm{\Φ}}_{j,j}\left(k\right)=\underset{l=1}{\overset{L}{\mathrm{\Σ}}}{\left|x_{\mathrm{lj}}\right|}^2;$

步骤3）：根据所述的步骤2）得到的η(k)从最大的开始依次选出q个相干函数模的平方；

步骤4）：根据所述的步骤3）选出的q个相干函数模的平方对应的声音数据X_l(k)进行波达方向估计，其中，k∈[k₁,...,k_q]，k_q表示频带序号，且[k₁，...,k_q]∈[1,...,N/2]。

2.根据权利要求1所述的针对宽带声源的波达方向估计方法，其特征在于，所述的步骤1）还包括：保存所述的传声器阵列到当前为止采集的连续L帧声音数据；并保存来自所述传声器阵列最新的1帧声音数据后，将保存的连续L帧声音数据中最早的1帧声音数据删除。

3.根据权利要求1所述的针对宽带声源的波达方向估计方法，其特征在于，所述的步骤2）中，两个传声器为传声器阵列中空间距离最大的两个传声器。

4.根据权利要求2所述的针对宽带声源的波达方向估计方法，其特征在于，所述的步骤2）中，两个传声器之间连续（L-1）帧的数据作为数据缓存，根据当前获得第L帧的N个数据和缓存的N×（L-1）个数据来求取相干函数的模平方。

5.根据权利要求4所述的针对宽带声源的波达方向估计方法，其特征在于，所述的步骤4）中，根据所述的步骤3）选出的q个相干函数模的平方对应的第L帧声音数据X_L(k)进行波达方向估计，其中，k∈[k₁,...,k_q]，k_q表示频带序号，且[k₁,...,k_q]∈[1,...,N/2]。

6.根据权利要求1所述的针对宽带声源的波达方向估计方法，其特征在于，所述的步骤3）中，比较所述的步骤2）中得到的η(k)大小，从最大的开始依次挑选前q个η(k)值，即η(k₁)≥η(k₂)≥...≥η(k_q)≥η(k_m)，其中，k_q表示频带序号，[k₁，...,k_q∈[1，2，....,N/2]，

7.一种针对宽带声源的波达方向估计装置，其特征在于，该装置包括：傅里叶变换模块、相干函数的模平方模块和波达方向估计模块；

所述的傅里叶变换模块，包括：数据选择单元、数据分帧单元和傅里叶变换单元；所述的数据选择单元从传声器阵列采集的声音数据中，选取一段声音数据X(n)=[x₁(n)，...,x_P(n)]^T，所述的数据分帧单元将所述的数据选择单元选择的声音数据X(n)均匀分成L帧，所述的傅里叶变换单元对每一帧声音数据X_l(n)做N点快速傅里叶变换，得到声音数据的频域表示X_l(k)=[x_l1(k),...,x_lP(k)]^T；并将所述频域声音数据输出至所述的相干函数的模平方模块；其中，P表示传声器阵列中传声器的个数，P≥2；每帧的数据X_l(n)的长度为N，l=1,...,L，X(n)的数据长度为N×L；n∈Z^*；k=1,2,...,N/2；

所述的相干函数的模平方模块，用于根据所述的傅里叶变换模块得到的X_l(k)按照式（2）计算传声器阵列中的两个传声器间相干函数的模平方η(k)；

$\mathrm{\η}\left(k\right)={\left|{\mathrm{\γ}}_{i,j\left(k\right)}\right|}^2=\frac{{\left|{\mathrm{\Φ}}_{i,j}\left(k\right)\right|}^2}{{\mathrm{\Φ}}_{i,i}\left(k\right){\mathrm{\Φ}}_{j,j}\left(k\right)}---\left(2\right)$

其中，k=1,2,...,N/2；i,j∈[1,2,...,P]；Φ_i,j(k)为第i个传声器和第j个传声器的互功率谱，

Φ_i,i(k)为第i个传声器的自功率谱，

${\mathrm{\Φ}}_{i,i}\left(k\right)=\underset{l=1}{\overset{L}{\mathrm{\Σ}}}|{{X_{\mathrm{lj}}\left(k\right)|}_{;}}^2$ Φ_j,j(k)为第j个传声器的自功率谱， ${\mathrm{\Φ}}_{j,j}\left(k\right)=\underset{l=1}{\overset{L}{\mathrm{\Σ}}}{\left|X_{\mathrm{lj}}\left(k\right)\right|}^2;$

所述的波达方向估计模块，用于根据从所述的相干函数的模平方模块得到的相干函数的模平方中从最大的开始依次选出q个相干函数模的平方对应的声音数据X_l(k)进行波达方向估计，其中，k∈[k₁,...,k_q]，k_q表示频带序号，且[k₁，...,k_q]∈[1,...,N/2]。

8.根据权利要求7所述的针对宽带声源的波达方向估计装置，其特征在于，所述的相干函数的模平方模块中，两个传声器为传声器阵列中空间距离最大的两个传声器。

9.根据权利要求7所述的针对宽带声源的波达方向估计装置，其特征在于，所述的装置还包括声音数据缓存模块；

所述的声音数据缓存模块，用于保存传声器阵列到当前为止采集的连续L帧声音数据；并保存接收到来自传声器阵列最新的1帧声音数据后，将保存的连续L帧声音数据中最早的1帧声音数据删除；并将缓存的声音数据输出至所述的傅里叶变换模块。

10.根据权利要求9所述的针对宽带声源的波达方向估计装置，其特征在于，所述的相干函数的模平方模块将选取的两个传声器之间连续（L-1）帧的数据作为数据缓存，根据当前获得第L帧的N个数据和缓存的N×（L-1）个数据来求取相干函数的模平方。

11.根据权利要求10所述的针对宽带声源的波达方向估计装置，其特征在于，所述的波达方向估计模块根据比较模块得到的q个相干函数模的平方对应的第L帧声音数据X_L(k)进行波达方向估计，其中，k∈[k₁,...,k_q]。

12.根据权利要求7所述的针对宽带声源的波达方向估计装置，其特征在于，该装置还包括一比较模块，用于根据比较所述的相干函数的模平方模块得到的η(k)大小，从最大的开始依次挑选前q个η(k)值，即η(k₁)≥η(k₂)≥...≥η(k_q)≥η(k_m)，其中，k_q表示频带序号，[k₁,...,k_q]∈[1,2,....,N/2]，

将挑选出的相干函数的模平方输出至所述的波达方向估计模块。

Images