CN105044675A - 一种srp声源定位的快速实现方法 - Google Patents

Abstract

本发明公开了一种SRP声源定位的快速实现方法，首先求得Farrow结构波束形成器的最优权值；然后计算属于Farrow结构波束形成器SRP-PHAT的相位加权函数；接着计算属于Farrow结构波束形成器SRP-PHAT所有麦克风两两麦克风的广义互相关函数之和并进行锁存；再次将锁存后的输出信号与Farrow结构波束形成器扫描角度唯一有关的调向参数D进行乘积相加运算，得到相应扫描点的瞬时功率，搜索最大峰值得到对应的定位角度θ_rˊ；最后进行消除Farrow结构波束形成器的系统误差的校准，得最终θ_r。本发明相比原SRP-PHAT运算量降低；且对混响和噪声具有更高的鲁棒性、方位角估计精度也较高。

Description

一种SRP声源定位的快速实现方法

技术领域

本发明涉及声源定位技术，特别是一种SRP声源定位的方法。

背景技术

声源定位技术主要分为3大类，即基于最大可控响应功率的波束形成方法、采用高分辨率谱估计的方法和基于时延估计的两步定位方法。2002年，J.H.DiBiase(见文献[1]J.H.DiBiase.Ahigh-accuracy,low-latencytechniquefortalkerlocalizationinreverberantenvironments[D].ProvidenceRI,USA:BrownUniversity,2000和文献[2]谭颖,殷福亮,李细林.改进的SRP-PHAT声源定位方法[J].电子与信息学报,2006,28(7)：1223～1227.见文献[3]JacekP.Dmochowski.AGeneralizedSteeredResponsePowerMethodforComputationallyViableSourceLocalization.IEEEtransactionsonaudio,speech,andlanguageprocessing,VOL.15,NO.8,2007)提出了联合可控响应功率和相位变换(SRP-PHAT)的声源定位方法，该方法将可控响应功率方法固有的鲁棒性、短时分析特征与时延估计中相位变换方法对信号周围环境的不敏感性相结合，使声源定位系统具有一定的抗噪性、抗混响性。但是,SRP-PHAT方法运算量庞大，并且在恶劣环境下(噪声干扰大、混响影响严重)性能急剧下降，因此在实时处理和实际应用时变得困难。

发明内容

发明目的：为了克服现有技术中存在的不足，本发明提供了一种SRP声源定位的快速实现方法，用于解决现有的SRP-PHAT方法运算量打且恶劣环境下性能差的技术问题。

技术方案：为实现上述目的，本发明采用的技术方案为：

一种SRP声源定位的快速实现方法，包括如下步骤：

步骤1、利用基于凸优化的波束指向可调的宽带稳健远场频率不变波束形成器算法得到Farrow结构波束形成器的最优权值，从而设计完成Farrow结构波束形成器；

步骤2、利用Farrow结构波束形成器在调整波束指向时无需重新计算权值这个结构特征，将各麦克风阵列的输入信号进行FFT运算，根据不同麦克风以及输入信号的阶数将FFT运算后的输入信号进行两两相乘、求得相位加权函数；

步骤3、将各麦克风阵列的输入信号输入至Farrow结构波束形成器中进行空域滤波，空域滤波后的输出信号进行FFT运算，根据不同麦克风以及输出信号的阶数将FFT运算后的输出信号进行两两相乘求得相位加权函数对应量，；

步骤4、根据两个麦克风以及各自输入信号的阶数所确定的相位加权函数以及相位加权函数对应量为一组进行对应，将属于同一组的相位加权函数与相位加权函数对应量进行相乘，然后将相乘的结果进行IFFT运算之后并累加求和得到所有的麦克风中的两两麦克风的广义互相关函数之和并进行锁存；

步骤5、将两两麦克风的广义互相关函数之和与扫描角度唯一有关的调向参数D进行乘积相加运算，得到相应扫描点的瞬时功率，搜索瞬时功率的最大峰值得到对应的指向角度θ_r'；

步骤6、对Farrow结构波束形成器的系统误差进行校准，并将步骤5中获得的指向角度θ_r'校准为最终的定位角度θ_r。

进一步的，在本发明中，所述步骤1的具体步骤为：基于凸优化的思想设计稳健波束指向可调的频率不变波束形成器，稳健的最大误差最小化的优化设计的代价函数：

$\underset{w}{\min }\underset{(\mathrm{\Φ},\mathrm{\θ},f)\∈Z}{\max }\underset{\left|\right|\mathrm{\Δ}\stackrel{\‾}{g}\left|\right|\≤\mathrm{\ϵ}}{\max }\underset{{\mathrm{\Φ}}_{d_i}\∈\mathrm{\Ψ}}{\Σ}\underset{f_k\∈\mathrm{\Ω}}{\Σ}\underset{{\mathrm{\θ}}_n\∈{{\mathrm{\Θ}}_p}^{{\mathrm{\Φ}}_{d_i}}}{\Σ}F(\mathrm{\Φ},\mathrm{\θ},f)|w^T\[g({\mathrm{\Φ}}_{d_i},f_k,{\mathrm{\θ}}_n)+\mathrm{\Δ}g({\mathrm{\Φ}}_{d_i},f_k,{\mathrm{\θ}}_n)\]-w^{T}g({\mathrm{\Φ}}_{d_i},f_r,{\mathrm{\θ}}_n)|$

$s.t.\{\begin{array}{c}\underset{{\mathrm{\&Phi;}}_{d_i}\&Element;\mathrm{\&Psi;}}{\&Sigma;}\underset{f_k\&Element;\mathrm{\&Omega;}}{\&Sigma;}\underset{{\mathrm{\&theta;}}_n\&Element;{{\mathrm{\&Theta;}}_s}^{{\mathrm{\&Phi;}}_{d_i}}}{\&Sigma;}F(\mathrm{\&Phi;},\mathrm{\&theta;},f)|w^T\&lsqb;g({\mathrm{\&Phi;}}_{d_i},f_k,{\mathrm{\&theta;}}_n)+\mathrm{\&Delta;}g({\mathrm{\&Phi;}}_{d_i},f_k,{\mathrm{\&theta;}}_n)\&rsqb;|<\mathrm{\&eta;}\\ w^{T}g({\mathrm{\&Phi;}}_{d_i},f_r,{\mathrm{\&Phi;}}_{d_i})\&GreaterEqual;1\end{array}---\left(4\right)$

上面代价函数(4)式中包含从上往下的3个式子：第一个式子表示的物理含义是在主瓣范围内保证频率不变；第二个式子表示在旁瓣范围内需要对实际响应值小于一个上限值，达到使旁瓣范围内的能量谱尽量小的目的；第三个式子的目的在于保证在指定波束指向方向参考频率f_r以及参考方向上的波束响应为1。

代价函数(4)式中：

g是导向量；Δg是导向矢量的误差；f是频率点；θ是波达方向；φ是波束指向角度；Ψ是预定义的可调方向范围；f_k、θ_n分别表示离散的可调方向角度点、频率点以及波达方向角度点；f_r是参考频率；Ω是感兴趣的频带范围；Ζ是预设的参数范围；Θ_p是主瓣角度范围；Θ_s是旁瓣角度范围；ε表示将导向矢量误差Δg的模值控制在ε值的界内；同理，η则表示其左侧的表达式的模值被认为控制在η值得界内；F是加权系数；w^T是权值的转置；

利用SeDumi软件对上式进行求解得出Farrow结构波束形成器的最优权值w。

基于Farrow结构波束形成器始终只需一组固定的权值，通过调整调向参数D这一参数，即可实现波束指向实时调整，无需重新求解权值，这正是我们设计的Farrow结构波束形成器用于SRP声源定位法能够简化计算量的关键原因所在。

所述步骤3中编号为j的麦克风阵元的输入信号经过Farrow结构波束形成器的空域滤波作用后输出信号y_j为：

$y_j=\underset{m=0}{\overset{M-1}{\&Sigma;}}{A}_m\&CenterDot;D_i^{M-1-m}---\left(6\right)$

(6)式中：

M是Farrow结构波束形成器中的滤波器的阶数；

并且有：

$\begin{array}{cc}{A}_m=\underset{k=0}{\overset{K-1}{\&Sigma;}}{x}_k\left(t\right)*w_{k,m}& m=0,1,\mathrm{...},M-1---\left(7\right)\end{array}$

(7)式中：

x_k(t)表示第k个麦克风阵元的输入信号，这里的n表示离散形式；w_k,m表示最优权值w中的第k个麦克风阵元的第m阶滤波器的权值；K表示麦克风阵元的个数；

根据Farrow波束形成器的结构，在这里A_m是只要通过一次计算之后就可以直接锁存起来的数据。即输出信号y_j只取决于调向参数D，即A_m可被直接锁存，从而达到实时调整波束指向的目的。

所述步骤2中，由Farrow结构波束形成器的SRP-PHAT方法的波束输出功率P(τ₀,...,τ_K-1)定义为：

(8)式中：

这里P(τ₀,...,τ_K-1)中的τ_k表示第k个阵元相对于参考阵元的时延差；为第k个麦克风接收到的输入信号x_k(t)的傅里叶变换；

为第l个麦克风接收到的输入信号x_l(t)的傅里叶变换；是的共轭转置；

且有：

D_i＝(θ_i-90°)/90°(11)

这里i表示可调角度离散形式的每个离散点；

多通道形式的PHAT的相位加权函数为：

(12)式表述的含义：将每个麦克风阵元输入信号进行FFT运算之后再两两进行相乘、取模、求倒数一系列运算之后得相位加权函数。

所述步骤4运算量降低的具体表现原理：

令：

(13)式中：

$A_{k,m_1}=\underset{n=0}{\overset{N-1}{\&Sigma;}}{x}_{k,m_1}\left(t\right)*w_{k,m_1,n}---\left(14\right)$

$A_{l,m_2}=\underset{n=0}{\overset{N-1}{\&Sigma;}}{x}_{l,m_2}\left(t\right)*w_{l,m_2,n}---\left(15\right)$

(14)式子表示的含义是：第k个麦克风阵元、第m₁阶FIR滤波器的输入信号经过空域滤波之后得到的输出信号将被锁存在的数据里。同理可解释(15)式中的

这样很容易理解(13)式所表达的含义是：将锁存起来的数据即输入信号与权值进行卷积再进行FFT后两两相乘，再与相应的加权函数进行相乘，IFFT运算之后累加求和，可得所有麦克风两两麦克风的广义互相关函数之和并进行锁存在H里。这一系列的过程随着权值的形成一步完成并锁存下来，无需重新计算。所述步骤4中利用Farrow结构波束形成器特殊的结构优势，降低计算量的原理过程。

所述步骤5中，按照Farrow结构波束形成器的结构图写出其瞬时功率：

$P({\mathrm{\&tau;}}_0,\mathrm{...},{\mathrm{\&tau;}}_{K-1})=H\&CenterDot;\underset{m_1=0}{\overset{M-1}{\&Sigma;}}\underset{m_2=0}{\overset{M-1}{\&Sigma;}}{D}_i^{2M-2-(m_1+m_2)}---\left(16\right)$

所述步骤5中按照Farrow结构波束形成器的结构图将调向参数D与已经锁存的输出数据结合起来，达到调向的目的。

所述步骤5中，传统的SRP-PHAT与本发明的Farrow-SRP-PHAT运算量对比：

假设传统的SRP-PHAT中的常规波束形成器有K个阵元、每个阵元有N个抽头、扫描点数为n′、信号序列长度为L点。

复乘次数：n′*[N*K*L+K²*[(3/2)*L*log₂L+3L]]次

复加次数：n′*[(N-1)*K*L+K²*[3*L*log₂L]+(K²-1)*L]次

假设本发明中的Farrow结构波束形成器有K个阵元、每个阵元有M阶N个抽头FIR滤波器、扫描点数为n′、信号序列长度为L点。

复乘次数：

L*M*N*K+n′*[K²*(M-1)²*L]+K²*M²*[(3/2)*L*log₂L+3*L]次

复加次数：

L*M*(N-1)*K+n′*[K²*(M-1)²*L]+K²*M²*(3*L*log₂L)次

所述步骤6中，进行校准工作的原因在于Farrow结构波束形成器本身是存在系统误差的。为了消除系统误差，从而将定位的精度提高，我们在算法的最后添加了一个校准工作。

有益效果：

本发明的一种SRP声源定位的快速实现方法，相比现有技术，本发明大大降低了运算量，简化了SRP-PHAT方法的实现复杂度，并且对混响和噪声具有更高的鲁棒性，方位估计精度也较高。具体为：

(1)SRP-PHAT方法运算量庞大，因此在实时处理变得困难，本发明利用Farrow结构波束形成器在调整波束指向时无需重新计算权值的结构优势，从而使得整体运算量大大降低，大大提高了实时性。

(2)SRP-PHAT方法在恶劣环境下(噪声干扰大、混响影响严重)性能急剧下降，所以使得实际应用变得困难，本发明利用了Farrow结构波束形成器鲁棒性较好的特性，使得本发明设计具有较高的抗混响能力和空间噪声抑制能力。

(3)本发明算法在解决了Farrow结构波束形成器本身具有系统误差这个问题的前提下，相比现有的方法提高了SRP-PHAT在恶劣环境下的定位精度。

附图说明

图1为Farrow结构波束形成器的结构图：

图2为本发明方法原理流程图；

图3为本发明基于凸优化束指向可调频率不变(FIB)波束指向可调指向60°的波束图的三维图；

图4为本发明基于凸优化束指向可调频率不变(FIB)波束指向可调指向60°的波束图的正视图；

图5为基于凸优化束指向可调最小二乘(LS)波束指向可调指向60°的波束图的三维图；

图6为基于凸优化束指向可调最小二乘(LS)波束指向可调指向60°的波束图的正视图；

图7为基于凸优化波束指向可调频率不变(FIB)算法在上述条件下设计的波束形成器指向60°时方向图的三维图；

图8为基于凸优化波束指向可调频率不变(FIB)算法在上述条件下设计的波束形成器指向60°时方向图的正视图；

图9为信噪比分别为为10dB、5dB，混响分别300ms到800ms以50ms为间隔的混响噪声环境下仿真的性能对比图；

图10为混响时间分别为250ms、300ms，信噪比分别-20dB到20dB以5dB为间隔的混响噪声环境下仿真的性能对比图；

图11表示信噪比10dB、0dB，混响分别为300ms到650ms以50ms为间隔的混响噪声环境下仿真的性能对比图；

图12为混响时间分别是600ms和750ms，信噪比分别为3dB到21dB以3dB为间隔的混响噪声环境下仿真的性能对比图；

图13为阵元数对本发明方法与SRP-PHAT算法在运算量的影响的对比图；

图14为扫描点数对本发明方法与SRP-PHAT算法在运算量的影响的对比图；

图15为抽头数对本发明方法与SRP-PHAT算法在运算量的影响对比图；

图16为信号序列的点数对本发明方法与SRP-PHAT算法在运算量的影响的对比图。

具体实施方式

下面结合附图对本发明作更进一步的说明。

本发明方法是在如图1所示的麦克风阵列的基础上，按照图2所示的原理流程图进行研究的。

步骤1、对Farrow结构波束形成器进行设计。所用到的算法是基于凸优化的波束指向可调的宽带稳健远场频率不变波束形成器算法来得到波束形成器的最优权值。

这里采用稳健的最大误差最小化的优化设计的代价函数：

$\underset{w}{\min }\underset{(\mathrm{\&Phi;},\mathrm{\&theta;},f)\&Element;Z}{\max }\underset{\left|\right|\mathrm{\&Delta;}\stackrel{\&OverBar;}{g}\left|\right|\&le;\mathrm{\&epsiv;}}{\max }\underset{{\mathrm{\&Phi;}}_{d_i}\&Element;\mathrm{\&Psi;}}{\&Sigma;}\underset{f_k\&Element;\mathrm{\&Omega;}}{\&Sigma;}\underset{{\mathrm{\&theta;}}_n\&Element;{{\mathrm{\&Theta;}}_p}^{{\mathrm{\&Phi;}}_{d_i}}}{\&Sigma;}F(\mathrm{\&Phi;},\mathrm{\&theta;},f)|w^T\&lsqb;g({\mathrm{\&Phi;}}_{d_i},f_k,{\mathrm{\&theta;}}_n)+\mathrm{\&Delta;}g({\mathrm{\&Phi;}}_{d_i},f_k,{\mathrm{\&theta;}}_n)\&rsqb;-w^{T}g({\mathrm{\&Phi;}}_{d_i},f_r,{\mathrm{\&theta;}}_n)|$

$s.t.\{\begin{array}{c}\underset{{\mathrm{\&Phi;}}_{d_i}\&Element;\mathrm{\&Psi;}}{\&Sigma;}\underset{f_k\&Element;\mathrm{\&Omega;}}{\&Sigma;}\underset{{\mathrm{\&theta;}}_n\&Element;{{\mathrm{\&Theta;}}_s}^{{\mathrm{\&Phi;}}_{d_i}}}{\&Sigma;}F(\mathrm{\&Phi;},\mathrm{\&theta;},f)|w^T\&lsqb;g({\mathrm{\&Phi;}}_{d_i},f_k,{\mathrm{\&theta;}}_n)+\mathrm{\&Delta;}g({\mathrm{\&Phi;}}_{d_i},f_k,{\mathrm{\&theta;}}_n)\&rsqb;|<\mathrm{\&eta;}\\ w^{T}g({\mathrm{\&Phi;}}_{d_i},f_r,{\mathrm{\&Phi;}}_{d_i})\&GreaterEqual;1\end{array}---\left(17\right)$

上式中，Δg(Φ,θ,f)表示由幅度和相位、位置等误差都将引起阵列的导向矢量g(Φ,θ,f)的畸变，假设g(Φ,θ,f)的畸变是在一定范围内，即：

按照SeDumi软件的求解规律，将(18)式化简为对偶形式再解代价函数，求得最优权值。

步骤2、使用本发明的设计方法设计好Farrow结构波束形成器之后，接下来我们用Farrow结构波束形成器在调整波束指向时无需重新计算权值这个结构特征，将输入信号进行FFT运算在两两进行相乘、取模、求倒数一系列运算之后得相位加权函数并进行锁存。

编号为j的麦克风阵元的输入信号经过Farrow结构波束形成器的空域滤波作用后输出：

$y_j=\underset{m=0}{\overset{M-1}{\&Sigma;}}{A}_m\&CenterDot;D_i^{M-1-m}---\left(19\right)$

(19)式中：

M是Farrow结构波束形成器中的滤波器的阶数；其中：

$\begin{array}{cc}{A}_m=\underset{k=0}{\overset{K-1}{\&Sigma;}}{x}_k\left(t\right)*w_{k,m}& m=0,1,\mathrm{...},M-1---\left(20\right)\end{array}$

其中x_k(t)是第k个阵元的输入信号。A_m是输入信号与权值卷积累加求和作用的结果。根据Farrow波束形成器的结构，在这里A_m是只要通过一次计算之后就可以直接锁存起来的数据。

设扫描角度是θ_i，这里的i表示离散点，在这里我们假设可调范围是[40°,140°]，扫描间隔为2°，共有离散点个数为(140-40)/2+1＝51个；

D_i＝(θ_i-90°)/90°i＝1,2,...,51(21)

已知Farrow-SRP-PHAT方法的波束输出功率定义为：

多通道形式的PHAT的相位加权函数为：

对输入每个麦克风的每一阶FIR滤波器的原输入信号分别做一次傅里叶变换，再将经过变换之后的两两信号经过相乘、取模、求倒数等一系列运算之后可得相位加权函数；

其中为第k个麦克风接收到的信号x_k(t)的傅里叶变换；为第l个麦克风接收到的信号x_l(t)的傅里叶变换；

步骤3、将各麦克风阵列的输入信号输入至Farrow结构波束形成器中进行空域滤波，空域滤波后的输出信号进行FFT运算后再两两相乘求得与式(23)中已获得的相位加权函数下标相对应的量；

令

它们分别为延时补偿后的信号做FFT变换之后的输出信号；

步骤4、将步骤2中获得的相位加权函数与步骤3中获得的与式(25)中已获得的相位加权函数下标相对应的量对应进行相乘，然后将相乘的结果进行IFFT运算之后并累加求和得到所有的麦克风中的两两麦克风的广义互相关函数之和并进行锁存；

结合(22)，(23)，(24)，(25)，(26)，(27)，整理(22)式，得：

步骤5、将两两麦克风的广义互相关函数之和与调向参数D做乘积并累加乘积相加运算，得到相应扫描点的瞬时功率，搜索瞬时功率的最大峰值得到对应的指向角度θ_r'；

令x_k,m(t)为输入第k个麦克风，第m阶FIR滤波器的的输入信号经过频率不变波束指向可调的波束形成器空域滤波作用之后的输出信号，x_l,m(t)为输入第l个麦克风、第m阶FIR滤波器的的输入信号同样经过空域滤波之后的输出信号，则：

$\begin{array}{cc}{\stackrel{\&OverBar;}{x}}_{k,m}\left(t\right)=\underset{n=0}{\overset{N-1}{\&Sigma;}}{x}_k\left(t\right)*w_{k,m,n}\&CenterDot;D_i^{M-1-m}& k=0,1,\mathrm{...},K-1\end{array}---\left(29\right)$

$\begin{array}{cc}{\stackrel{\&OverBar;}{x}}_{l,m}\left(t\right)=\underset{n=0}{\overset{N-1}{\&Sigma;}}{x}_l\left(t\right)*w_{1,m,n}\&CenterDot;D_i^{M-1-m}& l=0,1,\mathrm{...},K-1---\left(30\right)\end{array}$

D_i＝(θ_i-90°)/90°(31)

结合(26)(27)和(29)(30)式可知：

我们知道，调向参数D其实是和指向角度有关的一个常数，根据FFT的性质，式(32)、(33)可以写成

令:

$A_{k,m}=\underset{n=0}{\overset{N-1}{\&Sigma;}}{x}_{k,m}\left(t\right)*w_{k,m,n}---\left(36\right)$

$A_{l,m}=\underset{n=0}{\overset{N-1}{\&Sigma;}}{x}_{l,m}\left(t\right)*w_{l,m,n}---\left(37\right)$

所以，(34)、(35)可以改写成:

(38)的式子表示的含义是即第k个麦克风，第m₁阶输入信号经过空域滤波之后的输出信号经过FFT变换可以等价为输入信号与第k个麦克风第m₁阶滤波器的N个抽头系数卷积之后做FFT变换再与调向参数D做乘积运算。这个过程只需要做一次运算即可锁存，可重复使用。同理解释(39)式。

整理(28)，将(38)(39)带入(28)中：

令

式(28)可化简为：

$P({\mathrm{\&tau;}}_0,\mathrm{...},{\mathrm{\&tau;}}_{K-1})=H\&CenterDot;\underset{m_1=0}{\overset{M-1}{\&Sigma;}}\underset{m_2=0}{\overset{M-1}{\&Sigma;}}{D}_i^{2M-2-(m_1+m_2)}---\left(42\right)$

H所表达的含义是将锁存起来的数据即输入信号与权值进行卷积再进行FFT后两两相乘，再与相应的加权函数进行相乘，IFFT运算之后累加求和可得所有麦克风两两麦克风的广义互相关函数之和并进行锁存在H里。这一系列的过程随着权值的形成一步完成并锁存下来，无需重新计算。P(τ₀,...,τ_K-1)实际就是调向参数D根据扫描角度的不同，与H数据进行乘积累加求和运算。从而大大降低了运算量。

步骤6、对Farrow结构波束形成器的系统误差进行校准，并将步骤5中获得的指向角度θ_r'校准为最终的定位角度θ_r。

下面结合一些具体实施例以及对比例对本发明做进一步说明。

实施例1

本发明方法波束形成器的设计方法是基于凸优化的波束指向可调的宽带稳健远场频率不变波束形成器算法。

现有的波束形成器设计技术其中包括基于凸优化的最小二乘波束指向可调的宽带稳健远场波束形成器算法，本发明的算法在处理语音信号的频段分布较宽时，性能优势较为明显。为此，举例说明相同环境下的两种算法的性能对比。考虑Farrow结构波束形成器的阵型为：阵元数为10、阶数为5、抽头数为20，信号频段分布在[800Hz,3600Hz],阵元间距为5cm，可调方向范围为40°到140°，参考频率为8000Hz，过渡带为40°。图3和图4为本发明基于凸优化束指向可调频率不变(FIB)波束指向可调指向60°的波束图的三维图和正视图；图5和图6为基于凸优化束指向可调最小二乘(LS)波束指向可调指向60°的波束图的三维图和正视图。通过以上仿真分析表明，本发明方法在信号频率分布在[800Hz,3600Hz]的环境要求下，指向性能优势明显。

实施例2

本发明方法说明了Farrow结构波束形成器本身是存在系统误差的。校准工作在本发明中起着很大的作用。这个步骤的根本原因在于Farrow结构波束形成器本身是具有系统误差的。我们通过具体实例证明这个理论。

假设存在均匀分布的麦克风线列，维度是10×5×20，即有10个麦克风，每个麦克风接5阶FIR滤波器，每个FIR的抽头为20。相邻的麦克风间隔5cm，并且采样频率为f_s＝8000Hz。波束可调方向范围Ψ设为感兴趣的频率范围Ω为[1000-3000]Hz，参考频率f_r＝2600。在每个波束指向中，通带带宽(PW)设为40°。在我们定义角度范围内，每个波束指向下的左右两个阻带分别[0°,φ_d-PW/2-20°]，[φ_d+PW/2+20°,180°]，其中φ_d是指向角度，白噪声增益为0.001。我们按照上述条件来设计我们的波束形成器。

图7和图8分别给出了基于凸优化波束指向可调频率不变(FIB)算法在上述条件下设计的波束形成器指向60°时方向图的三维图和正视图。由正视图可以看出，方向图的最高点并不是指向60°，实际指向为66.5306°，实际值与理论值相差了6.5306°。这就是我们所说的系统误差。

实施例3

本发明方法说明校准工作对Farrow-SRP-PHAT算法性能的影响。校准工作是为了消除Farrow结构波束形成器系统误差。这个步骤在本发明方法中的起着很大的作用。

说话人位于坐标(2.6，2.8)m位置，麦克风阵列中心点位于(3，0)m；

考虑波束可调方向范围从40°到140°，以2°为间隔，进行10次蒙特卡洛实验。对这10次结果求取均值和方差，利用均方根误差作为评价指标。

图9为信噪比分别为为10dB、5dB，混响分别300ms到800ms以50ms为间隔的混响噪声环境下仿真的性能对比图；图10为混响时间分别为250ms、300ms，信噪比分别-20dB到20dB以5dB为间隔的混响噪声环境下仿真的性能对比图。通过以上仿真分析表明，本发明方法校准之后，定位精度有显著提高；并且抗混响性能以及抗噪声性能也有明显的提升。

实施例4

SRP-PHAT算法以及本发明方法仿真以及性能比较。

说话人位于坐标(2.6，2.8)m位置，麦克风阵列中心点位于(3，0)m；

考虑波束可调方向范围从40°到140°，以2°为间隔，进行10次蒙特卡洛实验。对这10次结果求取均值和方差，利用均方根误差作为评价指标。

图11表示信噪比10dB、0dB，混响分别为300ms到650ms以50ms为间隔的混响噪声环境下仿真的性能对比图；图12为混响时间分别是600ms和750ms，信噪比分别为3dB到21dB以3dB为间隔的混响噪声环境下仿真的性能对比图。从仿真结果来看，本发明方法的定位精度要高于SRP-PHAT方法，且抗噪声以及抗混响性能也是优于SRP-PHAT方法的。

实施例5

通过对各参数分别对SRP-PHAT以及本发明算法运算量的影响，比较本发明算法和SRP-PHAT算法运算量的对比。图13为阵元数对二者运算量的影响，图14为扫描点数对二者运算量的影响，图15为抽头数对二者运算量的影响，图16为信号序列的点数对二者运算量的影响。通过各参数对二者运算量的影响可以看出：Farrow-SRP-PHAT算法的的运算量远远小于SRP-PHAT算法。

以上所述仅是本发明的优选实施方式，应当指出：对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。

Claims (3)

1.一种SRP声源定位的快速实现方法，其特征在于：包括如下步骤：

步骤1、利用基于凸优化的波束指向可调的宽带稳健远场频率不变波束形成器算法得到Farrow结构波束形成器的最优权值，从而完成Farrow结构波束形成器的设计；

步骤2、将各麦克风阵列的输入信号进行FFT运算，根据不同麦克风以及输入信号的阶数将FFT运算后的输入信号进行两两相乘，求得相位加权函数；

步骤3、将各麦克风阵列的输入信号输入至Farrow结构波束形成器中进行空域滤波，空域滤波后的输出信号进行FFT运算，根据不同麦克风以及输出信号的阶数将FFT运算后的输出信号进行两两相乘求得相位加权函数对应量；

步骤4、根据两个麦克风以及各自输入信号的阶数所确定的相位加权函数以及相位加权函数对应量为一组进行对应，将属于同一组的相位加权函数与相位加权函数对应量进行相乘，然后将相乘的结果进行IFFT运算之后并累加求和得到所有的麦克风中的两两麦克风的广义互相关函数之和并进行锁存；

步骤5、将两两麦克风的广义互相关函数之和与调向参数D做乘积并累加运算，得到相应扫描点的瞬时功率，搜索瞬时功率的最大峰值得到对应的指向角度θ_r'；

步骤6、对Farrow结构波束形成器的系统误差进行校准，并将步骤5中获得的指向角度θ_r'校准为最终的定位角度θ_r。

2.如权利要求1所述的SRP声源定位的快速实现方法，其特征在于：步骤1的设计过程中，采用如下稳健的最大误差最小化的优化设计的代价函数：

$\begin{array}{c}\underset{w}{\min }\underset{(\mathrm{\&Phi;},\mathrm{\&theta;},f)\&Element;Z}{\max }\underset{\left|\right|\mathrm{\&Delta;}\stackrel{\&OverBar;}{g}\left|\right|\&le;\mathrm{\&epsiv;}}{\max }\underset{{\mathrm{\&Phi;}}_{d_i}\&Element;\mathrm{\&Psi;}}{\mathrm{\&Sigma;}}\underset{f_k\&Element;\mathrm{\&Omega;}}{\mathrm{\&Sigma;}}\underset{{\mathrm{\&theta;}}_n\&Element;{{\mathrm{\&Theta;}}_p}^{{\mathrm{\&Phi;}}_{d_i}}}{\mathrm{\&Sigma;}}F(\mathrm{\&Phi;},\mathrm{\&theta;},f)|w^T\&lsqb;g({\mathrm{\&Phi;}}_{d_i},f_k,{\mathrm{\&theta;}}_n)+\mathrm{\&Delta;}g({\mathrm{\&Phi;}}_{d_i},f_k,{\mathrm{\&theta;}}_n)\&rsqb;-w^{T}g({\mathrm{\&Phi;}}_{d_i},f_k,{\mathrm{\&theta;}}_n)|\\ s.t.\left\{\begin{array}{c}\underset{{\mathrm{\&Phi;}}_{d_i}\&Element;\mathrm{\&Psi;}}{\mathrm{\&Sigma;}}\underset{f_k\&Element;\mathrm{\&Omega;}}{\mathrm{\&Sigma;}}\underset{{\mathrm{\&theta;}}_n\&Element;{{\mathrm{\&Theta;}}_p}^{{\mathrm{\&Phi;}}_{d_i}}}{\mathrm{\&Sigma;}}F(\mathrm{\&Phi;},f,\mathrm{\&theta;})|w^T\&lsqb;g({\mathrm{\&Phi;}}_{d_i},f_k,{\mathrm{\&theta;}}_n)+\mathrm{\&Delta;}g({\mathrm{\&Phi;}}_{d_i},f_k,{\mathrm{\&theta;}}_n)\&rsqb;|<\mathrm{\&eta;}\\ w^{T}g({\mathrm{\&Phi;}}_{d_i},f_r,{\mathrm{\&Phi;}}_{d_i})\&le;1\end{array}\right.\end{array}---\left(1\right)$

(1)式中：

g是导向量；Δg是导向矢量的误差；f是频率点；θ是波达方向角度点；φ是波束指向角度；Ψ是预定义的可调方向范围；f_k、θ_n分别表示离散的可调方向角度点、频率点以及波达方向角度点；f_r是参考频率；Ω是感兴趣的频带范围；Ζ是预设的参数范围；Θ_p是主瓣角度范围；Θ_s是旁瓣角度范围；ε表示将导向矢量误差Δg的模值控制在ε值的界内；同理，η则表示其左侧的表达式的模值被认为控制在η值得界内；F是加权系数；w^T是权值的转置；

利用SeDumi软件对上式进行求解得出Farrow结构波束形成器的最优权值w。

3.如权利要求1所述的SRP声源定位的快速实现方法，其特征在于：所述步骤3步骤中，编号为j的麦克风阵元的输入信号经过Farrow结构波束形成器进行空域滤波后输出信号为：

$y_j=\underset{m=0}{\overset{M-1}{\&Sigma;}}{A}_m\&CenterDot;D_i^{M-1-m}---\left(2\right)$

(2)式中：

M是Farrow结构波束形成器中的滤波器的阶数；

并且有：

$A_m=\underset{k=0}{\overset{K-1}{\&Sigma;}}{x}_k\left(t\right)*w_{k,m},m=0,1,...,M-1---\left(3\right)$

(3)式中：

x_k(t)表示第k个麦克风阵元的输入信号；w_k,m表示最优权值中的第k个麦克风阵元的第m阶滤波器的权值；K表示麦克风阵元的个数。