简介面向信号处理与机器学习研究者的复数域独立成分分析实现资源重点基于快速ICA算法处理包含幅度与相位信息的复数混合信号适用于通信、雷达、生物电信号等需要盲源分离的工程场景。资源包内共3个文件总大小仅5KB核心为MATLAB函数源码.m文件附带MATLAB自动保存文件.asv与说明文本.txt结构清晰便于快速导入MATLAB环境运行与二次开发。目前已有489人浏览学习对于需要掌握复数ICA原理或直接复用MATLAB代码的研究生、工程师具有较高参考价值。源码覆盖数据预处理、白化、复数高阶矩计算、非线性函数映射与解混矩阵迭代求解等关键步骤阅读代码与说明文本可理解复数FASTICA如何同时利用幅度和相位信息避开常见初始化与收敛问题并可在此基础上扩展到自有数据实验。1. 复数 fastICA 到底在解决什么问题处理雷达阵列信号、OFDM 通信信号或者脑磁图MEG数据时你拿到的观测值不是实数序列而是 I/Q 两路基带采样本质上是复数向量。很多做盲源分离的工程师第一反应是把实部虚部分开拼成一个两倍长度的实数矩阵再丢给 scikit-learn 的 FastICA。这个做法能跑但分离结果往往是错的——因为实部虚部共享同一个混合矩阵拆开处理等于人为破坏了这个约束。复数 fastICA 解决的正是这个问题直接在复数域做独立成分分析输入输出都是复数混合矩阵、白化矩阵、迭代更新全部在复数域内完成。它适合处理阵列信号处理、无线通信干扰消除、脑电/脑磁信号分析这类源头就是复数模型的场景也是本文要拆解清楚的核心内容。2. 从实数 fastICA 到复数域数学上改了什么2.1 为什么不能直接套用实数算法实数 fastICA 的核心假设是观测信号 x As其中 A 是混合矩阵s 是统计独立的源信号。算法通过最大化非高斯性来估计分离矩阵 W使得 y Wx 尽可能接近 s。这个过程依赖峭度kurtosis或负熵negentropy作为目标函数而对复数信号这些统计量的定义方式完全不同。一个最直接的障碍是解析性。实数域里梯度上升和定点迭代都建立在实值函数对实变量可导的基础上。但复数域中任何一个实值函数比如负熵对复变量都不是解析的——它不满足柯西-黎曼条件直接求导会得到矛盾结果。这意味着你不能简单地把w^T x换成w^H x就完事迭代公式需要基于 Wirtinger 微积分重新推导。另一个问题是统计量本身。复数的二阶统计量不是协方差矩阵E[xx^T]而是共轭协方差矩阵E[xx^H]。如果信号是循环对称的circularly symmetricE[xx^T] 0这时实部虚部的互信息天然为零拆开处理也许碰巧可行但实际通信信号和雷达回波往往不具备循环对称性实部虚部之间存在二阶相关性拆开就丢了这部分信息。2.2 复数域的代价函数构造复数 fastICA 的常见做法是保留负熵作为代价函数但对复数随机变量重新定义。记 y w^H x其中 w 是复数权向量x 是复数观测向量。目标函数取为J(w) E[ G(|w^H x|^2) ]这里 G 是一个光滑偶函数常见的选择有G1(u) sqrt(a u)和G2(u) log(a u)其中 a 是防止数值溢出的小常数一般取 0.1 左右。注意 G 的输入是|y|^2不是 y 本身——这样 G 就是实数到实数的函数可以绕开复解析性的坑同时保留了 y 的相位信息。对 J(w) 求 Wirtinger 梯度可以得到更新方向。具体来说定义 g(u) dG(u)/du那么梯度上升方向为w ← E[ x (w^H x)^* g(|w^H x|^2) ] - E[ g(|w^H x|^2) |w^H x|^2 g(|w^H x|^2) ] w这个公式看起来复杂但它的结构和实数 fastICA 的定点迭代是平行的第一项是加权相关矩阵乘 w第二项是缩放修正。在代码实现中我们只需要在每个批次上计算这两个期望值。2.3 白化从协方差改成共轭协方差无论实数还是复数ICA 的第一步通常是白化目的是把混合矩阵 A 转化为酉矩阵正交矩阵的复数版从而降低后续优化的自由度。实数情形做特征分解E[xx^T] U D U^T复数情形则分解E[xx^H] U D U^H然后取白化矩阵 V D^(-1/2) U^H。这里有一个关键区别白化后的信号 z Vx 满足E[zz^H] I但E[zz^T]不一定为零。如果源信号是循环对称的那么E[zz^T] 0自动成立分离矩阵可以限制为酉矩阵如果不是比如某些调制方式E[zz^T]会携带额外信息但这不影响 fastICA 的收敛性只影响分离矩阵的约束形式。另一个工程上的坑是特征分解的排序。np.linalg.eigh返回的特征值按升序排列如果直接取前 M 个会拿到噪声子空间。正确做法是取最大的 M 个特征值对应的特征向量或者根据特征值大小截断。好在白化过程中我们只关心满秩情况一般直接全部保留。2.4 实数与复数算法的迭代对比下面用一个并排对比来展示两者的区别。假设我们已经有了白化后的数据 z实数 fastICA 对每个独立分量做w w / np.linalg.norm(w) # 归一化 w_new np.mean(z * g(w z), axis1) - np.mean(g_deriv(w z)) * w而复数版本在代码上只多了共轭和模长运算w w / np.sqrt(np.sum(np.abs(w)**2)) y np.conj(w) z # 复数域投影 w_new np.mean(z * np.conj(y) * g(np.abs(y)**2), axis1) - np.mean(g_deriv(np.abs(y)**2) np.abs(y)**2 * g2_deriv(np.abs(y)**2)) * w两者描述的是同一个算法只是统计矩的计算从w^T x换成了w^H x并把G(y^2)替换为G(|y|^2)。理解了这一层后面自己写实现或者读懂别人代码都不费劲。3. 用 Python 实现复数 fastICA最小可运行代码3.1 设计思路与算法流程为便于复现这里直接采用梯度上升加 Gram-Schmidt 正交化的实现路线。整体流程如下生成或读取复数观测矩阵 X形状为 (M, N)M 是观测通道数N 是采样点数。对 X 去均值X X - mean(X, axis1, keepdimsTrue)。用共轭协方差矩阵做白化z V X。对每个独立分量随机初始化单位复数权向量 w。迭代更新 w每次更新后做正交化消除已提取分量的影响。收敛后得到分离信号 Y W X。为什么不直接调用sklearn.decomposition.FastICA因为它的实现不支复数据传入complex128数组会直接报类型错误。而tensorflow或pytorch的 ICA 实现也很少见自己手写一个干净版本大约 60 行可控性更强。3.2 生成复数混合信号先用 QPSK 信号和调频信号构造一组已知源这样后面可以定量验证分离效果。QPSK 是典型的非循环对称复数信号用它做测试能同时验证算法对相位信息的利用。import numpy as np np.random.seed(42) n_samples 5000 n_sources 3 # 源1: QPSK 调制信号 data np.random.randint(0, 4, n_samples) s1 (2 * (data // 2) - 1) 1j * (2 * (data % 2) - 1) # 源2: 复数调频信号 t np.linspace(0, 10, n_samples) s2 np.exp(1j * (2 * np.pi * 0.5 * t 0.2 * np.sin(2 * np.pi * 3 * t))) # 源3: 复数噪声 正弦 s3 0.5 * np.exp(1j * 2 * np.pi * 0.1 * t) 0.3 * (np.random.randn(n_samples) 1j * np.random.randn(n_samples)) S np.vstack([s1, s2, s3]) # 随机复数混合矩阵 A np.random.randn(n_sources, n_sources) 1j * np.random.randn(n_sources, n_sources) X A SQPSK 源信号的实部和虚部都是 ±1峭度高是 fastICA 最喜欢的源类型。调频信号的包络恒定相位随时间变化测的是算法对相位结构的敏感度。第三个源加了高斯噪声模拟低信噪比场景。3.3 白化实现白化用特征分解注意一定要用eigh而不是eig因为共轭协方差矩阵是 Hermitian 矩阵eigh更稳定且更高效。def whiten(X): 复数白化: 使输出满足 E[zz^H] I n, m X.shape cov (X X.conj().T) / m D, U np.linalg.eigh(cov) # eigh 返回升序特征值取最大的 n 个注意特征向量也随之重排 idx np.argsort(D)[::-1] D, U D[idx], U[:, idx] V np.diag(1.0 / np.sqrt(D)) U.conj().T return V X, Veigh返回的特征向量每一列都是单位正交的满足 Hermitian 矩阵的要求。对特征值开方取倒数时加一个小的eps1e-12防止除零这在通道数多于源数时会遇到。白化后数据维度不变但各通道功率归一。3.4 核心迭代单分量提取非线性函数选用G(u) sqrt(0.1 u)它的一阶和二阶导数都是有理函数计算开销小且在 u 接近 0 时不会像log那样产生大梯度。def g(u): return 1.0 / (2 * np.sqrt(0.1 u)) def g_deriv(u): return -1.0 / (4 * np.power(0.1 u, 1.5)) def g2_deriv(u): return 3.0 / (8 * np.power(0.1 u, 2.5)) def complex_fastica_one_unit(z, w_init, tol1e-6, max_iter300): 提取一个复数独立分量 w w_init / np.sqrt(np.sum(np.abs(w_init)**2)) n z.shape[1] for i in range(max_iter): y w.conj() z # (n_sources,) 复数投影 u np.abs(y)**2 # 计算梯度上升方向 term1 (z * y.conj() * g(u)).mean(axis1) term2 (g_deriv(u) u * g2_deriv(u)).mean() w_new term1 - term2 * w # 归一化 w_new w_new / np.sqrt(np.sum(np.abs(w_new)**2)) # 收敛判断: 权向量变化小于阈值 if np.abs(np.abs(w_new.conj() w) - 1) tol: return w_new w w_new return w收敛条件用的是前后两次权向量的内积模长接近 1这等价于方向不再变化。注意np.abs(w_new.conj() w)可能大于 1所以不能直接比较差值要用绝对值减 1 的绝对值。3.5 多分量提取与正交化提取多个分量时每提取一个都要把所有已提取的权向量投影出去。标准的 Gram-Schmidt 过程在复数域的写法是def complex_fastica(X, n_componentsNone, tol1e-6, max_iter300): X X - X.mean(axis1, keepdimsTrue) Z, V whiten(X) if n_components is None: n_components Z.shape[0] W np.zeros((n_components, Z.shape[0]), dtypecomplex) for k in range(n_components): # 随机初始化单位复数向量 w np.random.randn(Z.shape[0]) 1j * np.random.randn(Z.shape[0]) w w / np.sqrt(np.sum(np.abs(w)**2)) for i in range(max_iter): y w.conj() Z u np.abs(y)**2 term1 (Z * y.conj() * g(u)).mean(axis1) term2 (g_deriv(u) u * g2_deriv(u)).mean() w_new term1 - term2 * w # 正交化: 减去与已提取分量的重叠部分 for j in range(k): w_new w_new - (W[j].conj() w_new) * W[j] w_new w_new / np.sqrt(np.sum(np.abs(w_new)**2)) if np.abs(np.abs(w_new.conj() w) - 1) tol: w w_new break w w_new W[k] w S_est W X return W, S_est, V正交化必须在归一化之前做否则投影的幅度不准确。如果先归一化再减投影会导致后续分量幅度偏小多次迭代后累积误差会很明显。另一个细节是每次迭代都在投影而不是只在初始化时投影——因为 w 在更新过程中会重新引入与已提取分量的相关性。3.6 运行与结果解释W, S_est, V complex_fastica(X, n_components3) # 输出分离信号的幅度范围检查是否合理 print(np.abs(S_est).mean(axis1))提取出的信号S_est与原始源S的顺序和尺度都不一样这是 ICA 固有的排列不确定性和尺度不确定性导致的。需要先做归一化再做相关性分析才能确认分离效果。下一章我们会做这个定量验证。4. 参数选型与工程中的坑4.1 非线性函数 G 的选择不同 G 函数对应不同的统计假设。本文实现的G(u) sqrt(0.1 u)适合超高斯源峭度大于零也就是大部分通信调制信号。G(u) log(0.1 u)对亚高斯源更鲁棒但梯度在 u 接近 0 时会变得很大收敛慢且容易震荡。G(u)g(u)适用源类型收敛速度稳定性sqrt(au)1/(2*sqrt(au))超高斯QPSK/16QAM快好log(au)1/(au)亚高斯/混合慢一般u^22u超高斯最快差对离群值敏感实际工程中如果不知道源的类型可以先用np.mean(np.abs(S)**4) - 2估计峭度峭度为正选 sqrt峭度为负选 log。峭度计算本身就基于复数信号的模长不需要拆实部虚部。4.2 收敛阈值和学习率本文的实现用的是定点迭代思想没有显式学习率所以不存在调学习率的烦恼。如果改成梯度上升的变体学习率建议从 0.1 开始每轮检查目标函数是否上升不升就减半。收敛阈值tol取 1e-6 已经足够严格取到 1e-8 不会有实质提升反而增加迭代次数。关于max_iter300 轮对大多数信号足够了。如果超过 300 轮仍未收敛优先检查白化是否正确np.max(np.abs(np.cov(Z) - np.eye(3)))应该在 1e-10 量级。白化不彻底的最常见原因是数据里有 NaN 或 Inf这会导致协方差矩阵奇异。4.3 复数数据预处理中心化和尺度去均值这一步在复数域同样必要但要注意减去的是复数均值不是分别减实部均值、虚部均值。两者在数学上等价但用X - X.mean(axis1, keepdimsTrue)一次性做更不容易出错。尺度方面fastICA 对幅度的绝对大小不敏感因为白化已经做了归一化但对不同通道之间的相对幅度敏感所以混合前各通道的增益差异会被白化自动校正。4.4 和 sklearn FastICA 拆实部虚部的区别常见的错误做法是X_real np.concatenate([X.real, X.imag], axis0) ica FastICA(n_components3) S_est ica.fit_transform(X_real.T)这样处理的问题有三个。第一分离矩阵的维度翻倍源数不变解得的是一个 6×3 的实数矩阵无法还原成复数域的 3×3 复矩阵。第二实部虚部的混合不是独立的拆开后等于无视了 I/Q 两路之间的耦合关系。第三即使强行把正负频率分量当作独立源分离结果的物理意义也不对应原始源信号。如果一定只能用实数工具正确做法是构建统计量完全的 2×2 增广矩阵X_aug np.vstack([X.real, X.imag]) # 对 X_aug 做实数 ICA再把结果的前半和后半合成复数但这样得到的源并不等于复数 fastICA 的源只在循环对称信号的特殊情况下才近似成立。与其绕路不如直接用复数实现。4.5 源数未知时的处理实际信号源数往往未知而 fastICA 要求输出维度不超过输入维度。通常做法是先对白化后的协方差矩阵做特征值分解看特征值之间存在明显落差的位置。特征值下降曲线从陡峭变平缓的拐点就是源数的估计值。在复数域同样适用D, _ np.linalg.eigh((Z Z.conj().T) / Z.shape[1]) D D[::-1] # 降序 # 找相邻特征值比值最大处 ratio D[:-1] / D[1:] n_est np.argmax(ratio) 1这个启发式方法在信噪比高于 10 dB 时基本可靠。信噪比更低时建议用更稳健的 MDL最小描述长度准则或随机矩阵理论的门限公式不过这些方法对复数协方差矩阵需要重新推导。5. 用归一化相关系数定量验证分离效果盲源分离的验证不能靠“看波形像不像”因为排列和尺度是任意的。学术界通用的做法是计算分离信号与源信号的相关系数矩阵再检查每行每列是否只有一个元素接近 1。受排列不确定性影响这个矩阵是置换矩阵在相位模糊下的推广归一化相关系数定义如下def complex_corr(s, y): 复数信号归一化相关系数, 消除尺度和相位影响 s s - np.mean(s) y y - np.mean(y) num np.abs(np.sum(s * np.conj(y))) den np.sqrt(np.sum(np.abs(s)**2) * np.sum(np.abs(y)**2)) return num / den corr_matrix np.zeros((n_sources, n_sources)) for i in range(n_sources): for j in range(n_sources): corr_matrix[i, j] complex_corr(S[i], S_est[j]) print(np.round(corr_matrix, 3))取模的目的是消除 ICA 固有的相位旋转——复数混合矩阵 A 引入的相位偏移会同时作用于源信号分离后每个源可能整体旋转一个任意角度这不算错误。相关系数的幅值对相位不敏感所以能正确反映匹配程度。理想情况下这个矩阵应该是置换矩阵比如[[0.98, 0.02, 0.01], [0.01, 0.03, 0.97], [0.02, 0.99, 0.02]]如果出现了某一行有两个元素都是 0.7 左右的数值说明两个源没有完全分离可能原因是非线性函数选错、迭代未收敛或者源之间存在微弱的相关性。进一步可以计算 Amari 误差指标这是 ICA 文献中最常用的定量评价标准。定义全局矩阵 C W A这里 W 包含了白化矩阵和分离矩阵置换矩阵 P 满足每行每列只有一个非零元素。Amari 误差为def amari_error(W_combined, A): W_combined 为综合分离矩阵, 需要乘以白化矩阵 V 一起计算 C W_combined A C np.abs(C) ** 2 # 每行最大值除以行和, 再转置求每列 row_sum C.sum(axis1, keepdimsTrue) col_sum C.sum(axis0, keepdimsTrue) row_term np.sum(C / row_sum, axis1) col_term np.sum(C / col_sum, axis0) return (np.sum(row_term / np.max(row_term)) np.sum(col_term / np.max(col_term))) / (2 * C.shape[0]) - 1注意前面的complex_fastica函数返回的 W 是作用在白化后数据 Z 上的因此综合分离矩阵应该是 W V。把 A、V、W 组合后算出 Amari 误差小于 0.1 说明分离效果良好0.01 量级就是非常理想了。多说一句验证时的取值范围随机复混合矩阵 A 的条件数最好不要太大条件数超过 100 时分离问题本身就病态算法再完善也救不回来。可以在生成 A 后检查np.linalg.cond(A)超过 100 就重新生成。这样你的验证脚本才是在测试算法而不是在测试 A 的随机性。本文还有配套的精品资源点击获取 复数FastICA：原理、Python实现与工程验证
简介面向信号处理与机器学习研究者的复数域独立成分分析实现资源重点基于快速ICA算法处理包含幅度与相位信息的复数混合信号适用于通信、雷达、生物电信号等需要盲源分离的工程场景。资源包内共3个文件总大小仅5KB核心为MATLAB函数源码.m文件附带MATLAB自动保存文件.asv与说明文本.txt结构清晰便于快速导入MATLAB环境运行与二次开发。目前已有489人浏览学习对于需要掌握复数ICA原理或直接复用MATLAB代码的研究生、工程师具有较高参考价值。源码覆盖数据预处理、白化、复数高阶矩计算、非线性函数映射与解混矩阵迭代求解等关键步骤阅读代码与说明文本可理解复数FASTICA如何同时利用幅度和相位信息避开常见初始化与收敛问题并可在此基础上扩展到自有数据实验。1. 复数 fastICA 到底在解决什么问题处理雷达阵列信号、OFDM 通信信号或者脑磁图MEG数据时你拿到的观测值不是实数序列而是 I/Q 两路基带采样本质上是复数向量。很多做盲源分离的工程师第一反应是把实部虚部分开拼成一个两倍长度的实数矩阵再丢给 scikit-learn 的 FastICA。这个做法能跑但分离结果往往是错的——因为实部虚部共享同一个混合矩阵拆开处理等于人为破坏了这个约束。复数 fastICA 解决的正是这个问题直接在复数域做独立成分分析输入输出都是复数混合矩阵、白化矩阵、迭代更新全部在复数域内完成。它适合处理阵列信号处理、无线通信干扰消除、脑电/脑磁信号分析这类源头就是复数模型的场景也是本文要拆解清楚的核心内容。2. 从实数 fastICA 到复数域数学上改了什么2.1 为什么不能直接套用实数算法实数 fastICA 的核心假设是观测信号 x As其中 A 是混合矩阵s 是统计独立的源信号。算法通过最大化非高斯性来估计分离矩阵 W使得 y Wx 尽可能接近 s。这个过程依赖峭度kurtosis或负熵negentropy作为目标函数而对复数信号这些统计量的定义方式完全不同。一个最直接的障碍是解析性。实数域里梯度上升和定点迭代都建立在实值函数对实变量可导的基础上。但复数域中任何一个实值函数比如负熵对复变量都不是解析的——它不满足柯西-黎曼条件直接求导会得到矛盾结果。这意味着你不能简单地把w^T x换成w^H x就完事迭代公式需要基于 Wirtinger 微积分重新推导。另一个问题是统计量本身。复数的二阶统计量不是协方差矩阵E[xx^T]而是共轭协方差矩阵E[xx^H]。如果信号是循环对称的circularly symmetricE[xx^T] 0这时实部虚部的互信息天然为零拆开处理也许碰巧可行但实际通信信号和雷达回波往往不具备循环对称性实部虚部之间存在二阶相关性拆开就丢了这部分信息。2.2 复数域的代价函数构造复数 fastICA 的常见做法是保留负熵作为代价函数但对复数随机变量重新定义。记 y w^H x其中 w 是复数权向量x 是复数观测向量。目标函数取为J(w) E[ G(|w^H x|^2) ]这里 G 是一个光滑偶函数常见的选择有G1(u) sqrt(a u)和G2(u) log(a u)其中 a 是防止数值溢出的小常数一般取 0.1 左右。注意 G 的输入是|y|^2不是 y 本身——这样 G 就是实数到实数的函数可以绕开复解析性的坑同时保留了 y 的相位信息。对 J(w) 求 Wirtinger 梯度可以得到更新方向。具体来说定义 g(u) dG(u)/du那么梯度上升方向为w ← E[ x (w^H x)^* g(|w^H x|^2) ] - E[ g(|w^H x|^2) |w^H x|^2 g(|w^H x|^2) ] w这个公式看起来复杂但它的结构和实数 fastICA 的定点迭代是平行的第一项是加权相关矩阵乘 w第二项是缩放修正。在代码实现中我们只需要在每个批次上计算这两个期望值。2.3 白化从协方差改成共轭协方差无论实数还是复数ICA 的第一步通常是白化目的是把混合矩阵 A 转化为酉矩阵正交矩阵的复数版从而降低后续优化的自由度。实数情形做特征分解E[xx^T] U D U^T复数情形则分解E[xx^H] U D U^H然后取白化矩阵 V D^(-1/2) U^H。这里有一个关键区别白化后的信号 z Vx 满足E[zz^H] I但E[zz^T]不一定为零。如果源信号是循环对称的那么E[zz^T] 0自动成立分离矩阵可以限制为酉矩阵如果不是比如某些调制方式E[zz^T]会携带额外信息但这不影响 fastICA 的收敛性只影响分离矩阵的约束形式。另一个工程上的坑是特征分解的排序。np.linalg.eigh返回的特征值按升序排列如果直接取前 M 个会拿到噪声子空间。正确做法是取最大的 M 个特征值对应的特征向量或者根据特征值大小截断。好在白化过程中我们只关心满秩情况一般直接全部保留。2.4 实数与复数算法的迭代对比下面用一个并排对比来展示两者的区别。假设我们已经有了白化后的数据 z实数 fastICA 对每个独立分量做w w / np.linalg.norm(w) # 归一化 w_new np.mean(z * g(w z), axis1) - np.mean(g_deriv(w z)) * w而复数版本在代码上只多了共轭和模长运算w w / np.sqrt(np.sum(np.abs(w)**2)) y np.conj(w) z # 复数域投影 w_new np.mean(z * np.conj(y) * g(np.abs(y)**2), axis1) - np.mean(g_deriv(np.abs(y)**2) np.abs(y)**2 * g2_deriv(np.abs(y)**2)) * w两者描述的是同一个算法只是统计矩的计算从w^T x换成了w^H x并把G(y^2)替换为G(|y|^2)。理解了这一层后面自己写实现或者读懂别人代码都不费劲。3. 用 Python 实现复数 fastICA最小可运行代码3.1 设计思路与算法流程为便于复现这里直接采用梯度上升加 Gram-Schmidt 正交化的实现路线。整体流程如下生成或读取复数观测矩阵 X形状为 (M, N)M 是观测通道数N 是采样点数。对 X 去均值X X - mean(X, axis1, keepdimsTrue)。用共轭协方差矩阵做白化z V X。对每个独立分量随机初始化单位复数权向量 w。迭代更新 w每次更新后做正交化消除已提取分量的影响。收敛后得到分离信号 Y W X。为什么不直接调用sklearn.decomposition.FastICA因为它的实现不支复数据传入complex128数组会直接报类型错误。而tensorflow或pytorch的 ICA 实现也很少见自己手写一个干净版本大约 60 行可控性更强。3.2 生成复数混合信号先用 QPSK 信号和调频信号构造一组已知源这样后面可以定量验证分离效果。QPSK 是典型的非循环对称复数信号用它做测试能同时验证算法对相位信息的利用。import numpy as np np.random.seed(42) n_samples 5000 n_sources 3 # 源1: QPSK 调制信号 data np.random.randint(0, 4, n_samples) s1 (2 * (data // 2) - 1) 1j * (2 * (data % 2) - 1) # 源2: 复数调频信号 t np.linspace(0, 10, n_samples) s2 np.exp(1j * (2 * np.pi * 0.5 * t 0.2 * np.sin(2 * np.pi * 3 * t))) # 源3: 复数噪声 正弦 s3 0.5 * np.exp(1j * 2 * np.pi * 0.1 * t) 0.3 * (np.random.randn(n_samples) 1j * np.random.randn(n_samples)) S np.vstack([s1, s2, s3]) # 随机复数混合矩阵 A np.random.randn(n_sources, n_sources) 1j * np.random.randn(n_sources, n_sources) X A SQPSK 源信号的实部和虚部都是 ±1峭度高是 fastICA 最喜欢的源类型。调频信号的包络恒定相位随时间变化测的是算法对相位结构的敏感度。第三个源加了高斯噪声模拟低信噪比场景。3.3 白化实现白化用特征分解注意一定要用eigh而不是eig因为共轭协方差矩阵是 Hermitian 矩阵eigh更稳定且更高效。def whiten(X): 复数白化: 使输出满足 E[zz^H] I n, m X.shape cov (X X.conj().T) / m D, U np.linalg.eigh(cov) # eigh 返回升序特征值取最大的 n 个注意特征向量也随之重排 idx np.argsort(D)[::-1] D, U D[idx], U[:, idx] V np.diag(1.0 / np.sqrt(D)) U.conj().T return V X, Veigh返回的特征向量每一列都是单位正交的满足 Hermitian 矩阵的要求。对特征值开方取倒数时加一个小的eps1e-12防止除零这在通道数多于源数时会遇到。白化后数据维度不变但各通道功率归一。3.4 核心迭代单分量提取非线性函数选用G(u) sqrt(0.1 u)它的一阶和二阶导数都是有理函数计算开销小且在 u 接近 0 时不会像log那样产生大梯度。def g(u): return 1.0 / (2 * np.sqrt(0.1 u)) def g_deriv(u): return -1.0 / (4 * np.power(0.1 u, 1.5)) def g2_deriv(u): return 3.0 / (8 * np.power(0.1 u, 2.5)) def complex_fastica_one_unit(z, w_init, tol1e-6, max_iter300): 提取一个复数独立分量 w w_init / np.sqrt(np.sum(np.abs(w_init)**2)) n z.shape[1] for i in range(max_iter): y w.conj() z # (n_sources,) 复数投影 u np.abs(y)**2 # 计算梯度上升方向 term1 (z * y.conj() * g(u)).mean(axis1) term2 (g_deriv(u) u * g2_deriv(u)).mean() w_new term1 - term2 * w # 归一化 w_new w_new / np.sqrt(np.sum(np.abs(w_new)**2)) # 收敛判断: 权向量变化小于阈值 if np.abs(np.abs(w_new.conj() w) - 1) tol: return w_new w w_new return w收敛条件用的是前后两次权向量的内积模长接近 1这等价于方向不再变化。注意np.abs(w_new.conj() w)可能大于 1所以不能直接比较差值要用绝对值减 1 的绝对值。3.5 多分量提取与正交化提取多个分量时每提取一个都要把所有已提取的权向量投影出去。标准的 Gram-Schmidt 过程在复数域的写法是def complex_fastica(X, n_componentsNone, tol1e-6, max_iter300): X X - X.mean(axis1, keepdimsTrue) Z, V whiten(X) if n_components is None: n_components Z.shape[0] W np.zeros((n_components, Z.shape[0]), dtypecomplex) for k in range(n_components): # 随机初始化单位复数向量 w np.random.randn(Z.shape[0]) 1j * np.random.randn(Z.shape[0]) w w / np.sqrt(np.sum(np.abs(w)**2)) for i in range(max_iter): y w.conj() Z u np.abs(y)**2 term1 (Z * y.conj() * g(u)).mean(axis1) term2 (g_deriv(u) u * g2_deriv(u)).mean() w_new term1 - term2 * w # 正交化: 减去与已提取分量的重叠部分 for j in range(k): w_new w_new - (W[j].conj() w_new) * W[j] w_new w_new / np.sqrt(np.sum(np.abs(w_new)**2)) if np.abs(np.abs(w_new.conj() w) - 1) tol: w w_new break w w_new W[k] w S_est W X return W, S_est, V正交化必须在归一化之前做否则投影的幅度不准确。如果先归一化再减投影会导致后续分量幅度偏小多次迭代后累积误差会很明显。另一个细节是每次迭代都在投影而不是只在初始化时投影——因为 w 在更新过程中会重新引入与已提取分量的相关性。3.6 运行与结果解释W, S_est, V complex_fastica(X, n_components3) # 输出分离信号的幅度范围检查是否合理 print(np.abs(S_est).mean(axis1))提取出的信号S_est与原始源S的顺序和尺度都不一样这是 ICA 固有的排列不确定性和尺度不确定性导致的。需要先做归一化再做相关性分析才能确认分离效果。下一章我们会做这个定量验证。4. 参数选型与工程中的坑4.1 非线性函数 G 的选择不同 G 函数对应不同的统计假设。本文实现的G(u) sqrt(0.1 u)适合超高斯源峭度大于零也就是大部分通信调制信号。G(u) log(0.1 u)对亚高斯源更鲁棒但梯度在 u 接近 0 时会变得很大收敛慢且容易震荡。G(u)g(u)适用源类型收敛速度稳定性sqrt(au)1/(2*sqrt(au))超高斯QPSK/16QAM快好log(au)1/(au)亚高斯/混合慢一般u^22u超高斯最快差对离群值敏感实际工程中如果不知道源的类型可以先用np.mean(np.abs(S)**4) - 2估计峭度峭度为正选 sqrt峭度为负选 log。峭度计算本身就基于复数信号的模长不需要拆实部虚部。4.2 收敛阈值和学习率本文的实现用的是定点迭代思想没有显式学习率所以不存在调学习率的烦恼。如果改成梯度上升的变体学习率建议从 0.1 开始每轮检查目标函数是否上升不升就减半。收敛阈值tol取 1e-6 已经足够严格取到 1e-8 不会有实质提升反而增加迭代次数。关于max_iter300 轮对大多数信号足够了。如果超过 300 轮仍未收敛优先检查白化是否正确np.max(np.abs(np.cov(Z) - np.eye(3)))应该在 1e-10 量级。白化不彻底的最常见原因是数据里有 NaN 或 Inf这会导致协方差矩阵奇异。4.3 复数数据预处理中心化和尺度去均值这一步在复数域同样必要但要注意减去的是复数均值不是分别减实部均值、虚部均值。两者在数学上等价但用X - X.mean(axis1, keepdimsTrue)一次性做更不容易出错。尺度方面fastICA 对幅度的绝对大小不敏感因为白化已经做了归一化但对不同通道之间的相对幅度敏感所以混合前各通道的增益差异会被白化自动校正。4.4 和 sklearn FastICA 拆实部虚部的区别常见的错误做法是X_real np.concatenate([X.real, X.imag], axis0) ica FastICA(n_components3) S_est ica.fit_transform(X_real.T)这样处理的问题有三个。第一分离矩阵的维度翻倍源数不变解得的是一个 6×3 的实数矩阵无法还原成复数域的 3×3 复矩阵。第二实部虚部的混合不是独立的拆开后等于无视了 I/Q 两路之间的耦合关系。第三即使强行把正负频率分量当作独立源分离结果的物理意义也不对应原始源信号。如果一定只能用实数工具正确做法是构建统计量完全的 2×2 增广矩阵X_aug np.vstack([X.real, X.imag]) # 对 X_aug 做实数 ICA再把结果的前半和后半合成复数但这样得到的源并不等于复数 fastICA 的源只在循环对称信号的特殊情况下才近似成立。与其绕路不如直接用复数实现。4.5 源数未知时的处理实际信号源数往往未知而 fastICA 要求输出维度不超过输入维度。通常做法是先对白化后的协方差矩阵做特征值分解看特征值之间存在明显落差的位置。特征值下降曲线从陡峭变平缓的拐点就是源数的估计值。在复数域同样适用D, _ np.linalg.eigh((Z Z.conj().T) / Z.shape[1]) D D[::-1] # 降序 # 找相邻特征值比值最大处 ratio D[:-1] / D[1:] n_est np.argmax(ratio) 1这个启发式方法在信噪比高于 10 dB 时基本可靠。信噪比更低时建议用更稳健的 MDL最小描述长度准则或随机矩阵理论的门限公式不过这些方法对复数协方差矩阵需要重新推导。5. 用归一化相关系数定量验证分离效果盲源分离的验证不能靠“看波形像不像”因为排列和尺度是任意的。学术界通用的做法是计算分离信号与源信号的相关系数矩阵再检查每行每列是否只有一个元素接近 1。受排列不确定性影响这个矩阵是置换矩阵在相位模糊下的推广归一化相关系数定义如下def complex_corr(s, y): 复数信号归一化相关系数, 消除尺度和相位影响 s s - np.mean(s) y y - np.mean(y) num np.abs(np.sum(s * np.conj(y))) den np.sqrt(np.sum(np.abs(s)**2) * np.sum(np.abs(y)**2)) return num / den corr_matrix np.zeros((n_sources, n_sources)) for i in range(n_sources): for j in range(n_sources): corr_matrix[i, j] complex_corr(S[i], S_est[j]) print(np.round(corr_matrix, 3))取模的目的是消除 ICA 固有的相位旋转——复数混合矩阵 A 引入的相位偏移会同时作用于源信号分离后每个源可能整体旋转一个任意角度这不算错误。相关系数的幅值对相位不敏感所以能正确反映匹配程度。理想情况下这个矩阵应该是置换矩阵比如[[0.98, 0.02, 0.01], [0.01, 0.03, 0.97], [0.02, 0.99, 0.02]]如果出现了某一行有两个元素都是 0.7 左右的数值说明两个源没有完全分离可能原因是非线性函数选错、迭代未收敛或者源之间存在微弱的相关性。进一步可以计算 Amari 误差指标这是 ICA 文献中最常用的定量评价标准。定义全局矩阵 C W A这里 W 包含了白化矩阵和分离矩阵置换矩阵 P 满足每行每列只有一个非零元素。Amari 误差为def amari_error(W_combined, A): W_combined 为综合分离矩阵, 需要乘以白化矩阵 V 一起计算 C W_combined A C np.abs(C) ** 2 # 每行最大值除以行和, 再转置求每列 row_sum C.sum(axis1, keepdimsTrue) col_sum C.sum(axis0, keepdimsTrue) row_term np.sum(C / row_sum, axis1) col_term np.sum(C / col_sum, axis0) return (np.sum(row_term / np.max(row_term)) np.sum(col_term / np.max(col_term))) / (2 * C.shape[0]) - 1注意前面的complex_fastica函数返回的 W 是作用在白化后数据 Z 上的因此综合分离矩阵应该是 W V。把 A、V、W 组合后算出 Amari 误差小于 0.1 说明分离效果良好0.01 量级就是非常理想了。多说一句验证时的取值范围随机复混合矩阵 A 的条件数最好不要太大条件数超过 100 时分离问题本身就病态算法再完善也救不回来。可以在生成 A 后检查np.linalg.cond(A)超过 100 就重新生成。这样你的验证脚本才是在测试算法而不是在测试 A 的随机性。本文还有配套的精品资源点击获取