CN108108666B - 一种基于小波分析和时频单源检测的混合矩阵估计方法 - Google Patents

Abstract

本发明属于阵列信号处理技术领域，公开的一种基于小波分析和时频单源检测的混合矩阵估计方法，通过小波分析方法对原观测信号进行小波分解、重构实现原始观测信号的扩维，然后对扩维后的观测信号矩阵进行奇异值分解，基于计算得到的奇异值向量，得到时频单源点检测的判断阈值计算方法，然后计算得到时频单源点集合，基于时频单源点集合应用K均值聚类法进行混合矩阵估计；将多次估计得到的混合矩阵组合成为一个矩阵，选择出现最频繁的n个列向量组合成为新的混合矩阵，作为最终估计得到的混合矩阵。本发明能够较好地克服K均值聚类法在混合矩阵估计中易受初始初始聚类中心选择影响从而估计精度不高的问题，得到混合矩阵估计的优化结果。

Description

一种基于小波分析和时频单源检测的混合矩阵估计方法

技术领域

本发明属于阵列信号处理技术领域，涉及一种基于小波分析和时频单源检测的混合矩阵估计方法。

背景技术

盲源分离(Blind Source Separation,BSS)是指在辐射源信号和传播信道信息未知的情况下，仅根据观测信号恢复得到辐射源信号的技术^[1]。鉴于其独有的优势，盲源分离技术引起了国内外信号处理界学者的广泛关注^[2]，在语音分离、通信、生物医学、图像恢复等领域具有重要的应用价值。

当天线阵元数量少于辐射源信号数量时，BSS问题称为欠定盲源分离^[3]。目前，稀疏分量分析是解决欠定盲源分离问题的主要方法，它通过利用源信号的稀疏性使分离后的源信号尽可能唯一。绝大多数欠定盲源分离算法都采用“两步法”的思路，即先估计出混合矩阵，然后结合估得的混合矩阵实现源信号的分离，因而混合矩阵估计是最终实现信号分离的前提和基础。其中信号稀疏性：在绝大多数采样点上为零或者较小，在少数采样点幅值较大的信号，非零值采样点越少，信号越稀疏^[1]。

与本发明相关的现有技术，现阶段，混合矩阵估计方法大致可以分为三类^[4-7]：一是基于超完备稀疏表示的估计方法，需要假设源信号满足特定的概率分布，通过自适应迭代把混合矩阵和源信号同时估计出来。这类方法计算复杂，并且当源信号为非充分稀疏信号时，会收敛到局部极值点；二是基于张量分解的估计算法，适用于源信号相互独立并且非高斯的情况，这些方法不需要设置初始参数，适应能力强，估计精度高。但是，目前已有的方法计算复杂，且只利用了源信号的独立性，对于假设条件要求较高；三是基于源信号稀疏性和聚类的混合矩阵估计方法，计算相对简单，当源信号能满足算法对应的稀疏性条件时，混合矩阵的估计精度较高，是混合矩阵估计的主要方法。

基于源信号的稀疏性，S.G.Kim等^[3]提出了一种基于时频单源点检测(Single-Source-Detection,SSD)的混合矩阵估计算法，但是没有考虑混合信号中叠加加性噪声的一般情况下时频单源点判断阈值的确定方法，影响了算法的应用范围和混合矩阵估计精度。实际上，混叠加性噪声的混合信号盲源分离问题与实际更具有一致性，因而成为该领域的研究热点和难点之一^[8-9]。算法检测得到时频单源点集合之后，应用K均值聚类法估计混合矩阵，而K均值聚类法应用效果严重依赖于初始聚类中心的选择，导致混合矩阵估计精度不稳定。参考文献如下：

[1]何昭水,谢胜利,傅予力.信号的稀疏性分析[J].自然科学进展,2006,16(9):1167-1173.

[2]张贤达.现代信号处理(第三版)[M].北京:清华大学出版社,2015.

[3]S.G.Kim and C.D.Yoo.Underdetermined Blind Source Separation Basedon Subspace Representation[J].IEEE Trans.Signal Process.,2009,57(7):2604-2614.

[4]陆凤波.复杂电磁环境下的欠定盲源分离技术研究[D].长沙:国防科学技术大学博士论文,2011.F.B.Lu.Research on Underdetermined Blind Source Separation inComplex Electromagnetic Environment[D].Changsha:doctoral thesis of NationalUniversity of Defense Technology,2011.

[5]D.Z.Peng andY.Xiang.Underdetermined Blind Source Separation Basedon Relaxed Sparsity Condition of Sources[J].IEEE Trans.Signal Process.,2009,57(2):809-814.

[6]Y.Zhang,H.L.Wang,W.W.Wang,et al.K-Plane Clustering Algorithm forAnalysis Dictionary Learning[J].2013 IEEE International Workshop on MachineLearning for Signal Processing,Sept.22-25,2013,Southampton,UK.

[7]谭北海,赵敏,谢胜利.基于无监督学习的源数估计及盲分离算法[J].系统工程与电子技术,2009,31(8):1790-1794.

[8]吴微,彭华,周正康.一种改进的FastICA算法及其在含噪盲源分离中的应用[J].信息工程大学学报,2013,14(6):708-712.

[9]张晗博,殷奕,殷奎喜.基于小波变换的非平稳信号分析与处理[J].南京师范大学学报,2014,14(1):63-69。

发明内容

为克服现有技术的不足，本发明提供一种基于小波分析和时频单源检测的混合矩阵估计方法，在观测信号叠加加性噪声的普遍情况下，基于小波分析方法提出时频单源点检测判断阈值确定方法，针对K均值聚类法存在的问题，提出基于列向量夹角检测的混合矩阵优化方法，使得时频单源点检测精度和效率得到提高，混合矩阵估计效果改善。

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

一种基于小波分析和时频单源检测的混合矩阵估计方法，在观测信号混叠加性噪声的情况下，通过小波分析方法对原观测信号进行小波分解、重构实现原始观测信号的扩维，然后对扩维后的观测信号矩阵进行奇异值分解，基于计算得到的奇异值向量，得到时频单源点检测的判断阈值计算方法，然后计算得到时频单源点集合，基于时频单源点集合应用K均值聚类法进行混合矩阵估计；重复前述混合矩阵估计步骤，得到多次混合矩阵估计结果，将多次估计得到的混合矩阵组合成为一个矩阵，对该矩阵进行列向量夹角检测，设辐射源数为n，选择出现最频繁的n个列向量组合成为新的混合矩阵，作为最终估计得到的混合矩阵，具体实现步骤如下：

(1)基于小波分析的时频单源点判断阈值确定方法

为了解决在观测信号混叠加性噪声情况下，时频单源点检测中判断阈值选取缺乏理论指导的问题，本发明基于小波分析原理，得到一种阈值的选取方法，步骤如下：

1)对原始观测数据进行白化处理；

2)对白化后的观测信号X(t)的小波分解，设小波分解层数为I，x₁(t)分解得到的多路小波系数分别为a_1j，d_1j，j＝1,2…I，上述小波系数的下标中，第一个“1”表示第一个原始观测信号x₁(t)，其中a_1j表示第1个观测信号的低频小波系数，d_1j表示第1个观测信号的高频小波系数；其他观测信号小波分解方法按照上述方法进行操作；

3)小波系数重构，将m个观测信号分解得到的小波系数序列信号分别进行小波重构，以第一个观测信号x₁(t)分解得到的多路小波系数a_1j，d_1j，j＝1,2…I为例，将低频小波系数重构为x₁₀(t)，每一层高频小波系数重构为x_1j(t)，j＝1,2…I，其他观测信号分解得到的小波系数重构方法及表示方法与此相同；

4)将白化后的观测信号和其小波重构信号组合成为新的多维观测信号，以第一路观测信号x₁(t)为例说明其组合方法。选取x₁(t)，x₁₀(t)，x_1j(t)，j＝1,2…I作为第一路观测信号扩维后的信号，x₁₀(t)为由x₁(t)进行I层小波系数分解后得到的低频系数重构的信号，x_1j(t)为第j层高频系数重构后的信号。其他观测信号扩维的方法与上述相同。经过扩维后的观测信号表示为：

x_w(t)＝[x₁(t),x₁₀(t),x_1j(t),x₂(t),x₂₀(t),x_2j(t),…,x_m(t),x_m0(t),x_mj(t),]^T,j＝1,2…I (1)

5)新的多维观测信号X_w(t)的去均值化处理，表示为

6)计算多维观测信号

的协方差矩阵R_w，

T₀为观测信号采样点数。

7)协方差矩阵R_w的奇异值分解。计算矩阵R_w的奇异值分解

得到其奇异值Λ_s＝{λ₁,λ₂, … λ_M}，其中M＝m×(I+2)为升维后的观测信号的维数。

8)得到时频单源点检测判断阈值。时频单源点判断阈值ε的取值范围为ε＝[λ_m+1,λ_m+2]。

以上即为基于小波分析的时频单源点检测阈值确定方法，该方法考虑到观测信号中叠加的加性噪声；

(2)基于时频单源点检测的混合矩阵估计

对于第i个辐射源信号s_i(n)进行短时傅里叶变换(Short-time Fouriertransform,STFT)：

其中，h(n)是窗函数，τ是STFT的加窗函数中心，k是频率点；令X(τ,k)＝[X₁(τ,k),X₂(τ,k),…X_m(τ,k)]^T，S(τ,k)＝[S₁(τ,k),S₂(τ,k),…S_m(τ,k)]^T分别表示观测信号X和源信号S的STFT变换；

单源点检测法基于观测信号时频变换比率寻找一个时频点集合S，在该集合中每个时频点只有一个辐射源信号是活跃的，基于时频单源点检测的混合矩阵估计算法核心是：对于只有单源活跃的时频点，观测信号的时频域比率是实数，而对于多源活跃的时频点，观测信号的时频域比率为复数；基于观测信号时频变换的比率，检测只有单源活跃的时频点，即可估计得到混合矩阵，算法步骤如下：

1)给定一个时频单源点检测判断阈值ε＞0，检测得到只有单源活跃的时频点集合S_s，其中观测信号STFT变换的实部和虚部都有足够的能量：

在式(3)中，Im(x)表示x的虚部，应用其他观测信号的STFT代替X₁(τ,k)作为分母；

在时频点(τ_s,k_s)，如果满足式(4)，则(τ_s,k_s)∈S_s，

其中，a_li为混合矩阵A的第(l,i)个元素。由式(4)可判断在时频点(τ_s,k_s)处只存在第i个源信号。相反，如果第i个和第j个源信号S_i(τ,k)和S_j(τ,k)在时频点(τ_n,k_n)处同时存在，则时频单源点集合S_s中将不包括时频点(τ_n,k_n)，因为：

式(5)表示S_i(τ,k)和S_j(τ,k)为复数，不能如式(4)一样从分子和分母中消去；

2)应用聚类算法，基于式(6)所示比率向量，将时频单源点集合S_s聚成m类

其中，Re(x)表示x的实部；应用k均值聚类算法对S_s进行聚类；在S_s中，若第i个源信号的时频点活跃，则有以下比率向量：

基于上述比率向量的时频单源点集合表示为

应用上述算法，检测得到每一个辐射源信号对应的单源点区域；

3)基于所有辐射源信号的时频单源点集合估计得到混合矩阵A，在混合矩阵估计之前，所有观测信号的实部和虚部都通过l₂范数进行归一化；

基于第i个源信号活跃的时频单源点集合

通过式(8)估计得到混合矩阵的第i列

(3)基于向量夹角检测的混合矩阵优化方法

检测得到时频单源点集合以后，应用K均值聚类算法进行混合矩阵估计，K均值聚类算法具有计算量小，收敛快的优势，但是聚类效果严重依赖于初始聚类中心的选择；对于相同的观测信号，即便时频单源点检测精度很高，由于K均值聚类法每次自动选择的初始聚类中心是不确定的，导致每次估计得到的混合矩阵相应地是变化的，每次估计得到的混合矩阵可能只有一部分或全部列向量与真正的混合矩阵相似；由此分析可知，基于时频单源点检测的混合矩阵估计效果受K均值聚类算法影响很大；

为了提高混合矩阵估计的精度，采用基于列向量夹角检测的混合矩阵优化算法，通过多次估计混合矩阵，对每次估计的矩阵进行列向量的相关性检测，出现次数频繁的列向量，是真正的混合矩阵的列向量；混合矩阵列向量顺序的改变，以及每一个列向量正负符号的改变，对混合矩阵的正确性都是不影响的，具体步骤如下：

1)对相同的观测信号，进行时频变换，设定时频单源点检测阈值，检测得到时频单源点集合；

2)多次估计混合矩阵，每次估计的结果写为A_e1，A_e2， … ，A_ek，k为设定的混合矩阵估计的最大次数。设混合矩阵的维数为m×n，则估计得到的混合矩阵维数也为m×n，将A_e1，A_e2， … ，A_ek组合成为一个新的矩阵并转置，记为W，则W的维数为N×m，N＝n×k；

3)计算矩阵W的行向量的夹角，设计算结果表示成矩阵W_e，显然W_e是N×N维的矩阵，该矩阵沿主对角线，右上半与左下半是对称的，包含的信息一样；可以令左下半元素等于0，只考虑右上半矩阵，设矩阵为W_e1；

4)令W_e1主对角线右上部元素大于10的元素等于0，矩阵记为W_e2。根据W_e2对同一或相似，这里设定列向量夹角小于10的列向量出现的次数进行统计，根据给定的辐射源个数，选择出现次数频繁的列向量作为预估混合矩阵的列向量，组成一个新的估计得到的混合矩阵A_e，其维数为m×n。

由于采用如上所述的技术方案，本发明具有如下优越性：

一种基于小波分析的时频单源点检测阈值确定方法，适用于观测信号叠加加性噪声的一般情况。能够较好地克服K均值聚类法在混合矩阵估计中易受初始初始聚类中心选择影响从而估计精度不高的问题，得到混合矩阵估计的优化结果。

附图说明

图1为阈值计算方法流程图；

图2(a)混合信号在时频域上的单源邻域、多源邻域以及没有源信号的邻域的信号图；图2(b)表示选择那些实部和虚部都具有较大能量的观测信号图；图2(c)表示计算图(b)中被选择信号的时频点比率图。图2(d)表示通过判断虚部与阈值ε的大小关系，检测得到时频单源点信号图。

具体实施方法

如图1、2所示，一种基于小波分析和时频单源检测的混合矩阵估计方法，在观测信号混叠加性噪声的情况下，通过小波分析方法对原观测信号进行小波分解、重构实现原始观测信号的扩维，然后对扩维后的观测信号矩阵进行奇异值分解，基于计算得到的奇异值向量，得到时频单源点检测的判断阈值计算方法，然后计算得到时频单源点集合，基于时频单源点集合应用K均值聚类法进行混合矩阵估计；重复前述混合矩阵估计步骤，得到多次混合矩阵估计结果，将多次估计得到的混合矩阵组合成为一个矩阵，对该矩阵进行列向量夹角检测，设辐射源数为n，选择出现最频繁的n个列向量组合成为新的混合矩阵，作为最终估计得到的混合矩阵，具体实现步骤如下：

(1)基于小波分析的时频单源点判断阈值确定方法

基于小波分析的时频单源点检测阈值确定方法步骤如下：

1)对原始观测数据进行白化处理；

2)对白化后的观测信号X(t)的小波分解，设小波分解层数为I，x₁(t)分解得到的多路小波系数分别为a_1j，d_1j，j＝1,2…I，上述小波系数的下标中，第一个“1”表示第一个原始观测信号x₁(t)，其中a_1j表示第1个观测信号的低频小波系数，d_1j表示第1个观测信号的高频小波系数；其他观测信号小波分解方法按照上述方法进行操作；

3)小波系数重构，将m个观测信号分解得到的小波系数序列信号分别进行小波重构，以第一个观测信号x₁(t)分解得到的多路小波系数a_1j，d_1j，j＝1,2…I为例，将低频小波系数重构为x₁₀(t)，每一层高频小波系数重构为x_1j(t)，j＝1,2…I，其他观测信号分解得到的小波系数重构方法及表示方法与此相同；

4)将白化后的观测信号和其小波重构信号组合成为新的多维观测信号，以第一路观测信号x₁(t)为例说明其组合方法。选取x₁(t)，x₁₀(t)，x_1j(t)，j＝1,2…I作为第一路观测信号扩维后的信号，x₁₀(t)为由x₁(t)进行I层小波系数分解后得到的低频系数重构的信号，x_1j(t)为第j层高频系数重构后的信号。其他观测信号扩维的方法与上述相同。经过扩维后的观测信号表示为：

x_w(t)＝[x₁(t),x₁₀(t),x_1j(t),x₂(t),x₂₀(t),x_2j(t),…,x_m(t),x_m0(t),x_mj(t),]^T,j＝1,2…I (9)

5)新的多维观测信号Xw(t)的去均值化处理，表示为

6)计算多维观测信号

的协方差矩阵R_w，

T₀为观测信号采样点数。

7)协方差矩阵R_w的奇异值分解。计算矩阵R_w的奇异值分解

得到其奇异值Λ_s＝{λ₁,λ₂, … λ_M}，其中M＝m×(I+2)为升维后的观测信号的维数。

8)得到时频单源点检测判断阈值。时频单源点判断阈值ε的取值范围为ε＝[λ_m+1,λ_m+2]。

以上即为基于小波分析的时频单源点检测阈值确定方法，该方法考虑到观测信号中叠加的加性噪声，拓展了文献[3]方法的应用范围，提高了计算效率和精度。阈值估计算法流程如图1所示：

需要说明的是，在小波分析理论中，小波基函数及分解层数的选择是小波分析方法应用的关键步骤。symN系列小波基函数和dbN系列小波基函数具有优良的特性，在信号的小波分析中效果比较突出，因此通常被选为小波基函数。分解层数对于信号的降噪，是一个重要的问题，一般分解2、3层，就会取得较好的应用效果。

(2)基于时频单源点检测的混合矩阵估计

对于第i个辐射源信号s_i(n)进行短时傅里叶变换(Short-time Fouriertransform,STFT)：

其中，h(n)是窗函数，τ是STFT的加窗函数中心，k是频率点。令X(τ,k)＝[X₁(τ,k),X₂(τ,k),…X_m(τ,k)]^T，S(τ,k)＝[S₁(τ,k),S₂(τ,k),…S_m(τ,k)]^T分别表示观测信号X和源信号S的STFT变换。

单源点检测法基于观测信号时频变换比率寻找一个时频点集合S，在该集合中每个时频点只有一个辐射源信号是活跃的。基于时频单源点检测的混合矩阵估计算法核心思想是：对于只有单源活跃的时频点，观测信号的时频域比率是实数，而对于多源活跃的时频点，观测信号的时频域比率为复数。基于观测信号时频变换的比率，检测只有单源活跃的时频点，即可估计得到混合矩阵，算法主要步骤如下：

1)给定一个时频单源点检测判断阈值ε＞0，检测得到只有单源活跃的时频点集合S_s，其中观测信号STFT变换的实部和虚部都有足够的能量：

在式(3)中，Im(x)表示x的虚部。可以应用其他观测信号的STFT代替X₁(τ,k)作为分母。

在时频点(τ_s,k_s)，如果满足式(12)，则(τ_s,k_s)∈S_s，

其中，a_li为混合矩阵A的第(l,i)个元素。由式(4)可判断在时频点(τ_s,k_s)处只存在第i个源信号。相反，如果第i个和第j个源信号S_i(τ,k)和S_j(τ,k)在时频点(τ_n,k_n)处同时存在，则时频单源点集合S_s中将不包括时频点(τ_n,k_n)，因为：

式(13)表示S_i(τ,k)和S_j(τ,k)为复数，不能如式(12)一样从分子和分母中消去。

图2为两路观测信号X₁和X₂的情况下，基于时频单源点检测的混合矩阵估计算法的主要步。示意图上的矩形区域分别表示不同的时频邻域，图(a)中不同的颜色分别表示混合信号在时频域上的单源邻域、多源邻域以及没有源信号的邻域。图(b)表示选择那些实部和虚部都具有较大能量的观测信号。图(c)表示计算图(b)中被选择信号的时频点比率，在SSO区域比率为实数，在MSO区域比率为复数。图(d)表示通过判断虚部与阈值ε的大小关系，检测得到时频单源点。

图2两路观测信号情况下SSD算法的示意图^[3](a)白色、轻度灰色和深度灰色区域分别表示该时频点没有源信号、只有一个源信号和多个源信号；(b)选择混合信号STFT的实部和虚部的能量都较大的时频点；(c)计算图(b)中所选择观测信号时频点的比率；(d)通过判断比率的虚部是否足够小，检测得到时频单源点

2)应用聚类算法，基于式(14)所示比率向量，将时频单源点集合S_s聚成m类

其中，Re(x)表示x的实部。本发明应用k均值聚类算法对S_s进行聚类。在S_s中，若第i个源信号的时频点活跃，则有以下比率向量：

基于上述比率向量的时频单源点集合表示为

应用上述算法，可以检测得到每一个辐射源信号对应的单源点区域。

3)基于所有辐射源信号的时频单源点集合估计得到混合矩阵A。在混合矩阵估计之前，所有观测信号的实部和虚部都通过l₂范数进行归一化。

基于第i个源信号活跃的时频单源点集合

通过式(16)估计得到混合矩阵的第i列

(2)基于向量夹角检测的混合矩阵优化方法

基于列向量夹角检测的混合矩阵优化方法是：通过多次估计混合矩阵，对每次估计的矩阵进行列向量的相关性检测，出现次数频繁的列向量，一般就是真正的混合矩阵的列向量。相当于概率论中的大概率事件。这里需要说明的是，混合矩阵列向量顺序的改变，以及每一个列向量正负符号的改变，对混合矩阵的正确性都是不影响的。方法的具体步骤如下：

1)对相同的观测信号，进行时频变换，设定时频单源点检测阈值，检测得到时频单源点集合；

2)多次估计混合矩阵，每次估计的结果写为A_e1，A_e2， … ，A_ek，k为设定的混合矩阵估计的最大次数。设混合矩阵的维数为m×n，则估计得到的混合矩阵维数也为m×n，将A_e1，A_e2， … ，A_ek组合成为一个新的矩阵并转置，记为W，则W的维数为N×m，N＝n×k；

3)利用式(17)，计算矩阵W的行向量的夹角，设计算结果表示成矩阵W_e，显然W_e是N×N维的矩阵，该矩阵沿主对角线，右上半与左下半是对称的，包含的信息一样。可以令左下半元素等于0，只考虑右上半矩阵，设矩阵为W_e1；

预估矩阵与原混合矩阵相应列矢量间的夹角计算表达式为：

式中，a—原始混合矩阵A的列矢量；

—估计的混合矩阵A_e的列矢量。

值越小表示混合矩阵估计精度越高。

4)令W_e1主对角线右上部元素大于10的元素等于0，矩阵记为W_e2。根据W_e2对同一或相似(这里设定列向量夹角小于10)列向量出现的次数进行统计，根据给定的辐射源个数(应用其他方法估计得到)，选择出现次数频繁的列向量作为预估混合矩阵的列向量，组成一个新的估计得到的混合矩阵A_e，其维数为m×n。

仿真实验

本发明的优点可以通过以下仿真试验进一步说明，源信号是文献[3]中应用的四个语音信号，数据长度是40000，将四路源信号混合成三路观测信号，即m＝3，n＝4，混合矩阵为：

设小波分解层数为3，小波基函数为db4，在不同信噪比(Signal to Noise Ratio,SNR)情况下，根据第3.2节方法计算得到的阈值范围如表1所示。

表1单源点判断阈值计算结果

在混合信号信噪比为20dB的情况下，取不同阈值时，估计得到的混合矩阵及向量夹角计算如下：

(1)阈值取值：0.0417

设混合矩阵估计次数为10，每次估计的混合矩阵为：

将上述估计得到的混合矩阵与真正的混合矩阵比较可以发现，估计得到的混合矩阵有的与原混合矩阵十分相似，而有的只有部分列向量与原混合矩阵相似。根据第3.3节方法，将上述矩阵组合成为W，基于W进行行向量夹角检测，得到矩阵W_e，令W_e主对角线左下部元素为0，主对角线右上部元素大于10的元素为0，得到矩阵W_e2，由于数据量庞大，这里不再列出。下面给出出现次数较多的W的行向量及其夹角计算结果：

ang(W₂,W₁₁)＝5.1615，ang(W₂,W₁₈)＝1.4263，ang(W₂,W₂₂)＝1.2539，ang(W₂,W₂₅)＝2.3668，ang(W₂,W₃₄)＝1.1231，可认为W的第2行出现5次；

ang(W₃,W₅)＝0.7574，ang(W₃,W₁₀)＝0.1825，ang(W₃,W₁₄)＝0.4078，ang(W₃,W₁₇)＝0.7497，ang(W₃,W₂₃)＝0.3641，ang(W₃,W₂₈)＝0.9987，ang(W₃,W₃₀)＝0.6175，ang(W₃,W₃₃)＝0.5555，ang(W₃,W₃₉)＝0.1020，可认为W的第3行出现9次；

ang(W₄,W₈)＝1.9218，ang(W₄,W₉)＝2.9391，ang(W₄,W₁₉)＝0.5493，ang(W₄,W₂₄)＝5.3649，ang(W₄,W₂₇)＝2.9136，ang(W₄,W₃₆)＝2.8871，可认为W的第4行出现6次；

ang(W₁₂,W₁₆)＝0.5829，ang(W₁₂,W₂₀)＝0.4805，ang(W₁₂,W₂₁)＝1.2405，ang(W₁₂,W₂₆)＝1.1877，ang(W₁₂,W₃₁)＝1.3788，ang(W₁₂,W₃₅)＝0.7874，ang(W₁₂,W₃₇)＝9.9659，ang(W₁₂,W₄₀)＝1.4246，可认为W的第12行出现8次；

ang(W₂₉,W₃₈)＝4.3206，可认为W的第29行出现1次。

选择W的第2、3、4、12行组成预估混合矩阵

计算其与原混合矩阵的列向量夹角，结果如下：

由此可见，估计的混合矩阵与原混合矩阵的列向量夹角都在5°以下，精度很高，表示很好地得到了原混合矩阵的估计。

(2)阈值取值：0.03

设混合矩阵估计次数为10，每次估计的混合矩阵为：

将上述矩阵组合成为W，基于W进行行向量夹角检测，得到矩阵W_e，令W_e主对角线左下部元素为0，主对角线右上部元素大于10的元素为0，得到矩阵W_e2，由于数据量庞大，这里不再列出。下面给出出现次数较多的W的行向量及其夹角计算结果：

ang(W₂,W₁₃)＝7.6701，ang(W₂,W₂₄)＝2.5049，ang(W₂,W₃₀)＝3.7955，ang(W₂,W₄₀)＝3.9075，可认为W的第2行出现4次；

ang(W₃,W₁₆)＝2.9258，ang(W₃,W₂₂)＝1.8824，ang(W₃,W₃₈)＝1.0815，可认为W的第3行出现3次；

ang(W₅,W₁₀)＝0.5221，ang(W₅,W₁₅)＝0.4693，ang(W₅,W₁₈)＝0.7092，ang(W₅,W₁₉)＝8.7969，ang(W₅,W₂₃)＝0.2832，ang(W₅,W₂₈)＝0.7491，ang(W₅,W₂₉)＝0.7912，ang(W₅,W₃₅)＝0.6204，ang(W₅,W₃₇)＝0.7821，可认为W的第5行出现9次；

ang(W₆,W₉)＝1.1910，ang(W₆,W₁₄)＝0.6675，ang(W₆,W₁₇)＝1.2382，ang(W₆,W₂₁)＝1.8564，ang(W₆,W₂₇)＝2.2968，ang(W₆,W₃₃)＝1.5670，ang(W₆,W₃₉)＝1.5724，可认为W的第6行出现7次；

ang(W₈,W₂₀)＝3.9659，ang(W₈,W₂₆)＝7.1850，ang(W₈,W₃₄)＝6.3789，可认为W的第8行出现3次；

选择W的第2、3、5、6行组成预估混合矩阵：

它与原混合矩阵的列向量夹角分别为：

第2、5、6、8行组成预估混合矩阵：

它与原混合矩阵的列向量夹角分别为：

由此可看出原混合矩阵中的4列被正确地估计出了3列。

出现上述得到两个混合矩阵的原因是最初只进行了10次混合矩阵估计，样本数较少，若增大估计次数，将能更准确地得到混合矩阵的唯一解。

(3)阈值取值：0.0227

运行20次，将出现次数最多的4列组成混合矩阵，结果如下：

估计的混合矩阵与原混合矩阵的列向量夹角：

可见应用本发明方法估计结果精度很好。

实验结果分析

针对实际的语音信号，开展了欠定条件下基于小波分析和时频单源点检测的混合矩阵估计仿真研究，从仿真实验结果可以看出，估计得到的混合矩阵与原混合矩阵列向量夹角偏差绝大部分都在10°以下，表明混合矩阵估计的精度较好。需要说明的是，时频单源点检测阈值取值范围并非唯一的，但是本发明方法使得时频单源点检测方法的关键参数—单源点检测阈值的设置一定程度上能够有章可循，在观测信号混叠加性噪声的一般情况下，本发明方法能够给出一个较可靠的阈值取值范围，提高时频单源点检测的效率和精度，结合基于向量夹角检测的混合矩阵优化方法，拓展了基于时频单源点检测的混合矩阵估计方法的应用范围，提高了混合矩阵的估计精度。

Claims (1)

1.一种基于小波分析和时频单源检测的混合矩阵估计方法，其特征是：在观测信号混叠加性噪声的情况下，通过小波分析方法对原观测信号进行小波分解、重构实现原始观测信号的扩维，然后对扩维后的观测信号矩阵进行奇异值分解，基于计算得到的奇异值向量，得到时频单源点检测的判断阈值计算方法，然后计算得到时频单源点集合，基于时频单源点集合应用K均值聚类法进行混合矩阵估计；重复前述混合矩阵估计步骤，得到多次混合矩阵估计结果，将多次估计得到的混合矩阵组合成为一个矩阵，对该矩阵进行行向量夹角检测，辐射源数为n，选择出现频繁的n个行向量，组合成为新的混合矩阵，作为最终估计得到的混合矩阵，其步骤如下：

(1)基于小波分析的时频单源点判断阈值确定方法

为了解决在观测信号混叠加性噪声情况下，时频单源点检测中判断阈值选取缺乏理论指导的问题，基于小波分析原理，得到一种阈值的选取方法，具体步骤如下：

1)对原始观测数据进行白化处理；

2)对白化后的观测信号X(t)的小波分解，小波分解层数为I，x₁(t)分解得到的多路小波系数分别为a_1-j，d_1-j，j＝1，2…I，上述小波系数的下标中，第一个“1”表示第一个原始观测信号x₁(t)，其中a_1-j表示第1个观测信号的低频小波系数，d_1-j表示第1个观测信号的高频小波系数；其他观测信号小波分解方法按照步骤2)方法进行操作；

3)小波系数重构，将m个观测信号分解得到的小波系数序列信号分别进行小波重构，以第一个观测信号x₁(t)分解得到的多路小波系数a_1-j，d_1-j，j＝1，2…I为例，将低频小波系数重构为x_1-0(t)，每一层高频小波系数重构为x_1-j(t)，j＝1，2，…，I；其他观测信号分解得到的小波系数重构方法及表示方法与本步骤3)相同；

4)将白化后的观测信号和其小波重构信号组合成为新的多维观测信号，以第一路观测信号x₁(t)为例说明其组合方法，选取x₁(t)，x_1-0(t)，x_1-j(t)，j＝1，2…I作为第一路观测信号扩维后的信号，x_1-0(t)为由x₁(t)进行I层小波系数分解后得到的低频系数重构的信号，x_1-j(t)为第j层高频系数重构后的信号；其他观测信号扩维的方法与本步骤4)相同，经过扩维后的观测信号表示为：

X_w(t)＝[x₁(t)，x_1-0(t)，x_1-j(t)，x₂(t)，x_2-0(t)，x_2-j(t)，…，x_m(t)，x_m-0(t)x_m-j(t)]^T，j＝1，2…I (1)

5)新的多维观测信号X_w(t)的去均值化处理，表示为

6)计算多维观测信号

的协方差矩阵R_w，

T₀为观测信号采样点数；

7)协方差矩阵R_w的奇异值分解：协方差矩阵R_w的奇异值分解表示为R_w＝V·Λ·V^T，得到其奇异值Λ＝{λ₁，λ₂，…，λ_M}，其中M＝m×(I+2)为扩维后的观测信号X_w(t)的维数；

8)得到时频单源点检测判断阈值；时频单源点判断阈值ε的取值范围为ε＝[λ_m+1，λ_m+2]；

以上即为基于小波分析的时频单源点检测阈值确定方法，该方法考虑到观测信号中叠加的加性噪声；

(2)基于时频单源点检测的混合矩阵估计

对于第i个辐射源信号s_i(t)进行短时傅里叶变换(Short-timeFouriertransform，STFT)：

其中，h(t)是窗函数，τ是STFT的加窗函数中心，k是频率点；令X(τ，k)＝[X₁(τ，k)，X₂(τ，k)，…X_m(τ，k)]^T，s(τ，k)＝[S₁(τ，k)，S₂(τ，k)，…S_n(τ，k)]^T分别表示观测信号X和源信号S的STFT变换；

单源点检测法基于观测信号时频变换比率寻找一个时频点集合S_s，在该集合中每个时频点只有一个辐射源信号是活跃的，基于时频单源点检测的混合矩阵估计算法核心是：对于只有单源活跃的时频点，观测信号的时频域比率是实数，而对于多源活跃的时频点，观测信号的时频域比率为复数；基于观测信号时频变换的比率，检测只有单源活跃的时频点，即可估计得到混合矩阵，算法步骤如下：

1)给定一个时频单源点检测判断阈值ε＞0，检测得到只有单源活跃的时频点集合S_s，其中观测信号STFT变换的实部和虚部都有足够的能量：

在式(3)中，Im(x)表示x的虚部，应用其他观测信号的STFT代替X₁(τ，k)作为分母；

在时频点(τ_s，k_s)，如果满足式(4)，则(τ_s，k_s)∈S_s，

其中，a_li为混合矩阵A的第(l，i)个元素；由式(4)可判断在时频点(τ_s，k_s)处只存在第i个源信号；相反，如果第i个和第j个源信号S_i(τ，k)和S_j(τ，k)在时频点(τ_p，k_p)处同时存在，则时频单源点集合S_s中将不包括时频点(τ_p，k_p)，因为：

式(5)表示S_i(τ，k)和S_j(τ，k)为复数，不能如式(4)一样从分子和分母中消去；

2)应用聚类算法，基于式(6)所示比率向量，将时频单源点集合S_s聚n类

其中，Re(x)表示x的实部；应用K均值聚类算法对S_s进行聚类；在S_s中，若第i个源信号的时频点活跃，则有以下比率向量：

基于上述比率向量的时频单源点集合表示为

应用上述算法，检测得到每一个辐射源信号对应的单源点区域；

3)基于所有辐射源信号的时频单源点集合估计得到混合矩阵A，在混合矩阵估计之前，所有观测信号的实部和虚部都通过l₂范数进行归一化；

基于第i个源信号活跃的时频单源点集合

通过式(8)估计得到混合矩阵的第i列

该式表示实部加虚部的求和操作；

(3)基于向量夹角检测的混合矩阵优化方法

检测得到时频单源点集合以后，应用K均值聚类算法进行混合矩阵估计，K均值聚类算法具有计算量小，收敛快的优势，但是聚类效果严重依赖于初始聚类中心的选择；对于相同的观测信号，即便时频单源点检测精度很高，由于K均值聚类法每次自动选择的初始聚类中心是不确定的，导致每次估计得到的混合矩阵相应地是变化的，每次估计得到的混合矩阵可能只有一部分或全部列向量与真正的混合矩阵相似；由此分析可知，基于时频单源点检测的混合矩阵估计效果受K均值聚类算法影响很大；

为了提高混合矩阵估计的精度，采用基于向量夹角检测的混合矩阵优化算法，通过多次估计混合矩阵，对每次估计的混合矩阵进行列向量的相关性检测，出现次数频繁的列向量，是真正的混合矩阵的列向量；混合矩阵列向量顺序的改变，以及每一个列向量正负符号的改变，对混合矩阵的正确性都是不影响的，具体步骤如下：

1)对相同的观测信号，进行时频变换，根据时频单源点检测阈值，检测得到时频单源点集合；

2)多次估计混合矩阵，每次估计的结果写为A_e1，A_e2，…，A_eL，L为混合矩阵估计的最大次数；混合矩阵的维数为m×n，则估计得到的混合矩阵维数也为m×n，将A_e1，A_e2，…，A_ek组合成为一个新的矩阵并转置，记为W，则W的维数为N×m，N＝n×L；

3)计算矩阵W的行向量的夹角，计算结果表示成矩阵W_e，显然W_e是N×N维的矩阵，该矩阵沿主对角线，右上半与左下半是对称的，包含的信息一样；可以令左下半元素等于0，只考虑右上半矩阵，该矩阵为W_e1；

4)令W_e1主对角线右上部元素大于10的元素等于10，矩阵记为W_e2，根据W_e2对矩阵W的同一或相似行向量出现的次数进行统计，同一或相似指行向量夹角小于10，根据辐射源个数n，选择出现次数频繁的n个行向量，作为预估混合矩阵的列向量，组合成为最终估计的混合矩阵A_e，其维数为m×n。

Images