MATLAB JADE盲源分离：四阶累积量与联合对角化拆解混合信号

简介：面向信号处理、语音分析与算法复现的学习者，这份资源围绕JADE盲源分离展开，适合已有一定MATLAB基础、希望理解盲信号分离原理并动手实现的研究生与工程技术人员。内容从盲源分离的基本问题切入：在源信号与混合方式均未知的条件下，利用观测信号的统计独立性恢复原始信号，并重点说明JADE依赖高阶累积量（常用四阶，也有三阶方案）与源信号非高斯性的思路。资源还提示了方法边界：除统计独立外，源信号中最多只能有一个高斯信号，算法主要适用于统计独立的非高斯信号，对非白、非平稳特性未作专门考虑。压缩包约426KB，文件总数与类型明细暂未提供，从标题看应包含原理说明及MATLAB程序，便于对照公式理解累积量计算、分离矩阵求解与结果验证。已有953人学习，适合作为盲源分离入门与课程实验的参考材料。

1. 三路信号混在一起时，这套 JADE 的 matlab 算法能拆到什么程度

实验室里那台三通道采集设备，采回来的信号永远是混在一起的：三路传感器同时录，同一个冲击事件在三路里都留下痕迹，谁的贡献占多少说不清。想拆开，深度学习的路子要先有标注、要调 GPU、要反复试网络结构；而一套完整的 matlab 算法——JADE（Joint Approximate Diagonalization of Eigenmatrices，特征矩阵联合近似对角化）——核心代码不到两百行，跑一次几秒钟，就能把线性瞬时混合的多路信号按统计独立性拆回源信号。它押注的是四阶累积量：高斯噪声的四阶累积量为零，所以现场噪声大一点，它反而不慌。适合做课程设计、做传感器阵列预处理、想给 ICA 找一个确定性基准的人。下面按能跑通的顺序拆：原理、代码结构、参数、踩坑。

2. 四阶累积量凭什么能定位那个旋转：从白化到联合对角化

2.1 线性瞬时混合与两个模糊：为什么二阶统计量不够

先把模型钉死，后面的代码才有依据。JADE 面对的是线性瞬时混合：

x(t) = A s(t) + n(t)

x 是 n 维观测，s 是 m 维源（这里假设 m = n，源数靠主成分谱自己判断），A 是 n×m 的未知混合矩阵，n(t) 是观测噪声。要求只有三条：源之间统计独立、源里最多一个高斯分布、A 列满秩。这三条不满足，后面无论怎么调参数都是白费。

只用二阶统计量能做到什么程度？对 x 求协方差 R = E[xx^T] = A E[ss^T] A^T，特征分解后能得到一个白化矩阵 W，使 z = Wx 满足 E[zz^T] = I。但白化只消除了二阶相关性，剩下的自由度是一个正交旋转：z = U s'，U 是任意正交矩阵，你怎么转，二阶统计量都看不出差别。这就是 ICA 里常说的「白化之后只剩一个旋转」，这个旋转必须靠更高阶的统计信息来定，JADE 用的是四阶累积量。

另一个绕不开的点是两个模糊：分离出来的信号幅度不确定（源乘任意标量，混合矩阵除以该标量，观测不变），顺序也不确定（源的排列可以任意换）。所以判断分离好坏，不能拿分离信号和原始信号逐点相减看误差，必须用带排列和缩放不变性的指标，这一点在第 3 章的 PI 函数里会落实到代码。

2.2 白化：把混合矩阵压成一个正交旋转

白化的实现有讲究，不是一句zscore能替代的。标准做法是对协方差矩阵做特征分解，取特征值的负平方根做缩放。这里有两个工程细节：一是协方差必须强制对称，eig在病态矩阵上算出来的 R 可能有微小的非对称，导致特征向量出现复数；二是特征值要设下限，当某一路传感器几乎没信号时，最小特征值接近零，取负平方根会直接炸掉，分离结果全是 NaN。

function [Z, W] = whiten(X, r) % X: n x N 观测矩阵，每行一个通道 % r: 保留的主分量个数，缺省等于通道数 n [n, N] = size(X); if nargin < 2 || isempty(r), r = n; end Xc = X - mean(X, 2); % 逐通道去均值，累积量对均值敏感 R = (Xc * Xc') / N; % 协方差 [E, D] = eig((R + R') / 2); % 强制对称，避免复数特征向量 [d, idx] = sort(diag(D), 'descend'); E = E(:, idx); d = d(1:r); floor_val = max(d) * 1e-10; % 特征值下限，防止病态放大 d(d < floor_val) = floor_val; W = diag(d.^(-1/2)) * E'; % 白化矩阵 Z = W * Xc; % 白化后 E[zz^T] ≈ I

r这个参数决定保留几个主分量，也就间接决定了后面认为有几个源。常见错误是直接取r = n，明明只有两路有效信号，硬塞三路进去，多出来的那一维全是噪声，累积量矩阵束会被噪声主导，分离结果看起来「拆开了但每一路都很脏」。判断方法在 4.2 节给：看特征值谱有没有断崖。

白化还有一个容易忽略的前提——各通道的量纲要大致可比。三路传感器一个是毫伏级、一个是伏级，协方差矩阵的条件数会非常大，白化相当于把弱通道放大一千倍，同时把它的噪声也放大一千倍。常见做法是白化前对每个通道做一次标准差归一化，让协方差矩阵的对角线都在 1 附近，这一步在实测数据上经常是决定性的。

2.3 四阶累积量矩阵束：为什么是 Q(M) 而不是 n 的四次方张量

四阶累积量完整形式是一个 n×n×n×n 的张量，元素个数是 n 的四次方。三通道是 81 个，八通道就是 4096 个，直接对张量做联合对角化，计算量和存储都不划算。JADE 的做法是把张量沿着若干方向投影成一组 n×n 的矩阵，这叫矩阵束（matrix pencil）。

对白化后的数据 z，定义累积量矩阵：

Q(M) = E[ (z^T M z) · z z^T ] − 2M − tr(M)·I

M 是任意 n×n 的权重矩阵。这个式子里最关键的性质是：如果 z 服从高斯分布，Q(M) 恒等于零矩阵。推导不复杂，高斯分布的四阶累积量为零，而白化后的二阶项刚好被 −2M − tr(M)I 抵消掉。这意味着观测噪声只要是高斯的，无论多大方差，对 Q(M) 的贡献都是零，JADE 的抗噪能力就来自这里，而不是来自什么滤波预处理。

权重矩阵 M 取多少组？理论上取 n² 个线性无关的基就够了，实际实现中利用 Q(M) 的对称性，取 n(n+1)/2 个对称基即可覆盖。三通道是 6 个矩阵，八通道是 36 个，规模可控。

function Q = cum4mat(Z, M) % Z: n x N，已经白化且去均值 % M: n x n 权重矩阵 [n, N] = size(Z); MZ = M * Z; % n x N u = sum(Z .* MZ, 1); % 1 x N，即每个时刻的 z^T M z Q = (Z .* u) * Z' / N - 2*M - trace(M)*eye(n);

这段代码里Z .* u是隐式扩展，把第 j 个时刻的标量 u(j) 乘到 z(:,j) 上，再乘 Z' 求平均，得到的就是 E[(z^T M z) z z^T]。老版本 MATLAB 不支持隐式扩展，要写成bsxfun(@times, Z, u)，这个兼容性坑在第 5 章还会提一次。

2.4 联合对角化的旋转角闭式解

现在手上有一组矩阵 Q₁…Q_K，性质是：当分离矩阵 U 取对了，U^T Q_k U 对每个 k 都是对角阵。原因很直接：源之间独立，四阶互累积量为零，只有自累积量留在对角线上。于是问题变成找一个正交矩阵 U，让这组矩阵同时尽可能对角，目标函数写成非对角元素平方和：

J(U) = Σ_k off(U^T Q_k U)²

求解用 Jacobi 旋转：每次挑两个坐标轴 (p, q)，做一个 2×2 的 Givens 旋转，把所有矩阵的 (p,q) 和 (q,p) 元素往零推，逐对扫完一遍算一次迭代，反复扫到收敛。关键是这个旋转角有闭式解，不用做一维搜索。

记 u_k = Q_k(p,p) − Q_k(q,q)，v_k = Q_k(p,q) + Q_k(q,p)，累加 C = Σu_k²、D = Σu_k v_k、E = Σv_k²，则最优旋转角满足：

θ = 0.25 · atan2(2D, C − E)

这个式子来自「最大化旋转后对角元素平方和」的推导，和最小化非对角平方和是等价的。验证一下：单个矩阵 [[1, 0.5],[0.5, 1]]，u = 0、v = 1，C = 0、D = 0、E = 1，atan2(0, −1) = π，θ = π/4，正好把它转成对角阵 [[0.5, 0],[0, 1.5]]。旋转周期是 π/2，θ 落在 (−π/4, π/4] 恰好覆盖一个周期，不会漏解。

提示：这里的 θ 定义和教科书上某些 Givens 旋转差一个符号，取决于旋转矩阵写成 G^T A G 还是 G A G^T。抄公式的时候必须和代码里的更新方式对齐，否则迭代会发散而不是收敛。

3. 代码怎么组织：五个环节的接口与回代顺序

3.1 主脚本：先造一组自己知道答案的混合信号

验证任何盲分离程序，第一步都是造一个「知道真值」的数据集。不要上来就拿实测数据调，那样分离错了你根本不知道是算法问题还是数据问题。下面的主脚本生成三路特性不同的源：纯正弦（二阶统计量强、四阶弱）、方波（四阶很强）、调幅信号（时变包络，考验算法对非平稳性的鲁棒性），再用一个条件数适中的矩阵混合，最后叠一点高斯噪声。

rng(7); % 固定随机种子，保证结果可复现 fs = 500; T = 4; N = fs * T; t = (0:N-1) / fs; S_true = [ sin(2*pi*7*t); % 源1：纯正弦 sign(sin(2*pi*13*t)); % 源2：方波，四阶累积量强 (1 + 0.6*sin(2*pi*2.5*t)) .* sin(2*pi*41*t) ]; % 源3：调幅 A_true = [1.0 0.6 0.3; 0.4 1.0 0.5; 0.2 0.5 1.0]; % 对角线占优，条件数约 2，数值友好 X = A_true * S_true; X = X + 0.05 * randn(size(X)); % 高斯观测噪声，JADE 对它不敏感 [S_est, W, A_est] = jade(X, 3); % r=3，认为有三个源 PI = perf_index(W, A_true); fprintf('分离矩阵与真实混合矩阵的 PI = %.4f\n', PI);

这里的rng(7)不是可有可无的装饰。四阶累积量的估计方差比二阶大得多，换一个随机种子，PI 可能从 0.02 跳到 0.15，尤其是样本量在两三千点的时候。做对比实验必须固定种子，否则你调的「参数优化」很可能只是随机波动。

A_true我特意选成对角占优，条件数约 2。如果混合矩阵接近奇异，比如两列几乎平行，那本质上就是两路信号在观测空间里分不开，任何盲分离算法都会翻车，这属于数据问题而不是算法问题。

3.2 jade 主函数：白化、累积量、回代三步的接口约定

主函数只做三件事：白化、构造矩阵束、联合对角化，然后把分离矩阵回代到原始观测上。接口设计上有个细节值得说：返回分离矩阵 W，而不是只返回分离信号 S。因为做性能评估、做后续的通道匹配、做在线分离，都需要 W 这个线性算子；只给 S 的话，新来一批数据你还得重新跑一遍整个算法。

function [S, W, A_est] = jade(X, r) % X: n x N 观测 % r: 源数目（保留的主分量个数） [n, N] = size(X); if nargin < 2 || isempty(r), r = n; end Xc = X - mean(X, 2); % 去均值 [Z, Wz] = whiten(Xc, r); % 白化，Wz 是 r x n Qset = cum4_set(Z); % n x n x K 的矩阵束 [U, ~] = joint_diag(Qset, 1e-8, 200); % U 是 r x r 正交阵 W = U' * Wz; % 总分离矩阵 r x n S = W * Xc; % 分离信号 r x N A_est = pinv(W); % 估计的混合矩阵 n x r

矩阵束的生成函数把 n(n+1)/2 个基矩阵遍历一遍，注意非对角基要对称化，写成 E(p,q) = E(q,p) = 1，这样 Q(M) 才是实对称矩阵，后面的 Jacobi 旋转才能用实数的闭式解。如果只写 E(p,q) = 1，得到的 Q 是非对称的，联合对角化就得处理复数，代码复杂度翻倍不说，收敛也慢。

function Qset = cum4_set(Z) n = size(Z, 1); K = n * (n + 1) / 2; Qset = zeros(n, n, K); k = 0; for p = 1:n for q = p:n E = zeros(n); if p == q E(p,p) = 1; % 对角基 else E(p,q) = 1; E(q,p) = 1; % 对称非对角基 end k = k + 1; Qset(:,:,k) = cum4mat(Z, E); end end

接口一览：

本文还有配套的精品资源，点击获取