基于混合信号局部峰值方差检测的盲源分离方法 【技术领域】
本发明涉及盲源分离技术领域，具体涉及时频分量分析中对DUET盲源分离方法局部信号峰值检测方法的改进。
背景技术
盲源分离(Blind Source Separation，BSS)是指在不知道源信号和传输信道的先验信息的情况下，根据输入源信号的统计特性，仅由传感器观测到的混叠信号来恢复或分离出源信号的过程。由于其特殊的性质，在无线通信、阵列信号处理、生物医学信号处理、语音信号处理及图像信号处理等许多领域都具有广阔的实际应用前景。
目前独立分量分析(Independent Component Analysis，ICA)方法被应用于大多数盲源分离问题，但是这种方法不能解决观测信号个数少于源信号个数的欠定盲源分离问题(Underdetermined Blind Source Separation，UBSS)。随着研究的深入，出现了稀疏分量分析以及时频分量分析方法来解决欠定盲分离问题，而Scott Rickard等提出的DUET(Degenerate Unmixing Estimation Technique，DUET)盲源分离方法是时频分量分析中具有代表性的方法，它可以在时频域中利用信号的W-错位正交性(W-disjoint orthogonal，W-DO)并通过时频掩码对信号进行分离，对于分离实际信号具有重要意义。
如图1所示，DUET盲源分离方法所处理的盲源分离的基本模型如下：
公式二x1(t)=Σj=1Psj(t)x2(t)=Σj=1Pajsj(t-δj)]]>
其中，t＝1，2，...，T表示离散时刻，P为源信号的个数，x1(t)和x2(t)分别表示两个传感器接收到的观测信号，sj(t)表示第j个源信号，xi(t)和vi(t)表示第i个传感器的观测信号和噪声信号，DUET盲源分离方法进行盲源分离时，将噪声信号忽略，故公式二的表达式中未列出噪声信号，aj∈R+和δj∈R分别为第j个源信号到第二个传感器的相对衰减和相对延迟，表示第j个源信号的估计信号。
DUET盲源分离方法步骤如下：
1、方法假设
DUET盲源分离方法里面假设信号之间是满足W-错位正交性，即各个源信号在时频域中表现为每个时频点至多有一个源信号占优。信号从时域到时频域往往通过短时傅里叶变换(Short-Time FouriesTransform，STFT)实现。
信号的短时傅里叶变换定义为：
SW(ω,τ)=∫-∞∞W(t-τ)s(t)e-iωtdt]]>
其中W(t)为窗函数，公式二在时频域中可以写成如下形式：
X1(ω)X2(ω)=1...1a1e-iωδ1...aPe-iωδPS1(ω)...SP(ω)]]>
所以在时频域中W-错位正交性可以表示为：
SiW(ω,τ)SjW(ω,τ)=0,∀i≠∀j,∀ω,τ]]>
从上式中可以看出，在任意给定的一个时频点，至多只有一个信号是占优的。于是，在某些时频点，有下式成立：
X1(ω)X2(ω)=1aie-iωδiSi(ω)]]>
2、混合矩阵参数估计
利用上式，我们可以通过下式计算一个信号的相对衰减和延迟参数：
(ai,δi)=(||X2(ω)X1(ω)||,Im(log(X2(ω)X1(ω)))/ω)]]>
如果W是一个有限长的窗函数，那么上式可以写成：
(ai(ω,τ),δi(ω,τ))=(||X2W(ω,τ)X1W(ω,τ)||,Im(log(X2W(ω,τ)X1W(ω,τ)))/ω)]]>
其中τ为短时傅里叶变换的窗长度。因此，通过上式我们可以计算(a(ω，τ)，δ(ω，τ))在所有的(ω，τ)，这样，可以画出信号衰减-延迟的直方图，如果信号之间满足W-错位正交性，那么，在直方图上，衰减和延迟将会表现为峰值的形式，直方图上的峰值则为各观测信号短时傅里叶变换模的乘积和频率的加权值，如图1所示。
3、源信号估计
通过在直方图上识别出各个峰值，检测到的横纵坐标即为混合矩阵参数。然后利用二进制时频掩码，在时频域中提取出各个源信号，最后反变换到时域得到需要的源信号。
在DUET盲源分离方法的峰值检测部分，比较有代表性的是标准K-means方法，加权K-means方法以及人眼识别来寻找。标准K-means方法没有考虑到各个聚类变量对聚类结果的贡献程度的差异，因此聚类结果精度较差；加权K-means方法是对标准K-means方法的一个改进，但是对于实录语音信号，时频域上的聚类点并不明显，聚类效果也不精确；通过人眼在直方图上识别峰值，显然，人眼精度也不够，且对于盲源分离方法的智能性有严重的影响，对于一个实际的语音盲分离系统来说是不可行的，实用价值小。
【发明内容】
为了解决目前DUET盲源分离方法不能有效地检测信号衰减-延迟直方图中局部峰值、对混合矩阵参数估计不精确，导致源信号估计精度不高的问题，提出了一种基于混合信号局部峰值方差检测的盲源分离方法(Sub-areavariance peak detecting，SV-peakdet)。
基于混合信号局部峰值方差检测的盲源分离方法是基于以下系统实现的：盲源分离系统由混合系统和解混系统组成，源信号sj(t)与噪声信号vi(t)在混合系统中混合，通过混合系统的观测信号为xi(t)，观测信号xi(t)经过解混系统分离出源信号的估计信号为
基于混合信号局部峰值方差检测的盲源分离方法的具体步骤如下：
步骤A、通过观测信号xi(t)的相对衰减和延迟参数画出信号衰减-延迟直方图；
步骤B、在步骤A获得的信号衰减-延迟直方图平面区域Ω内遍历x轴和y轴，找出所有子区域Ωi，j∈Ω，所述子区域Ωi，j为公式一所示的N×N的网格区域，N为奇数，其中a(i+m，j+n)为观测信号xi(t)短时傅里叶变换模的乘积与频率的加权值，也为信号衰减-延迟直方图的幅值，其中其中I、J分别为平面区域Ω的横坐标最大值和纵坐标最大值；
公式一a(i-N-12,j+N-12)...a(i-1,j+N-12)a(i,j+N-12)a(i+1,j+N-12)...a(i+N-12,j+N-12)...............a(i-N-12,j+1)...a(i-1,j+1)a(i,j+1)a(i+1,j+1)...a(i+N-12,j+1)a(i-N-12,j)...a(i-1,j)a(i,j)a(i+1,j)...a(i+N-12,j)a(i-N-12,j-1)...a(i-1,j-1)a(i,j-1)a(i+1,j-1)...a(i+N-12,j-1)...............a(i-N-12,j-N-12)...a(i-1,j-N-12)a(i,j-N-12)a(i+1,j-N-12)...a(i+N-12,j-N-12)]]>
步骤C、从步骤B获得地所有子区域Ωi，j中选出满足a(i，j)≥a(i+m，j+n)的多个子区域作为峰值子区域其中a(i，j)是子区域Ωi，j中心点的幅值，
m,n∈{-N-12,...,-1,0,1,...,N-12};]]>
步骤D、分别计算获得步骤C中每一个峰值子区域的方差VARk，每个峰值子区域的方差VARk的计算过程为：求取峰值子区域中所有数据点Xk的三维坐标的平均值和所述平均值的三维坐标确定的点为平均点
X‾=1N2Σk=1N2Xk]]>
依次求得峰值子区域内的每个数据点Xk到平均点的距离
d2(Xk,X‾)=|Xk-X‾|=(xk-x‾)2+(yk-y‾)2+(ak-a‾)2]]>
其中Xk＝(xk，yk，ak)，k∈{1，2，…，N2}；
计算峰值子区域内的所有数据点Xk三维坐标与平均点三维坐标的方差，即获得该峰值子区域的方差VARk
VARk=1N2Σk=1N2d2(Xk-X‾)]]>
步骤E、对步骤D中取得的所有峰值子区域的方差VARk进行从大到小排序，提取前P个方差VARk，P为源信号的个数，也是所需要的峰值的个数；
步骤F、将步骤E获得的P个方差VARk对应的横、纵坐标分别呈递给衰减矩阵和延迟矩阵，然后通过二进制时频掩码来提取峰值，在时频域上分离源信号，最后反变换到时域得到最终分离出的源信号的估计信号
本发明解决了目前DUET盲源分离方法中不能有效地进行信号衰减-延迟直方图中局部峰值的检测问题，本发明是根据信号衰减-延迟直方图的特点，首先进行局部峰值检测，然后计算每个子区域的方差大小，最后根据源信号的数目确定出需要的子区域，子区域峰值点的横纵坐标就是衰减-延迟参数。本发明能够有效而准确的估计所有峰值点，并且筛选出需要的峰值点，在DUET盲源分离方法的峰值检测部分有很好的估计能力，且对于最后源信号估计精度的提高有一定作用。本发明适用于DUET盲源分离方法的峰值检测，同样适用于一般的峰值检验。
【附图说明】
图1为盲源分离的基本模型示意图。图2为源信号衰减-延迟直方图。图3为语音盲分离的实验系统原理框图。图4为实录语音实验环境示意图。
【具体实施方式】
具体实施方式一：结合图1说明本实施方式，基于混合信号局部峰值方差检测的盲源分离方法是基于以下系统实现的：盲源分离系统由混合系统和解混系统组成，源信号sj(t)与噪声信号vi(t)在混合系统中混合，通过混合系统的观测信号为xi(t)，观测信号xi(t)经过解混系统分离出源信号的估计信号为
基于混合信号局部峰值方差检测的盲源分离方法的具体步骤如下：
步骤A、通过观测信号xi(t)的相对衰减和延迟参数画出信号衰减-延迟直方图；
步骤B、在步骤A获得的信号衰减-延迟直方图平面区域Ω内遍历x轴和y轴，找出所有子区域Ωi，j∈Ω，所述子区域Ωi，j为公式一所示的N×N的网格区域，N为奇数，其中a(i+m，j+n)为观测信号xi(t)短时傅里叶变换模的乘积与频率的加权值，也为信号衰减-延迟直方图的幅值，其中其中I、J分别为平面区域Ω的横坐标最大值和纵坐标最大值；
公式一a(i-N-12,j+N-12)...a(i-1,j+N-12)a(i,j+N-12)a(i+1,j+N-12)...a(i+N-12,j+N-12)...............a(i-N-12,j+1)...a(i-1,j+1)a(i,j+1)a(i+1,j+1)...a(i+N-12,j+1)a(i-N-12,j)...a(i-1,j)a(i,j)a(i+1,j)...a(i+N-12,j)a(i-N-12,j-1)...a(i-1,j-1)a(i,j-1)a(i+1,j-1)...a(i+N-12,j-1)...............a(i-N-12,j-N-12)...a(i-1,j-N-12)a(i,j-N-12)a(i+1,j-N-12)...a(i+N-12,j-N-12)]]>
步骤C、从步骤B获得的所有子区域Ωi，j中选出满足a(i，j)≥a(i+m，j+n)的多个子区域作为峰值子区域其中a(i，j)是子区域Ωi，j中心点的幅值，
m,n∈{-N-12,...,-1,0,1,...,N-12};]]>
步骤D、分别计算获得步骤C中每一个峰值子区域的方差VARk，每个峰值子区域的方差VARk的计算过程为：求取峰值子区域中所有数据点Xk的三维坐标的平均值和所述平均值的三维坐标确定的点为平均点
X‾=1N2Σk=1N2Xk]]>
依次求得峰值子区域内的每个数据点Xk到平均点的距离
d2(Xk,X‾)=|Xk-X‾|=(xk-x‾)2+(yk-y‾)2+(ak-a‾)2]]>
其中Xk＝(xk，yk，ak)，k∈{1，2，…，N2}；
计算峰值子区域内的所有数据点Xk三维坐标与平均点三维坐标的方差，即获得该峰值子区域的方差VARk
VARk=1N2Σk=1N2d2(Xk-X‾)]]>
步骤E、对步骤D中取得的所有峰值子区域的方差VARk进行从大到小排序，提取前P个方差VARk，P为源信号的个数，也是所需要的峰值的个数；
步骤F、将步骤E获得的P个方差VARk对应的横、纵坐标分别呈递给衰减矩阵和延迟矩阵，然后通过二进制时频掩码来提取峰值，在时频域上分离源信号，最后反变换到时域得到最终分离出的源信号的估计信号
步骤B中的划分子区域Ωi，j的大小可以根据实际情况调整，对于子区域Ωi，j的寻找，通过i、j两重循环实现，遍历x轴和y轴。
当步骤B中的N＝5时，则所述Ωi，j为公式六所示的5×5的网格区域。
公式六a(i-2,j+2)a(i-1,j+2)a(i,j+2)a(i+1,j+2)a(i+2,j+2)a(i-2,j+1)a(i-1,j+1)a(i,j+1)a(i+1,j+1)a(i+2,j+1)a(i-2,j)a(i-1,j)a(i,j)a(i+1,j)a(i+2,j)a(i-2,j-1)a(i-1,j-1)a(i,j-1)a(i+1,j-1)a(i+2,j-1)a(i-2,j-2)a(i-1,j-2)a(i,j-2)a(i+1,j-2)a(i+2j-2)]]>
步骤F属于原DUET盲源分离方法内容，在这里陈述是为体现本发明与DUET盲源分离方法衔接的完整性。
具体实施方式二：本实施方式对具体实施方式一中步骤C的进一步限定，在步骤C中设定一个阈值门限δ，对峰值子区域进一步筛选，选取中心点的值a(i，j)大于阈值门限δ的峰值子区域参与后续计算，即满足a(i，j)≥δ的峰值子区域参与后续计算。
本实施方式中最大值的阈值门限δ是根据信号的实际需求由技术人员人为调整设定的，一般该阈值门限δ与信号衰减-延迟直方图中的最大幅值M成正比。例如，本实施方式中的最大值的阈值门限δ可以选择
本实施方式中，根据设置的阈值门限δ，对峰值子区域做了进一步筛选，只有中心点的值a(i，j)是本子区域中的最大值、并且大于阈值门限δ的子区域才能够作为峰值子区域参与后续计算，其余的子区域舍去，减少了后续的计算量，提高了本发明的运算效率。
具体实施方式三：本实施方式是对具体实施方式一中步骤E的进一步限定，步骤E中对所有峰值子区域的方差VARk进行排序的方法采用“冒泡法”。
具体实施方式四：结合图2、图3和图4说明本实施方式，本实施方式是将本发明方法和现有方法分别应用于语音盲源分离的仿真实验和实录实验，并将本发明方法和现有方法的平均信号失真比进行对比。
实验环境参见图4所示，是一个长13m、宽5m、高3m的立方体空旷房间，定义所述房间中的一个宽5m的侧墙为基准墙面，定义同时垂直于所述基准墙面和房间地面、并穿过所述基准墙面的对称中心线的平面为房间对称面。实验设备包括两个麦克1、三个声源2，所述实验设备放置于房间内，并且以所述对称面为中心左右对称，所述两个麦克1、三个声源2的中心位于同一个平面内，其中两个麦克1间隔4cm分别放在对称面的两侧，麦克1距离基准墙壁4m，一个声源2放置在对称面上，该声源2的发声面朝向麦克1，该声源2距离基准墙壁5m，另两个声源2分别位于对称面的两侧，所述两个声源2的发声面朝向基准墙壁，其中一个声源2与两个麦克1的拾音器连线的中点的连线与对称面成45°夹角，并且所述三个声源2的中心都位于以两个麦克1的拾音器连线的中点为中心、半径为1m的圆弧上。
图2为本实施方式采用的源信号的衰减-延迟直方图。
实验过程中的信号处理流程参见图3所示：三个声源2发出的混合语音信号经过麦克进行采集，然后麦克将采集到的音频信号发送给音频采集卡进行识别，所述音频采集卡将转换后的信号传递给计算机进行分析处理，获得分离出的信号输出。其中，计算机中的信号分析处理分别采用现有的基于标准K-means的DUET盲源分离方法(以下简称标准K-means方法)、基于加权K-means的DUET盲源分离方法(以下简称加权K-means方法)和本发明的方法。
本实验过程中，估计精度用平均信号失真比来衡量。由于峰值检测是对混合参数的估计，直接影响到最后源信号的估计，故用指标能够衡量本发明的优劣。
表1给出了在上述实验环境中，分别采用现有的标准K-means方法、加权K-means方法和本发明的方法在语音盲源分离的仿真实验中源信号的估计精度，其中每种方法运行20次。表中用“KmPs”表示一种混合情况，例如，“2m3s”表示3个源信号、2个观测信号的情况。由表1可看出，本发明比以往的方法对源信号的估计精度有所提高。
表1各盲源分离方法分离仿真实验结果
表2给出了在上述实验环境中，分别采用现有的标准K-means方法、加权K-means方法和本发明的方法在分离实录语音实验中源信号的估计精度，其中每种方法运行20次。从表2中可以看出，对于实录语音的分离，本发明方法要明显优于以往的方法，这是由于实录语音时存在环境回声及噪声的影响，在衰减-延迟图上数据聚类并不像仿真实验那么明显，基于聚类的盲源分离方法在估计聚类中心时与实际的衰减-延迟矩阵会有所偏差，所以最后本发明方法得到的都要优于标准K-means方法和加权K-means方法，且稳定性很好，基于聚类的盲源分离方法在估计聚类中心时和初始值有关，所以稳定性方面不如本发明方法。
表2各盲源分离方法分离实录语音结果