简介针对欠定盲源分离UBSS问题的一份MATLAB实现工具包面向信号处理、通信与机器学习方向的研究者及学生。当观测通道数少于源信号数时UBSS需要利用稀疏性、独立性等先验从混合信号中恢复独立源SCAN相关算法提供了可行的分离方案。压缩包共5个文件以4个.m脚本为主配合1个txt说明文档整体体积仅约5KB代码轻量、易于阅读与二次修改。.m文件包含SCAN演示脚本、DBSE算法函数以及经典的SOBI等分离方法readme.txt则对原理与调用方式作了简要说明。已有254人学习下载适合希望在欠定盲源分离方向快速上手或开展实验对比的研究者。通过学习这些代码可直观理解欠定情景下的稀疏分解、矩阵联合对角化等关键步骤并以此为基底进行算法改进与应用扩展。1. 欠定盲源分离的 SCAN 路线为什么两个传感器救不回三个源在振动监测或者多麦克风录音现场最常遇到的一个尴尬是想分离的源有 3 个手上却只有 2 路采集通道。直接跑 ICA 往往会得到一堆听不出内容的噪声因为传统盲源分离的前提是“观察数不少于源数”。对不少做现场调试的工程师来说scan 这个词首先让人想到 Modbus 扫描、设备地址排查这类活儿但在信号处理里SCAN 代表另一条路线Sparse Component Analysis 的工程化变体。它靠“稀疏性”把欠定问题重新变成可解问题。这篇按工程落地方式来拆解 SCAN 做欠定盲源分离的完整路径原理、最小可跑代码、参数坑和验证手段。适合正在做传感器阵列、语音分离、故障诊断的从业者。2. 欠定盲源分离为什么让 ICA 翻车稀疏性假设才是 SCAN 的立身之本盲源分离的标准模型写出来并不复杂x A s v。x 是 m 路观测s 是 n 个源A 是 m×n 的混合矩阵v 是噪声。传统 ICA 能工作的前提是 m ≥ n也就是传感器数量不少于源数量。可现实里往往是反过来的轴承故障、齿轮啮合、电机噪声同时存在麦克风或加速度计却只有两只。m n 的时候问题从“求一个可逆分离矩阵”变成了“解一个欠定方程组”数学上立刻变得麻烦起来。很多第一次接触欠定盲源分离的人会想既然 A 不知道那把 A 和 s 一起估计不就行了理论上可以但实践里几乎必翻车。欠定问题里 A 不是一个方阵A 的列数比行数多它能张成的空间维度最多只有 m。给定一次采样 x方程 As x 里未知数个数 n 比方程个数 m 多解空间是 n-m 维的。也就是说光是“满足观测方程”这个条件就能找出无穷多个 s。如果不加额外约束算法输出的结果里有很大一部分其实是随便猜的跟真实源没有任何关系。2.1 从 xAs 说起欠定问题到底难在哪欠定的“难”首先难在信息量不够。以 2 个传感器、3 个源为例某一个采样时刻上两个方程三个未知数初中生都能看出解不唯一。ICA 之所以能解决传统盲源分离是因为它假设混合矩阵可逆然后通过统计独立性找到一个可逆的分离矩阵 W让 y Wx 的各分量尽可能独立。这个思路在 m ≥ n 时是自洽的因为 W 有足够的自由度把每个源分别“揪”出来。可当 m nW 再怎么做也只有 m 个输出通道不可能线性表示出 n 个相互独立的源。所以“ICA 在欠定场景翻车”不是实现细节问题而是前提失效。很多刚转过来的人还抱着“多试几个 ICA 变体”的念头实际效果要么是几个输出通道里混着多个源要么是算法不收敛。做欠定盲源分离的第一课就是接受这个现实不能再用“求逆”的思路必须引入额外信息。这个额外信息工程上最容易落地的是稀疏性。2.2 稀疏性救场为什么时频域比时域更适合做 SCAN稀疏性的意思不难理解不是所有源在所有时刻都同时活跃也不是所有频率上都同时有能量。语音、机械振动、生物电信号这类天然信号在时域上往往叠成一团但在时频域里能量却会集中在少数区域。拿语音来说一个时频点某个时间窗、某个频率通常只由说话人的一个谐波分量主导拿齿轮箱振动来说某一阶啮合频率也往往只对应一个故障源。这个特性直接改变了欠定问题的难度。假设某一个时频点 (f, t) 只有一个源 k 起主导作用那么观测的 STFT 向量 z(f,t) ≈ a_k · s_k(f,t)。也就是说这个时频点的观测向量在方向上和混合矩阵的第 k 列一致。把所有这样的时频点放到一起看它们会在方向空间里聚成若干个簇每个簇对应 A 的一列。于是“估计 A”这个难题被转化成了“在方向空间里做聚类”这个成熟问题。时域里为什么不行因为任意一个采样时刻多个源几乎必然同时发生时域波形看不出谁是主导。时频域之所以行是因为频率分辨率把信号摊开了每个点的活跃源数目大幅下降。实际调试时判断数据适不适合 SCAN最粗暴的办法就是看 STFT 幅值分布如果大部分能量集中在一小部分时频点上这条路就值得走如果时频图上一片均匀的噪声底稀疏性不成立后面所有操作都是白搭。2.3 SCAN 的两步解耦先估混合矩阵再恢复源SCAN 在工程上不是一个封装好的黑匣子它更像一套固定流程的统称我一般按两步来落地。第一步是混合矩阵估计把观测信号做 STFT提取每个时频点的方向特征用聚类算法把 A 的列方向找出来。第二步是源恢复已知 A 之后对每一个时频点求解一个稀疏约束的欠定线性问题再把时频结果逆变换回时域。两步解耦的最大好处是出错能定位。如果恢复出来的源有串音可能是 A 估计偏了如果源本身断续、有金属声可能是时频掩码的边界处理问题。每一步都有独立的验证手段不用把一堆非线性迭代搅在一起。相比之下联合估计 A 和 s 的目标函数是非凸的对初始值極度敏感参数调起来像是在摸黑。工程上我几乎只在论文对比实验里才会用联合优化的方法实际项目里两步法已经足够稳而且每一步都可以被替换成更高效的算法。SCAN 的完整处理流程可以概括成这五步对观测做 STFT在每个时频点上提取方向特征用聚类估计混合矩阵 A已知 A 后用掩码或 l1 最小化恢复源用 ISTFT 回到时域。下面两章就照着这五步把代码写出来。3. 用 SCAN 估混合矩阵STFT 加方向聚类的最小可跑代码要验证 SCAN 能不能用我一般先不上真实数据。真实数据的问题是你永远不知道标准答案分离结果对不对只能靠耳朵或谱图猜。合成数据是先把“3 个源、2 个传感器、已知混合矩阵”这些条件定死让算法去还原这样 A 估得好不好、源恢复得准不准一目了然。下面这套代码就是干这个事的。3.1 用合成信号和混合矩阵构造欠定观测import numpy as np from scipy.signal import stft from sklearn.cluster import KMeans np.random.seed(0) # 固定随机种子保证结果可复现 fs 8000 # 采样率8 kHz 够演示 dur 4 # 4 秒信号 t np.arange(dur * fs) / fs # 三个源不同基频的调幅信号包络不完全重叠 s1 np.sin(2 * np.pi * 230 * t) * (0.6 0.4 * np.sin(2 * np.pi * 1.3 * t)) s2 np.sin(2 * np.pi * 510 * t) * (0.5 0.5 * np.sin(2 * np.pi * 2.1 * t)) s3 np.sin(2 * np.pi * 920 * t) * (0.7 0.3 * np.sin(2 * np.pi * 0.7 * t)) S np.vstack([s1, s2, s3]) # 3 x N三个源 # 2 个传感器、3 个源m n这就是欠定场景 A_true np.array([[0.80, 0.50, 0.35], [0.30, 0.70, 0.90]]) X A_true S 0.002 * np.random.randn(2, len(t)) # 2 x N两路观测代码里三个源选在不同的基频段目的是让它们在时频域里稀疏性足够好。调幅包络则是为了制造非平稳性否则纯正弦信号在 STFT 后能量太集中聚类会过于简单掩盖真实数据里的问题。A_true 的三列方向差异要拉开0.8/0.3、0.5/0.7、0.35/0.9这样三个源在方向空间里分得开如果两列几乎平行那再好的算法也难分那不是算法问题是传感器布局问题。噪声幅度 0.002 乘以 randn 是相对信号很低的水平放到真实采集里大概相当于中等偏上的信噪比。你可以把它调大到 0.01 或者 0.05观察聚类什么时候开始崩这对理解算法边界很有帮助。3.2 特征提取把二维 STFT 变成方向点云# 每通道 STFT窗长 256、半窗重叠hann 窗满足 COLA 重建条件 f, t_stft, Z0 stft(X[0], fsfs, nperseg256, noverlap128, windowhann) _, _, Z1 stft(X[1], fsfs, nperseg256, noverlap128, windowhann) # 只保留能量较大的时频点降低噪声和混叠点的干扰 mag np.abs(Z0) np.abs(Z1) strong mag np.percentile(mag, 85) # 方向特征第一通道能量占比范围 0~1 ratio np.abs(Z0[strong]) / (np.abs(Z0[strong]) np.abs(Z1[strong]) 1e-10) feat ratio[:, None]STFT 的窗长选了 256时域分辨率 32 ms频率分辨率 31.25 Hz对这三个基频在 230/510/920 Hz 的信号来说足够区分。如果你处理的是低频机械振动窗长可以放到 1024 或更大但代价是时间分辨率下降后面会专门讲这个坑。方向特征我直接用了“第一通道能量占比”。对某个单源主导的时频点这个占比近似等于 a0k² / (a0k² a1k²)不同类型的源会落在不同的数值上。这本质上是在一维方向空间里做聚类简单到不会出幺蛾子。如果你做的是麦克风阵列还可以加入通道间的相位差构成二维特征这里先不展开。strong这个掩码很关键。它把能量低于 85 分位的时频点全部滤掉。这些低能量点里噪声占比高方向信息被噪声污染聚进去反而会带来大量离群点。3.3 KMeans 聚类估计 A 的列方向# KMeans 聚类簇数设为源数 3 km KMeans(n_clusters3, initk-means, n_init20, random_state0).fit(feat) centers km.cluster_centers_ # 每簇的典型能量占比 A_hat np.stack([centers[:, 0], 1 - centers[:, 0]], axis1) # 3 x 2 print(A_hat)KMeans 在这种一维特征上表现稳定但要注意两点。第一n_init20是必须的不然 KMeans 的随机初始中心可能让其中一簇被拆成两半这个现象在高维特征里更常见。第二簇的顺序是随机的第一簇未必对应真实源 1后续恢复源的时候要按频率或能量排序去对或者干脆只看分离效果不管顺序。代码里A_hat每一行是对应混合矩阵某列的归一化方向把能量占比还原成两通道幅值比。因为 A 的列绝对幅度在盲源分离里本来就无法确定所以这里不追求还原 A 的绝对尺度。跑完可以打印一下看看三个 center 应该分别接近 0.73、0.42、0.28 附近对应的 A_hat 行就是 (0.73, 0.27)、(0.42, 0.58)、(0.28, 0.72)和真实的 A_true 列方向是吻合的。如果你用真实采集数据这一步的 A_hat 就是后续所有分离的基础一定要先画图确认簇是清晰的再往下走。提示源数 K 在实测数据里往往未知。可以先画方向特征的直方图数峰也可以用 DBSCAN 做一次探测没必要一上来就 KMeans 硬指定。4. 已知混合矩阵之后时频掩码恢复源的最短实现A 估计出来了源恢复看起来就是个普通反问题但实际没这么简单。已知 A 的情况下每个时频点依然要解 z A s 这样一个 m 个方程、n 个未知数的欠定方程组。直接最小二乘会把误差摊到所有源上结果每一个源都带着别人的残影。恢复这一步还是要靠稀疏性来兜底。4.1 恢复阶段为什么不能直接求逆如果 m2、n3伪逆解出来的 s 在数学上是“能量最小”的解不代表真实源。真实语音或振动信号的源往往只在少数时频点活跃真正合理的解应该是“这个点主要属于谁就把它给谁”而不是“每个源各分一点”。这就是时频掩码的思路。时频掩码的做法非常直接既然聚类阶段已经知道每个强时频点属于哪个方向簇那就把这些点整体分配给对应的源不属于它的点幅值置 0。硬掩码是 0/1 二值软掩码是 0~1 连续权重。硬掩码实现简单软掩码边界更干净实际项目里我一般从硬掩码起步确认流程通了再换软掩码。4.2 硬掩码恢复十几行代码的实用做法from scipy.signal import istft # 把标签映射回二维时频坐标 f_idx, t_idx np.where(strong) # 顺序与 Z0[strong] 的扁平顺序一致 rec np.zeros((3, Z0.shape[0], Z0.shape[1]), dtypecomplex) for k in range(3): sel np.zeros_like(strong) sel[f_idx[km.labels_ k], t_idx[km.labels_ k]] True rec[k] Z0 * sel # 只保留该簇的时频点 _, x_rec istft(rec[k], fsfs, nperseg256, noverlap128, windowhann) print(fsource {k} recovered, peak abs {np.max(np.abs(x_rec)):.4f})这里的sel是当前源对应的时频掩码把观测第一通道的 STFT 系数按掩码筛选后做 ISTFT。因为每个时频点只归属一个源ISTFT 之后得到的是这个源在第一通道上的贡献也就是 a0k · s_k 的近似。所以恢复出来的源幅值是被 A 的第一行系数缩放过的这是欠定盲源分离的尺度不定性不是 bug。如果一定需要真实幅值必须用标定信号去定标算法本身给不了这个信息。ISTFT 之所以能重建得比较干净是因为 hann 窗加 50% 重叠满足 COLA 条件也就是重叠相加的窗函数和为常数不会引入额外的幅度调制。如果你自己换了窗函数必须确认它满足这个约束否则 recovery 出来的信号会有周期性起伏。4.3 什么时候改用 l1 最小化软掩码与正则化硬掩码的适用场景是“每个时频点确实只有一个源主导”。当源之间有较多频谱重叠、或者噪声底偏高时硬掩码会把两个源的交叉点武断地分给其中一方产生梳状滤波听感。这时候有两个升级方向。第一个升级是软掩码。算出每个点时各个簇心的距离转成连续权重而不是 0/1。这样做的好处是边界时频点的能量被按比例分配听感上的金属声会明显降低。代价是串音比硬掩码稍高需要在掩码权重公式里调一个软化系数常见做法是加一个以距离为分母的软阈值系数取 0.05 到 0.2 之间。第二个升级是逐点做 l1 最小化。对单个时频点 z求 s 使 z A s 且 s 的 l1 范数最小import cvxpy as cp zf np.array([Z0[10, 20], Z1[10, 20]]) # 随便取一个强点 s_var cp.Variable(3, complexTrue) prob cp.Problem(cp.Minimize(cp.norm(s_var, 1)), [A_true s_var zf]) prob.solve()这个做法的理论意义是在源稀疏的前提下l1 最小化能给出更接近真实源的最优解。但实际工程里它有一个很大的问题——慢。4 秒信号有上千个强时频点逐点解优化问题要跑很久。我一般只在验证算法上限时用 l1工程落地还是以掩码为主如果有必要再对特定频段局部做 l1 精修。5. SCAN 落地避坑五个让结果面目全非的细节这一章全是参数和细节的教训每一条都是我至少踩过一次的真问题。按“现象、原因、解决”写方便你对照排查。5.1 聚类数 K 设错KMeans 把两个方向卷成一个现象估计出的混合矩阵少一列恢复之后有两个源严重串音另一个源反而特别干净。原因KMeans 必须指定 K。但方向点云里“应该有几个簇”取决于时频点里真正单源主导的比例不一定等于源数。源 2 和源 3 方向太接近、或者某个源能量太弱KMeans 就会把两簇合并。簇数设得比实际源数多又会把一簇硬拆成两半。解决先画方向特征的直方图数一下峰的个数再定 K。峰不够明显时用 DBSCAN 做一次方向聚类它不需要预定簇数。工程上更稳的做法是把 K 设大一号聚类完成后再按簇心距离合并相近的簇避免一上来就被 K 绑架。5.2 强点筛选阈值把噪声点也当成方向点现象A 估出来之后簇心分散和前面对不上恢复源里全是毛刺。原因低幅值时频点的信噪比低。z a·s v 里 s 很弱v 主导方向信息是噪声的方向不是 A 列的方向。把这些点拿进聚类等于在方向点云里撒了一大把随机点KMeans 的簇心会被拽偏。解决用能量分位数筛掉弱点时85%~90% 是比较常见的区间。分位数太高会丢掉低频强点太低则噪声混入。更稳的做法是先估计噪声底比如用频谱最低 5% 的幅值均值做底再设一个“噪声底乘倍数”的绝对阈值。这个倍数一般取 5~10比你换不同分位数更容易保持参数稳定。5.3 STFT 窗长一改可分离度就大变现象窗长 256 时 A 估得挺准改成 128 或者 1024 之后方向特征图变成一团分离结果完全没法听。原因窗长直接决定了时频里“每个点由几个源主导”。窗太短频率分辨率差两个不同频率的源在同一个频率箱里打架单源占优假设失效窗太长时间分辨率差一个窗内可能横跨多个源的瞬变事件同样破坏假设。这个 tradeoff 不是玄学是短时傅里叶变换的基本性质。解决平稳信号电机稳态振动、持续音可以用长窗窗长 1024 甚至 2048 都行瞬变多的信号语音、冲击故障用短窗256 左右起步。最靠谱的方式是固定其他参数把窗长按 128、256、512、1024 各跑一遍看简化 SIR 分数或者方向直方图峰值对应的窗长就是这套数据的最佳值。注意采样率不同同样的窗长对应的实际物理时长完全不同跨采样率对比时应该按毫秒来选窗。5.4 硬掩码恢复出来的源有“闷罐子”听感现象分离结果在频谱上看着干净但听起来发闷、不连续像隔着棉被说话。原因硬掩码把每个时频点整块归给一个源边界上幅值跳变重构时在掩码边界产生伪影。另外重建只用了一个通道的 STFT其实浪费了第二通道的信息而且相位保留了混合观测的相位不是源本身的相位听觉上就会发闷。解决把二值掩码换成软掩码掩码在时频方向上做一次 3×3 的中值滤波边界过渡平滑很多。重建时不要只取第一通道可以用两通道做加权合成比如按各自簇心方向的比例把能量合并相当于做了一次极简波束形成。这样能多挽回一些信噪比。相位问题在欠定下基本无解能做的是在掩码阶段减少边界效应。5.5 尺度不定性别拿归一化当万能后悔药现象恢复出来的源时域幅值忽大忽小和原始源完全对不上归一化到 0~1 之后看着正常但一到定量分析就露馅。原因欠定盲源分离本身有尺度不定性即使 A 估计完全正确也只能确定 A 列的方向不能确定每一列的模长。源 s 和 A 的列可以互相补偿A 第 k 列乘 0.5、s_k 乘 2观测完全不变。所以分离结果天然缺少物理量纲。解决需要真实幅值时必须在系统里加标定环节。比如让已知幅值的源依次单独工作标出每个通道的增益系数最后乘回去。没有标定条件时SCAN 适合做故障检测、语音内容分离、模式识别这类相对比较的任务不适合直接做幅值计量。别在分离后强行把信号拼到某个幅值范围那只是自欺欺人掩盖不了一个被缩放了 10 倍的事实。6. 用方向直方图和 SIR 验证 SCAN 没做偏进阶检验技巧调试 SCAN 的时候我最常犯的错是盯着分离出的时域波形看半天也看不出是 A 估错还是恢复参数不对。后来养成两个习惯问题定位快很多。第一个习惯是先画方向直方图。把第 3 章里提的方向特征ratio放到直方图里如果能看到三个清晰的峰说明稀疏性前提成立A 有救如果是一团平滑的大包别急着调参数先去检查源本身是否稀疏。直方图比聚类结果更直观因为它不受 KMeans 初始化影响是一个纯数据视角。import matplotlib.pyplot as plt plt.hist(ratio, bins80, edgecolornone) plt.xlabel(channel 0 energy ratio) plt.ylabel(count) plt.show()第二个习惯是算简化 SIR 分数。在合成数据里每个源都是已知的分离结果可以直接和参考源对比。工程上我常用一段很短的函数先把时延对齐、幅值对齐再算能量比。分数超过 10 dB 基本可听15 dB 以上算清晰。def simple_sir(ref, est): c np.correlate(est, ref, modefull) lag np.argmax(np.abs(c)) - (len(ref) - 1) est_shift np.roll(est, -lag) alpha np.dot(ref, est_shift) / (np.dot(est_shift, est_shift) 1e-10) est_shift * alpha noise ref - est_shift return 10 * np.log10(np.dot(ref, ref) / (np.dot(noise, noise) 1e-10))算分数时要固定变量一次只扫一个参数。我一般固定源数和聚类数然后扫窗长和强点阈值看分数随参数的变化曲线。峰值所在的窗长就是当前数据的合理取值这个值在换数据后很可能会变需要重新扫。这些年做欠定分离踩过最深的坑是急着调恢复算法而没确认 A 估得准不准。现在我一般先把方向直方图打出来再跑一遍参数扫描确认前提成立之后才继续调掩码或 l1。这套验证习惯比任何参数技巧都管用。希望帮到你。本文还有配套的精品资源点击获取