
简介面向信号处理、机器学习与盲源分离领域的初学者和研究者资源包提供 JADE独立成分分析中的联合近似对角化算法的 Matlab 实现。JADE 利用四阶累积量最大化的思想可从多通道混合信号中恢复独立源且支持复数形式数据处理适用于音频信号分离、脑电图EEG信号分解、通信干扰抑制等场景。压缩包内仅含 1 个 m 文件大小约 4KB即核心函数 jade.m结构紧凑、免安装依赖便于直接运行、阅读测试和二次开发。目前已有 518 人学习浏览。文件中封装了去均值、白化、四阶累积阵构造、特征值分解和联合对角化等关键步骤无需混合先验信息即可实现盲分离读者可直接传入混合信号调用函数也可逐段剖析算法原理结合 ICA 经典论文深入理解非高斯性度量与高阶统计量在盲分离中的作用为课程设计、毕业设计或实际工程预研提供便捷起点。1. 复数 JADE 盲源分离从一份 jade.zip 开始看到jade.zip这个压缩包名不用打开也能猜到大半里面是 JADE 盲源分离的一套实现常见形态是主函数加演示数据甚至可能是 MATLAB 老代码。JADE 全称 Joint Approximate Diagonalization of Eigenmatrices核心动作是把多个矩阵联合近似对角化从而把混合信号中的独立源分离出来。它在复数信号上的价值尤其明显因为复数基带、频谱域处理、阵列接收这几类场景里源信号天然带相位信息直接套实值盲分离代码会漏掉相位这一维。下面的步骤从数学模型推到可运行代码参数怎么调、结果怎么验都会覆盖到。2. 复数盲源分离的统计依据四阶累积量矩阵2.1 为什么二阶统计量不够JADE 要盯四阶累积量盲源分离的观测模型写作 X A SX 是 M×T 的观测矩阵A 是 M×N 混合矩阵S 是 N×T 源信号矩阵。独立成分分析的基本矛盾在于只用协方差 R E[x x^H]最多能做到白化而白化后的信号仍然保留任意酉变换自由度这个自由度恰恰是无法从二阶统计量里读出来的独立性信息。独立性最终要落在四阶累积量上也就是形如Cum(z_i, z_j*, z_k, z_l*)的量。高斯源的四阶累积量恒为零所以 JADE 能分离的前提是源里最多只有一个高斯成分。实际通信信号基本都满足这条QPSK、OFDM、脉冲成形后的复数基带都有明显的非零四阶累积量。复数情形下四阶量同时刻画幅度与相位的联合分布分离效果比实值情形更依赖这个统计量的估计质量。选四阶而不是更高阶的原因也很现实四阶累积量张量有 N^4 个复数元素源数 N8 时是 4096 个N16 时是 65536 个六阶累积量规模到 N^6样本量需求跟着爆炸工程上很难喂饱。2.2 复数四阶累积量的定义与移植陷阱复数源的四阶累积量约定为C_{ijkl} Cum(z_i, z_j*, z_k, z_l*)展开形式是E[z_i z_j* z_k z_l*] - E[z_i z_j*] E[z_k z_l*] - E[z_i z_k*] E[z_j z_l*] - E[z_i z_l*] E[z_j z_k*]注意共轭位置必须固定。实值代码往复数移植时最常见的错误是某个E[z_i z_j]位置没加共轭导致相位信息被当成实部平均掉表现是分离矩阵抖动、谱峰扩散怎么调阈值都没用。有了四阶张量 C_{ijkl}JADE 需要的是一个矩阵集合{Q(M_r)}作用方式定义为Q(M_r)_{ij} Σ_{kl} C_{ijkl} M_{r,kl}每个 M_r 是一个 N×N 的复数矩阵相当于对四阶张量做一次加权投影。M_r 的选取直接决定联合对角化能不能在有限次扫掠里收敛这也是后面要重点看的一个参数。补充一个边界如果源是强非圆信号比如纯实幅调制那上述共轭位置的配置会损失部分统计量需要换非圆扩展版本普通 JADE 在这里会退化。2.3 JADE 的三步框架白化、累积量、联合对角化JADE 的流程可以压成三步。第一步去均值并白化得到零均值、协方差为单位阵的 Z V XV 是白化矩阵。第二步在 Z 上计算四阶累积量张量并从中选出一组有代表性的 M_r构造待对角化矩阵族。第三步用一系列复数 Givens 旋转做联合近似对角化找到一个酉矩阵 U使所有的U^H Q(M_r) U尽量对角。最终分离矩阵为W U^H V分离输出为Y W X。整个算法是批处理的没有学习率也没有非线性激活函数这是 JADE 和 FastICA 在工程习惯上最大的区别。下面用 Python 把这三步完整跑通并且把累积量、旋转参数都写在代码里方便照着改。3. 用 Python 从零实现复数 JADE 的最小可运行版本3.1 构造复数混合信号的最小实验先造三路统计特性差异明显的复数源一路 QPSK 符号一路带缓慢相位漂移的复指数一路稀疏脉冲。三路源互不相同且都非高斯适合做 JADE 的测试对象。import numpy as np rng np.random.default_rng(7) n, T 3, 8000 t np.arange(T) # 源1离散 QPSK幅度恒定、相位四值 s1 rng.choice(np.array([11j, 1-1j, -11j, -1-1j]), sizeT) # 源2复正弦加相位调制相位连续变化 s2 np.exp(1j * (2 * np.pi * 0.03 * t 0.2 * np.sin(2 * np.pi * 0.002 * t))) # 源360 个随机时刻的稀疏复脉冲 pos rng.choice(T, size60, replaceFalse) s3 np.zeros(T, dtypecomplex) s3[pos] rng.standard_normal(60) 1j * rng.standard_normal(60) S np.vstack([s1, s2, s3]) # 随机复数混合矩阵列方向归一化能稳定噪声尺度 A (rng.standard_normal((n, n)) 1j * rng.standard_normal((n, n))) / np.sqrt(2) X A S 0.02 * (rng.standard_normal((n, T)) 1j * rng.standard_normal((n, T)))这段代码里源信号的类型选择决定了 JADE 的辨识条件QPSK 的四阶累积量非零且与另两路完全不同复正弦的相位结构连续稀疏脉冲是高峰度信号。混合矩阵用1/sqrt(2)缩放是为了让实部和虚部方差一致避免某一路在预处理阶段就占据主导。噪声系数 0.02 对应约 34 dB 信噪比足够让累积量估计稳定。3.2 白化复数协方差特征分解与特征值截断白化在 JADE 里不只是预处理它还承担降维和噪声抑制的作用。复数白化必须使用共轭转置并且特征向量矩阵要取共轭转置后与对角阵相乘。def whiten(x, tol_eig1e-6): x x - x.mean(axis1, keepdimsTrue) R x x.conj().T / x.shape[1] evals, evecs np.linalg.eigh(R) keep evals tol_eig * evals.max() V np.diag(1.0 / np.sqrt(evals[keep])) evecs[:, keep].conj().T return V x, Vnp.linalg.eigh对复数厄米矩阵返回升序特征值evecs[:, keep]只保留特征值大于阈值的列。这里的tol_eig是白化截断阈值默认按最大特征值的 1e-6 截断。阈值设得太小会放大噪声维度设得太大则会把能量较弱的源直接丢掉这个参数在第 4 章会专门讲。返回值里V是白化矩阵V x是白化后的信号 ZZ 的协方差在保留子空间内近似为单位阵。3.3 累积量矩阵集合与复数 Jacobi 联合对角化白化后的数据 Z 已经去掉了二阶相关接下来要构造四阶累积量张量。这里直接按 2.2 的展开式用einsum计算避免写多层循环。def compute_cumulants(z): n, T z.shape E2 z z.conj().T / T M4 np.einsum(it,jt,kt,lt-ijkl, z, z.conj(), z, z.conj(), optimizeTrue) / T C4 M4 - np.einsum(ij,kl-ijkl, E2, E2) \ - np.einsum(ik,jl-ijkl, E2, E2) \ - np.einsum(il,jk-ijkl, E2, E2) return C4einsum的四个索引i,j,k,l对应z_i z_j* z_k z_l*的样本均值这个顺序和 2.2 里的定义一一对应。后两行的三个einsum分别减去三个二阶矩乘积项。白化后 E2 理论上就是单位阵但保留通用写法能避免小样本下数值误差累积。复杂度上这个函数最坏是 O(N^4 T)N 小于等于 10 时可以直接跑N 再大建议分块估计或换用累积量矩阵递推。特征矩阵集合的选取沿用 JADE 原版思路把四阶张量展开成 N^2×N^2 矩阵取特征值最大的前 num 个特征向量再还原成 N×N 矩阵每个矩阵做一次厄米对称化。def select_eigematrices(C4, num): n C4.shape[0] Rc C4.reshape(n * n, n * n) Rc (Rc Rc.conj().T) / 2 evals, evecs np.linalg.eigh(Rc) idx np.argsort(np.abs(evals))[::-1][:num] Ms [] for ix in idx: M evecs[:, ix].reshape(n, n) Ms.append((M M.conj().T) / 2) return Ms特征值最大的特征向量对应累积量能量最集中的方向前 num 个已经能覆盖主要统计结构。如果直接取 num N通常能获得串音最小的分离取 2N 到 3N 会提高冗余度但联合对角化耗时也会上升。联合对角化是 JADE 的核心迭代过程。复数 2×2 旋转比实值多一个相位参数每次处理一对坐标 (p, q) 时要解一个三维优化问题。标准做法是把旋转参数映射到三维向量再求一个 3×3 实对称矩阵的最大特征向量def joint_diag(Ms, tol1e-8, max_sweeps50): n Ms[0].shape[0] U np.eye(n, dtypecomplex) for _ in range(max_sweeps): change 0.0 for p in range(n - 1): for q in range(p 1, n): G np.zeros((3, 3)) for M in Ms: g np.array([ M[p, p] - M[q, q], M[p, q] M[q, p], 1j * (M[q, p] - M[p, q]) ]) G np.real(np.outer(g, g.conj())) _, vecs np.linalg.eigh(G) u vecs[:, -1] theta 0.5 * np.arctan2(np.hypot(u[1], u[2]), u[0]) phi np.arctan2(u[2], u[1]) c, s np.cos(theta), np.sin(theta) J np.array([ [c, s * np.exp(1j * phi)], [-s * np.exp(-1j * phi), c] ], dtypecomplex) for k, M in enumerate(Ms): Ms[k] J.conj().T M J U U J change np.abs(s) if change tol: break return U这里的g是复数三维向量第一维是两对角线元素之差第二维是两个非对角元素的和第三维捕获非对角元素的反转差。累加np.real(np.outer(g, g.conj()))得到实对称矩阵最大特征向量u的球坐标直接映射为旋转角theta和相位补偿角phi。相位参数phi是复数 JADE 区别于实值 Jacobi 的关键没有它非对角元素的虚部永远无法被压掉。主流程把三段串起来Z, V whiten(X, tol_eig1e-6) C4 compute_cumulants(Z) Ms select_eigematrices(C4, numn) U joint_diag(Ms, tol1e-8, max_sweeps50) W U.conj().T V Y W XW的作用是直接从观测 X 估计源信号所以最后一步用W X。输出的 Y 每一行对应一个源但行序和幅度、相位都是不确定的这是盲源分离的固有歧义后续要靠导频或调制特性来对齐。4. 复数 JADE 的参数设置与高频排错4.1 三个必调参数特征矩阵数、收敛阈值、白化截断JADE 用起来顺手是因为它的参数极少但每个参数失效时的表现差异很大。先看一张参数表再逐个展开。参数推荐范围失效时的表现调整方向特征矩阵数 numN 到 3N串音明显个别源混叠从 N 开始扫观察对角化目标函数对角化阈值 tol1e-8 到 1e-6迭代次数打满仍不退出放宽到 1e-6或检查白化白化截断 tol_eig1e-6 × 最大特征值源数变少或噪声被放大画特征值谱看肘部位置快拍数 T大于 50N通信建议大于 200N两块数据估计的 W 不一致增大 T或减少源数特征矩阵数是影响计算量最直接的开关。取 num N 时速度最快但遇到累积量估计噪声偏大的场景分离质量会明显下降取 2N 到 3N 相当于给联合对角化提供更多约束。判断方法是看每次 Jacobi 扫掠的change是否收敛到接近零不收敛时先加特征矩阵数不要急着放宽阈值。白化截断 tol_eig 的坑最隐蔽。复数混合矩阵的协方差特征值如果出现明显的台阶说明有效源数小于观测通道数此时保留所有特征值会让低能量噪声维度参与累积量计算四阶张量的信噪比被拉低。正确的做法是先打印evals / evals.max()看谱选出肘部位置再定阈值。4.2 用分块重采样检查分离矩阵稳定性调完参数后判断 JADE 是否可信的最好方法不是看单次分离波形而是用两段互不相交的数据分别估计分离矩阵比较两者的一致性。这是复数场景下最值得做的一步因为相位歧义会导致单次结果的波形看起来正常但换了数据段就完全变样。def estimate_w(X, n_src, blockNone, tol_eig1e-6): if block is None: block X.shape[1] idx rng.choice(X.shape[1], sizeblock, replaceFalse) Z, V whiten(X[:, idx], tol_eigtol_eig) C4 compute_cumulants(Z) Ms select_eigematrices(C4, numn_src) U joint_diag(Ms, tol1e-8, max_sweeps50) return U.conj().T V W1 estimate_w(X, n, block4000) W2 estimate_w(X, n, block4000) G np.abs(W1 np.linalg.pinv(W2))理想情况下W1 pinv(W2)应该是一个排列矩阵乘以对角相位取模之后每一行只有一个接近 1 的元素。用两个指标评价row_hit (np.argmax(G, axis1) 互不重复)给出排列是否一致G.max(axis1).mean()给出能量集中度。如果第二项低于 0.7第一反应不是调阈值而是回看源数和样本量。这个分块重采样检查对参数调优和自动化回归都有价值建议写进验证脚本而不是只在调试时手动跑。4.3 高斯源、样本量与累积量塌方累积量塌方是复数 JADE 最常见的失败模式快拍数不足时四阶矩的估计方差远大于二阶统计量导致 C4 中出现伪结构联合对角化会把噪声当成信号去对齐。表现是分离矩阵的奇异值分布平坦输出里每个通道都带背景噪声。经验上 T 至少要大于 50N通信符号流建议 200N 以上。如果源里混入了高斯噪声源JADE 的性能会随高斯成分增多而下降。两个高斯源的线性混合仍然是高斯四阶累积量无法感知这个自由度所以这类场景需要换基于二阶或高阶混合统计的方法。还有一种常见误用是把 JADE 直接用在强相关的源上比如两路来自同一发射机不同时延的多径信号它们的四阶累积量结构性重合JADE 只能分离第一路其余会残留在同一输出通道。5. 用两组分离矩阵校验相位歧义与排列一致性最后一层实操是把第 4 章的稳定性检查压成一个可回归的指标顺便把盲分离的相位歧义转成可量化的验收项。定义一个恢复指数函数def recovery_index(W1, W2): G np.abs(W1 np.linalg.pinv(W2)) n G.shape[0] row_pos np.argmax(G, axis1) perm_ok len(set(row_pos.tolist())) n energy np.mean(np.max(G, axis1)) return perm_ok, energy ok, energy recovery_index(W1, W2) print(ok, round(energy, 4))perm_ok为 True 且energy接近 1 时说明两次估计在排列意义上完全一致JADE 留下的不确定性只剩对角复增益。这个复增益在复基带系统里通常不是问题因为接收机后续还有载波同步和信道均衡。如果要把分离结果直接用于解调可以拿已知导频做一次最小二乘相位校正pilot ref_source[0, :] # 已知导频序列 gain (pilot y_ref.conj()) / (y_ref y_ref.conj()) y_corrected y_ref * gain.conj()注意这里gain是一个复数标量同时补偿幅度和相位计算时用共轭相乘保证复数相位对齐。这个技巧不改变 JADE 本身但能把分离出的复数信号快速接回后续的信号处理链适合写进验收文档的最后一个步骤。实际工程里建议把recovery_index和分块重采样一起固化到自动化测试里源数变化、通道数变化、信噪比变化时只需调 4.1 的参数表不需要改算法主体。本文还有配套的精品资源点击获取