t)]T，则观测信号与源信号的 关系为：
[0056] x(t)=As(t)+n(t)
[0057] 其中，si(t)，SN(t)分别表示W观测时间t为自变量的第1个源信号和第N个观测信 号，N表示观测信号总数，A表示未知混合矩阵，n(t)W观测时间t为自变量的表示传感器噪 声，T表示转置操作。
[0058] 步骤2,获得观测信号向量的时频分布矩阵。
[0059] 计算观测信号向量x(t)在时频点（ti，。)上的魏格纳-瑞利分布矩阵，将计算的结 果作为观测信号向量X(t)在该时频点（ti，。)上的时频分布矩阵Ru，其中，ti表示时频点中 的第i个时间点，fj表示时频点中的第j个频率点，i = l,. ..，1，j = l,. ..，J，I表示由观测信 号确定的时频点中的总时间点数，J表示由观测信号确定的时频点中的总频率点数。
[0060] 计算观测信号向量X (t)在时频点（ti，f j)上的魏格纳-瑞利分布矩阵的具体步骤如 下：
[0061] 第一步，对每个观测信号做希尔伯特变换，获得每个观测信号的解析信号。
[0062] 希尔伯特变换数学公式如下：
[0063]
[0064] 其中，氧,,倘表示第m个观测信号在观测时刻为t时的希尔伯特变换信号，t表示观 测时刻，表示第m个观测信号在延迟时间为T时观测信号，T表示观测信号相对t时刻 的延迟时间，m=l，...，M，M表示观测信号总数，表示积分操作。
[0065] 按照下式，获得观测信号的解析信号：
[0066] (t) = ,?；," (t) + jXm (t)
[0067] 其中，ym(t)表示第m个观测信号xm(t)的解析信号，m=l，. . .，M，M表示观测信号总 数，t表示观测时刻，Xm(t)表示第m个观测信号，.Mt)表示第m个观测信号Xm(t)的希尔伯特 变换信号，J表示虚数单位。
[0068] 第二步，按照下式，计算时频点(ti山)上的第m个解析信号的自魏格纳-瑞利分布：
[0069]
[0070] 其中，表示第m个解析信号在时频点上的自魏格纳瑞利分布，m = 1，...，M，M表示观测信号总数，ti表示时频点中的第i个时间点，表示时频点中的第j个频 率点，+ f)表示相对于观测时刻t超前了责时的第m个解析信号，τ表示观测信号相对于 时刻t超前的时间，t表示观测信号的观测时刻，f表示观测信号的频率，j表示虚数单位，e 表示自然常数，约为2.718，（-)^表示取共辆操作，1(-)表示积分操作。
[0071] 第Ξ步，按照下式，计算时频点（ti，。)上的第m个解析信号与第η个信号的互魏格 纳-瑞利分布；
[0072]
[0073] 其中，rmn(ti，fj)表示时频点上的第m个解析信号与第η个信号的互魏格纳-瑞利分布，m=l, . . .，Μ，η=1, . . .，M，m辛η，Μ表示观测信号总数，ti表示时频点中的第i个时 间点，。表示时频点中的第j个频率点，抓化牛荀表示相对于观测时亥Ijt超前了f时的第m个 解析信号，-号)表示相对于观测时刻t滞后了 Y时的第η个解析信号的共辆，T表示观测 信号相对于时刻t超前的时间，t表示观测信号的观测时刻，f表示观测信号的频率，？表示 虚数单位，6表示自然常数，约为2.718，（-)^表示取共辆操作，1(-)表示积分操作。
[0074] 第四步，将时频点（ti，。）上计算得到的自魏格纳-瑞利分布和互魏格纳-瑞利分 布，按如下方式排列成该时频点上的时频分布矩阵Rij:
[00对 Ru的第m个对角元素为时频点（ti，。）上的第m个解析信号的自魏格纳-瑞利分布 r?(ti山），Ru的第m行η列的元素为时频点（ti，。化的第m个解析信号与第η个解析信号的 互魏格纳-瑞利分布打in(ti,f j)，其中，m= 1,...,Μ，η=1,... ,M，m辛η，Μ表示观测信号的总 数。
[0076] 步骤3,构造待估计的目标矩阵集合。
[0077] 选取具有可零对角化结构的时频分布矩阵，组成待估计的目标矩阵集合化。图2是 步骤3构造待估计的目标矩阵集合的流程图，图2中所示的流程对应W下具体步骤：
[0078] 第1步，使用主分量分析方法估计噪声方差和白化矩阵S，其中，0表示噪声标准 差，（·）2表示取平方操作。
[0079] 第2步，将时频点中的时间点数i的初值设置为1。
[0080] 第3步，将时频点中的频率点数j的初值设置为1。
[0081] 第4步，按照下式，对时频点（ti，。)上的时频分布矩阵Ru进行去除噪声操作，得到 去噪后时频分布矩阵：
[0082]
[00削其中，择。，表示时频点上的去噪后时频分布矩阵，t康示时频点中的第i个 时间点，表示时频点中的第j个频率点，σ表示噪声标准差，（·）2表示取平方操作，I表示单 位矩阵，X表示相乘操作。
[0084] 第5步，按照下式，对时频点上的去噪后时频分布矩阵反,3进行白化处理，得 到白化后时频分布矩阵：
[0085]
[0086] 其中，瓦3表示时频点（ti，。)上的白化后时频分布矩阵，t康示时频点中的第i个 时间点，表示时频点中的第j个频率点，（·）H表示共辆转置操作，X表示相乘操作。
[0087] 第6步，计算时频点（ti，。)上的白化后时频分布矩阵反,,虚部阵的弗罗贝尼乌斯范 数^，其中，t康示时频点中的第i个时间点，f康示时频点中的第j个频率点。
[008引第7步，计算时频点（ti，f j化的白化后时频分布矩阵瓦,的迹的绝对值柏，其中，ti 表示时频点中的第i个时间点，fj表示时频点中的第j个频率点。
[0089] 第8步，判断〇^>6且01^< 丫是否成立，若是，则执行第9步，否则，执行第10步，其 中，曰1徒示时频点化的白化后时频分布矩阵私虚部阵的弗罗贝尼乌斯范数，βι徒 示白化后时频分布矩阵瓦的迹的绝对值，ti表示时频点中的第i个时间点，fj表示时频点中 的第j个频率点，ε表示白化后时频分布矩阵虚部阵的弗罗贝尼乌斯范数的阔值，γ表示时 频点上的白化后时频分布矩阵的迹的绝对值的阔值。
[0090] 第9步，将时频点(ti，。)上的白化后时频分布矩阵氏,加入待估计的目标矩阵集合 巧，.其中，ti表示时频点中的第i个时间点，fj表示时频点中的第j个频率点。
[0091 ]第10步，将时频点中的频率点数j增加1。
[0092] 第11步，判断当前j的值是否小于等于时频点中的总频率点数J，若是，执行第4步， 否则，执行第12步。
[0093] 第12步，将时频点中的时间点数i增加1。
[0094] 第13步，判断当前i的值是否小于等于时频点中的总时间点数I，若是，执行第4步， 否则结束循环，输出待估计的目标矩阵集合记。
[00M]步骤4,使用乘积型联合零对角化方法估计分离矩阵。
[0096] 图3是步骤4使用乘积型联合零对角化方法估计分离矩阵的流程图，图3中所示的 流程对应W下具体步骤：
[0097] (4a)判断观测信号的总数Μ是否大于源信号的总数N，若是，则执行步骤(4b)，否 贝1J，执行步骤(4c)。
[0098] (4b)使用子空间方法，获得降维矩阵P，具体步骤为：
[0099] 第一步，按照下式，计算目标矩阵集合取中的所有矩阵的平均值矩阵：
[0100]
[0101] 其中，化表示目标矩阵集合中的第k个目标矩阵，k=l，...，K，K表示目标矩阵集合 中矩阵的总数，Σ表示求和操作。
[0102] 第二步，按照下式，对平均值矩阵R进行奇异值分解，得到平均值矩阵R的非零奇异 值和左奇异向量矩阵U:
[0103] R = UXAXyH
[0104] 其中，Λ为对角矩阵，由平均值矩阵R的非零奇异值组成，U为非零奇异值对应的左 奇异向量矩阵，V为非零奇异值对应的右奇异向量矩阵，X表示相乘操作。
[0105] 第Ξ步，按照下式，计算降维矩阵：
[0106] p = ljH
[0107] 其中，P表示降维矩阵，U表示平均值矩阵R的左奇异向量矩阵，（·）H表示共辆转置 操作。
[0108] (4c)设置降维矩阵P为腺隹的单位矩阵，其中N表示源信号总数。
[0109] (4d)按照下式，对目标矩阵集合'反中的目标矩阵分别进行降维处理，得到降维后 目标矩阵集合G二{Gi，...，Ga-}:
[0110] Gk=PX 化 ΧρΗ
[0111] 其中，扣表示降维后目标矩阵集合?中的第k个矩阵，k = 1，. . .，Κ，Κ表示目标矩阵 集合中的矩阵总数，Ρ表示降维矩阵，化表示降维前的目标矩阵集合巧中的第k个矩阵，（· )H表示共辆转置操作，X表示相乘操作。
[0112] (4e)初始化更新矩阵W为腺隹的零对角矩阵，其中N表示源信号总数。
[0113] (4f)初始化降维后目标矩阵集合巧的零对角化矩阵B/为N维的单位矩阵，其中N表 示源信号总数。
[0114] (4g)使用最小二乘方法，估计更新矩阵W的每个非对角元素的最优解，具体实施步 骤为：
[0115] 第1步，将更新矩阵的行数P的初值设置为2。
[0116] 第2步，将更新矩阵的列数q的初值设置为1。
[0117] 第3步，按照下式，构造两个位置向量：
[0120] 其中，好表示与降维后目标矩阵集合9中的第k个矩阵Gk对应的第1个位置向