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

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

Info

Publication number
CN108108666A
CN108108666A CN201711236736.9A CN201711236736A CN108108666A CN 108108666 A CN108108666 A CN 108108666A CN 201711236736 A CN201711236736 A CN 201711236736A CN 108108666 A CN108108666 A CN 108108666A
Authority
CN
China
Prior art keywords
mrow
msub
time
matrix
mtd
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
CN201711236736.9A
Other languages
English (en)
Other versions
CN108108666B (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.)
STATE KEY LABORATORY OF COMPLEX ELECTROMAGNETIC ENVIRONMENTAL EFFECTS ON ELECTRONICS & INFORMATION SYSTEM
Original Assignee
STATE KEY LABORATORY OF COMPLEX ELECTROMAGNETIC ENVIRONMENTAL EFFECTS ON ELECTRONICS & INFORMATION SYSTEM
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 STATE KEY LABORATORY OF COMPLEX ELECTROMAGNETIC ENVIRONMENTAL EFFECTS ON ELECTRONICS & INFORMATION SYSTEM filed Critical STATE KEY LABORATORY OF COMPLEX ELECTROMAGNETIC ENVIRONMENTAL EFFECTS ON ELECTRONICS & INFORMATION SYSTEM
Priority to CN201711236736.9A priority Critical patent/CN108108666B/zh
Publication of CN108108666A publication Critical patent/CN108108666A/zh
Application granted granted Critical
Publication of CN108108666B publication Critical patent/CN108108666B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2218/00—Aspects of pattern recognition specially adapted for signal processing
    • G06F2218/02—Preprocessing
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00—Pattern recognition
    • G06F18/20—Analysing
    • G06F18/23—Clustering techniques
    • G06F18/232—Non-hierarchical techniques
    • G06F18/2321—Non-hierarchical techniques using statistics or function optimisation, e.g. modelling of probability density functions
    • G06F18/23213—Non-hierarchical techniques using statistics or function optimisation, e.g. modelling of probability density functions with fixed number of clusters, e.g. K-means clustering

Landscapes

  • Engineering & Computer Science (AREA)
  • Theoretical Computer Science (AREA)
  • Data Mining & Analysis (AREA)
  • Physics & Mathematics (AREA)
  • Artificial Intelligence (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • General Engineering & Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Probability & Statistics with Applications (AREA)
  • Signal Processing (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Evolutionary Biology (AREA)
  • Evolutionary Computation (AREA)
  • Other Investigation Or Analysis Of Materials By Electrical Means (AREA)
  • Complex Calculations (AREA)

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，x1(t)分解得到的多路小波系数分别为a1j，d1j，j＝1,2…I，上述小波系数的下标中，第一个“1”表示第一个原始观测信号x1(t)，其中a1j表示第1个观测信号的低频小波系数，d1j表示第1个观测信号的高频小波系数；其他观测信号小波分解方法按照上述方法进行操作；
3)小波系数重构，将m个观测信号分解得到的小波系数序列信号分别进行小波重构，以第一个观测信号x1(t)分解得到的多路小波系数a1j，d1j，j＝1,2…I为例，将低频小波系数重构为x10(t)，每一层高频小波系数重构为x1j(t)，j＝1,2…I，其他观测信号分解得到的小波系数重构方法及表示方法与此相同；
4)将白化后的观测信号和其小波重构信号组合成为新的多维观测信号，以第一路观测信号x1(t)为例说明其组合方法。选取x1(t)，x10(t)，x1j(t)，j＝1,2…I作为第一路观测信号扩维后的信号，x10(t)为由x1(t)进行I层小波系数分解后得到的低频系数重构的信号，x1j(t)为第j层高频系数重构后的信号。其他观测信号扩维的方法与上述相同。经过扩维后的观测信号表示为：
xw(t)＝[x1(t),x10(t),x1j(t),x2(t),x20(t),x2j(t),…,xm(t),xm0(t),xmj(t),]T,j＝1,2…I (1)
5)新的多维观测信号Xw(t)的去均值化处理，表示为
6)计算多维观测信号的协方差矩阵Rw，T0为观测信号采样点数。
7)协方差矩阵Rw的奇异值分解。计算矩阵Rw的奇异值分解得到其奇异值Λs＝{λ1,λ2, … λM}，其中M＝m×(I+2)为升维后的观测信号的维数。
8)得到时频单源点检测判断阈值。时频单源点判断阈值ε的取值范围为ε＝[λm+1,λm+2]。
以上即为基于小波分析的时频单源点检测阈值确定方法，该方法考虑到观测信号中叠加的加性噪声；
(2)基于时频单源点检测的混合矩阵估计
对于第i个辐射源信号si(n)进行短时傅里叶变换(Short-time Fouriertransform,STFT)：
其中，h(n)是窗函数，τ是STFT的加窗函数中心，k是频率点；令X(τ,k)＝[X1(τ,k),X2(τ,k),…Xm(τ,k)]T，S(τ,k)＝[S1(τ,k),S2(τ,k),…Sm(τ,k)]T分别表示观测信号X和源信号S的STFT变换；
单源点检测法基于观测信号时频变换比率寻找一个时频点集合S，在该集合中每个时频点只有一个辐射源信号是活跃的，基于时频单源点检测的混合矩阵估计算法核心是：对于只有单源活跃的时频点，观测信号的时频域比率是实数，而对于多源活跃的时频点，观测信号的时频域比率为复数；基于观测信号时频变换的比率，检测只有单源活跃的时频点，即可估计得到混合矩阵，算法步骤如下：
1)给定一个时频单源点检测判断阈值ε＞0，检测得到只有单源活跃的时频点集合Ss，其中观测信号STFT变换的实部和虚部都有足够的能量：
在式(3)中，Im(x)表示x的虚部，应用其他观测信号的STFT代替X1(τ,k)作为分母；
在时频点(τs,ks)，如果满足式(4)，则(τs,ks)∈Ss，
其中，ali为混合矩阵A的第(l,i)个元素。由式(4)可判断在时频点(τs,ks)处只存在第i个源信号。相反，如果第i个和第j个源信号Si(τ,k)和Sj(τ,k)在时频点(τn,kn)处同时存在，则时频单源点集合Ss中将不包括时频点(τn,kn)，因为：
式(5)表示Si(τ,k)和Sj(τ,k)为复数，不能如式(4)一样从分子和分母中消去；
2)应用聚类算法，基于式(6)所示比率向量，将时频单源点集合Ss聚成m类
其中，Re(x)表示x的实部；应用k均值聚类算法对Ss进行聚类；在Ss中，若第i个源信号的时频点活跃，则有以下比率向量：
基于上述比率向量的时频单源点集合表示为应用上述算法，检测得到每一个辐射源信号对应的单源点区域；
3)基于所有辐射源信号的时频单源点集合估计得到混合矩阵A，在混合矩阵估计之前，所有观测信号的实部和虚部都通过l2范数进行归一化；
基于第i个源信号活跃的时频单源点集合通过式(8)估计得到混合矩阵的第i列
(3)基于向量夹角检测的混合矩阵优化方法
检测得到时频单源点集合以后，应用K均值聚类算法进行混合矩阵估计，K均值聚类算法具有计算量小，收敛快的优势，但是聚类效果严重依赖于初始聚类中心的选择；对于相同的观测信号，即便时频单源点检测精度很高，由于K均值聚类法每次自动选择的初始聚类中心是不确定的，导致每次估计得到的混合矩阵相应地是变化的，每次估计得到的混合矩阵可能只有一部分或全部列向量与真正的混合矩阵相似；由此分析可知，基于时频单源点检测的混合矩阵估计效果受K均值聚类算法影响很大；
为了提高混合矩阵估计的精度，采用基于列向量夹角检测的混合矩阵优化算法，通过多次估计混合矩阵，对每次估计的矩阵进行列向量的相关性检测，出现次数频繁的列向量，是真正的混合矩阵的列向量；混合矩阵列向量顺序的改变，以及每一个列向量正负符号的改变，对混合矩阵的正确性都是不影响的，具体步骤如下：
1)对相同的观测信号，进行时频变换，设定时频单源点检测阈值，检测得到时频单源点集合；
2)多次估计混合矩阵，每次估计的结果写为Ae1，Ae2， … ，Aek，k为设定的混合矩阵估计的最大次数。设混合矩阵的维数为m×n，则估计得到的混合矩阵维数也为m×n，将Ae1，Ae2， … ，Aek组合成为一个新的矩阵并转置，记为W，则W的维数为N×m，N＝n×k；
3)计算矩阵W的行向量的夹角，设计算结果表示成矩阵We，显然We是N×N维的矩阵，该矩阵沿主对角线，右上半与左下半是对称的，包含的信息一样；可以令左下半元素等于0，只考虑右上半矩阵，设矩阵为We1；
4)令We1主对角线右上部元素大于10的元素等于0，矩阵记为We2。根据We2对同一或相似，这里设定列向量夹角小于10的列向量出现的次数进行统计，根据给定的辐射源个数，选择出现次数频繁的列向量作为预估混合矩阵的列向量，组成一个新的估计得到的混合矩阵Ae，其维数为m×n。
由于采用如上所述的技术方案，本发明具有如下优越性：
一种基于小波分析的时频单源点检测阈值确定方法，适用于观测信号叠加加性噪声的一般情况。能够较好地克服K均值聚类法在混合矩阵估计中易受初始初始聚类中心选择影响从而估计精度不高的问题，得到混合矩阵估计的优化结果。
附图说明
图1为阈值计算方法流程图；
图2(a)混合信号在时频域上的单源邻域、多源邻域以及没有源信号的邻域的信号图；图2(b)表示选择那些实部和虚部都具有较大能量的观测信号图；图2(c)表示计算图(b)中被选择信号的时频点比率图。图2(d)表示通过判断虚部与阈值ε的大小关系，检测得到时频单源点信号图。
具体实施方法
如图1、2所示，一种基于小波分析和时频单源检测的混合矩阵估计方法，在观测信号混叠加性噪声的情况下，通过小波分析方法对原观测信号进行小波分解、重构实现原始观测信号的扩维，然后对扩维后的观测信号矩阵进行奇异值分解，基于计算得到的奇异值向量，得到时频单源点检测的判断阈值计算方法，然后计算得到时频单源点集合，基于时频单源点集合应用K均值聚类法进行混合矩阵估计；重复前述混合矩阵估计步骤，得到多次混合矩阵估计结果，将多次估计得到的混合矩阵组合成为一个矩阵，对该矩阵进行列向量夹角检测，设辐射源数为n，选择出现最频繁的n个列向量组合成为新的混合矩阵，作为最终估计得到的混合矩阵，具体实现步骤如下：
(1)基于小波分析的时频单源点判断阈值确定方法
基于小波分析的时频单源点检测阈值确定方法步骤如下：
1)对原始观测数据进行白化处理；
2)对白化后的观测信号X(t)的小波分解，设小波分解层数为I，x1(t)分解得到的多路小波系数分别为a1j，d1j，j＝1,2…I，上述小波系数的下标中，第一个“1”表示第一个原始观测信号x1(t)，其中a1j表示第1个观测信号的低频小波系数，d1j表示第1个观测信号的高频小波系数；其他观测信号小波分解方法按照上述方法进行操作；
3)小波系数重构，将m个观测信号分解得到的小波系数序列信号分别进行小波重构，以第一个观测信号x1(t)分解得到的多路小波系数a1j，d1j，j＝1,2…I为例，将低频小波系数重构为x10(t)，每一层高频小波系数重构为x1j(t)，j＝1,2…I，其他观测信号分解得到的小波系数重构方法及表示方法与此相同；
4)将白化后的观测信号和其小波重构信号组合成为新的多维观测信号，以第一路观测信号x1(t)为例说明其组合方法。选取x1(t)，x10(t)，x1j(t)，j＝1,2…I作为第一路观测信号扩维后的信号，x10(t)为由x1(t)进行I层小波系数分解后得到的低频系数重构的信号，x1j(t)为第j层高频系数重构后的信号。其他观测信号扩维的方法与上述相同。经过扩维后的观测信号表示为：
xw(t)＝[x1(t),x10(t),x1j(t),x2(t),x20(t),x2j(t),…,xm(t),xm0(t),xmj(t),]T,j＝1,2…I (9)
5)新的多维观测信号Xw(t)的去均值化处理，表示为
6)计算多维观测信号的协方差矩阵Rw，T0为观测信号采样点数。
7)协方差矩阵Rw的奇异值分解。计算矩阵Rw的奇异值分解得到其奇异值Λs＝{λ1,λ2, … λM}，其中M＝m×(I+2)为升维后的观测信号的维数。
8)得到时频单源点检测判断阈值。时频单源点判断阈值ε的取值范围为ε＝[λm+1,λm+2]。
以上即为基于小波分析的时频单源点检测阈值确定方法，该方法考虑到观测信号中叠加的加性噪声，拓展了文献[3]方法的应用范围，提高了计算效率和精度。阈值估计算法流程如图1所示：
需要说明的是，在小波分析理论中，小波基函数及分解层数的选择是小波分析方法应用的关键步骤。symN系列小波基函数和dbN系列小波基函数具有优良的特性，在信号的小波分析中效果比较突出，因此通常被选为小波基函数。分解层数对于信号的降噪，是一个重要的问题，一般分解2、3层，就会取得较好的应用效果。
(2)基于时频单源点检测的混合矩阵估计
对于第i个辐射源信号si(n)进行短时傅里叶变换(Short-time Fouriertransform,STFT)：
其中，h(n)是窗函数，τ是STFT的加窗函数中心，k是频率点。令X(τ,k)＝[X1(τ,k),X2(τ,k),…Xm(τ,k)]T，S(τ,k)＝[S1(τ,k),S2(τ,k),…Sm(τ,k)]T分别表示观测信号X和源信号S的STFT变换。
单源点检测法基于观测信号时频变换比率寻找一个时频点集合S，在该集合中每个时频点只有一个辐射源信号是活跃的。基于时频单源点检测的混合矩阵估计算法核心思想是：对于只有单源活跃的时频点，观测信号的时频域比率是实数，而对于多源活跃的时频点，观测信号的时频域比率为复数。基于观测信号时频变换的比率，检测只有单源活跃的时频点，即可估计得到混合矩阵，算法主要步骤如下：
1)给定一个时频单源点检测判断阈值ε＞0，检测得到只有单源活跃的时频点集合Ss，其中观测信号STFT变换的实部和虚部都有足够的能量：
在式(3)中，Im(x)表示x的虚部。可以应用其他观测信号的STFT代替X1(τ,k)作为分母。
在时频点(τs,ks)，如果满足式(12)，则(τs,ks)∈Ss，
其中，ali为混合矩阵A的第(l,i)个元素。由式(4)可判断在时频点(τs,ks)处只存在第i个源信号。相反，如果第i个和第j个源信号Si(τ,k)和Sj(τ,k)在时频点(τn,kn)处同时存在，则时频单源点集合Ss中将不包括时频点(τn,kn)，因为：
式(13)表示Si(τ,k)和Sj(τ,k)为复数，不能如式(12)一样从分子和分母中消去。
图2为两路观测信号X1和X2的情况下，基于时频单源点检测的混合矩阵估计算法的主要步。示意图上的矩形区域分别表示不同的时频邻域，图(a)中不同的颜色分别表示混合信号在时频域上的单源邻域、多源邻域以及没有源信号的邻域。图(b)表示选择那些实部和虚部都具有较大能量的观测信号。图(c)表示计算图(b)中被选择信号的时频点比率，在SSO区域比率为实数，在MSO区域比率为复数。图(d)表示通过判断虚部与阈值ε的大小关系，检测得到时频单源点。
图2两路观测信号情况下SSD算法的示意图[3](a)白色、轻度灰色和深度灰色区域分别表示该时频点没有源信号、只有一个源信号和多个源信号；(b)选择混合信号STFT的实部和虚部的能量都较大的时频点；(c)计算图(b)中所选择观测信号时频点的比率；(d)通过判断比率的虚部是否足够小，检测得到时频单源点
2)应用聚类算法，基于式(14)所示比率向量，将时频单源点集合Ss聚成m类
其中，Re(x)表示x的实部。本发明应用k均值聚类算法对Ss进行聚类。在Ss中，若第i个源信号的时频点活跃，则有以下比率向量：
基于上述比率向量的时频单源点集合表示为应用上述算法，可以检测得到每一个辐射源信号对应的单源点区域。
3)基于所有辐射源信号的时频单源点集合估计得到混合矩阵A。在混合矩阵估计之前，所有观测信号的实部和虚部都通过l2范数进行归一化。
基于第i个源信号活跃的时频单源点集合通过式(16)估计得到混合矩阵的第i列
(2)基于向量夹角检测的混合矩阵优化方法
基于列向量夹角检测的混合矩阵优化方法是：通过多次估计混合矩阵，对每次估计的矩阵进行列向量的相关性检测，出现次数频繁的列向量，一般就是真正的混合矩阵的列向量。相当于概率论中的大概率事件。这里需要说明的是，混合矩阵列向量顺序的改变，以及每一个列向量正负符号的改变，对混合矩阵的正确性都是不影响的。方法的具体步骤如下：
1)对相同的观测信号，进行时频变换，设定时频单源点检测阈值，检测得到时频单源点集合；
2)多次估计混合矩阵，每次估计的结果写为Ae1，Ae2， … ，Aek，k为设定的混合矩阵估计的最大次数。设混合矩阵的维数为m×n，则估计得到的混合矩阵维数也为m×n，将Ae1，Ae2， … ，Aek组合成为一个新的矩阵并转置，记为W，则W的维数为N×m，N＝n×k；
3)利用式(17)，计算矩阵W的行向量的夹角，设计算结果表示成矩阵We，显然We是N×N维的矩阵，该矩阵沿主对角线，右上半与左下半是对称的，包含的信息一样。可以令左下半元素等于0，只考虑右上半矩阵，设矩阵为We1；
预估矩阵与原混合矩阵相应列矢量间的夹角计算表达式为：
式中，a—原始混合矩阵A的列矢量；—估计的混合矩阵Ae的列矢量。
值越小表示混合矩阵估计精度越高。
4)令We1主对角线右上部元素大于10的元素等于0，矩阵记为We2。根据We2对同一或相似(这里设定列向量夹角小于10)列向量出现的次数进行统计，根据给定的辐射源个数(应用其他方法估计得到)，选择出现次数频繁的列向量作为预估混合矩阵的列向量，组成一个新的估计得到的混合矩阵Ae，其维数为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进行行向量夹角检测，得到矩阵We，令We主对角线左下部元素为0，主对角线右上部元素大于10的元素为0，得到矩阵We2，由于数据量庞大，这里不再列出。下面给出出现次数较多的W的行向量及其夹角计算结果：
ang(W2,W11)＝5.1615，ang(W2,W18)＝1.4263，ang(W2,W22)＝1.2539，ang(W2,W25)＝2.3668，ang(W2,W34)＝1.1231，可认为W的第2行出现5次；
ang(W3,W5)＝0.7574，ang(W3,W10)＝0.1825，ang(W3,W14)＝0.4078，ang(W3,W17)＝0.7497，ang(W3,W23)＝0.3641，ang(W3,W28)＝0.9987，ang(W3,W30)＝0.6175，ang(W3,W33)＝0.5555，ang(W3,W39)＝0.1020，可认为W的第3行出现9次；
ang(W4,W8)＝1.9218，ang(W4,W9)＝2.9391，ang(W4,W19)＝0.5493，ang(W4,W24)＝5.3649，ang(W4,W27)＝2.9136，ang(W4,W36)＝2.8871，可认为W的第4行出现6次；
ang(W12,W16)＝0.5829，ang(W12,W20)＝0.4805，ang(W12,W21)＝1.2405，ang(W12,W26)＝1.1877，ang(W12,W31)＝1.3788，ang(W12,W35)＝0.7874，ang(W12,W37)＝9.9659，ang(W12,W40)＝1.4246，可认为W的第12行出现8次；
ang(W29,W38)＝4.3206，可认为W的第29行出现1次。
选择W的第2、3、4、12行组成预估混合矩阵
计算其与原混合矩阵的列向量夹角，结果如下：
由此可见，估计的混合矩阵与原混合矩阵的列向量夹角都在5°以下，精度很高，表示很好地得到了原混合矩阵的估计。
(2)阈值取值：0.03
设混合矩阵估计次数为10，每次估计的混合矩阵为：
将上述矩阵组合成为W，基于W进行行向量夹角检测，得到矩阵We，令We主对角线左下部元素为0，主对角线右上部元素大于10的元素为0，得到矩阵We2，由于数据量庞大，这里不再列出。下面给出出现次数较多的W的行向量及其夹角计算结果：
ang(W2,W13)＝7.6701，ang(W2,W24)＝2.5049，ang(W2,W30)＝3.7955，ang(W2,W40)＝3.9075，可认为W的第2行出现4次；
ang(W3,W16)＝2.9258，ang(W3,W22)＝1.8824，ang(W3,W38)＝1.0815，可认为W的第3行出现3次；
ang(W5,W10)＝0.5221，ang(W5,W15)＝0.4693，ang(W5,W18)＝0.7092，ang(W5,W19)＝8.7969，ang(W5,W23)＝0.2832，ang(W5,W28)＝0.7491，ang(W5,W29)＝0.7912，ang(W5,W35)＝0.6204，ang(W5,W37)＝0.7821，可认为W的第5行出现9次；
ang(W6,W9)＝1.1910，ang(W6,W14)＝0.6675，ang(W6,W17)＝1.2382，ang(W6,W21)＝1.8564，ang(W6,W27)＝2.2968，ang(W6,W33)＝1.5670，ang(W6,W39)＝1.5724，可认为W的第6行出现7次；
ang(W8,W20)＝3.9659，ang(W8,W26)＝7.1850，ang(W8,W34)＝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，x1(t)分解得到的多路小波系数分别为a1j，d1j，j＝1,2…I，上述小波系数的下标中，第一个“1”表示第一个原始观测信号x1(t)，其中a1j表示第1个观测信号的低频小波系数，d1j表示第1个观测信号的高频小波系数；其他观测信号小波分解方法按照上述方法进行操作；
3)小波系数重构，将m个观测信号分解得到的小波系数序列信号分别进行小波重构，以第一个观测信号x1(t)分解得到的多路小波系数a1j，d1j，j＝1,2…I为例，将低频小波系数重构为x10(t)，每一层高频小波系数重构为x1j(t)，j＝1,2…I，其他观测信号分解得到的小波系数重构方法及表示方法与此相同；
4)将白化后的观测信号和其小波重构信号组合成为新的多维观测信号，以第一路观测信号x1(t)为例说明其组合方法，选取x1(t)，x10(t)，x1j(t)，j＝1,2…I作为第一路观测信号扩维后的信号，x10(t)为由x1(t)进行I层小波系数分解后得到的低频系数重构的信号，x1j(t)为第j层高频系数重构后的信号。其他观测信号扩维的方法与上述相同，经过扩维后的观测信号表示为：
xw(t)＝[x1(t),x10(t),x1j(t),x2(t),x20(t),x2j(t),…,xm(t),xm0(t),xmj(t),]T,j＝1,2…I (1)
5)新的多维观测信号Xw(t)的去均值化处理，表示为
6)计算多维观测信号的协方差矩阵Rw，T0为观测信号采样点数；
7)协方差矩阵Rw的奇异值分解；计算矩阵Rw的奇异值分解得到其奇异值Λs＝{λ1,λ2,…λM}，其中M＝m×(I+2)为升维后的观测信号的维数；
8)得到时频单源点检测判断阈值。时频单源点判断阈值ε的取值范围为ε＝[λm+1,λm+2]；
以上即为基于小波分析的时频单源点检测阈值确定方法，该方法考虑到观测信号中叠加的加性噪声；
(2)基于时频单源点检测的混合矩阵估计
对于第i个辐射源信号si(n)进行短时傅里叶变换(Short-time Fouriertransform,STFT)：
<mrow> <msub> <mi>S</mi> <mi>i</mi> </msub> <mrow> <mo>(</mo> <mi>&amp;tau;</mi> <mo>,</mo> <mi>k</mi> <mo>)</mo> </mrow> <mo>=</mo> <munderover> <mi>&amp;Sigma;</mi> <mrow> <mi>n</mi> <mo>=</mo> <mo>-</mo> <mi>&amp;infin;</mi> </mrow> <mi>&amp;infin;</mi> </munderover> <mi>h</mi> <mrow> <mo>(</mo> <mi>n</mi> <mo>-</mo> <mi>&amp;tau;</mi> <mo>)</mo> </mrow> <msub> <mi>s</mi> <mi>i</mi> </msub> <mrow> <mo>(</mo> <mi>n</mi> <mo>)</mo> </mrow> <msup> <mi>e</mi> <mrow> <mo>-</mo> <mi>j</mi> <mi>k</mi> <mi>n</mi> </mrow> </msup> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>2</mn> <mo>)</mo> </mrow> </mrow>
其中，h(n)是窗函数，τ是STFT的加窗函数中心，k是频率点；令X(τ,k)＝[X1(τ,k),X2(τ,k),…Xm(τ,k)]T，S(τ,k)＝[S1(τ,k),S2(τ,k),…Sm(τ,k)]T分别表示观测信号X和源信号S的STFT变换；
单源点检测法基于观测信号时频变换比率寻找一个时频点集合Ss，在该集合中每个时频点只有一个辐射源信号是活跃的，基于时频单源点检测的混合矩阵估计算法核心是：对于只有单源活跃的时频点，观测信号的时频域比率是实数，而对于多源活跃的时频点，观测信号的时频域比率为复数；基于观测信号时频变换的比率，检测只有单源活跃的时频点，即可估计得到混合矩阵，算法步骤如下：
1)给定一个时频单源点检测判断阈值ε＞0，检测得到只有单源活跃的时频点集合Ss，其中观测信号STFT变换的实部和虚部都有足够的能量：
<mrow> <msub> <mi>S</mi> <mi>s</mi> </msub> <mo>=</mo> <mo>{</mo> <mrow> <mo>(</mo> <mi>&amp;tau;</mi> <mo>,</mo> <mi>k</mi> <mo>)</mo> </mrow> <mo>|</mo> <mo>|</mo> <mo>|</mo> <mi>Im</mi> <mrow> <mo>(</mo> <mfrac> <mrow> <mi>X</mi> <mrow> <mo>(</mo> <mi>&amp;tau;</mi> <mo>,</mo> <mi>k</mi> <mo>)</mo> </mrow> </mrow> <mrow> <msub> <mi>X</mi> <mn>1</mn> </msub> <mrow> <mo>(</mo> <mi>&amp;tau;</mi> <mo>,</mo> <mi>k</mi> <mo>)</mo> </mrow> </mrow> </mfrac> <mo>)</mo> </mrow> <mo>|</mo> <mo>|</mo> <mo>&lt;</mo> <mi>&amp;epsiv;</mi> <mo>,</mo> <msub> <mi>X</mi> <mn>1</mn> </msub> <mrow> <mo>(</mo> <mi>&amp;tau;</mi> <mo>,</mo> <mi>k</mi> <mo>)</mo> </mrow> <mo>&amp;NotEqual;</mo> <mn>0</mn> <mo>}</mo> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>3</mn> <mo>)</mo> </mrow> </mrow>
在式(3)中，Im(x)表示x的虚部，应用其他观测信号的STFT代替X1(τ,k)作为分母；
在时频点(τs,ks)，如果满足式(4)，则(τs,ks)∈Ss，
<mrow> <mo>|</mo> <mo>|</mo> <mi>Im</mi> <mrow> <mo>(</mo> <mfrac> <mrow> <mi>X</mi> <mrow> <mo>(</mo> <mi>&amp;tau;</mi> <mo>,</mo> <mi>k</mi> <mo>)</mo> </mrow> </mrow> <mrow> <msub> <mi>X</mi> <mn>1</mn> </msub> <mrow> <mo>(</mo> <mi>&amp;tau;</mi> <mo>,</mo> <mi>k</mi> <mo>)</mo> </mrow> </mrow> </mfrac> <mo>)</mo> </mrow> <mo>|</mo> <mo>|</mo> <mo>=</mo> <mo>|</mo> <mo>|</mo> <mi>Im</mi> <mrow> <mo>(</mo> <msup> <mfenced open = "[" close = "]"> <mtable> <mtr> <mtd> <mfrac> <mrow> <msub> <mi>a</mi> <mrow> <mn>1</mn> <mi>i</mi> </mrow> </msub> <msub> <mi>S</mi> <mi>i</mi> </msub> <mrow> <mo>(</mo> <msub> <mi>&amp;tau;</mi> <mrow> <mi>s</mi> <mo>,</mo> </mrow> </msub> <msub> <mi>k</mi> <mi>s</mi> </msub> <mo>)</mo> </mrow> </mrow> <mrow> <msub> <mi>a</mi> <mrow> <mn>1</mn> <mi>i</mi> </mrow> </msub> <msub> <mi>S</mi> <mi>i</mi> </msub> <mrow> <mo>(</mo> <msub> <mi>&amp;tau;</mi> <mrow> <mi>s</mi> <mo>,</mo> </mrow> </msub> <msub> <mi>k</mi> <mi>s</mi> </msub> <mo>)</mo> </mrow> </mrow> </mfrac> <mfrac> <mrow> <msub> <mi>a</mi> <mrow> <mn>2</mn> <mi>i</mi> </mrow> </msub> <msub> <mi>S</mi> <mi>i</mi> </msub> <mrow> <mo>(</mo> <msub> <mi>&amp;tau;</mi> <mrow> <mi>s</mi> <mo>,</mo> </mrow> </msub> <msub> <mi>k</mi> <mi>s</mi> </msub> <mo>)</mo> </mrow> </mrow> <mrow> <msub> <mi>a</mi> <mrow> <mn>1</mn> <mi>i</mi> </mrow> </msub> <msub> <mi>S</mi> <mi>i</mi> </msub> <mrow> <mo>(</mo> <msub> <mi>&amp;tau;</mi> <mrow> <mi>s</mi> <mo>,</mo> </mrow> </msub> <msub> <mi>k</mi> <mi>s</mi> </msub> <mo>)</mo> </mrow> </mrow> </mfrac> </mtd> </mtr> <mtr> <mtd> <mo>...</mo> <mfrac> <mrow> <msub> <mi>a</mi> <mrow> <mi>m</mi> <mi>i</mi> </mrow> </msub> <msub> <mi>S</mi> <mi>i</mi> </msub> <mrow> <mo>(</mo> <msub> <mi>&amp;tau;</mi> <mrow> <mi>s</mi> <mo>,</mo> </mrow> </msub> <msub> <mi>k</mi> <mi>s</mi> </msub> <mo>)</mo> </mrow> </mrow> <mrow> <msub> <mi>a</mi> <mrow> <mn>1</mn> <mi>i</mi> </mrow> </msub> <msub> <mi>S</mi> <mi>i</mi> </msub> <mrow> <mo>(</mo> <msub> <mi>&amp;tau;</mi> <mrow> <mi>s</mi> <mo>,</mo> </mrow> </msub> <msub> <mi>k</mi> <mi>s</mi> </msub> <mo>)</mo> </mrow> </mrow> </mfrac> </mtd> </mtr> </mtable> </mfenced> <mi>T</mi> </msup> <mo>)</mo> </mrow> <mo>|</mo> <mo>|</mo> <mo>=</mo> <mo>|</mo> <mo>|</mo> <mi>Im</mi> <mrow> <mo>(</mo> <msup> <mfenced open = "[" close = "]"> <mtable> <mtr> <mtd> <mn>1</mn> </mtd> <mtd> <mfrac> <msub> <mi>a</mi> <mrow> <mn>2</mn> <mi>i</mi> </mrow> </msub> <msub> <mi>a</mi> <mrow> <mn>1</mn> <mi>i</mi> </mrow> </msub> </mfrac> </mtd> <mtd> <mo>...</mo> </mtd> <mtd> <mfrac> <msub> <mi>a</mi> <mrow> <mi>m</mi> <mi>i</mi> </mrow> </msub> <msub> <mi>a</mi> <mrow> <mn>1</mn> <mi>i</mi> </mrow> </msub> </mfrac> </mtd> </mtr> </mtable> </mfenced> <mi>T</mi> </msup> <mo>)</mo> </mrow> <mo>|</mo> <mo>|</mo> <mo>=</mo> <mn>0</mn> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>4</mn> <mo>)</mo> </mrow> </mrow>
其中，ali为混合矩阵A的第(l,i)个元素。由式(4)可判断在时频点(τs,ks)处只存在第i个源信号。相反，如果第i个和第j个源信号Si(τ,k)和Sj(τ,k)在时频点(τn,kn)处同时存在，则时频单源点集合Ss中将不包括时频点(τn,kn)，因为：
<mrow> <mo>|</mo> <mo>|</mo> <mi>Im</mi> <mrow> <mo>(</mo> <mfrac> <mrow> <mi>X</mi> <mrow> <mo>(</mo> <msub> <mi>&amp;tau;</mi> <mi>n</mi> </msub> <mo>,</mo> <msub> <mi>k</mi> <mi>n</mi> </msub> <mo>)</mo> </mrow> </mrow> <mrow> <msub> <mi>X</mi> <mn>1</mn> </msub> <mrow> <mo>(</mo> <msub> <mi>&amp;tau;</mi> <mi>n</mi> </msub> <mo>,</mo> <msub> <mi>k</mi> <mi>n</mi> </msub> <mo>)</mo> </mrow> </mrow> </mfrac> <mo>)</mo> </mrow> <mo>|</mo> <mo>|</mo> <mo>=</mo> <mo>|</mo> <mo>|</mo> <mi>Im</mi> <mo>(</mo> <msup> <mfenced open = "[" close = "]"> <mtable> <mtr> <mtd> <mtable> <mtr> <mtd> <mn>1</mn> </mtd> <mtd> <mfrac> <mrow> <msub> <mi>a</mi> <mrow> <mn>2</mn> <mi>i</mi> </mrow> </msub> <msub> <mi>S</mi> <mi>i</mi> </msub> <mrow> <mo>(</mo> <msub> <mi>&amp;tau;</mi> <mi>n</mi> </msub> <mo>,</mo> <msub> <mi>k</mi> <mi>n</mi> </msub> <mo>)</mo> </mrow> <mo>+</mo> <msub> <mi>a</mi> <mrow> <mn>2</mn> <mi>j</mi> </mrow> </msub> <msub> <mi>S</mi> <mi>j</mi> </msub> <mrow> <mo>(</mo> <msub> <mi>&amp;tau;</mi> <mi>n</mi> </msub> <mo>,</mo> <msub> <mi>k</mi> <mi>n</mi> </msub> <mo>)</mo> </mrow> </mrow> <mrow> <msub> <mi>a</mi> <mrow> <mn>1</mn> <mi>i</mi> </mrow> </msub> <msub> <mi>S</mi> <mi>i</mi> </msub> <mrow> <mo>(</mo> <msub> <mi>&amp;tau;</mi> <mi>n</mi> </msub> <mo>,</mo> <msub> <mi>k</mi> <mi>n</mi> </msub> <mo>)</mo> </mrow> <mo>+</mo> <msub> <mi>a</mi> <mrow> <mn>1</mn> <mi>j</mi> </mrow> </msub> <msub> <mi>S</mi> <mi>j</mi> </msub> <mrow> <mo>(</mo> <msub> <mi>&amp;tau;</mi> <mi>n</mi> </msub> <mo>,</mo> <msub> <mi>k</mi> <mi>n</mi> </msub> <mo>)</mo> </mrow> </mrow> </mfrac> </mtd> </mtr> </mtable> </mtd> </mtr> <mtr> <mtd> <mo>...</mo> <mfrac> <mrow> <msub> <mi>a</mi> <mrow> <mi>m</mi> <mi>i</mi> </mrow> </msub> <msub> <mi>S</mi> <mi>i</mi> </msub> <mrow> <mo>(</mo> <msub> <mi>&amp;tau;</mi> <mi>n</mi> </msub> <mo>,</mo> <msub> <mi>k</mi> <mi>n</mi> </msub> <mo>)</mo> </mrow> <mo>+</mo> <msub> <mi>a</mi> <mrow> <mi>m</mi> <mi>j</mi> </mrow> </msub> <msub> <mi>S</mi> <mi>j</mi> </msub> <mrow> <mo>(</mo> <msub> <mi>&amp;tau;</mi> <mi>n</mi> </msub> <mo>,</mo> <msub> <mi>k</mi> <mi>n</mi> </msub> <mo>)</mo> </mrow> </mrow> <mrow> <msub> <mi>a</mi> <mrow> <mn>1</mn> <mi>i</mi> </mrow> </msub> <msub> <mi>S</mi> <mi>i</mi> </msub> <mrow> <mo>(</mo> <msub> <mi>&amp;tau;</mi> <mi>n</mi> </msub> <mo>,</mo> <msub> <mi>k</mi> <mi>n</mi> </msub> <mo>)</mo> </mrow> <mo>+</mo> <msub> <mi>a</mi> <mrow> <mn>1</mn> <mi>j</mi> </mrow> </msub> <msub> <mi>S</mi> <mi>j</mi> </msub> <mrow> <mo>(</mo> <msub> <mi>&amp;tau;</mi> <mi>n</mi> </msub> <mo>,</mo> <msub> <mi>k</mi> <mi>n</mi> </msub> <mo>)</mo> </mrow> </mrow> </mfrac> </mtd> </mtr> </mtable> </mfenced> <mi>T</mi> </msup> <mo>)</mo> <mo>|</mo> <mo>|</mo> <mo>&gt;</mo> <mi>&amp;epsiv;</mi> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>5</mn> <mo>)</mo> </mrow> </mrow>
式(5)表示Si(τ,k)和Sj(τ,k)为复数，不能如式(4)一样从分子和分母中消去；
2)应用聚类算法，基于式(6)所示比率向量，将时频单源点集合Ss聚成m类
<mrow> <mi>Re</mi> <mrow> <mo>(</mo> <mfrac> <mrow> <mi>X</mi> <mrow> <mo>(</mo> <mi>&amp;tau;</mi> <mo>,</mo> <mi>k</mi> <mo>)</mo> </mrow> </mrow> <mrow> <msub> <mi>X</mi> <mn>1</mn> </msub> <mrow> <mo>(</mo> <mi>&amp;tau;</mi> <mo>,</mo> <mi>k</mi> <mo>)</mo> </mrow> </mrow> </mfrac> <mo>)</mo> </mrow> <mo>,</mo> <mo>&amp;ForAll;</mo> <mrow> <mo>(</mo> <mi>&amp;tau;</mi> <mo>,</mo> <mi>k</mi> <mo>)</mo> </mrow> <mo>&amp;Element;</mo> <msub> <mi>S</mi> <mi>s</mi> </msub> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>6</mn> <mo>)</mo> </mrow> </mrow>
其中，Re(x)表示x的实部；应用k均值聚类算法对Ss进行聚类；在Ss中，若第i个源信号的时频点活跃，则有以下比率向量：
<mrow> <mi>Re</mi> <mrow> <mo>(</mo> <mfrac> <mrow> <mi>X</mi> <mrow> <mo>(</mo> <mi>&amp;tau;</mi> <mo>,</mo> <mi>k</mi> <mo>)</mo> </mrow> </mrow> <mrow> <msub> <mi>X</mi> <mn>1</mn> </msub> <mrow> <mo>(</mo> <mi>&amp;tau;</mi> <mo>,</mo> <mi>k</mi> <mo>)</mo> </mrow> </mrow> </mfrac> <mo>)</mo> </mrow> <mo>=</mo> <mi>Re</mi> <mrow> <mo>(</mo> <msup> <mfenced open = "[" close = "]"> <mtable> <mtr> <mtd> <mn>1</mn> </mtd> <mtd> <mfrac> <msub> <mi>a</mi> <mrow> <mn>2</mn> <mi>i</mi> </mrow> </msub> <msub> <mi>a</mi> <mrow> <mn>1</mn> <mi>i</mi> </mrow> </msub> </mfrac> </mtd> <mtd> <mo>...</mo> </mtd> <mtd> <mfrac> <msub> <mi>a</mi> <mrow> <mi>m</mi> <mi>i</mi> </mrow> </msub> <msub> <mi>a</mi> <mrow> <mn>1</mn> <mi>i</mi> </mrow> </msub> </mfrac> </mtd> </mtr> </mtable> </mfenced> <mi>T</mi> </msup> <mo>)</mo> </mrow> <mo>=</mo> <msup> <mfenced open = "[" close = "]"> <mtable> <mtr> <mtd> <mn>1</mn> </mtd> <mtd> <mfrac> <msub> <mi>a</mi> <mrow> <mn>2</mn> <mi>i</mi> </mrow> </msub> <msub> <mi>a</mi> <mrow> <mn>1</mn> <mi>i</mi> </mrow> </msub> </mfrac> </mtd> <mtd> <mo>...</mo> </mtd> <mtd> <mfrac> <msub> <mi>a</mi> <mrow> <mi>m</mi> <mi>i</mi> </mrow> </msub> <msub> <mi>a</mi> <mrow> <mn>1</mn> <mi>i</mi> </mrow> </msub> </mfrac> </mtd> </mtr> </mtable> </mfenced> <mi>T</mi> </msup> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>7</mn> <mo>)</mo> </mrow> </mrow>
基于上述比率向量的时频单源点集合表示为i＝1,2,…,m；应用上述算法，检测得到每一个辐射源信号对应的单源点区域；
3)基于所有辐射源信号的时频单源点集合估计得到混合矩阵A，在混合矩阵估计之前，所有观测信号的实部和虚部都通过l2范数进行归一化；
基于第i个源信号活跃的时频单源点集合i＝1，2，…，m，通过式(8)估计得到混合矩阵的第i列
<mrow> <msub> <mover> <mi>a</mi> <mo>^</mo> </mover> <mi>i</mi> </msub> <mo>=</mo> <mfrac> <mn>1</mn> <mrow> <mn>2</mn> <mrow> <mo>|</mo> <msub> <mi>S</mi> <msub> <mi>c</mi> <mi>i</mi> </msub> </msub> <mo>|</mo> </mrow> </mrow> </mfrac> <munder> <mo>&amp;Sigma;</mo> <mrow> <mo>(</mo> <mi>&amp;tau;</mi> <mo>,</mo> <mi>k</mi> <mo>)</mo> <mo>&amp;Element;</mo> <msub> <mi>S</mi> <msub> <mi>c</mi> <mi>i</mi> </msub> </msub> </mrow> </munder> <mi>Re</mi> <mo>&amp;lsqb;</mo> <mi>X</mi> <mrow> <mo>(</mo> <mi>&amp;tau;</mi> <mo>,</mo> <mi>k</mi> <mo>)</mo> </mrow> <mo>&amp;rsqb;</mo> <mo>+</mo> <mi>Im</mi> <mo>&amp;lsqb;</mo> <mi>X</mi> <mrow> <mo>(</mo> <mi>&amp;tau;</mi> <mo>,</mo> <mi>k</mi> <mo>)</mo> </mrow> <mo>&amp;rsqb;</mo> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>8</mn> <mo>)</mo> </mrow> </mrow>
(3)基于向量夹角检测的混合矩阵优化方法
检测得到时频单源点集合以后，应用K均值聚类算法进行混合矩阵估计，K均值聚类算法具有计算量小，收敛快的优势，但是聚类效果严重依赖于初始聚类中心的选择；对于相同的观测信号，即便时频单源点检测精度很高，由于K均值聚类法每次自动选择的初始聚类中心是不确定的，导致每次估计得到的混合矩阵相应地是变化的，每次估计得到的混合矩阵可能只有一部分或全部列向量与真正的混合矩阵相似；由此分析可知，基于时频单源点检测的混合矩阵估计效果受K均值聚类算法影响很大；
为了提高混合矩阵估计的精度，采用基于列向量夹角检测的混合矩阵优化算法，通过多次估计混合矩阵，对每次估计的矩阵进行列向量的相关性检测，出现次数频繁的列向量，是真正的混合矩阵的列向量；混合矩阵列向量顺序的改变，以及每一个列向量正负符号的改变，对混合矩阵的正确性都是不影响的，具体步骤如下：
1)对相同的观测信号，进行时频变换，设定时频单源点检测阈值，检测得到时频单源点集合；
2)多次估计混合矩阵，每次估计的结果写为Ae1，Ae2，…，Aek，k为设定的混合矩阵估计的最大次数。设混合矩阵的维数为m×n，则估计得到的混合矩阵维数也为m×n，将Ae1，Ae2，…，Aek组合成为一个新的矩阵并转置，记为W，则W的维数为N×m，N＝n×k；
3)计算矩阵W的行向量的夹角，设计算结果表示成矩阵We，显然We是N×N维的矩阵，该矩阵沿主对角线，右上半与左下半是对称的，包含的信息一样；可以令左下半元素等于0，只考虑右上半矩阵，设矩阵为We1；
4)令We1主对角线右上部元素大于10的元素等于0，矩阵记为We2。根据We2对同一或相似，这里设定列向量夹角小于10的列向量出现的次数进行统计，根据给定的辐射源个数，选择出现次数频繁的列向量作为预估混合矩阵的列向量，组成一个新的估计得到的混合矩阵Ae，其维数为m×n。
CN201711236736.9A 2017-11-30 2017-11-30 一种基于小波分析和时频单源检测的混合矩阵估计方法 Active CN108108666B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201711236736.9A CN108108666B (zh) 2017-11-30 2017-11-30 一种基于小波分析和时频单源检测的混合矩阵估计方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201711236736.9A CN108108666B (zh) 2017-11-30 2017-11-30 一种基于小波分析和时频单源检测的混合矩阵估计方法

Publications (2)

Publication Number Publication Date
CN108108666A true CN108108666A (zh) 2018-06-01
CN108108666B CN108108666B (zh) 2020-09-04

Family

ID=62208035

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201711236736.9A Active CN108108666B (zh) 2017-11-30 2017-11-30 一种基于小波分析和时频单源检测的混合矩阵估计方法

Country Status (1)

Country Link
CN (1) CN108108666B (zh)

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109344754A (zh) * 2018-09-21 2019-02-15 电子信息系统复杂电磁环境效应国家重点实验室 一种改进式最短路径欠定源信号恢复方法
CN110060698A (zh) * 2019-04-11 2019-07-26 哈尔滨工程大学 一种基于改进势函数的语音信号混合矩阵估计方法
CN114925726A (zh) * 2022-05-13 2022-08-19 华侨大学 亚采样工作模态参数的识别方法、装置、设备和存储介质
CN115363586A (zh) * 2022-09-08 2022-11-22 山东大学 一种基于脉搏波信号的心理压力等级评估系统及方法

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104332161A (zh) * 2014-09-28 2015-02-04 武汉理工大学 一种基于接收先验和单源点检测的欠定盲辨识方法
US8958750B1 (en) * 2013-09-12 2015-02-17 King Fahd University Of Petroleum And Minerals Peak detection method using blind source separation
CN106371070A (zh) * 2016-08-30 2017-02-01 电子信息系统复杂电磁环境效应国家重点实验室 一种改进的基于小波分析的欠定盲源分离源数估计方法
CN106778001A (zh) * 2016-12-26 2017-05-31 西安电子科技大学 基于改进时频单源区的欠定混合矩阵盲估计方法

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US8958750B1 (en) * 2013-09-12 2015-02-17 King Fahd University Of Petroleum And Minerals Peak detection method using blind source separation
CN104332161A (zh) * 2014-09-28 2015-02-04 武汉理工大学 一种基于接收先验和单源点检测的欠定盲辨识方法
CN106371070A (zh) * 2016-08-30 2017-02-01 电子信息系统复杂电磁环境效应国家重点实验室 一种改进的基于小波分析的欠定盲源分离源数估计方法
CN106778001A (zh) * 2016-12-26 2017-05-31 西安电子科技大学 基于改进时频单源区的欠定混合矩阵盲估计方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
M.Y.ABBASS: "Blind source separation with wavelet based ICA technique using kurtosis", 《ICCTA 2013》 *
陈永强: "欠定信号分离及其在语音处理中的应用", 《万方学位论文数据库》 *

Cited By (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109344754A (zh) * 2018-09-21 2019-02-15 电子信息系统复杂电磁环境效应国家重点实验室 一种改进式最短路径欠定源信号恢复方法
CN110060698A (zh) * 2019-04-11 2019-07-26 哈尔滨工程大学 一种基于改进势函数的语音信号混合矩阵估计方法
CN114925726A (zh) * 2022-05-13 2022-08-19 华侨大学 亚采样工作模态参数的识别方法、装置、设备和存储介质
CN114925726B (zh) * 2022-05-13 2024-07-26 华侨大学 亚采样工作模态参数的识别方法、装置、设备和存储介质
CN115363586A (zh) * 2022-09-08 2022-11-22 山东大学 一种基于脉搏波信号的心理压力等级评估系统及方法

Also Published As

Publication number Publication date
CN108108666B (zh) 2020-09-04

Similar Documents

Publication Publication Date Title
CN108108666A (zh) 一种基于小波分析和时频单源检测的混合矩阵估计方法
CN103728551B (zh) 一种基于级联集成分类器的模拟电路故障诊断方法
CN105488466B (zh) 一种深层神经网络和水声目标声纹特征提取方法
CN103245907B (zh) 一种模拟电路故障诊断方法
CN104375976B (zh) 基于张量正则分解的欠定盲源分离中的混合矩阵识别方法
CN109870729B (zh) 基于离散余弦变换的深度神经网络磁共振信号消噪方法
CN101587186B (zh) 一种雷达脉内调制信号的特征提取方法
CN103854660B (zh) 一种基于独立成分分析的四麦克语音增强方法
CN110084148A (zh) 一种高压断路器机械故障诊断方法
CN104459668A (zh) 基于深度学习网络的雷达目标识别方法
CN112329914B (zh) 地埋式变电站的故障诊断方法、装置及电子设备
CN109471074A (zh) 基于奇异值分解与一维cnn网络的雷达辐射源识别方法
CN102279358A (zh) 一种基于mcskpca的神经网络模拟电路故障诊断方法
Qin et al. Radar waveform recognition based on deep residual network
CN104793124A (zh) 基于小波变换和ica特征提取的开关电路故障诊断方法
CN112613493A (zh) 一种基于多特征提取和woa-elm的滚动轴承故障诊断方法
CN104504403A (zh) 一种基于散射变换的旋转机械故障预测方法
CN114428234A (zh) 基于gan和自注意力的雷达高分辨距离像降噪识别方法
CN114595732A (zh) 基于深度聚类的雷达辐射源分选方法
Chen et al. Automatic modulation classification of radar signals utilizing X-net
CN110716532A (zh) 一种基于小波包能量与fft的水下机器人推进器弱故障辨识方法
CN113704688A (zh) 基于变分贝叶斯平行因子分解的缺失振动信号的恢复方法
CN107132518B (zh) 一种基于稀疏表示和时频特征的距离扩展目标检测方法
Shi et al. Intelligent fault diagnosis of bearings based on feature model and Alexnet neural network
CN116337449A (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