[go: up one dir, main page]

CN114035157A - 一种基于期望最大化算法的分频带时延估计方法及其系统 - Google Patents

一种基于期望最大化算法的分频带时延估计方法及其系统 Download PDF

Info

Publication number
CN114035157A
CN114035157A CN202111274630.4A CN202111274630A CN114035157A CN 114035157 A CN114035157 A CN 114035157A CN 202111274630 A CN202111274630 A CN 202111274630A CN 114035157 A CN114035157 A CN 114035157A
Authority
CN
China
Prior art keywords
band
cross
full
signal
sub
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.)
Granted
Application number
CN202111274630.4A
Other languages
English (en)
Other versions
CN114035157B (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.)
Institute of Acoustics CAS
Original Assignee
Institute of Acoustics CAS
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 Institute of Acoustics CAS filed Critical Institute of Acoustics CAS
Priority to CN202111274630.4A priority Critical patent/CN114035157B/zh
Publication of CN114035157A publication Critical patent/CN114035157A/zh
Application granted granted Critical
Publication of CN114035157B publication Critical patent/CN114035157B/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)
  • Measurement Of Velocity Or Position Using Acoustic Or Ultrasonic Waves (AREA)

Abstract

本发明属于通信技术领域，具体涉及一种基于期望最大化算法的分频带时延估计方法，包括：传感器阵列接收任意的两个声源信号，分别对第一声源信号和第二声源信号进行离散时间傅里叶变换，得到对应的第一声信号和第二声信号；采用广义互相关时延估计算法，得到第一声信号和第二声信号之间广义互相关函数；提取两个声信号的互功率谱，将其拆分成多个子频带，进而建立全频带互相关观测数据，并得到全频带矩阵；建立全频带矩阵与全频带互相关观测数据之间的映射关系，进而建立全频带矩阵与由均值和协方差矩阵组成的观测数据之间的新映射关系；采用期望最大化算法，计算每个子频带的估计时延，并进行多次迭代，并将该稳定估计时延作为最终的时延。

Description

一种基于期望最大化算法的分频带时延估计方法及其系统
技术领域
本发明属于通信技术领域，具体地说，涉及一种基于期望最大化算法的分频带时延估计方法及其系统。
背景技术
基于传感器阵列声学测量模型，时延估计(Time Delay Estimatin，简称TDE)是准确测量声源信息的关键技术之一，主要是利用各个传感器节点与参考传感器节点的信号到达时间差(TDOA)。传统上，时延估计TDE在许多定位系统中发挥了重要作用，包括雷达、声纳、无线系统或地震学。在声学信号处理中，时延估计TDE对于定位和跟踪声源至关重要。
时延估计中，应用最广的算法为广义互相关函数(Generalized Cross-Correlation,GCC)。该算法是由Knapp和Carter在1976年使用最大似然估计器提出。然而，在实际应用环境中，有色噪声、脉冲噪声和混响、延迟估计性能的影响，广义互相关算法性能急剧衰减。为此，出现了多种方法来提高互相关函数加权，锐化互相关函数峰值。根据加权形式和准则的不同，如ROTH-GCC(GCC with Roth transform)、SCOT-GCC(GCC withsmoothed coherence transform)、PHAT-GCC(GCC with phase transform)等时延估计方法。Cobos et.al提出了频率滑动广义互相关(FS-GCC)，其旨在解决两个传感器子带时间延迟差异估计问题。FS-GCC利用滑动窗口方法求得频谱相位，获得一组子带GCC，对不同频段进行加权，以准确估计时延差。这些方法本质上利用频域滤波获得高信噪比频带下的时延估计，需要预先知道声源信号所在的高信噪比频带。
发明内容
为解决现有技术存在的上述缺陷，本发明提出了一种基于期望最大化算法的分频带时延估计方法，本发明则对各个子带进行概率建模，利用EM算法优化收敛到高信噪比频带，获得准确的时延估计结果。该方法包括：
传感器阵列接收任意的两个声源信号，记为第一声源信号和第二声源信号；分别对第一声源信号和第二声源信号进行离散时间傅里叶变换，得到对应的第一声信号和第二声信号；
根据第一声信号和第二声信号，采用广义互相关时延估计算法，得到第一声信号和第二声信号之间广义互相关函数；
根据该互相关函数，提取两个声信号的互功率谱，并将其拆分成多个子频带，进而建立全频带互相关观测数据，并基于建立的全频带互相关观测数据，得到全频带矩阵；
建立全频带矩阵与全频带互相关观测数据之间的映射关系，根据该映射关系，确定均值和协方差矩阵，进而建立全频带矩阵与由均值和协方差矩阵组成的观测数据之间的新映射关系；
将该新映射关系带入条件期望函数，采用期望最大化算法，计算每个子频带的估计时延，并进行多次迭代，直至收敛于某一个稳定不变的稳定估计时延，且该稳定估计时延不再增加，则停止迭代；
并将该稳定估计时延作为最终的时延，完成分频带时延估计。
作为上述技术方案的改进之一，所述传感器阵列接收任意的两个声源信号，记为第一声源信号和第二声源信号；分别对第一声源信号和第二声源信号进行离散时间傅里叶变换，得到对应的第一声信号和第二声信号；其具体过程为：
传感器阵列接收任意的两个声源信号，记为第一声源信号x1(t)和第二声源信号x2(t)：
Figure BDA0003329044100000021
其中，a1是第一振幅衰减因子；a2是第二振幅衰减因子；s(t)是声源发射的信号，n1(t)是第一加性白噪声；n2(t)是第二加性白噪声；τ1为第一声源信号的时延；τ2为第二声源信号的时延；
分别对第一声源信号x1(t)和第二声源信号x2(t)进行离散时间傅里叶变换，得到对应的第一声信号X1(ω)和第二声信号X2(ω)：
Figure BDA0003329044100000022
其中，S(ω)为声源发射信号s(t)的傅里叶频谱；W1(ω)为第一加性白噪声n1(t)的傅里叶频谱；W2(ω)为第二加性白噪声n2(t)的傅里叶频谱。
作为上述技术方案的改进之一，所述根据第一声信号和第二声信号，采用广义互相关时延估计算法，得到第一声信号和第二声信号之间广义互相关函数；其具体过程为：
根据第一声信号X1(ω)和第二声信号X2(ω)，采用广义互相关时延估计算法，得到第一声信号和第二声信号之间广义互相关函数R(τ)；
Figure BDA0003329044100000031
其中，F-1为傅里叶逆变换；
Figure BDA0003329044100000032
为第一声源信号x1(t)和第二声源信号x2(t)的互功率谱；
Figure BDA0003329044100000033
为PHAT加权滤波后的第一声源信号x1(t)和第二声源信号x2(t)的互功率谱；Ψ(ω)为采用PHAT进行加权滤波后的信号；
其中，
Figure BDA0003329044100000034
其中，H1(ω)为第一通道传递函数；H2(ω)为第二通道传递函数；*为共轭。
作为上述技术方案的改进之一，所述根据该互相关函数，提取两个声信号的互功率谱，并将其拆分成多个子频带，进而建立全频带互相关观测数据，并基于建立的全频带互相关观测数据，得到全频带矩阵；其具体过程为：
根据该互相关函数R(τ)，提取两个声信号的互功率谱
Figure BDA0003329044100000035
并将其拆分成多个子频带rk(n)，
Figure BDA0003329044100000036
其中，
Figure BDA0003329044100000037
为理想的、不具有噪声的第k子频带的互相关函数；wk(n)为第k子频带的噪声；
其中，K为分频带数，则：
Figure BDA0003329044100000038
其中，ωα为分频带宽；T为全频带带宽；
进而建立全频带互相关观测数据R(n)；
Figure BDA0003329044100000039
其中，w(n)为加性噪声；r(n)为理想下全频带的广义互相关函数；
并基于建立的全频带互相关观测数据R(n)，得到全频带矩阵L(n)；
Figure BDA00033290441000000310
其中，lK(n)为第k个子频带的广义互相关函数；
Figure BDA0003329044100000041
其中，wK(n)为相互独立的零均值高斯变量。
作为上述技术方案的改进之一，所述建立全频带矩阵与全频带互相关观测数据之间的映射关系，根据该映射关系，确定均值和协方差矩阵，进而建立全频带矩阵与由均值和协方差矩阵组成的观测数据之间的新映射关系；其具体过程为：
根据得到的全频带互相关观测数据R(n)和全频带矩阵L(n)，建立全频带矩阵与全频带互相关观测数据之间的映射关系：
R(n)＝f[L(n)]＝[1，...，1]L(n) (8)根据该映射关系，确定均值μ(n)和协方差矩阵Q(n)；
Figure BDA0003329044100000042
Figure BDA0003329044100000043
其中，
Figure BDA0003329044100000044
为理想的、不具有噪声的第k子频带的互相关函数；
Figure BDA0003329044100000045
为方差；
Figure BDA0003329044100000046
Figure BDA0003329044100000047
其中，σ2为加性噪声w(n)的方差；
进而建立全频带矩阵与由均值和协方差矩阵组成的观测数据之间的新映射关系；
L(n)～N(μ(n))，Q(n))，n＝1，...，N (13)
其中，n为第n个离散采样点。
作为上述技术方案的改进之一，所述将该新映射关系带入条件期望函数，采用期望最大化算法，计算每个子频带的估计时延，并进行多次迭代，直至收敛于某一个稳定不变的稳定估计时延，且该稳定估计时延不再增加，则停止迭代；并将该稳定估计时延作为最终的时延，完成分频带时延估计；其具体过程为：
建立全频带矩阵与由均值和协方差矩阵组成的观测数据之间的新映射关系；
L(n)～N(μ(n))，Q(n))，n＝1，...，N (13)
其中，n为第n个离散采样点；
将建立的映射关系带入到期望最大化算法(即EM算法)的条件期望函数U(θ，θ′)中：
Figure BDA0003329044100000051
其中，θ为待估计时延参数，θ＝τ；E[]为期望最大化算法的条件期望的期望运算；lnp()为自然对数运算；θ′为最大化期望算法迭代过程中θ的估计值；R为全频带互相关估计结果；e为时延估计误差；
Figure BDA0003329044100000052
为第k个子带的lK(n)估计结果；
从上式，U(θ，θ′)要取得最大值，需满足：
Figure BDA0003329044100000053
对每个子频带分别进行E步和M步为：
E步：对于k＝1，...，K；并利用前面迭代的时延估计结果
Figure BDA0003329044100000054
代替θ′，从而构造各个子频带的
Figure BDA0003329044100000055
是最大化期望算法第m次迭代过程中第k个子带的互相关函数：
Figure BDA0003329044100000056
其中，
Figure BDA0003329044100000057
为第m次迭代计算时的第k个子带的lK(n)估计结果；
M步：根据上一步得到的
Figure BDA0003329044100000058
得到第(m+1)次迭代的时延估计结果
Figure BDA0003329044100000059
Figure BDA00033290441000000510
其中，在无回响与噪声的条件下
Figure BDA00033290441000000511
其中，Ψ1(ω)为在无回响与噪声的条件下，采用广义互相关的加权函数进行加权滤波后的信号；a1为第一信号的振幅因子；a2为第二信号的振幅银子；S(ω)为声源发射信号的傅里叶频谱；ω为角频率；τ0为真实时延估计结果；j＝-1；
因此，
Figure BDA00033290441000000512
其中，φk为包含位移窗口响应的N个样本向量；
Figure BDA00033290441000000513
为为理想的、不具有噪声的第k子频带的互相关函数；
Figure BDA00033290441000000514
为第k个子频带的频率范围；
Figure BDA0003329044100000061
其中，ωα为分频带宽，K为分频带数；
因此，上述对每个子频带分别进行E步和M步进一步化简为：
E步：
Figure BDA0003329044100000062
其中，
Figure BDA0003329044100000063
为第m次迭代计算时的时延估计结果；φk为包含位移窗口响应的N个样本向量；
M步：
Figure BDA0003329044100000064
其中，
Figure BDA0003329044100000065
f0为信号的频率，fs为采样频率，*代表共轭；
因此，得到第(m+1)次迭代的时延估计
Figure BDA0003329044100000066
在上述迭代过程中，每次迭代都会增大每个子频带的最大似然函数值，当收敛于似然函数的某一个稳定点时，如果继续迭代，参数的估计值都不再发生变化，似然函数也不再变化，则判断迭代终止；
并将该稳定估计时延作为最终的时延，完成分频带时延估计。
本发明还提供了一种基于期望最大化算法的分频带时延估计系统，该系统包括：
声信号获取模块，用于传感器阵列接收任意的两个声源信号，记为第一声源信号和第二声源信号；分别对第一声源信号和第二声源信号进行离散时间傅里叶变换，得到对应的第一声信号和第二声信号；
函数获取模块，用于根据第一声信号和第二声信号，采用广义互相关时延估计算法，得到第一声信号和第二声信号之间广义互相关函数；
全频带矩阵建立模块，用于根据该互相关函数，提取两个声信号的互功率谱，并将其拆分成多个子频带，进而建立全频带互相关观测数据，并基于建立的全频带互相关观测数据，得到全频带矩阵；
映射模块，用于建立全频带矩阵与全频带互相关观测数据之间的映射关系，根据该映射关系，确定均值和协方差矩阵，进而建立全频带矩阵与由均值和协方差矩阵组成的观测数据之间的新映射关系；和
估计模块，用于将该新映射关系带入条件期望函数，采用期望最大化算法，计算每个子频带的估计时延，并进行多次迭代，直至收敛于某一个稳定不变的稳定估计时延，且该稳定估计时延不再增加，则停止迭代；
并将该稳定估计时延作为最终的时延，完成分频带时延估计。
本发明还提供了一种计算机设备，包括存储器、处理器及存储在所述存储器上并可在所述处理器上运行的计算机程序，所述处理器执行所述计算机程序时实现所述的方法。
本发明还提供了一种计算机可读存储介质，所述计算机可读存储介质存储有计算机程序，所述计算机程序当被处理器执行时使所述处理器执行所述的方法。
本发明与现有技术相比的有益效果是：
本发明的方法通过对各子带信号进行概率建模，利用EM算法迭代优化收敛获得高信噪比频带下的时延估计结果，从而提高宽带信号在低信噪比下的时延估计精度。
附图说明
图1是本发明的一种基于期望最大化算法的分频带时延估计方法的流程图；
图2(a)是采用EM-FSGCC算法，时间t与幅度之间的曲线示意图；
图2(b)是采用GCC-PHAT算法，时间t与幅度之间的曲线示意图；
图2(c)是采用WSVD FS-GCC算法，时间t与幅度之间的曲线示意图；
图3(a)是采用GCC-EM算法、GCC-PHAT算法、WSVD FS-GCC算法，信噪比SNR与MAE之间的曲线示意图；
图3(b)采用GCC-EM算法、GCC-PHAT算法、WSVD FS-GCC算法，信噪比SNR与P之间的曲线示意图；
图3(c)采用GCC-EM算法、GCC-PHAT算法、WSVD FS-GCC算法，信噪比SNR与SDAE之间的曲线示意图。
具体实施方式
现结合附图和实例对本发明作进一步的描述。
如图1所示，本发明提供了一种基于期望最大化算法的分频带时延估计方法，该方法将GCC-PHAT(GCC with phase transform)分成多个无重叠频带，然后利用最大似然分别对每个频带估计时延；该方法对每个频带反复迭代，使用当前时延估计，用于分解GCC-PHAT，从而改进下一次时延估计。在正则性条件下，该方法收敛于似然函数一个平稳点，其中每个迭代周期都增加了估计时延的似然性。
该方法具体包括：
传感器阵列接收任意的两个声源信号，记为第一声源信号和第二声源信号；分别对第一声源信号和第二声源信号进行离散时间傅里叶变换，得到对应的第一声信号和第二声信号；
具体地，传感器阵列接收任意的两个声源信号，记为第一声源信号x1(t)和第二声源信号x2(t)：
Figure BDA0003329044100000081
其中，a1是第一振幅衰减因子；a2是第二振幅衰减因子；s(t)是声源发射的信号，n1(t)是第一加性白噪声；n2(t)是第二加性白噪声；τ1为第一声源信号的时延；τ2为第二声源信号的时延；
分别对第一声源信号x1(t)和第二声源信号x2(t)进行离散时间傅里叶变换，得到对应的第一声信号X1(ω)和第二声信号X2(ω)：
Figure BDA0003329044100000082
其中，S(ω)为声源发射信号s(t)的傅里叶频谱；W1(ω)为第一加性白噪声n1(t)的傅里叶频谱；W2(ω)为第二加性白噪声n2(t)的傅里叶频谱。
根据第一声信号和第二声信号，采用广义互相关时延估计算法，得到第一声信号和第二声信号之间广义互相关函数；
具体地，根据第一声信号X1(ω)和第二声信号X2(ω)，采用广义互相关时延估计算法，得到第一声信号和第二声信号之间广义互相关函数R(τ)；
Figure BDA0003329044100000083
其中，F-1为傅里叶逆变换；
Figure BDA0003329044100000084
为第一声源信号x1(t)和第二声源信号x2(t)的互功率谱；
Figure BDA0003329044100000085
为PHAT加权滤波后的第一声源信号x1(t)和第二声源信号x2(t)的互功率谱；Ψ(ω)为采用PHAT进行加权滤波后的信号；
其中，
Figure BDA0003329044100000086
其中，H1(ω)为第一通道传递函数；H2(ω)为第二通道传递函数；*为共轭。
根据该互相关函数，提取两个声信号的互功率谱，并将其拆分成多个子频带，进而建立全频带互相关观测数据，并基于建立的全频带互相关观测数据，得到全频带矩阵；
具体地，根据该互相关函数R(τ)，提取两个声信号的互功率谱
Figure BDA0003329044100000091
并将其拆分成多个子频带rk(n)，
Figure BDA0003329044100000092
其中，
Figure BDA0003329044100000093
为理想的、不具有噪声的第k子频带的互相关函数；wk(n)为第k子频带的噪声；
其中，K为分频带数，则：
Figure BDA0003329044100000094
其中，ωα为分频带宽；T为全频带带宽；
进而建立全频带互相关观测数据R(n)；
Figure BDA0003329044100000095
其中，w(n)为加性噪声；r(n)为理想情况下全频带的广义互相关函数；
并基于建立的全频带互相关观测数据R(n)，得到全频带矩阵L(n)；
Figure BDA0003329044100000096
其中，lK(n)为第K个子频带的广义互相关函数；
Figure BDA0003329044100000097
其中，wK(n)为相互独立的零均值高斯变量。
建立全频带矩阵与全频带互相关观测数据之间的映射关系，根据该映射关系，确定均值和协方差矩阵，进而建立全频带矩阵与由均值和协方差矩阵组成的观测数据之间的新映射关系；
具体地，根据得到的全频带互相关观测数据R(n)和全频带矩阵L(n)，建立全频带矩阵与全频带互相关观测数据之间的映射关系：
R(n)＝f[L(n)]＝[1，...，1]L(n) (8)根据该映射关系，确定均值μ(n)和协方差矩阵Q(n)；
Figure BDA0003329044100000101
Figure BDA0003329044100000102
其中，
Figure BDA0003329044100000103
为理想的、不具有噪声的第k子频带的互相关函数；
Figure BDA0003329044100000104
为方差；
Figure BDA0003329044100000105
Figure BDA0003329044100000106
其中，σ2为加性噪声w(n)的方差；
进而建立全频带矩阵与由均值和协方差矩阵组成的观测数据之间的新映射关系；
L(n)～N(μ(n))，Q(n))，n＝1，...，N (13)
其中，n为第n个离散采样点。
将该新映射关系带入条件期望函数，采用期望最大化算法，计算每个子频带的估计时延，并进行多次迭代，直至收敛于某一个稳定不变的稳定估计时延，且该稳定估计时延不再增加，则停止迭代；
并将该稳定估计时延作为最终的时延，完成分频带时延估计。
具体地，建立全频带矩阵与由均值和协方差矩阵组成的观测数据之间的新映射关系；
L(n)～N(μ(n))，Q(n))，n＝1，...，N (13)
其中，n为第n个离散采样点；
将建立的映射关系带入到期望最大化算法(即EM算法)的条件期望函数U(θ，θ′)中：
Figure BDA0003329044100000107
其中，θ为待估计时延参数，θ＝τ；E[]为期望最大化算法的条件期望的期望运算；lnp()为自然对数运算；θ′为最大化期望算法迭代过程中θ的估计值；R为全频带互相关估计结果；e为时延估计误差；
Figure BDA0003329044100000108
为第k个子带的lK(n)估计结果；
从上式，U(θ，θ′)要取得最大值，需满足：
Figure BDA0003329044100000111
对每个子频带分别进行E步和M步为：
E步：对于k＝1，...，K；并利用前面迭代的时延估计结果
Figure BDA0003329044100000112
代替θ′，从而构造各个子频带的
Figure BDA0003329044100000113
是最大化期望算法第m次迭代过程中第k个子带的互相关函数：
Figure BDA0003329044100000114
其中，
Figure BDA0003329044100000115
为第m次迭代计算时的第k个子带的lK(n)估计结果；
M步：根据上一步得到的
Figure BDA0003329044100000116
得到第(m+1)次迭代的时延估计结果
Figure BDA0003329044100000117
Figure BDA0003329044100000118
其中，在无回响与噪声的条件下
Figure BDA0003329044100000119
其中，Ψ1(ω)为在无回响与噪声的条件下，采用广义互相关的加权函数进行加权滤波后的信号；a1为第一信号的振幅因子；a2为第二信号的振幅银子；S(ω)为声源发射信号的傅里叶频谱；ω为角频率；τ0为真实时延估计结果；j＝-1；
因此，
Figure BDA00033290441000001110
其中，φk为包含位移窗口响应的N个样本向量；
Figure BDA00033290441000001111
为为理想的、不具有噪声的第k子频带的互相关函数；
Figure BDA00033290441000001112
为第k个子频带的频率范围；
Figure BDA00033290441000001113
其中，ωα为分频带宽，K为分频带数；
因此，上述对每个子频带分别进行E步和M步进一步化简为：
E步：
Figure BDA00033290441000001114
其中，
Figure BDA00033290441000001115
为第m次迭代计算时的时延估计结果；φk为包含位移窗口响应的N个样本向量；
M步：
Figure BDA0003329044100000121
其中，
Figure BDA0003329044100000122
f0为信号的频率，fs为采样频率，*代表共轭；
因此，得到第(m+1)次迭代的时延估计
Figure BDA0003329044100000123
在上述迭代过程中，每次迭代都会增大每个子频带的最大似然函数值，当收敛于似然函数的某一个稳定点时，如果继续迭代，参数的估计值都不再发生变化，似然函数也不再变化，则判断迭代终止；
并将该稳定估计时延作为最终的时延，完成分频带时延估计。
本发明还提供了一种基于期望最大化算法的分频带时延估计系统，该系统包括：
声信号获取模块，用于传感器阵列接收任意的两个声源信号，记为第一声源信号和第二声源信号；分别对第一声源信号和第二声源信号进行离散时间傅里叶变换，得到对应的第一声信号和第二声信号；
函数获取模块，用于根据第一声信号和第二声信号，采用广义互相关时延估计算法，得到第一声信号和第二声信号之间广义互相关函数；
全频带矩阵建立模块，用于根据该互相关函数，提取两个声信号的互功率谱，并将其拆分成多个子频带，进而建立全频带互相关观测数据，并基于建立的全频带互相关观测数据，得到全频带矩阵；
映射模块，用于建立全频带矩阵与全频带互相关观测数据之间的映射关系，根据该映射关系，确定均值和协方差矩阵，进而建立全频带矩阵与由均值和协方差矩阵组成的观测数据之间的新映射关系；和
估计模块，用于将该新映射关系带入条件期望函数，采用期望最大化算法，计算每个子频带的估计时延，并进行多次迭代，直至收敛于某一个稳定不变的稳定估计时延，且该稳定估计时延不再增加，则停止迭代；
并将该稳定估计时延作为最终的时延，完成分频带时延估计。
本发明还提供了一种计算机设备，包括存储器、处理器及存储在所述存储器上并可在所述处理器上运行的计算机程序，所述处理器执行所述计算机程序时实现上述方法。
本发明还提供了一种计算机可读存储介质，所述计算机可读存储介质存储有计算机程序，所述计算机程序当被处理器执行时使所述处理器执行上述方法。
实施例1.
1.1性能准则
根据其绝对误差
Figure BDA0003329044100000131
将其分类为异常估计或者非异常估计；其中，τ0是真实时延，
Figure BDA0003329044100000132
是时延估计。
如果e＞Tc/2，则将其分为异常估计，其中，Tc是信号相关的时间。
对于仿真的特定信号源，Tc计算为自相关函数主瓣宽度(取SNR＝-3db)。
TDE(time delay estimation)性能包括：
定义异常估计百分比
Figure BDA0003329044100000133
绝对误差的平均值
Figure BDA0003329044100000134
标准偏差
Figure BDA0003329044100000135
(针对非异常估计的子集)等指标对算法的时延估计性能进行评估，这些指标定义为：
Figure BDA0003329044100000136
Figure BDA0003329044100000137
Figure BDA0003329044100000138
其中，NT表示估计的总数；Na被确定为异常值的估计数；FSPR被定义为最大GCC峰值相对于第二个较大峰值的平均增益(针对非异常估计的子集)。
1.2仿真设置及算法参数
1.3外场实验
考虑了在一个单源场景中用图像源方法模拟的矩形房间。在房间内设置传感器阵列的位置和方向，以及每个麦克风配置的随机声源位置，一对分离的传感器产生合成脉冲响应。对每个混响条件都重复进行了模拟。使用了以下参数：
1)房间尺寸：6×7×3米(长×宽×高)。
2)声源位置：平面上的随机位置(x，y，z＝1.25)。
3)麦克风位置：两个麦克风阵列，传感器间距为0.5m，在x-y平面(z＝1.25)上有随机位置和方向。
4)SNR：在-15dB和10dB之间变化。每个SNR条件生成不同的噪声实现。为了控制SNR，相互独立的高斯白噪声被适当地缩放并添加到每个麦克风信号中。
5)源信号：2秒的男性语音信号，以44.1kHz的16位分辨率数字化。
对于cobos等人的现有的方法，将采用4096个样本的帧长度和具有75％重叠的Hann窗口，所以B＝128(频谱窗口)，M＝32(跳频)；而对于本发明所提出的方法，则能够避免子频带之间的重叠，因此，B＝M＝32。因此，在每个信噪比和混响条件下，用于评估每种方法的估计值的总数为NT＝500。
1.3TDE结果
如图2(a)、2(b)和2(c)所示，展示出了，在信噪比较低时，不同算法所展现的效果图；其中，采用EM-FSGCC算法，在时间t＝91处，其对应的幅度是最大的，最大值为0.196；采用GCC-PHAT算法，在时间t＝-186处，其对应的幅度是最大的，最大值为0.058；采用WSVDFS-GCC算法，在时间t＝89处，其对应的幅度是最大的，最大值为0.056；
随信噪比变化的结果，如图3(a)、3(b)、3(c)所示。与传统的GCC-PHAT相比，本发明提出的方法，要优于WSVD-GCC。从图中可以看出，在低信噪比下，观察到了最大的改进，其中EM-FSGCC和传统GCC之间差异接近60个点，与WSVD-FSGCC相差40个点，随着信噪比的提升，所有算法都趋于一致。
最后所应说明的是，以上实施例仅用以说明本发明的技术方案而非限制。尽管参照实施例对本发明进行了详细说明，本领域的普通技术人员应当理解，对本发明的技术方案进行修改或者等同替换，都不脱离本发明技术方案的精神和范围，其均应涵盖在本发明的权利要求范围当中。

Claims (9)

1.一种基于期望最大化算法的分频带时延估计方法，该方法包括：
传感器阵列接收任意的两个声源信号，记为第一声源信号和第二声源信号；分别对第一声源信号和第二声源信号进行离散时间傅里叶变换，得到对应的第一声信号和第二声信号；
根据第一声信号和第二声信号，采用广义互相关时延估计算法，得到第一声信号和第二声信号之间广义互相关函数；
根据该互相关函数，提取两个声信号的互功率谱，并将其拆分成多个子频带，进而建立全频带互相关观测数据，并基于建立的全频带互相关观测数据，得到全频带矩阵；
建立全频带矩阵与全频带互相关观测数据之间的映射关系，根据该映射关系，确定均值和协方差矩阵，进而建立全频带矩阵与由均值和协方差矩阵组成的观测数据之间的新映射关系；
将该新映射关系带入条件期望函数，采用期望最大化算法，计算每个子频带的估计时延，并进行多次迭代，直至收敛于某一个稳定不变的稳定估计时延，且该稳定估计时延不再增加，则停止迭代；
并将该稳定估计时延作为最终的时延，完成分频带时延估计。
2.根据权利要求1所述的基于期望最大化算法的分频带时延估计方法，其特征在于，所述传感器阵列接收任意的两个声源信号，记为第一声源信号和第二声源信号；分别对第一声源信号和第二声源信号进行离散时间傅里叶变换，得到对应的第一声信号和第二声信号；其具体过程为：
传感器阵列接收任意的两个声源信号，记为第一声源信号x1(t)和第二声源信号x2(t)：
Figure FDA0003329044090000011
其中，a1是第一振幅衰减因子；a2是第二振幅衰减因子；s(t)是声源发射的信号，n1(t)是第一加性白噪声；n2(t)是第二加性白噪声；τ1为第一声源信号的时延；τ2为第二声源信号的时延；
分别对第一声源信号x1(t)和第二声源信号x2(t)进行离散时间傅里叶变换，得到对应的第一声信号X1(ω)和第二声信号X2(ω)：
Figure FDA0003329044090000012
其中，S(ω)为声源发射信号s(t)的傅里叶频谱；W1(ω)为第一加性白噪声n1(t)的傅里叶频谱；W2(ω)为第二加性白噪声n2(t)的傅里叶频谱。
3.根据权利要求1所述的基于期望最大化算法的分频带时延估计方法，其特征在于，所述根据第一声信号和第二声信号，采用广义互相关时延估计算法，得到第一声信号和第二声信号之间广义互相关函数；其具体过程为：
根据第一声信号X1(ω)和第二声信号X2(ω)，采用广义互相关时延估计算法，得到第一声信号和第二声信号之间广义互相关函数R(τ)；
Figure FDA0003329044090000021
其中，F-1为傅里叶逆变换；
Figure FDA0003329044090000022
为第一声源信号x1(t)和第二声源信号x2(t)的互功率谱；
Figure FDA0003329044090000023
为PHAT加权滤波后的第一声源信号x1(t)和第二声源信号x2(t)的互功率谱；Ψ(ω)为采用PHAT进行加权滤波后的信号；
其中，
Figure FDA0003329044090000024
其中，H1(ω)为第一通道传递函数；H2(ω)为第二通道传递函数；*为共轭。
4.根据权利要求3所述的基于期望最大化算法的分频带时延估计方法，其特征在于，所述根据该互相关函数，提取两个声信号的互功率谱，并将其拆分成多个子频带，进而建立全频带互相关观测数据，并基于建立的全频带互相关观测数据，得到全频带矩阵；其具体过程为：
根据该互相关函数R(τ)，提取两个声信号的互功率谱
Figure FDA0003329044090000025
并将其拆分成多个子频带rk(n)，
Figure FDA0003329044090000026
其中，
Figure FDA0003329044090000027
为理想的、不具有噪声的第k子频带的互相关函数；wk(n)为第k子频带的噪声；
其中，K为分频带数，则：
Figure FDA0003329044090000028
其中，ωα为分频带宽；T为全频带带宽；
进而建立全频带互相关观测数据R(n)；
Figure FDA0003329044090000029
其中，w(n)为加性噪声；r(n)为理想下全频带的广义互相关函数；
并基于建立的全频带互相关观测数据R(n)，得到全频带矩阵L(n)；
Figure FDA0003329044090000031
其中，lK(n)为第k个子频带的广义互相关函数；
Figure FDA0003329044090000032
其中，wK(n)为相互独立的零均值高斯变量。
5.根据权利要求4所述的基于期望最大化算法的分频带时延估计方法，其特征在于，所述建立全频带矩阵与全频带互相关观测数据之间的映射关系，根据该映射关系，确定均值和协方差矩阵，进而建立全频带矩阵与由均值和协方差矩阵组成的观测数据之间的新映射关系；其具体过程为：
根据得到的全频带互相关观测数据R(n)和全频带矩阵L(n)，建立全频带矩阵与全频带互相关观测数据之间的映射关系：
R(n)＝f[L(n)]＝[1，...，1]L(n) (8)
根据该映射关系，确定均值μ(n)和协方差矩阵Q(n)；
Figure FDA0003329044090000033
Figure FDA0003329044090000034
其中，
Figure FDA0003329044090000035
为理想的、不具有噪声的第k子频带的互相关函数；
Figure FDA0003329044090000036
为方差；
Figure FDA0003329044090000037
Figure FDA0003329044090000038
其中，σ2为加性噪声w(n)的方差；
进而建立全频带矩阵与由均值和协方差矩阵组成的观测数据之间的新映射关系；
L(n)～N(μ(n))，Q(n))，n＝1，...，N (13)
其中，n为第n个离散采样点。
6.根据权利要求5所述的基于期望最大化算法的分频带时延估计方法，其特征在于，所述将该新映射关系带入条件期望函数，采用期望最大化算法，计算每个子频带的估计时延，并进行多次迭代，直至收敛于某一个稳定不变的稳定估计时延，且该稳定估计时延不再增加，则停止迭代；并将该稳定估计时延作为最终的时延，完成分频带时延估计；其具体过程为：
建立全频带矩阵与由均值和协方差矩阵组成的观测数据之间的新映射关系；
L(n)～N(μ(n))，Q(n))，n＝1，...，N (13)
其中，n为第n个离散采样点；
将建立的映射关系带入到期望最大化算法(即EM算法)的条件期望函数U(θ，θ′)中：
Figure FDA0003329044090000041
其中，θ为待估计时延参数，θ＝τ；E[]为期望最大化算法的条件期望的期望运算；lnp()为自然对数运算；θ′为最大化期望算法迭代过程中θ的估计值；R为全频带互相关估计结果；e为时延估计误差；
Figure FDA0003329044090000042
为第k个子带的lK(n)估计结果；
从上式，U(θ，θ′)要取得最大值，需满足：
Figure FDA0003329044090000043
对每个子频带分别进行E步和M步为：
E步：对于k＝1，...，K；并利用前面迭代的时延估计结果
Figure FDA0003329044090000044
代替θ′，从而构造各个子频带的
Figure FDA0003329044090000045
是最大化期望算法第m次迭代过程中第k个子带的互相关函数：
Figure FDA0003329044090000046
其中，
Figure FDA0003329044090000047
为第m次迭代计算时的第k个子带的lK(n)估计结果；
M步：根据上一步得到的
Figure FDA0003329044090000048
得到第(m+1)次迭代的时延估计结果
Figure FDA0003329044090000049
Figure FDA00033290440900000410
其中，在无回响与噪声的条件下
Figure FDA00033290440900000411
其中，Ψ1(ω)为在无回响与噪声的条件下，采用广义互相关的加权函数进行加权滤波后的信号；a1为第一信号的振幅因子；a2为第二信号的振幅银子；S(ω)为声源发射信号的傅里叶频谱；ω为角频率；τ0为真实时延估计结果；j＝-1；
因此，
Figure FDA0003329044090000051
其中，φk为包含位移窗口响应的N个样本向量；
Figure FDA0003329044090000052
为为理想的、不具有噪声的第k子频带的互相关函数；
Figure FDA0003329044090000053
为第k个子频带的频率范围；
Figure FDA0003329044090000054
其中，ωα为分频带宽，K为分频带数；
因此，上述对每个子频带分别进行E步和M步进一步化简为：
E步：
Figure FDA0003329044090000055
其中，
Figure FDA0003329044090000056
为第m次迭代计算时的时延估计结果；φk为包含位移窗口响应的N个样本向量；
M步：
Figure FDA0003329044090000057
其中，
Figure FDA0003329044090000058
f0为信号的频率，fs为采样频率，*代表共轭；
因此，得到第(m+1)次迭代的时延估计
Figure FDA0003329044090000059
在上述迭代过程中，每次迭代都会增大每个子频带的最大似然函数值，当收敛于似然函数的某一个稳定点时，如果继续迭代，参数的估计值都不再发生变化，似然函数也不再变化，则判断迭代终止；
并将该稳定估计时延作为最终的时延，完成分频带时延估计。
7.一种基于期望最大化算法的分频带时延估计系统，其特征在于，该系统包括：
声信号获取模块，用于传感器阵列接收任意的两个声源信号，记为第一声源信号和第二声源信号；分别对第一声源信号和第二声源信号进行离散时间傅里叶变换，得到对应的第一声信号和第二声信号；
函数获取模块，用于根据第一声信号和第二声信号，采用广义互相关时延估计算法，得到第一声信号和第二声信号之间广义互相关函数；
全频带矩阵建立模块，用于根据该互相关函数，提取两个声信号的互功率谱，并将其拆分成多个子频带，进而建立全频带互相关观测数据，并基于建立的全频带互相关观测数据，得到全频带矩阵；
映射模块，用于建立全频带矩阵与全频带互相关观测数据之间的映射关系，根据该映射关系，确定均值和协方差矩阵，进而建立全频带矩阵与由均值和协方差矩阵组成的观测数据之间的新映射关系；和
估计模块，用于将该新映射关系带入条件期望函数，采用期望最大化算法，计算每个子频带的估计时延，并进行多次迭代，直至收敛于某一个稳定不变的稳定估计时延，且该稳定估计时延不再增加，则停止迭代；
并将该稳定估计时延作为最终的时延，完成分频带时延估计。
8.一种计算机设备，包括存储器、处理器及存储在所述存储器上并可在所述处理器上运行的计算机程序，其特征在于，所述处理器执行所述计算机程序时实现权利要求1-6中任一项所述的方法。
9.一种计算机可读存储介质，其特征在于，所述计算机可读存储介质存储有计算机程序，所述计算机程序当被处理器执行时使所述处理器执行权利要求1-6中任一项所述的方法。
CN202111274630.4A 2021-10-29 2021-10-29 一种基于期望最大化算法的分频带时延估计方法及其系统 Active CN114035157B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202111274630.4A CN114035157B (zh) 2021-10-29 2021-10-29 一种基于期望最大化算法的分频带时延估计方法及其系统

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202111274630.4A CN114035157B (zh) 2021-10-29 2021-10-29 一种基于期望最大化算法的分频带时延估计方法及其系统

Publications (2)

Publication Number Publication Date
CN114035157A true CN114035157A (zh) 2022-02-11
CN114035157B CN114035157B (zh) 2022-06-14

Family

ID=80142405

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202111274630.4A Active CN114035157B (zh) 2021-10-29 2021-10-29 一种基于期望最大化算法的分频带时延估计方法及其系统

Country Status (1)

Country Link
CN (1) CN114035157B (zh)

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114814728A (zh) * 2022-04-22 2022-07-29 安徽大学 一种声源定位方法、系统、电子设备及介质
CN116660829A (zh) * 2022-11-23 2023-08-29 江南大学 一种声源定位方法及系统
CN117892098A (zh) * 2024-03-15 2024-04-16 江苏翼电科技有限公司 一种应用于高压输电线路作业机器人的时延估计算法
CN121027982A (zh) * 2025-10-27 2025-11-28 杭州电子科技大学 一种波达方向估计方法、系统、设备及介质

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN107479030A (zh) * 2017-07-14 2017-12-15 重庆邮电大学 基于分频和改进的广义互相关双耳时延估计方法
US9989633B1 (en) * 2017-03-15 2018-06-05 Cypress Semiconductor Corporation Estimating angle measurements for source tracking using a phased array system
CN108957403A (zh) * 2018-06-09 2018-12-07 西安电子科技大学 一种基于广义互相关的高斯拟合包络时延估计方法及系统
CN112565119A (zh) * 2020-11-30 2021-03-26 西北工业大学 一种基于时变混合信号盲分离的宽带doa估计方法

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US9989633B1 (en) * 2017-03-15 2018-06-05 Cypress Semiconductor Corporation Estimating angle measurements for source tracking using a phased array system
CN107479030A (zh) * 2017-07-14 2017-12-15 重庆邮电大学 基于分频和改进的广义互相关双耳时延估计方法
CN108957403A (zh) * 2018-06-09 2018-12-07 西安电子科技大学 一种基于广义互相关的高斯拟合包络时延估计方法及系统
CN112565119A (zh) * 2020-11-30 2021-03-26 西北工业大学 一种基于时变混合信号盲分离的宽带doa估计方法

Non-Patent Citations (5)

* Cited by examiner, † Cited by third party
Title
JUN ZHANG等: "High-Resolution DOA Estimation Algorithm for a Single Acoustic Vector Sensor at Low SNR", 《IEEE TRANSACTIONS ON SIGNAL PROCESSING》 *
ZHIHUA LU等: "EM-Based TDOA Estimation of a Speech Source via Gaussian Mixture Models in Noisy and Anechoic Environments", 《 IEEE ACCESS》 *
张亚斌等: "基于高阶累积量的高精度时延估计算法", 《振动与冲击》 *
张鹏等: "加权最大似然波达方向估计算法及其应用研究", 《声学学报》 *
黄澎等: "用于被动声探测的时间延迟估计技术", 《电脑开发与应用》 *

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114814728A (zh) * 2022-04-22 2022-07-29 安徽大学 一种声源定位方法、系统、电子设备及介质
CN116660829A (zh) * 2022-11-23 2023-08-29 江南大学 一种声源定位方法及系统
CN117892098A (zh) * 2024-03-15 2024-04-16 江苏翼电科技有限公司 一种应用于高压输电线路作业机器人的时延估计算法
CN121027982A (zh) * 2025-10-27 2025-11-28 杭州电子科技大学 一种波达方向估计方法、系统、设备及介质

Also Published As

Publication number Publication date
CN114035157B (zh) 2022-06-14

Similar Documents

Publication Publication Date Title
CN114035157B (zh) 一种基于期望最大化算法的分频带时延估计方法及其系统
JP7158806B2 (ja) オーディオ認識方法、ターゲットオーディオを位置決める方法、それらの装置、およびデバイスとコンピュータプログラム
CN106251877B (zh) 语音声源方向估计方法及装置
Erdogan et al. Improved MVDR beamforming using single-channel mask prediction networks.
US9984702B2 (en) Extraction of reverberant sound using microphone arrays
CN114242104B (zh) 语音降噪的方法、装置、设备及存储介质
JP2019503107A (ja) 音響信号を向上させるための音響信号処理装置および方法
CN107703486A (zh) 一种基于卷积神经网络cnn的声源定位算法
Schmalenstroeer et al. Multi-stage coherence drift based sampling rate synchronization for acoustic beamforming
CN110111802B (zh) 基于卡尔曼滤波的自适应去混响方法
CN107290717B (zh) 针对非圆信号的多目标直接定位方法
CN115097378B (zh) 一种基于卷积神经网络的非相干散射源检测与定位方法
CN109188362A (zh) 一种麦克风阵列声源定位信号处理方法
CN111798869A (zh) 一种基于双麦克风阵列的声源定位方法
CN117932319A (zh) 一种适用于宽带线性调频信号的快速无网格波达方向估计方法
Hsiao et al. Super-resolution time-of-arrival estimation using neural networks
CN109658944B (zh) 直升机声信号增强方法及装置
CN105429720B (zh) 基于emd重构的相关时延估计方法
Nesta et al. Enhanced multidimensional spatial functions for unambiguous localization of multiple sparse acoustic sources
Oliinyk et al. Time delay estimation for noise-like signals embedded in non-Gaussian noise using pre-filtering in channels
Fu et al. Blind DOA estimation in a reverberant environment based on hybrid initialized multichannel deep 2-D convolutional NMF with feedback mechanism
CN111968671B (zh) 基于多维特征空间的低空声目标综合识别方法及装置
Li et al. Low complex accurate multi-source RTF estimation
Raikar et al. Effect of Microphone Position Measurement Error on RIR and its Impact on Speech Intelligibility and Quality.
CN114114353A (zh) 卫星信号定位及误差校正方法、系统、定位终端及程序

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