CN105354594B - 一种针对欠定盲源分离的混合矩阵估计方法 - Google Patents

Abstract

本发明涉及一种针对欠定盲源分离的混合矩阵估计方法，对接收到的两路观测信号分别进行短时傅里叶变换得到两路观测信号的短时傅里叶系数X₁(t,f)和X₂(t,f)，形成多个散点，去除低能量的散点；求取剩余散点的比值，根据比值将散点进行分类，得到散点数目最多的M类，计算每一类中散点比值的均值，然后将均值转化成直线的斜率角，再转换成列向量，根据方向接近的列向量得到旋转矩阵T，进行旋转变换，得到新的两路观测信号X′₁(t,f)和X′₂(t,f)；针对两路观测信号X′₁(t,f)和X′₂(t,f)，再次根据两路观测信号的比值将散点进行分类，得到散点数目最多的S类，计算每一类中散点比值的均值，然后得到相近的斜率角，通过这些斜率角得到对应的列向量，最后得到整个混合矩阵。

Description

一种针对欠定盲源分离的混合矩阵估计方法

技术领域

本发明涉及盲信号处理技术领域，尤其涉及一种针对欠定盲源分离的混合矩阵估计方法。

背景技术

如何从传感器接收到的信息中提取出各个单一信息称为盲源分离问题，相应的解决方法被称为盲源分离技术。随着国内外专家和学者的不断研究，盲源分离技术逐渐发展，并且被广泛应用到机械故障检测、通信信号处理、语音信号处理、图像信号处理、生物医学工程等领域。根据源信号的数目和观测信号的数目，盲源分离问题又可以被分为两种：欠定盲源分离问题和非欠定盲源分离问题。如果源信号的数目小于等于观测信号的数目，这时的问题被称为非欠定盲源分离问题；如果源信号的数目大于观测信号的数目，这时的问题就被称为欠定盲源分离问题。在实际环境中，由于受实际条件的限制，传感器数目有限而潜在的源信号的数目是未知的，截获到的观测信号中源信号的数目往往大于观测信号数目，因此，进一步研究欠定盲源分离问题具有更重要的现实意义。

现有技术中，针对欠定盲源分离问题，基于稀疏分量分析的两步法已成为解决该问题的重要手段。两步法具体实施过程为：首先根据观测信号，利用聚类算法估计出混合矩阵，然后根据混合矩阵估计出源信号。稀疏分量分析假设源信号为稀疏信号。稀疏信号在绝大多数时刻取值为零或者较小，在少数时刻取值远离零。当源信号均为稀疏信号时，在大多数采样时刻仅有一个源信号起主导作用，传感器接收到的观测信号会呈现线性特征。当源信号不为稀疏信号时，这种线性特征会不明显。为了得到更好的稀疏特性和更明显的线性特征，通常采用短时傅里叶变换或者小波变换将观测信号变换到变换域。现实存在的各种源信号很难全为稀疏信号，因此需要将观测信号变换到变换域进行处理。根据线性特征，采用聚类算法估计混合矩阵，然后根据混合矩阵和观测信号恢复源信号。因此，混合矩阵的估计在两步法中尤为关键，其估计精度将直接影响源信号的估计精度。

国内外的学者针对欠定盲源分离中的混合矩阵估计问题提出了解决方法。Li Y在Neural computation期刊第16卷第6期上发表的《Analysis of sparse representationand blind source separation》一文中采用了K均值统计聚类方法。该方法复杂度低，容易实现，但性能易受初始值影响,且需给出源信号个数，而实际中源信号个数可能是未知的。张烨在信号处理期刊第25卷第11期上发表的《基于拉普拉斯势函数的欠定盲分离中源数的估计》一文中提出势函数方法，该方法可在源信号个数未知时估计出混合矩阵，但该方法缺乏一定的理论依据，主观经验性太强，且仅适用于二维空间。同时，无论是上述哪种方法，当混合矩阵任意两列向量的方向非常接近，即由于源信号作用导致观测信号形成的直线中有两条直线的斜率十分接近且较小时会导致混合矩阵的估计产生很大的误差。

发明内容

本发明的目的在于提出一种欠定盲源分离的混合矩阵估计方法，能够克服混合矩阵中列向量方向接近的情况下混合矩阵估计精度出现下降的不足，实现在此种情况下精确估计混合矩阵。

实现本发明目的的技术方案：

一种针对欠定盲源分离的混合矩阵估计方法，其特征在于：

步骤1：对接收到的两路观测信号分别进行短时傅里叶变换得到两路观测信号的短时傅里叶系数X₁(t,f)和X₂(t,f)，将X₁(t,f)和X₂(t,f)对应的值分别作横坐标和纵坐标，形成多个散点，t表示观测时刻，f表示频率；

步骤2：针对步骤1形成的多个散点，去除低能量的散点；

步骤3：求取剩余散点的比值X₂(t,f)/X₁(t,f)，根据比值将散点进行分类，得到散点数目最多的M类，M为源信号个数，计算每一类中散点比值的均值R_i(i＝1,...M)，然后利用θ_i＝arctanR_i，arctan表示反正切函数，将均值转化成直线的斜率角，再转换成列向量[cosθ_i,sinθ_i]^T，cos和sin分别表示余弦函数和正弦函数；

步骤4：存储方向不接近的列向量；根据方向接近的列向量得到旋转矩阵T，按照[X'₁(t,f),X'₂(t,f)]^T＝T[X₁(t,f),X₂(t,f)]^T进行旋转变换，得到新的两路观测信号X'₁(t,f)和X'₂(t,f)；

步骤5：针对步骤4获得的两路观测信号X'₁(t,f)和X'₂(t,f)，按照步骤1和步骤2所述方法进行处理，然后，对剩余散点用tanθ₁＜Re{X'₂(t,f)}/Re{X′₁(t,f)}＜tanθ₂去除偏离目标方向的散点，θ₁和θ₂为设定的角度上限值和下限值，Re{}表示取实部；

再次根据两路观测信号的比值将散点进行分类，得到散点数目最多的S类，S为方向不接近的列向量的个数，计算每一类中散点比值的均值Q_i(i＝1,...S)，然后得到相近的斜率角，通过这些斜率角得到对应的列向量，最后得到整个混合矩阵。

步骤3中，对除低能量点以外的散点求取比值，得到比值中的最大值和最小值，再根据最大值和最小值将这些散点分成L类，L>M，得到散点数目最多的M类，M为源信号个数。

步骤4中，通过旋转矩阵T将散点形成的直线旋转到斜率角接近90度的位置。

本发明具有的有益效果：

本发明首先粗略估计出混合矩阵中各个列向量,其次根据拥有相近方向的列向量得到旋转矩阵，然后对原来两路观测信号进行旋转变换，旋转之后再次估计拥有相近方向的列向量，最终估计出整个混合矩阵。本发明通过旋转矩阵T将散点形成的直线旋转到斜率角接近90度的位置，即通过旋转变换提高了方向相近的列向量的估计精度，从而提高了整个混合矩阵的估计精度。

本发明针对混合矩阵中两列或者多列向量方向接近的情况下会导致矩阵估计出现较大误差的事实，提出了一种基于旋转变换的混合矩阵估计方法，适用于混合矩阵中两列(或者多列向量)方向接近的情况，混合矩阵的估计精度得到有效提高。

附图说明

图1是本发明针对欠定盲源分离的混合矩阵估计方法的流程图；

图2是两路观测信号经过短时傅里叶变换后的散点图；

图3是去除低能量散点后两路经过短时傅里叶变换的观测信号散点图；

图4是旋转变换及去除低能量和偏离目标方向散点后的观测信号散点图。

具体实施方式

如图1所示，本发明针对欠定盲源分离的混合矩阵估计方法包括以下步骤：

步骤1：对接收到的两路观测信号分别进行短时傅里叶变换得到两路观测信号的短时傅里叶系数X₁(t,f)和X₂(t,f)(t表示观测时刻，f表示频率)，将X₁(t,f)和X₂(t,f)对应的值分别作横坐标和纵坐标，形成多个散点。

线性瞬时混合情况下的欠定盲分离可以用以下数学模型表示：

x(t)＝As(t) (1)

式中，s(t)＝[s₁(t),s₂(t),......,s_M(t)]^T为M维源信号矢量，x(t)＝[x₁(t),x₂(t),......,x_N(t)]^T为N维观测信号矢量，其中M＞N。t为观测时刻，t＝1,......,n。A为N×M的混合矩阵，可以表示为A＝[a₁,a₂,...,a_M]，其中a_i(i＝1,2,......,M)为A的列向量，因此x(t)还可以表示为

x(t)＝a₁s₁(t)+a₂s₂(t)+...+a_Ms_M(t) (2)

如果在每个观测时刻，最多只有一个源信号的取值较大，这样就可以认为源信号是充分稀疏的，这时观测信号具有线性特征，分布在a_i(i＝1,......,M)的方向上。只有源信号s_i取值较大时，式(2)可以简化为

x(t)＝a_is_i(t) (3)

通常情况下，源信号不够稀疏。为了提高稀疏性，对信号进行变换。本发明针对的是一维信号，采用短时傅里叶变换。对式(2)进行短时傅里叶变换可以得到：

X(t,f)＝a₁(t)S₁(t,f)+a₂(t)S₂(t,f)+...+a_M(t)S_M(t,f) (4)

式中，X(t,f)表示观测信号的傅里叶系数，S_i(t,f)(i＝1,2,......,M)表示第i个源信号的傅里叶系数。假设在某个时频点(t_p,f_p)处只有一个源信号s_i存在，则式(4)可以化简为

X(t_p,f_p)＝a_iS_i(t_p,f_p) (5)

式中，X(t_p,f_p)和S₁(t_p,f_p)都是复数形式。但由式(5)可得

Re{X(t_p,f_p)}＝a_iRe{S_i(t_p,f_p)} (6)

Im{X(t_p,f_p)}＝a_iIm{S_i(t_p,f_p)} (7)

式中，Re{X(t_p,f_p)}和Im{X(t_p,f_p)}分别表示X(t_p,f_p)的实部和虚部。Re{S_i(t_p,f_p)}和Im{S_i(t_p,f_p)}分别表示S_i(t_p,f_p)的实部和虚部。

假设有M个源信号，其中M＞2，通过短时傅里叶变换，可以得到

其中，STFT()表示短时傅里叶变换。当变换后的信号足够稀疏时，在时频点(t_p,f_p)处会得到以下式子

两路观测信号通过短时傅里叶变换之后对应的散点，如图2所示。

步骤2：针对步骤1形成的多个散点，

为了去除低能量的散点，设定一个阈值λ，满足下面条件的点被认为是低能量点从而被去除

λ的值越大时，留下的点就越少；λ的值越小时，留下的点就越多。λ的具体取值需根据实际情况而定。去除相应的低能量点后，得到的两路短时傅里叶变换系数形成的散点，如图3所示。

步骤3：求取剩余散点的比值X₂(t,f)/X₁(t,f)，根据这些比值将散点进行分类，得到散点数目最多的M类(M为源信号个数)，计算每一类中散点比值的均值R_i(i＝1,...M)，然后利用θ_i＝arctanR_i(arctan表示反正切函数)将均值转化成直线的斜率角，再转换成列向量[cosθ_i,sinθ_i]^T(cos和sin分别表示余弦函数和正弦函数)；

定义以下的式子：

D(t,f)为X₂(t,f)和X₁(t,f)的比值。根据式(11)求出这些比值的最大值和最小值

D_max＝max{D(t,f)} (12)

D_min＝min{D(t,f)} (13)

又因为

R_i＝tanθ_i(i＝1,2,...,n) (15)

这样，列向量[cosθ_i,sinθ_i]^T就可以被得到。

步骤4：存储不相近斜率角所形成的列向量，即存储方向不接近的列向量，根据方向接近的列向量得到旋转矩阵T，按照[X'₁(t,f),X'₂(t,f)]^T＝T[X₁(t,f),X₂(t,f)]^T进行旋转变换，得到新的两路观测信号X'₁(t,f)和X'₂(t,f)。

假设θ_i中有两个斜率角θ_x和θ_y接近且满足θ_x＜θ_y＜90°，首先得到不接近的θ_i的估计值但是上述估计值中不包括θ_x和θ_y的估计值和即i≠x,i≠y。

对于θ_x和θ_y的估计，定义旋转角度为

其中，INT()为取整函数。由旋转角度得到旋转矩阵

通过旋转变换得到新的观测信号矢量X'(t,f)为

然后得到：

由于

并且因为可以得到

原本由于θ_x和θ_y比较接近，导致tanθ_x和tanθ_y的差值十分接近，旋转变换以后，和的差值增大。

步骤5：针对步骤4获得的两路观测信号X'₁(t,f)和X'₂(t,f)，按照步骤1和步骤2所述方法进行处理，然后，对剩余散点用tanθ₁＜Re{X'₂(t,f)}/Re{X'₁(t,f)}＜tanθ₂(θ₁和θ₂为设定的角度上限值和下限值，Re{}表示取实部)去除偏离目标方向的散点；再次根据两路观测信号的比值将散点进行分类，得到散点数目最多的S类(S为方向不接近的列向量的个数)，计算每一类中散点比值的均值Q_i(i＝1,...S)，然后得到相近的斜率角，通过这些斜率角得到对应的列向量，最后得到整个混合矩阵。

按照步骤1和步骤2进行处理，设剩下的散点为(X₁(t_i,f_i),X₂(t_i,f_i))(i＝1,...,r.r＜n)，由于经过旋转变换以后原来斜率角接近的散点形成的直线旋转到了接近90^°的位置，这样可以设定阈值，来去除偏离目标方向的点，判定规则为

其中，θ₁和θ₂的设定满足以下两个条件

将满足以上判定规则的散点留下，剩余散点如图4所示。对剩余散点的比值进行分类，得到最大的S类(S为方向不接近的列向量的个数)，计算每一类中散点比值的均值，转换成斜率值，得到两个旋转之后的估计角θ'_x和θ'_y，进而得到θ_x和θ_y的估计值和

结合存储的方向不接近的列向量，可以得到最终的混合矩阵

下面结合具体实施例进一步说明本发明的有益效果。

仿真参数具体设置如下：四个线性调频信号分别被转成中频为1.2MHz，0.8MHz，0.6MHz和0.6MHz的信号，调频斜率分别为-30MHz/s，-25MHz/s，20MHz/s和25MHz/s。抽样频率为2.5MSps，抽样点数为20000，进行STFT变换时需要设置的参数如下：每次变换的频率点数为1024，重叠部分的点数为128，加窗函数选择汉宁窗，源信号数目M＝4，N＝2，两个相近的角分别为θ₁＝24°，θ₂＝25°，两个不相近的角分别为θ₃＝50°，θ₄＝70°。

为了评价本发明的有效性，定义均方误差(MSE)作为评价指标，公式如下：

式中，θ_j表示第j个角，表示第j个角第i次的估计结果，Num表示总共的试验次数。测试在不同信噪比下两个相近角在旋转前后1000次试验的MSE，计算结果如表1所示。

表1 两个相近角旋转前后的MSE

表1的计算结果表明：通过旋转变换以后的两个相近角的估计精度得到了很大的提高，并且在信噪比降低时，这种提高会更加明显。

Claims (4)

1.一种针对欠定盲源分离的混合矩阵估计方法，其特征在于：

步骤1：对接收到的两路观测信号分别进行短时傅里叶变换得到两路观测信号的短时傅里叶系数X₁(t,f)和X₂(t,f)，将X₁(t,f)和X₂(t,f)对应的值分别作横坐标和纵坐标，形成多个散点，t表示观测时刻，f表示频率；

步骤2：针对步骤1形成的多个散点，去除低能量的散点；

步骤3：求取剩余散点的比值X₂(t,f)/X₁(t,f)，根据比值将散点进行分类，得到散点数目最多的M类，M为源信号个数，计算每一类中散点比值的均值R_i,i＝1,...M，然后利用θ_i＝arctanR_i，arctan表示反正切函数，将均值转化成直线的斜率角，再转换成列向量[cosθ_i,sinθ_i]^T，cos和sin分别表示余弦函数和正弦函数；

步骤4：存储方向不接近的列向量；根据方向接近的列向量得到旋转矩阵T，按照[X₁'(t,f),X'₂(t,f)]^T＝T[X₁(t,f),X₂(t,f)]^T进行旋转变换，得到新的两路观测信号X₁'(t,f)和X'₂(t,f)；

步骤5：针对步骤4获得的两路观测信号X₁'(t,f)和X'₂(t,f)，按照步骤1和步骤2所述方法进行处理，然后，对剩余散点用tanθ₁＜Re{X'₂(t,f)}/Re{X₁'(t,f)}＜tanθ₂去除偏离目标方向的散点，θ₁和θ₂为设定的角度上限值和下限值，Re{}表示取实部；

再次根据两路观测信号的比值将散点进行分类，得到散点数目最多的S类，S为方向不接近的列向量的个数，计算每一类中散点比值的均值Q_i,i＝1,...S，然后得到相近的斜率角，通过这些斜率角得到对应的列向量，最后得到整个混合矩阵。

2.根据权利要求1所述的针对欠定盲源分离的混合矩阵估计方法，其特征在于：步骤3中，对除低能量点以外的散点求取比值，得到比值中的最大值和最小值，再根据最大值和最小值将这些散点分成L类，L>M，得到散点数目最多的M类，M为源信号个数。

3.根据权利要求2所述的针对欠定盲源分离的混合矩阵估计方法，其特征在于：步骤4中，通过旋转矩阵T将散点形成的直线旋转到斜率角接近90度的位置。

4.根据权利要求1至3任何一项所述的针对欠定盲源分离的混合矩阵估计方法，其特征在于：步骤2中，针对步骤1形成的多个散点，利用去除低能量的散点，λ为0与1之间的值，max{}表示取最大值。