CN108549052B - 一种时频-空域联合加权的圆谐域伪声强声源定位方法 - Google Patents

一种时频-空域联合加权的圆谐域伪声强声源定位方法 Download PDF

Info

Publication number
CN108549052B
CN108549052B CN201810228816.8A CN201810228816A CN108549052B CN 108549052 B CN108549052 B CN 108549052B CN 201810228816 A CN201810228816 A CN 201810228816A CN 108549052 B CN108549052 B CN 108549052B
Authority
CN
China
Prior art keywords
time
weighted
domain
phat
frequency domain
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN201810228816.8A
Other languages
English (en)
Other versions
CN108549052A (zh
Inventor
宋宫琨琨
陈华伟
吕中兴
张普芬
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Nanjing University of Aeronautics and Astronautics
Original Assignee
Nanjing University of Aeronautics and Astronautics
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Nanjing University of Aeronautics and Astronautics filed Critical Nanjing University of Aeronautics and Astronautics
Priority to CN201810228816.8A priority Critical patent/CN108549052B/zh
Publication of CN108549052A publication Critical patent/CN108549052A/zh
Application granted granted Critical
Publication of CN108549052B publication Critical patent/CN108549052B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S5/00—Position-fixing by co-ordinating two or more direction or position line determinations; Position-fixing by co-ordinating two or more distance determinations
    • G01S5/18—Position-fixing by co-ordinating two or more direction or position line determinations; Position-fixing by co-ordinating two or more distance determinations using ultrasonic, sonic, or infrasonic waves
    • G01S5/22—Position of source determined by co-ordinating a plurality of position lines defined by path-difference measurements

Landscapes

  • Physics & Mathematics (AREA)
  • Engineering & Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Remote Sensing (AREA)
  • Circuit For Audible Band Transducer (AREA)

Abstract

本发明公开了一种时频‑空域联合加权的圆谐域伪声强声源定位方法，在该方法中，用先设计的六元圆形麦克风阵列采集语音信号；考虑到语音信号具有短时平稳特性，将采集到的语音信号变换到时频域，利用时频域PHAT加权进行滤波处理，降低噪声和混响对定位性能的影响；利用归一化后零阶特征波束、一阶特征波束，并由时频域PHAT加权的圆谐伪声强声源定位方法进行粗估计角度；利用粗估计角度构造空域波束指向性函数，得到归一化后时频域PHAT加权联合指向性加权零阶特征波束，再利用时频‑空域联合加权的圆谐域伪声强方法求解得到精确声源估计角度。

Description

一种时频-空域联合加权的圆谐域伪声强声源定位方法
技术领域
本发明涉及声源定位技术领域，特别是一种时频-空域联合加权的圆谐域伪声强声源定位方法。
背景技术
在音频与语音信号处理中，采用麦克风阵列的声源定位技术是中的一个重要研究方向，线性麦克风阵列因其简单且易于理解、实现已经在声源定位中得到广泛的应用，如声呐(见文献：王燕,邹男,梁国龙.强多途环境下水听器阵列位置近场有源校正方法[J].物理学报,2015,64(2):024304 1-10)、视频电话会议(见文献：Barbara Rauch,FriedrichFaubel,Dietrich Klakow.An analysis of nonstationary variance estimates in themaximum negentropy beamformer[C].Joint Workshop on Hands-free SpeechCommunication and Microphone Arrays,2011:201-206)、人工智能(见文献：梁瑞宇，周健，王青云，奚吉，赵力.仿人耳听觉的助听器双耳声源定位算法.声学学报，2015；40(3):446-454)、地震研究(见文献：吴晓平，顾治华，舒红波，冯海林.一种线性最小二乘法的声源目标精确定位方法.声学学报，2016；41(1):87-93)、声源定位与追踪(见文献：DespoinaPavlidi,Anthony Griffin,Matthieu Puigt,Athanasios Mouchtaris.Real-TimeMultiple Sound Source Localization and Counting Using a Circular MicrophoneArray.IEEE Transactions on Audio,Speech,and Language Processing,2013；21(10):2193-2206)、监控系统(见文献：林志斌，徐柏龄.基于传声器阵列的声源定位.电声技术，2004；28(5):19-23)等。
目前常用的声源定位的方法主要有三类：基于可控波束形成器的声源定位法，主要是将各阵元采集到的信号加权求和，通过调控权值(权值取决于阵元信号的相位延迟，主要利用波达方向等方法进行调控)，使阵列输出信号功率最大，从而进行声源定位；基于高分辨率谱估计的声源定位法，主要运用最小方差估计法、子空间法，如MUSIC法、ESPRIT法等进行声源定位；基于到达时间差(TDOA)的声源定位法，主要结合自适应、广义互相关等进行声源定位。(见文献：居太亮.基于麦克风阵列声源定位算法研究[D].博士学位论文(成都：电子科技大学),2006)。
现有定位技术中，主要是利用一维线性阵列、二维十字阵列进行声源定位，与圆形麦克风阵列相比，这些阵列由于自身结构所限制，在二维空间中的定位角度范围只能在0～180°，然而圆形麦克风阵列因其自身优势其定位角度范围在0～360°，能进行全方位声场分析。另外，在噪声与混响同时干扰下，现有声源定位方法的性能的提升往往以增加麦克风个数或阵列尺寸为代价，但是定位精确度提高甚微。
如何解决现有技术的不足已成为声源定位领域亟待解决的一大难题。
发明内容
本发明所要解决的技术问题是克服现有技术的不足，而提供一种时频-空域联合加权的圆谐域伪声强声源定位方法，本发明方法实时、有效地解决声源定位问题，在一定程度上降低了噪声、混响等对语音信号的干扰，提高了定位的精确性和鲁棒性。
本发明为解决上述技术问题采用以下技术方案：
根据本发明提出的一种适用于圆形麦克风阵列的时频-空域联合加权的圆谐域伪声强声源定位方法，具体如下：
构建六元圆形麦克风阵列，采用构建的六元圆形麦克风阵列采集语音信号；
将采集到的语音信号变换到时频域，利用时频域PHAT加权进行滤波处理；
利用归一化后零阶特征波束、一阶特征波束，并由时频域PHAT加权的圆谐伪声强声源定位方法进行粗估计角度；
最后用粗估计角度构造空域波束指向性函数，得到归一化后时频域PHAT加权联合指向性加权零阶特征波束，并利用时频-空域联合加权的圆谐域伪声强方法求解得到精确声源估计角度。
作为本发明所述的一种适用于圆形麦克风阵列的时频-空域联合加权的圆谐域伪声强声源定位方法进一步优化方案，包括如下步骤：
步骤一：采用Q个相同的全向性麦克风等间距的排列成半径为r的圆形麦克风阵列；
步骤二：对麦克风阵列采集到的声压信号Pq(t,ω)进行短时傅里叶变换，然后利用时频域PHAT加权对短时傅里叶变换后的声压信号进行预滤波，得到时频域PHAT加权后的n阶特征波束F′n(t,ω)；具体如下：
(201)、对声源发出的信号s(t)进行采样，得到第q个麦克风采集到的声压信号
Figure BDA0001602091820000021
其中，hq(t)表示声源到第q个麦克风间的脉冲响应，nq(t)表示第q个麦克风接收到的加性噪声信号，符号
Figure BDA0001602091820000022
表示线性卷积运算，q＝1,…,Q；
(202)、将(201)中采集到的声压信号pq(t)变换到时频域处理，经短时傅里叶变换得
Pq(t,ω)＝Hq(t,ω)·S(t,ω)+Nq(t,ω)
其中，Pq(t,ω)、S(t,ω)、Hq(t,ω)以及Nq(t,ω)分别表示麦克风接收声压信号pq(t)、声源信号s(t)、脉冲响应hq(t)和加性噪声信号nq(t)的短时傅里叶变换，(t,ω)表示时频率单元，t表示时间，ω表示频率；
(203)、利用PHAT加权对(202)中经短时傅里叶变换后的声压信号Pq(t,ω)进行预滤波，得到经过时频域PHAT加权后的n阶特征波束
Figure BDA0001602091820000031
其中，n表示阶数，
Figure BDA0001602091820000032
表示第1个麦克风与第q个麦克风之间按逆时针方向的夹角，
Figure BDA0001602091820000033
表示虚数单位，e表示自然指数；
步骤三：对步骤二中的n阶特征波束F′n(t,ω)进行归一化处理，得到归一化后时频域PHAT加权的零阶特征波束D′0(t,ω)、加权的一阶特征波束的两个正交分量D′x(t,ω)和D′y(t,ω)，然后用时频域PHAT加权的圆谐伪声强声源定位方法求解得到粗估计角度
Figure BDA0001602091820000034
具体如下：
(301)、对步骤二中的n阶特征波束F′n(t,ω)进行归一化处理，取阶数n＝0时得到归一化后时频域PHAT加权的零阶特征波束
Figure BDA0001602091820000035
其中，b0(t,ω)表示不同时频点的0阶贝塞尔函数；
取阶数n＝1时特征波束x轴和y轴的两个正交分量得到归一化后时频域PHAT加权的一阶特征波束的两个正交分量D′x(t,ω)和Dy′(t,ω)，
Figure BDA0001602091820000036
Figure BDA0001602091820000037
其中，b1(t,ω)表示不同时频点的1阶贝塞尔函数，γx,1＝ei·1·0表示一阶特征波束的x轴旋转系数，
Figure BDA0001602091820000038
表示一阶特征波束的y轴旋转系数，F′x(t,ω)表示1阶特征波束x轴分量，F′y(t,ω)表示1阶特征波束y轴分量；
(302)、将(301)中得到加权的零阶特征波束D′0(t,ω)、加权的一阶特征波束的两个正交分量D′x(t,ω)和D′y(t,ω)，采用时频域PHAT加权的圆谐伪声强声源定位方法求解得到粗估计角度
Figure BDA0001602091820000041
其中，tα和ωβ表示第α个时间单元和第β个频率单元，*表示取共轭运算，Re表示取实部运算，arctan表示反正切运算，I′x(tα,ωβ)和I′y(tα,ωβ)分别表示时频域PHAT加权后各时频点处的瞬时伪声强的x轴分量和y轴分量；
步骤四：利用步骤三中求得的粗估计角度
Figure BDA0001602091820000042
构造空域波束指向性函数
Figure BDA0001602091820000043
结合步骤三中求得的归一化后时频域PHAT加权零阶特征波束，得到归一化后时频域PHAT加权联合指向性加权零阶特征波束
Figure BDA0001602091820000044
最后利用时频-空域联合加权的圆谐域伪声强方法求解得到精确声源估计角度
Figure BDA0001602091820000045
具体如下：
(401)、利用步骤三中求得的粗估计角度
Figure BDA0001602091820000046
构造空域波束指向性函数
Figure BDA0001602091820000047
其中，n表示阶数，且n的最高阶数为N＝kr，
Figure BDA0001602091820000048
为波数，f为频率，c为声波传播的速度，bn(t,ω)表示不同时频点的n阶贝塞尔函数；
(402)、将(401)中的空域波束指向性函数的幅度
Figure BDA0001602091820000049
加权到步骤三中得到的归一化后时频域PHAT加权零阶特征波束上，得到归一化后时频域PHAT加权联合指向性加权零阶特征波束
Figure BDA00016020918200000410
其中，ψ(t,ω)表示
Figure BDA00016020918200000411
的相位；
(403)、利用时频域-空域联合加权的圆谐伪声强方法求解得到精确声源估计角度
Figure BDA0001602091820000051
作为本发明所述的一种适用于圆形麦克风阵列的时频-空域联合加权的圆谐域伪声强声源定位方法进一步优化方案，所述(302)中时频域PHAT加权的圆谐伪声强声源定位方法是采用有功伪声强进行瞬时方位估计，对各时频点伪声强进行平均相互补偿。
作为本发明所述的一种适用于圆形麦克风阵列的时频-空域联合加权的圆谐域伪声强声源定位方法进一步优化方案，所述(403)中时频域-空域联合加权的圆谐伪声强方是采用有功伪声强进行瞬时方位估计，对各时频点伪声强进行平均相互补偿。
作为本发明所述的一种适用于圆形麦克风阵列的时频-空域联合加权的圆谐域伪声强声源定位方法进一步优化方案，所述(402)中只保留步骤三中得到的归一化后时频域PHAT加权零阶特征波束D′0(t,ω)的相位ψ(t,ω)，幅度值替换为
Figure BDA0001602091820000052
得到归一化后时频域PHAT加权联合指向性加权零阶特征波束
Figure BDA0001602091820000053
本发明采用以上技术方案与现有技术相比，具有以下技术效果：
(1)本发明构建了六元圆形麦克风阵列，运用一种适用于圆形麦克风阵列的时频-空域联合加权的圆谐域伪声强声源定位方法，对室内远场单声源进行定位；首先用设计的六元圆形麦克风阵列采集语音信号；然后考虑到语音信号具有短时平稳特性，将采集到的语音信号变换到时频域，利用时频域PHAT加权进行滤波处理，降低噪声和混响对定位性能的影响；其次利用归一化后零阶特征波束、一阶特征波束，并由时频域PHAT加权的圆谐伪声强声源定位方法进行粗估计角度；最后用粗估计角度构造空域波束指向性函数，得到归一化后时频域PHAT加权联合指向性加权零阶特征波束，并利用时频-空域联合加权的圆谐域伪声强方法求解得到精确声源估计角度；
(2)本发明方法降低了噪声、混响等对语音信号的干扰，提高了定位的精确性和鲁棒性；
(3)本发明在仿真与实测实验中都能准确地确定声源方位，定位精度高与稳定性强，在语音信号处理领域，具有较强的实用性。
附图说明
图1是本发明设计的六元麦克风房间仿真模型。
图2是本发明在信噪比10dB、不同混响时间下，基本伪声强法、时频域PHAT加权法、谐波域PHAT加权法平均均方根误差的对比。
图3是本发明在混响时间300ms、不同信噪比下，基本伪声强法、时频域PHAT加权法、谐波域PHAT加权法平均均方根误差的对比。
图4是本发明在信噪比10dB、不同混响时间下，圆谐可控相应功率法、基本伪声强法、时频域PHAT加权法、时频-空域联合加权法的平均均方根误差比较对比。
图5是本发明在混响时间300ms、不同信噪比下，圆谐可控相应功率法、基本伪声强法、时频域PHAT加权法、时频-空域联合加权法的平均均方根误差比较对比。
图6是本发明实测实验进行声源定位结果的对比图。
图7是本发明的流程图。
具体实施方式
下面结合附图对本发明的技术方案做进一步的详细说明：
本发明是一种适用于圆形麦克风阵列的时频-空域联合加权的圆谐域伪声强声源定位方法，利用六元麦克风阵列，结合语音信号特性进行声源定位，图7所示是本发明的流程图，其具体实施步骤如下：
步骤一：采用Q个相同的全向性麦克风等间距的排列成半径为r的圆形麦克风阵列；
建立圆形麦克风阵列模型，如图1所示，由Q个全向性麦克风M1,...,MQ组成，选择阵列中心作为坐标原点O，麦克风按逆时针方向排列等间隔地分布在半径为r的圆周上，
Figure BDA0001602091820000061
表示第1个麦克风与第q个麦克风之间按逆时针方向的夹角。对于单声源远场情况，设声源s(t)入射方向与x轴正方向的夹角为φs∈[0°,360°)；
步骤二：对麦克风阵列采集到的声压信号Pq(t,ω)进行短时傅里叶变换，然后利用时频域PHAT(The Phase Transform，相位变换)加权对短时傅里叶变换后的声压信号进行预滤波，得到时频域PHAT加权后的n阶特征波束F′n(t,ω)；
(201)、对声源发出的信号s(t)进行采样，得到第q个麦克风采集到的声压信号
Figure BDA0001602091820000062
式中，hq(t)表示声源到第q个麦克风间的脉冲响应，nq(t)表示第q个麦克风接收到的加性噪声信号，符号
Figure BDA0001602091820000063
表示线性卷积运算；
在频域中，式(1)中麦克风接收到的声压信号可表示为
Pq(k)＝Hq(k)S(k)+Nq(k) (2)
式中，
Figure BDA0001602091820000071
表示波数，f表示频率，c≈340m/s表示声速，S(k)、Hq(k)以及Nq(k)分别表示声源信号s(t)、脉冲响应hq(t)和加性噪声信号nq(t)的傅里叶变换。
根据圆谐傅里叶变换可知，当声源的入射方位角为φs，则声源的声压信号P(kr,φs)可以由(201)中的声压信号
Figure BDA0001602091820000072
的n阶特征波束Fn(kr)表示，
Figure BDA0001602091820000073
Figure BDA0001602091820000074
其中，n表示阶数，
Figure BDA0001602091820000075
表示虚数单位，e表示自然指数，
Figure BDA0001602091820000076
表示n阶圆谐波，a(kr)表示声波幅度，bn(kr)表示模态系数。根据圆形孔径是否遮挡，模态系数bn(kr)有两种选取方式，即
Figure BDA0001602091820000077
式中，Jn(kr)表示n阶贝塞尔函数，Hn(kr)表示n阶汉克尔函数，J′n(kr)表示n阶贝塞尔函数的一阶导数，H′n(kr)表示n阶汉克尔函数的一阶导数。本发明中所考虑的圆形麦克风阵列为无遮挡的开放型圆阵。
由于圆形阵列实质上是对圆形孔径的空间采样，则根据式(4)可导出圆形阵列的n阶特征波束Fn(kr)近似表达式，即
Figure BDA0001602091820000078
在式(6)的近似过程中，有两个误差需要考虑，一是阶数的截断误差，二是阵列对孔径的采样误差。理想情况下，频域信号可以展开为无穷多个不相关的圆谐波，但是实际应用中，谐波数量必须截断，取一个最高阶数N＝kr。当N为最高阶数时，频域信号可以展开成2N+1个不相关的圆谐波。所以为了减少信息丢失，阵列的传声器数量必须满足条件：Q≥2N+1。
(202)、考虑到语音信号具有短时平稳特性，因此将(201)中采集到的声压信号pq(t)变换到时频域处理，经短时傅里叶变换得
Pq(t,ω)＝Hq(t,ω)·S(t,ω)+Nq(t,ω) (7)
式中，Pq(t,ω)、S(t,ω)、Hq(t,ω)以及Nq(t,ω)分别表示麦克风接收声压信号pq(t)、声源信号s(t)、脉冲响应hq(t)和加性噪声nq(t)的短时傅里叶变换，(t,ω)表示时频率单元，t表示时间，ω表示频率；
将(7)式代入(6)式可得时频域的n阶特征波束Fn(t,ω)为
Figure BDA0001602091820000081
(203)、利用PHAT加权对(202)中经短时傅里叶变换后的声压信号Pq(t,ω)进行预滤波，减少多径信道的畸变，从而降低混响对定位结果的影响，得到经过时频域PHAT加权后的n阶特征波束
Figure BDA0001602091820000082
步骤三：对步骤二中的n阶特征波束F′n(t,ω)进行归一化处理，得到归一化后时频域PHAT加权的零阶特征波束D′0(t,ω)、加权的一阶特征波束的两个正交分量D′x(t,ω)和D′y(t,ω)，然后用时频域PHAT加权的圆谐伪声强声源定位方法求解得到粗估计角度
Figure BDA0001602091820000083
(301)、对步骤二中的n阶特征波束F′n(t,ω)进行归一化处理，取阶数n＝0时得到归一化后时频域PHAT加权的零阶特征波束
Figure BDA0001602091820000084
式中，b0(t,ω)表示不同时频点的0阶贝塞尔函数；
取阶数n＝1时特征波束x轴和y轴的两个正交分量得到归一化后时频域PHAT加权的一阶特征波束的两个正交分量D′x(t,ω)和D′y(t,ω)，
Figure BDA0001602091820000085
Figure BDA0001602091820000086
式中，b1(t,ω)表示不同时频点的1阶贝塞尔函数，γx,1＝ei·1·0表示一阶特征波束的x轴旋转系数，
Figure BDA0001602091820000091
表示一阶特征波束的y轴旋转系数，F′x(t,ω)表示1阶特征波束x轴分量，F′y(t,ω)表示1阶特征波束y轴分量；
(302)、各时频点处的瞬时伪声强x轴分量Ix(t,ω)和y轴分量Iy(t,ω)可分别表示为
Figure BDA0001602091820000092
Figure BDA0001602091820000093
式中，D0(t,ω)表示归一化后的零阶特征波束，Dx(t,ω)和Dy(t,ω)表示一阶特征波束归一化后的x轴和y轴分量，*表示取共轭运算，Re表示取实部运算；
故由式(13)、(14)可得计算各时频点的瞬时方位
Figure BDA0001602091820000094
估计公式为
Figure BDA0001602091820000095
式中，arctan表示反正切运算。
从理论上讲，仅利用一阶特征波束的信息就可以估计出瞬时方位。但采用有功伪声强进行方位估计则可以提高抗噪声能力，因此这里采用有功伪声强进行瞬时方位估计。而且由于各时频点的瞬时方位估计结果差异性较大，对噪声鲁棒性不够，故对时频点伪声强进行平均，使得各时频点间相互补偿。可得计算声源的
Figure BDA0001602091820000096
估计公式为
Figure BDA0001602091820000097
式中，tα和ωβ表示第α个时间单元和第β个频率单元
根据式(13)、(14)、(15)、(16),利用(301)中得到归一化后时频域PHAT加权的零阶特征波束D′0(t,ω)、加权的一阶特征波束的两个正交分量D′x(t,ω)和D′y(t,ω)，用时频域PHAT加权的圆谐伪声强声源定位方法求解得到粗估计角度
Figure BDA0001602091820000101
式中，I′x(tα,ωβ)和I′y(tα,ωβ)分别表示时频域PHAT加权后各时频点处的瞬时伪声强的x轴分量和y轴分量。
步骤四：利用步骤三中求得的粗估计角度
Figure BDA0001602091820000102
构造空域波束指向性函数
Figure BDA0001602091820000103
结合步骤三中求得的归一化后时频域PHAT加权零阶特征波束，得到归一化后时频域PHAT加权联合指向性加权零阶特征波束
Figure BDA0001602091820000104
最后利用时频-空域联合加权的圆谐域伪声强方法求解得到精确声源估计角度
Figure BDA0001602091820000105
(401)、对于全向麦克风而言，并没有对特定角度的声源信号进行增强或削弱，这就使得在进行声源定位时往往受到其他方向的干扰。采用空域波束指向性加权的方法，可以只增强期望方向上的信号，削弱其他方向上的干扰信号，提高信干比，使得输出结果中特定方向上的信息能量增大。但是对于指向性加权方法，其定位性能又与加权函数的方位角估计值偏差大小紧密相关，若方位角估计值不准确，会导致其在偏离真实值的方向上形成波束，进而使得最终定位精度的偏差。如果，选择一种更加精准的方法进行方位角粗估计来得到指向性加权函数，最终的定位性能将会有提升。
本发明考虑用时频域PHAT加权的圆谐伪声强声源定位方法进行方位角的粗估计；然后用粗估计的结果构造波束指向性加权函数，将该加权函数作用于时频域PHAT加权后的各时频点的零阶特征波束信息；最后将加权后的时频点利用平均伪声强方法进行方位估计。
在圆谐域中，对于圆形麦克风阵列，可以通过对n阶特征波束进行加权组合，从而可以指向任意期望的方向。因此，空域波束形成器的输出响应表达式为
Figure BDA0001602091820000106
式中，δu表示任意方位角度。wn(t,ω,δu)表示n阶波束形成器系数，可表示为
Figure BDA0001602091820000107
式中，dn(t,ω)表示波束图调整参数常取为1，bn(t,ω)表示不同时频点的n阶贝塞尔函数。将式(19)代入式(18)，可将空域波束形成器的输出响应表达式化简为
Figure BDA0001602091820000111
且式(20)是对所有角度进行扫描后得到的每个角度的幅度响应。
在实际情况下，加权之前可以估计出大致的声源方位，只要计算出该声源方位上的指向性幅度响应即可继续进行后续的精确定位，这样在简化扫描过程的同时，减小计算量，并提高了定位精确性。
根据式(20)，利用步骤三中求得的粗估计角度
Figure BDA0001602091820000112
构造空域波束指向性函数
Figure BDA0001602091820000113
(402)、将(401)中的空域波束指向性函数的幅度
Figure BDA0001602091820000114
加权到步骤三中得到的归一化后时频域PHAT加权零阶特征波束上，得到归一化后时频域PHAT加权联合指向性加权零阶特征波束
Figure BDA0001602091820000115
其中，ψ(t,ω)表示
Figure BDA0001602091820000116
的相位；
(403)、利用时频-空域联合加权的圆谐域伪声强方法求解得到精确声源估计角度
Figure BDA0001602091820000117
仿真实验环境为7.8m×7.1m×3m的房间冲激响应模型，运用含有6个麦克风的圆形阵列，相邻麦克风间的夹角为60°，阵列尺寸直径为0.04m，声速c＝340m/s，阵列中心选取为房间中心位置，即坐标(3.9,3.55,1.5)m，声源到阵列中心的距离为2m，信号采样频率为16kHz，噪声为加性高斯白噪声，短时傅里叶变换的帧长为512点，帧移为256点。
实测实验环境为实验室，房间尺寸为9.7m×7.1m×3m，房间混响时间约为300ms，信噪比约为15dB，采用直径为0.04m的圆形传声器阵列为，麦克风选用MAP201型号全向性传声器，信号采集卡选用型号为USB 4432的采集卡，信号采样频率为16kHz。实验时，阵列中心置于房间中心位置，阵列距地面高度为1.5m，声源到阵列中心的距离均为2m，信号采集长度为1s，短时傅里叶变换帧长为512点，帧移为256点。
主要采用以下评价指标：均方根误差(Root Mean Square Error，RMSE)、平均均方根误差。其中，第η个方位角的均方根误差的计算公式表示为
Figure BDA0001602091820000121
式中，K是蒙特卡洛次数，
Figure BDA0001602091820000122
表示第η个方位角第k次蒙特卡洛实验的估计值，φ(η)表示第η个方位角的真实值。
平均均方根误差公式表示为
Figure BDA0001602091820000123
其中，L表示进行蒙特卡洛实验的方位角个数，这里选择L＝50。
分别从仿真实验和实测实验对不同信噪比、不同混响时间进行分析对比，计算圆谐域可控响应功率方法(Steered Response Power，SRP)、基本伪声强估计方法、时频域PHAT加权方法、谐波域PHAT加权方法和时频空域联合加权处理方法的平均均方根误差。为保证以上方法的公平比较，方法所用特征波束的最高阶数N＝1。
本发明中，仿真实验中声源方位角0°～360°以10°为间隔进行50次蒙特卡洛实验，观察平均均方根误差随混响时间的变化情况。实测实验中声源方位角由0°到330°以30°间隔逆时针转动，每个角度进行10次方位估计实验，以每个方位角的平均均方根误差作为评价指标。
图2是本发明在信噪比10dB、不同混响时间下，基本伪声强法、时频域PHAT加权法、谐波域PHAT加权法平均均方根误差的对比。图3是本发明在混响时间300ms、不同信噪比下，基本伪声强法、时频域PHAT加权法、谐波域PHAT加权法平均均方根误差的对比。图4是本发明在信噪比10dB、不同混响时间下，圆谐可控相应功率法、基本伪声强法、时频域PHAT加权法、时频-空域联合加权法的平均均方根误差比较对比。图5是本发明在混响时间300ms、不同信噪比下，圆谐可控相应功率法、基本伪声强法、时频域PHAT加权法、时频-空域联合加权法的平均均方根误差比较对比。图6是本发明实测实验进行声源定位结果的对比图。图2、图3、图4、图5、图6都是用来说明本发明定位效果。
仿真实验声源定位结果：
图2为信噪比10dB、不同混响时间下，基本伪声强估计方法、时频域PHAT加权方法、谐波域PHAT加权方法平均均方根误差的对比；图3为混响时间300ms、不同信噪比下平均均方根误差的对比。图4为信噪比10dB、不同混响时间下，各方法的平均均方根误差比较对比；图5为混响时间300ms、不同信噪比下各方法的平均均方根误差的对比。
由图2、图3可知，时频域PHAT加权法和谐波域PHAT加权法在混响时间较低、信噪比较高时，与基本伪声强方法性能相当。但从总体上看，时频域PHAT加权法整体性能优于谐波域PHAT加权法，且谐波域PHAT加权法性能整体差于基本伪声强法。在混响时间250ms时，时频域PHAT加权方法比谐波域PHAT加权方法平均均方根误差低0.4°，在混响时间600ms时，平均均方根误差低1.9°。
由图4的仿真结果可知，所有方法的平均均方根误差都随着混响时间的增加而增加。其中，圆谐SRP方法整体性能比较差；基本伪声强估计方法性能要优于圆谐SRP方法，这是由于圆谐SRP方法的定位性能依赖于特征波束的多少，当最高阶数越大圆谐SRP方法定位性能越好，而基本伪声强方法只需要零阶和一阶特征波束的信息就可以定位较为精确；由于PHAT加权具有一定的抑制干扰的能力，时频域PHAT加权方法要优于基本伪声强方法；本文所提出的时频-空域联合加权指向性方法的性能均优于圆谐域SRP方法、基本伪声强估计方法、时频域PHAT加权方法，其误差范围在3.5°～9.5°范围之内。
由图5的结果可知，随着信噪比的增大，各方法的平均均方根误差都有所减小。其中，圆谐SRP方法误差下降的幅度最大，即该方法对噪声的鲁棒性较差；基本伪声强估计方法相对圆谐SRP方法，定位误差有所减小；时频域PHAT加权方法整体性能都优于上述二种方法；而本文所提出的时频-空域联合加权指向性方法的定位误差相对更小。
因此从总体上看，时频-空域联合加权的圆谐域伪声强方法的定位性能最优。
实测实验声源定位结果：
图6是本发明实测实验进行声源定位结果的对比图。
由图6的实测结果可知，时频-空域联合加权的圆谐域伪声强方法的定位性能要优于其他方法，与仿真实验结果一致。
本发明方案所公开的技术手段不仅限于上述实施方式所公开的技术手段，还包括由以上技术特征任意组合所组成的技术方案。应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也视为本发明的保护范围。

Claims (4)

1.一种适用于圆形麦克风阵列的时频-空域联合加权的圆谐域伪声强声源定位方法，其特征在于，具体如下：
构建六元圆形麦克风阵列，采用构建的六元圆形麦克风阵列采集语音信号；
将采集到的语音信号变换到时频域，利用时频域PHAT加权进行滤波处理；
利用归一化后零阶特征波束、一阶特征波束，并由时频域PHAT加权的圆谐伪声强声源定位方法进行粗估计角度；
最后用粗估计角度构造空域波束指向性函数，得到归一化后时频域PHAT加权联合指向性加权零阶特征波束，并利用时频-空域联合加权的圆谐域伪声强方法求解得到精确声源估计角度；
圆谐域伪声强声源定位方法的具体步骤如下：
步骤一：采用Q个相同的全向性麦克风等间距的排列成半径为r的圆形麦克风阵列；
步骤二：对麦克风阵列采集到的声压信号Pq(t,ω)进行短时傅里叶变换，然后利用时频域PHAT加权对短时傅里叶变换后的声压信号进行预滤波，得到时频域PHAT加权后的n阶特征波束F′n(t,ω)；具体如下：
(201)、对声源发出的信号s(t)进行采样，得到第q个麦克风采集到的声压信号
Figure FDA0002929367970000011
其中，hq(t)表示声源到第q个麦克风间的脉冲响应，nq(t)表示第q个麦克风接收到的加性噪声信号，符号
Figure FDA0002929367970000012
表示线性卷积运算，q＝1,…,Q；
(202)、将(201)中采集到的声压信号pq(t)变换到时频域处理，经短时傅里叶变换得
Pq(t,ω)＝Hq(t,ω)·S(t,ω)+Nq(t,ω)
其中，Pq(t,ω)、S(t,ω)、Hq(t,ω)以及Nq(t,ω)分别表示麦克风接收声压信号pq(t)、声源信号s(t)、脉冲响应hq(t)和加性噪声信号nq(t)的短时傅里叶变换，(t,ω)表示时频率单元，t表示时间，ω表示频率；
(203)、利用PHAT加权对(202)中经短时傅里叶变换后的声压信号Pq(t,ω)进行预滤波，得到经过时频域PHAT加权后的n阶特征波束
Figure FDA0002929367970000013
其中，n表示阶数，
Figure FDA0002929367970000021
表示第1个麦克风与第q个麦克风之间按逆时针方向的夹角，
Figure FDA0002929367970000022
表示虚数单位，e表示自然指数；
步骤三：对步骤二中的n阶特征波束F′n(t,ω)进行归一化处理，得到归一化后时频域PHAT加权的零阶特征波束D′0(t,ω)、加权的一阶特征波束的两个正交分量D′x(t,ω)和D′y(t,ω)，然后用时频域PHAT加权的圆谐伪声强声源定位方法求解得到粗估计角度
Figure FDA0002929367970000023
具体如下：
(301)、对步骤二中的n阶特征波束F′n(t,ω)进行归一化处理，取阶数n＝0时得到归一化后时频域PHAT加权的零阶特征波束
Figure FDA0002929367970000024
其中，b0(t,ω)表示不同时频点的0阶贝塞尔函数；
取阶数n＝1时特征波束x轴和y轴的两个正交分量得到归一化后时频域PHAT加权的一阶特征波束的两个正交分量D′x(t,ω)和D′y(t,ω)，
Figure FDA0002929367970000025
Figure FDA0002929367970000026
其中，b1(t,ω)表示不同时频点的1阶贝塞尔函数，γx,1＝ei·1·0表示一阶特征波束的x轴旋转系数，
Figure FDA0002929367970000027
表示一阶特征波束的y轴旋转系数，F′x(t,ω)表示1阶特征波束x轴分量，F′y(t,ω)表示1阶特征波束y轴分量；
(302)、将(301)中得到加权的零阶特征波束D′0(t,ω)、加权的一阶特征波束的两个正交分量D′x(t,ω)和D′y(t,ω)，采用时频域PHAT加权的圆谐伪声强声源定位方法求解得到粗估计角度
Figure FDA0002929367970000031
其中，tα和ωβ表示第α个时间单元和第β个频率单元，*表示取共轭运算，Re表示取实部运算，arctan表示反正切运算，I′x(tα,ωβ)和I′y(tα,ωβ)分别表示时频域PHAT加权后各时频点处的瞬时伪声强的x轴分量和y轴分量；
步骤四：利用步骤三中求得的粗估计角度
Figure FDA0002929367970000032
构造空域波束指向性函数
Figure FDA0002929367970000033
结合步骤三中求得的归一化后时频域PHAT加权零阶特征波束，得到归一化后时频域PHAT加权联合指向性加权零阶特征波束
Figure FDA0002929367970000034
最后利用时频-空域联合加权的圆谐域伪声强方法求解得到精确声源估计角度
Figure FDA0002929367970000035
具体如下：
(401)、利用步骤三中求得的粗估计角度
Figure FDA0002929367970000036
构造空域波束指向性函数
Figure FDA0002929367970000037
其中，n表示阶数，且n的最高阶数为N＝kr，
Figure FDA0002929367970000038
为波数，f为频率，c为声波传播的速度，bn(t,ω)表示不同时频点的n阶贝塞尔函数；
(402)、将(401)中的空域波束指向性函数的幅度
Figure FDA0002929367970000039
加权到步骤三中得到的归一化后时频域PHAT加权零阶特征波束上，得到归一化后时频域PHAT加权联合指向性加权零阶特征波束
Figure FDA00029293679700000310
其中，ψ(t,ω)表示
Figure FDA00029293679700000311
的相位；
(403)、利用时频域-空域联合加权的圆谐伪声强方法求解得到精确声源估计角度
Figure FDA00029293679700000312
2.根据权利要求1所述的一种适用于圆形麦克风阵列的时频-空域联合加权的圆谐域伪声强声源定位方法，其特征在于，所述(302)中时频域PHAT加权的圆谐伪声强声源定位方法是采用有功伪声强进行瞬时方位估计，对各时频点伪声强进行平均相互补偿。
3.根据权利要求1所述的一种适用于圆形麦克风阵列的时频-空域联合加权的圆谐域伪声强声源定位方法，其特征在于，所述(403)中时频域-空域联合加权的圆谐伪声强方是采用有功伪声强进行瞬时方位估计，对各时频点伪声强进行平均相互补偿。
4.根据权利要求1所述的一种适用于圆形麦克风阵列的时频-空域联合加权的圆谐域伪声强声源定位方法，其特征在于，所述(402)中只保留步骤三中得到的归一化后时频域PHAT加权零阶特征波束D′0(t,ω)的相位ψ(t,ω)，幅度值替换为
Figure FDA0002929367970000041
得到归一化后时频域PHAT加权联合指向性加权零阶特征波束
Figure FDA0002929367970000042
CN201810228816.8A 2018-03-20 2018-03-20 一种时频-空域联合加权的圆谐域伪声强声源定位方法 Active CN108549052B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201810228816.8A CN108549052B (zh) 2018-03-20 2018-03-20 一种时频-空域联合加权的圆谐域伪声强声源定位方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201810228816.8A CN108549052B (zh) 2018-03-20 2018-03-20 一种时频-空域联合加权的圆谐域伪声强声源定位方法

Publications (2)

Publication Number Publication Date
CN108549052A CN108549052A (zh) 2018-09-18
CN108549052B true CN108549052B (zh) 2021-04-13

Family

ID=63516670

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201810228816.8A Active CN108549052B (zh) 2018-03-20 2018-03-20 一种时频-空域联合加权的圆谐域伪声强声源定位方法

Country Status (1)

Country Link
CN (1) CN108549052B (zh)

Families Citing this family (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109684723B (zh) * 2018-12-24 2022-05-10 哈尔滨工程大学 一种二维结构内部声学性能分析方法
CN109669160B (zh) * 2019-02-21 2022-07-05 哈尔滨工程大学 一种水下瞬态声信号的检测方法
CN109802717A (zh) * 2019-03-12 2019-05-24 西北工业大学 一种基于张量的极化-空域波束形成抗干扰方法
CN110596644B (zh) * 2019-09-24 2022-03-08 中国科学院声学研究所 一种使用移动环形传声器阵列的声源定位方法及系统
CN110954866B (zh) * 2019-11-22 2022-04-22 达闼机器人有限公司 声源定位方法、电子设备及存储介质
CN113655440B (zh) * 2021-08-09 2023-05-30 西南科技大学 一种自适应折中预白化的声源定位方法
CN115407270B (zh) * 2022-08-19 2023-11-17 苏州清听声学科技有限公司 一种分布式阵列的声源定位方法

Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101064942A (zh) * 2006-04-30 2007-10-31 北京信威通信技术股份有限公司 一种动态信道分配及空间特征提取和波束赋形方法
CN101566683A (zh) * 2009-03-24 2009-10-28 西北工业大学 基于相位差波束形成的目标方位估计方法
CN102854494A (zh) * 2012-08-08 2013-01-02 Tcl集团股份有限公司 一种声源定位方法及装置
CN104142492A (zh) * 2014-07-29 2014-11-12 佛山科学技术学院 一种srp-phat多源空间定位方法
CN104898086A (zh) * 2015-05-19 2015-09-09 南京航空航天大学 适用于微型麦克风阵列的声强估计声源定向方法
CN106093864A (zh) * 2016-06-03 2016-11-09 清华大学 一种麦克风阵列声源空间实时定位方法
CN106886010A (zh) * 2017-01-17 2017-06-23 南京航空航天大学 一种基于微型麦克风阵列的声源方位识别方法
CN107170441A (zh) * 2017-06-22 2017-09-15 西北工业大学 圆环阵最优频率不变响应超指向性波束形成方法

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US8289274B2 (en) * 2004-01-13 2012-10-16 Sliwa John W Microdroplet-based 3-D volumetric displays utilizing emitted and moving droplet projection screens

Patent Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101064942A (zh) * 2006-04-30 2007-10-31 北京信威通信技术股份有限公司 一种动态信道分配及空间特征提取和波束赋形方法
CN101566683A (zh) * 2009-03-24 2009-10-28 西北工业大学 基于相位差波束形成的目标方位估计方法
CN102854494A (zh) * 2012-08-08 2013-01-02 Tcl集团股份有限公司 一种声源定位方法及装置
CN104142492A (zh) * 2014-07-29 2014-11-12 佛山科学技术学院 一种srp-phat多源空间定位方法
CN104898086A (zh) * 2015-05-19 2015-09-09 南京航空航天大学 适用于微型麦克风阵列的声强估计声源定向方法
CN106093864A (zh) * 2016-06-03 2016-11-09 清华大学 一种麦克风阵列声源空间实时定位方法
CN106886010A (zh) * 2017-01-17 2017-06-23 南京航空航天大学 一种基于微型麦克风阵列的声源方位识别方法
CN107170441A (zh) * 2017-06-22 2017-09-15 西北工业大学 圆环阵最优频率不变响应超指向性波束形成方法

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
"3D SOURCE LOCALIZATION IN THE SPHERICAL HARMONIC DOMAIN USING A PSEUDOINTENSITY VECTOR";Daniel P.Jarrett etc.;《EUSIPCO-2010》;20100827;442-446 *
"Direction of Arrival Estimation in the Spherical Harmonic Domain Using Subspace Pseudointensity Vectors";Alastair H. Moore etc.;《IEEE/ACM TRANSACTIONS ON AUDIO, SPEECH, AND LANGUAGE PROCESSING》;20170131;第25卷(第1期);178-192 *
"基于差分麦克风阵列的声源定位方法研究";何赛娟;《中国优秀硕士学位论文全文数据库信息科技辑》;20170315;I136-601 *

Also Published As

Publication number Publication date
CN108549052A (zh) 2018-09-18

Similar Documents

Publication Publication Date Title
CN108549052B (zh) 一种时频-空域联合加权的圆谐域伪声强声源定位方法
Jarrett et al. 3D source localization in the spherical harmonic domain using a pseudointensity vector
Ma et al. Theoretical and practical solutions for high-order superdirectivity of circular sensor arrays
CN111123192B (zh) 一种基于圆形阵列和虚拟扩展的二维doa定位方法
US20050195988A1 (en) System and method for beamforming using a microphone array
US10455323B2 (en) Microphone probe, method, system and computer program product for audio signals processing
CN103181190A (zh) 用于远场多源追踪和分离的系统、方法、设备和计算机可读媒体
Huleihel et al. Spherical array processing for acoustic analysis using room impulse responses and time-domain smoothing
CN111798869B (zh) 一种基于双麦克风阵列的声源定位方法
Zhao et al. Open-lake experimental investigation of azimuth angle estimation using a single acoustic vector sensor
Xu et al. A modified differential beamforming and its application for DOA estimation of low frequency underwater signal
Diaz-Guerra et al. Direction of arrival estimation of sound sources using icosahedral CNNs
He et al. Closed-form DOA estimation using first-order differential microphone arrays via joint temporal-spectral-spatial processing
Imran et al. A methodology for sound source localization and tracking: Development of 3D microphone array for near-field and far-field applications
Tourbabin et al. Optimal real-weighted beamforming with application to linear and spherical arrays
Calmes et al. Azimuthal sound localization using coincidence of timing across frequency on a robotic platform
KR101354960B1 (ko) 영역 개념을 이용한 음파 입사 방향 추정 방법
Moore et al. 2D direction of arrival estimation of multiple moving sources using a spherical microphone array
Firoozabadi et al. Combination of nested microphone array and subband processing for multiple simultaneous speaker localization
CN113491137B (zh) 具有分数阶的灵活差分麦克风阵列
Swanson et al. Small-aperture array processing for passive multi-target angle of arrival estimation
Gur Modal beamforming for small circular arrays of particle velocity sensors
Berkun et al. A tunable beamformer for robust superdirective beamforming
Liu Spherical array superdirective beamforming based on spherical harmonic decomposition of the soundfield
Zhu et al. High-order Superdirectivity of Line Acoustic Vector Sensor Arrays

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant