量，k = 1，...，Κ，Κ表示目标矩阵集合中的矩阵总数，掠表示与降维后目标矩阵集合島中的第k个 矩阵Gk对应的第2个位置向量，热和疏分别表示降维后目标矩阵集合谷中的第k个矩阵Gk 的第P行第q列元素的实部和虚部，威和絲分别表示降维后目标矩阵集合e中的第k个矩阵 Gk的第q行第P列元素的实部和虚部，T表示转置操作。
[0121 ]第4步，按照下式，计算海森矩阵Q和两个中间向量：
[0125]其中，h表示第一个中间向量，V表示第二个中间向量，巧表示与降维后目标矩阵集 合g中的第k个矩阵Gk对应的第1个位置向量，k=l，. . .，Κ，Κ表示目标矩阵集合中的矩阵总 数，按表示与降维后目标矩阵集合巧中的第k个矩阵Gk对应的第2个位置向量，疏和疏分别 表示降维后目标矩阵集合白中的第k个矩阵Gk的第P个对角元素的实部和虚部，謀和热分 别表示降维后目标矩阵集合谷中的第k个矩阵Gk的第q个对角元素的实部和虚部，X表示相 乘操作，T表示转置操作，Σ表示求和操作。
[01 %]第5步，按照下式，计算两个元素向量：
[0127] Wpq = -Q-ixh [012 引 Wqp = -DXQ-iXv
[0129] 其中，wpq表示更新矩阵W的第P行第q列元素对应的元素向量，Q表示海森矩阵，h表 示第一个中间向量，wqp表示更新矩阵W的第q行第P列元素对应的元素向量，D表示对角矩阵，
I V表示第二个中间向量，（·尸表示求逆操作，X表示相乘操作。
[0130] 第6步，获取更新矩阵W的第P行第q列元素 wpq和W的第q行第P列元素 wqp的最优解：
[0131] Wpq的实部取更新矩阵W的第P行第q列元素对应的元素向量Wpq的第一个元素，Wpq的 虚部取更新矩阵W的第P行第q列元素对应的元素向量Wpq的第二个元素，Wqp的实部取更新矩 阵W的第q行第P列元素对应的元素向量Wqp的第一个元素，Wqp的虚部取更新矩阵W的第q行第 P列元素对应的元素向量Wqp的第二个元素。
[0132] 第7步，将更新矩阵的列数q增加1。
[0133] 第8步，判断当前的q>p-l是否成立，若是，则执行第9步，否则执行第3步，其中，P 表示更新矩阵的行数。
[0134] 第9步，将更新矩阵的行数P增加1。
[0135] 第10步，判断当前的P含N是否成立，若是，则执行第2步，否则终止迭代，输出更新 矩阵W，其中，N表示源信号总数。
[0136] (4h)计算更新矩阵W的弗罗贝尼乌斯范数沪。
[0137] (4i)判断更新矩阵W的弗罗贝尼乌斯范数是否大于1，若是，则执行步骤(4j)，否 贝1J，执行步骤(4k)。
[0138] (4 j)对更新矩阵W归一化。
[0139] (4k)按照下式，用更新矩阵W更新零对角化矩阵B/和降维后目标矩阵集合g中的 所有目标矩阵：
[0142]其中，B"表示更新之后的零对角化矩阵，I表示单位矩阵，W表示更新矩阵，B/表示 更新前的零对角化矩阵，G/k表示更新后的目标矩阵集合α中的第k个矩阵，Gk表示更新前的 目标矩阵集合0中的第4个矩阵^=1，...，1(，1(表示目标矩阵集合中的矩阵总数，（-卢表 示共辆转置操作。
[01创 （41)判断资 <《是否成立，若是，执行步骤(4m)，否则，执行步骤(4g)，其中，與表示 更新矩阵W的弗罗贝尼乌斯范数，δ表示终止迭代阔值，取值为10-3。
[0144] (4m)按照下式，计算分离矩阵Β:
[0145] B = B' XP
[0146] 其中，β/表示降维后目标矩阵集合径的零对角化矩阵，P表示降维矩阵，X表示相 乘操作。
[0147] (5)按照下式，从观测信号中分离源信号：
[014 引 y(t)=BXx(t)
[0149] 其中，y(t)表示W观测时间t为自变量的源信号向量，B表示分离矩阵，x(t)表示步 骤(1)获取的W观测时间t为自变量的观测信号向量，X表示相乘操作。
[0150] 下面结合附图对本发明的效果做进一步的描述。
[0151] 1.仿真条件：
[0152] 本发明的仿真运行系统为Intel(R)Core(TM)i7-4790 CPU@3.60GHz，64位Windows 操作系统，仿真软件采用MTLAB(R2014b)。
[0153] 2.仿真内容与结果分析：
[0154] 本发明的实施例是分别采集自3个电信号传感器的观测信号，运3个观测信号由3 个人工合成的源信号通过与混合矩阵相乘而得到，混合矩阵的元素均随机产生，且服从标 准正态分布。
[0155] 本发明仿真实验中采用的3个人工合成的电信号作为源信号，它们分别为线性调 频信号，正弦调频信号W及正弦信号。运3个源信号均采集128个样本点。由于样本点很少， 它们在时域几乎是统计相关的。图4是运Ξ个源信号的时域波形图W及时频域魏格纳-瑞利 分布图。其中，图4的左上、中上和右上方的子图分别是线性调频信号，正弦调频信号W及正 弦信号的时域波形图，运Ξ幅子图的横坐标表示时间，纵坐标表示信号幅值；图4的左下、中 下、右下方的子图分别是线性调频信号，正弦调频信号W及正弦信号的时频域魏格纳-瑞利 分布图，运3幅子图的横坐标表示时间，纵坐标表示频率。
[0156] 参照附图5中的3幅子图分别是信噪比为20地时3个观测信号的时频域魏格纳-瑞 利分布图，其中，横坐标表示时间，纵坐标表示频率。从附图5可W看出任意一个观测信号均 不同程度地包含有3个源信号成分，表明观测信号由源信号混合得到。
[0157] 现将本发明方法与其它两种时频域代数方法进行比较，其中方法一为 E.M.Fadaili ,Ν. Thir ion-Moreau , and E. Moreau 在文章 "Nonorthogonal joint diagonalization zero diagonalization for source separation based on time-frequency distributions"（IEEE Tr曰ns.Signal Process.,vol.55,no.5,pp.1673-1687, 2007)中公布的二次型优化非正交联合零对角化迭代方法，方法二为G.Chabriel, J . Barrere , N . Thir ion-Moreau , and E . Moreau 在文章"A1 gebrai c joint zero di曰gon曰liz曰tion 曰nd blind sources sep曰r曰tion"（IEEE Tr曰ns.Signsl Process., vol. 56，no. 3，pp.980-989，Mar. 2008)中公布的联合零对角化非迭代方法，运Ξ种方法在每 一个信噪比条件下均进行500次蒙特卡洛仿真。
[0158] 用干扰信号比ISR来衡量盲信号分离方法的性能，干扰信号比ISR的定义式为：
[0159]
[0160] 其中，N表示源信号总数，CU表示全局变换矩阵C的第i行第j列元素 ，C = BA，B表示 分离矩阵，A表示混合矩阵，cii表示全局变换矩阵的第i行第1列元素，CU表示全局变换矩阵 的第1行第i列元素，Σ( ·)表示求和操作，I · I表示取绝对值操作，（·）2表示取平方操作， maxi( ·)表示取所有行下标或列下标为1的元素中最大元素的操作。若ISR趋于ο则表示盲 分离性能较好。
[0161] 参照附图6所示的本发明分离出的3个源信号的时频域魏格纳-瑞利分布图。其中 横坐标表示时间，纵坐标表示频率。从图6与图4的对比可W看出，图6的3幅子图与图4中的3 幅子图形状相同，且不含有其它源信号成分。说明本发明可W有效地恢复所有源信号，而且 恢复出的源信号不再含有其它源信号的成分，盲信号分离效果较好。
[0162] 参照附图7所示的本发明方法与方法一和方法二在不同信噪比下的盲信号分离性 能对比曲线图。其中横坐标表示信噪比，纵坐标表示干扰信号比。从图中所示的干扰信号比 随信噪比变化的曲线可W看出，本发明方法在所有信噪比下的干扰信号比都比其它两种方 法更加接近0,说明本发明方法具有更好的盲信号分离性能。
【主权项】
1. 一种基于联合零对角化的时频域盲信号分离方法，包括如下步骤： (1) 获取观测信号： (la) 从电信号传感器中接收来自N个源信号的Μ个观测信号，其中，每一个传感器发出 一个观测信号，M2 Ν; (lb) 将观测信号按电信号传感器的编号顺序，依次排列成观测信号向量X(t) = [X1 (t)，...，XM(t)]T，其中，X1(t)，XM(t)分别表示以观测时间t为自变量的第1个观测信号和第 Μ个观测信号，Μ表示观测信号总数，T表示转置操作； (2) 获得观测信号向量的时频分布矩阵： 计算观测信号向量x(t)在时频点（tut)上的魏格纳-瑞利分布矩阵，将计算的结果作 为观测信号向量X (t)在该时频点(ti，f j)上的时频分布矩阵Rij，其中，ti表示时频点中的第i 个时间点，f j表示时频点中的第j个频率点，i = 1，. . .，I，j = 1，. . .，J，I表示由观测信号确 定的时频点中的总时间点数，J表示由观测信号确定的时频点中的总频率点数； (3) 构造待估计的目标矩阵集合： 选取具有可零对角化结构的时频分布矩阵，组成待估计的目标矩阵集合m; (4) 使用乘积型联合零对角化方法估计分离矩阵： (4a)判断观测信号的总数Μ是否大于源信号的总数N，若是，则执行步骤(4b)，否则，执