简介：这份资源聚焦欠定盲源分离（UBSS）这一信号处理方向，面向音频处理、医学成像、通信解码及金融数据分析等领域的学习者与研究人员，帮助理解观测通道少于源信号数量时的分离思路。压缩包共5个文件，以4个m脚本文件和1个txt说明文件为主，脚本涵盖SCAN、DBSE、sobi、vecsymm等算法实现，readme则提供使用说明，整体约5KB，体量轻便，适合快速上手与二次修改。资源围绕独立成分分析、稀疏表示、联合诊断学习及机器学习等主流方法展开，可配合音频、脑电等场景验证分离效果。目前已有254人学习下载，适合作为算法入门与实验对照的参考材料，便于读者在理解混合模型与线性系统的基础上，进一步探索噪声抑制、非高斯假设限制等挑战，提升盲信号分离的实践能力。

1. SCAN 欠定盲源分离：当观测通道比源信号还少时怎么把信号拆开

做设备故障诊断或者无线频谱监测的工程师迟早会撞上这个场景：手上只有两三个传感器，现场却混着四五个甚至更多独立信号源，想把它们一一拆出来。这就是欠定盲源分离要解决的问题——观测通道数少于源信号数，传统独立成分分析直接失效，因为混叠矩阵连方阵都不是，求逆这条路根本走不通。SCAN 这类基于稀疏分量分析的方案，核心思路是利用信号在某个变换域里的稀疏性，先估计混叠矩阵，再逐个恢复源信号。它适合做机械振动多故障源分离、脑电多通道伪迹剔除、通信信号盲解调前处理的从业者。下面把我实际落地时踩过的路讲清楚，从数据模型到代码复现到参数调优，能直接抄。

2. 欠定盲源分离的数学底子与 SCAN 的切入点

2.1 为什么观测通道少于源信号就变成了病态问题

先把这个问题的数学形式摆出来。标准线性瞬时混合模型是x(t) = A·s(t) + n(t)，其中 x(t) 是 M 维观测向量，s(t) 是 N 维源向量，A 是 M×N 的混叠矩阵。适定和超定情况下 M ≥ N，A 是方阵或高矩阵，可以用独立成分分析（ICA）通过最大化非高斯性来估计 A 的逆。但欠定情况下 M < N，A 是扁矩阵，它的零空间非平凡，意味着有无数个 s(t) 能产生同一组观测 x(t)。这就是欠定盲源分离的根本困难：解不唯一。

那为什么还能做？关键在于引入额外假设。最常用也最有效的假设是稀疏性——源信号在某个变换域（时域、频域、时频域）中，任意时刻只有少数几个源是活跃的。在这个假设下，虽然全局解不唯一，但每个时刻的局部混合可以看成是少数几个源列的线性组合，几何上表现为观测点聚集在由 A 的列向量张成的低维子空间里。SCAN（Sparse Component Analysis via clustering and Nullspace，也有文献叫 Sparse Component Analysis Network）正是利用这个几何结构：先通过聚类估计混叠矩阵 A 的列方向，再用最短路径或线性规划恢复源信号。

2.2 SCAN 方法的两阶段框架：混叠矩阵估计与源恢复

SCAN 的流程分两大阶段，每个阶段都有独立的算法选择和参数敏感性。

第一阶段：混叠矩阵估计。对观测信号做短时傅里叶变换（STFT），得到时频域表示 X(t,f)。在稀疏假设下，每个时频点上只有一个源主导，此时观测向量 X(t,f) 的方向就近似等于 A 中对应列向量的方向。把所有时频点的方向向量画在单位超球面上，它们会聚成 N 个簇，每个簇的中心方向就是 A 的一列。聚类算法常用 K-means 或者更鲁棒的 DBSCAN，后者不需要预先知道源数目，对噪声点也有更好的容忍度。

第二阶段：源信号恢复。有了 A 的估计值  之后，每个时频点上求解一个欠定方程 x = ·s。由于稀疏性，可以用 L1 范数最小化来求最稀疏解：min ||s||₁ subject to x = ·s。这本质上是一个线性规划问题，可以用基追踪（Basis Pursuit）或者匹配追踪（Matching Pursuit）系列算法求解。最后对每个源的时频表示做逆 STFT，得到时域波形。

整个框架听起来清晰，但实际落地时每一步都有坑。下面用 Python 代码把最小可复现版本跑通。

2.3 用 Python 跑通 SCAN 的最小复现流程

先安装依赖，核心是 numpy、scipy 和 scikit-learn（用于聚类）：

pip install numpy scipy scikit-learn matplotlib

然后生成模拟数据并执行完整的 SCAN 流程：

import numpy as np from scipy.signal import stft, istft from sklearn.cluster import KMeans np.random.seed(42) # ========== 1. 生成模拟源信号 ========== fs = 1000 # 采样率 Hz T = 2.0 # 时长秒 t = np.arange(0, T, 1/fs) N = 4 # 源信号数（欠定：源多于观测） # 四个源：两个正弦、一个方波、一个调频 s1 = np.sin(2 * np.pi * 50 * t) s2 = np.sin(2 * np.pi * 120 * t) s3 = np.sign(np.sin(2 * np.pi * 30 * t)) s4 = np.sin(2 * np.pi * (200 + 100 * t) * t) S = np.array([s1, s2, s3, s4]) # shape: (4, len(t)) # ========== 2. 构造欠定混合 ========== M = 2 # 观测通道数（M < N，欠定） A_true = np.array([[0.8, 0.3, 0.5, 0.2], [0.2, 0.7, 0.4, 0.9]]) # 2x4 混叠矩阵 X = A_true @ S + 0.01 * np.random.randn(M, len(t)) # 加少量噪声 # ========== 3. STFT 变换到稀疏域 ========== f, t_stft, Zxx = stft(X, fs=fs, nperseg=128, noverlap=96) # Zxx shape: (M, freq_bins, time_frames) # ========== 4. 提取方向向量并聚类估计混叠矩阵 ========== # 取能量足够大的时频点，避免噪声主导方向 power = np.sum(np.abs(Zxx)**2, axis=0) # (freq_bins, time_frames) threshold = np.percentile(power, 70) # 只保留能量前30%的点 mask = power > threshold # 提取方向向量（归一化到单位球面） directions = [] for fi in range(Zxx.shape[1]): for ti in range(Zxx.shape[2]): if mask[fi, ti]: v = Zxx[:, fi, ti] norm = np.linalg.norm(v) if norm > 1e-10: directions.append(v / norm) directions = np.array(directions) # shape: (num_points, M) # K-means 聚类，簇数设为源数 N kmeans = KMeans(n_clusters=N, n_init=10, random_state=42) labels = kmeans.fit_predict(directions) A_est = kmeans.cluster_centers_.T # shape: (M, N) # 修正列顺序和符号（与真实 A 对齐，仅用于评估） from scipy.optimize import linear_sum_assignment cost = np.abs(A_est.T @ A_true) row_ind, col_ind = linear_sum_assignment(-cost) A_est = A_est[:, col_ind] print("真实混叠矩阵 A:\n", A_true) print("估计混叠矩阵 Â:\n", np.round(A_est, 3))

这段代码的关键逻辑说明：

• STFT 参数选择：nperseg=128对应频率分辨率约 7.8 Hz，noverlap=96给出 75% 重叠，保证时频面上有足够密的采样点用于聚类。窗长太短频率分辨率不够，太长则时域稀疏性被破坏。
• 能量阈值：np.percentile(power, 70)只保留能量最高的 30% 时频点。这是欠定盲源分离里非常关键的一步——低能量点基本是噪声，方向向量随机散布，会严重干扰聚类。阈值太高则有效点太少，聚类不稳定。
• 方向向量归一化：把每个时频点的观测向量投影到单位球面上，这样聚类只关心方向不关心幅度，符合稀疏假设下“单源主导”的几何直觉。
• K-means 的局限：这里用了 K-means 是因为源数已知。实际场景中源数往往未知，需要用 DBSCAN 或者基于特征值衰减的估计方法。另外 K-means 对初始值敏感，n_init=10是必要的。

跑完这一步，你会看到估计的混叠矩阵和真实值已经比较接近。接下来是源恢复。

# ========== 5. L1 范数最小化恢复源信号 ========== from scipy.optimize import linprog def l1_min_recover(x_obs, A_est): """对单个时频点求解 min ||s||_1 s.t. A_est @ s = x_obs""" M, N = A_est.shape # linprog 标准形式: min c^T z, s.t. A_eq z = b_eq, bounds # 变量拆分: s = u - v, u,v >= 0, |s| = u + v c = np.ones(2 * N) A_eq = np.hstack([A_est, -A_est]) b_eq = x_obs.real # 先处理实部 res = linprog(c, A_eq=A_eq, b_eq=b_eq, bounds=[(0, None)] * (2 * N), method='highs') if res.success: s = res.x[:N] - res.x[N:] return s return np.zeros(N) # 对每个时频点恢复 S_est_stft = np.zeros((N, Zxx.shape[1], Zxx.shape[2]), dtype=complex) for fi in range(Zxx.shape[1]): for ti in range(Zxx.shape[2]): x_point = Zxx[:, fi, ti] s_real = l1_min_recover(x_point.real, A_est) s_imag = l1_min_recover(x_point.imag, A_est) S_est_stft[:, fi, ti] = s_real + 1j * s_imag # ========== 6. 逆 STFT 恢复时域波形 ========== S_est = np.zeros((N, len(t))) for i in range(N): _, s_rec = istft(S_est_stft[i], fs=fs, nperseg=128, noverlap=96) S_est[i, :len(t)] = s_rec[:len(t)] # 计算恢复信号与真实源的相关系数 for i in range(N): corr = np.corrcoef(S_est[i], S[i])[0, 1] print(f"源 {i+1} 恢复相关系数: {corr:.4f}")

这段恢复代码有几个实操要点：

• 逐点线性规划非常慢。上面这个双重循环在 2 秒数据上可能要跑几分钟。实际工程中会用匹配追踪（OMP）或者迭代硬阈值（IHT）来加速，牺牲一点精度换十倍以上的速度。
• 实部和虚部分开处理是一种简化。更严谨的做法是把复数线性规划转成实数形式统一求解，但代码复杂度会上升。
• method='highs'是 scipy 目前最快的线性规划求解器，比默认的 simplex 快很多。如果 scipy 版本较老没有这个选项，升级到 1.6 以上。
• 相关系数是评估恢复质量的直接指标。如果某个源的相关系数低于 0.8，通常意味着聚类阶段把两个源的方向混在了一起，或者该源在时频域不够稀疏。

3. 混叠矩阵估计的精度决定一切：聚类参数怎么调

3.1 源数目未知时用 DBSCAN 替代 K-means

上一章用了 K-means，前提是知道源有几个。实际项目中源数往往未知，这时候 DBSCAN 是更务实的选择。它不需要预设簇数，还能把噪声点标记为 -1 直接丢弃。代价是要调两个参数：邻域半径eps和最小样本数min_samples。

from sklearn.cluster import DBSCAN from sklearn.neighbors import NearestNeighbors # 用 k-距离图辅助确定 eps k = 5 nbrs = NearestNeighbors(n_neighbors=k).fit(directions) distances, _ = nbrs.kneighbors(directions) k_dist = np.sort(distances[:, -1]) # 找 k-距离曲线的拐点作为 eps # 实际使用时画出 k_dist 曲线肉眼找拐点，这里用差分近似 diff = np.diff(k_dist) eps_auto = k_dist[np.argmax(diff) + 1] print(f"自动估计的 eps: {eps_auto:.4f}") db = DBSCAN(eps=eps_auto, min_samples=10) labels_db = db.fit_predict(directions) n_clusters = len(set(labels_db)) - (1 if -1 in labels_db else 0) print(f"DBSCAN 发现 {n_clusters} 个簇") # 用各簇均值作为混叠矩阵列向量 A_est_db = np.zeros((M, n_clusters)) for i in range(n_clusters): cluster_points = directions[labels_db == i] center = cluster_points.mean(axis=0) A_est_db[:, i] = center / np.linalg.norm(center)

eps的选取是 DBSCAN 最玄学的地方。k-距离图是个靠谱的辅助工具：把每个点到第 k 近邻的距离排序画出来，曲线拐点对应的距离就是合适的eps。min_samples一般取观测维度的 2 到 3 倍，比如 M=2 时取 5 到 10。如果聚类结果里 -1 标签占比超过 30%，说明eps太小，大量有效点被当成噪声丢掉了。

3.2 时频点筛选：能量阈值和单源主导判据

方向向量聚类的前提是“每个时频点只有一个源主导”。但现实中这个条件并不总是满足，尤其在两个源频率重叠的区域。除了能量阈值，还可以加一个单源主导判据：计算该时频点的观测向量与最近簇中心的余弦相似度，低于某个阈值（比如 0.95）就丢弃。

# 在聚类完成后，用簇中心反筛时频点 cos_threshold = 0.95 valid_points = [] valid_indices = [] for idx, (d, label) in enumerate(zip(directions, labels)): if label == -1: continue center = kmeans.cluster_centers_[label] cos_sim = np.dot(d, center) / (np.linalg.norm(d) * np.linalg.norm(center)) if cos_sim > cos_threshold: valid_points.append(d) valid_indices.append(idx) print(f"原始有效点: {len(directions)}, 筛选后: {len(valid_points)}")

这个二次筛选能显著提升混叠矩阵估计精度，代价是可用于恢复的时频点变少，可能影响源恢复的完整性。实践中cos_threshold取 0.9 到 0.98 之间，根据源之间的频率重叠程度来定。重叠严重就放宽到 0.9，源之间频率分得很开可以严到 0.98。

3.3 混叠矩阵列顺序模糊与幅度歧义的工程处理

盲源分离有两个固有的模糊性：排列模糊和尺度模糊。排列模糊是指估计出来的 A 的列顺序和真实源顺序对不上，尺度模糊是指每列可以乘以任意非零常数而观测不变。这两个模糊性在数学上无法消除，但工程上可以处理。

排列模糊的常见解法是：如果你有每个源的先验频率范围或者波形特征，可以在恢复后对每个源做频谱分析，按特征匹配重新排序。尺度模糊则更简单——大多数应用只关心波形形状不关心绝对幅度，把每个恢复源归一化到单位方差即可。如果确实需要绝对幅度，可以用一个已知幅度的参考源做标定。

# 按主频排序解决排列模糊 from scipy.fft import fft, fftfreq def dominant_freq(signal, fs): """返回信号的主频""" spectrum = np.abs(fft(signal)) freqs = fftfreq(len(signal), 1/fs) positive = freqs > 0 return freqs[positive][np.argmax(spectrum[positive])] # 假设已知真实源的主频排序 true_freqs = [dominant_freq(s, fs) for s in S] est_freqs = [dominant_freq(S_est[i], fs) for i in range(N)] # 按主频匹配重排 order = np.argsort(est_freqs) S_est_sorted = S_est[order] print("重排后各源主频:", [f"{dominant_freq(S_est_sorted[i], fs):.1f}" for i in range(N)])

这段代码展示了用主频特征解决排列模糊的思路。实际项目中特征可以更丰富：谱峰位置、包络形状、循环平稳特征等。关键是找到对尺度变化不敏感、对源身份有区分度的特征。

4. 源恢复阶段的加速与精度权衡

4.1 匹配追踪替代线性规划：快十倍但精度掉多少

逐点线性规划在长数据上不可接受。匹配追踪（Matching Pursuit）系列是更实用的选择。正交匹配追踪（OMP）每次迭代选一个与残差最相关的原子，然后正交投影更新残差。在欠定盲源分离里，“原子”就是混叠矩阵的列向量。

def omp_recover(x_obs, A_est, sparsity=2): """正交匹配追踪恢复，sparsity 为预期活跃源数""" M, N = A_est.shape residual = x_obs.copy() indices = [] s_est = np.zeros(N, dtype=complex) for _ in range(sparsity): # 找与残差最相关的列 correlations = np.abs(A_est.conj().T @ residual) # 排除已选列 for idx in indices: correlations[idx] = 0 new_idx = np.argmax(correlations) indices.append(new_idx) # 最小二乘更新 A_sub = A_est[:, indices] s_sub = np.linalg.lstsq(A_sub, x_obs, rcond=None)[0] residual = x_obs - A_sub @ s_sub for i, idx in enumerate(indices): s_est[idx] = s_sub[i] return s_est # 用 OMP 替换线性规划，速度对比 import time start = time.time() S_omp = np.zeros((N, Zxx.shape[1], Zxx.shape[2]), dtype=complex) for fi in range(Zxx.shape[1]): for ti in range(Zxx.shape[2]): S_omp[:, fi, ti] = omp_recover(Zxx[:, fi, ti], A_est, sparsity=2) print(f"OMP 耗时: {time.time() - start:.2f} 秒")

sparsity参数控制每次恢复时假设有几个源同时活跃。在稀疏假设下这个值通常取 1 到 3。取太大会引入虚假源，取太小会丢失同时活跃的源。一个经验做法是先用 1 跑一遍，看恢复信号的相关系数，如果某些源质量差就增加到 2 再试。

OMP 的速度比线性规划快一个数量级以上，但精度在源之间相关性高时会下降。如果对精度要求极高且数据量不大，线性规划仍然是金标准。

4.2 时频掩膜法：更简单但依赖混叠矩阵精度

如果混叠矩阵估计得足够准，还有一个更快的路子：时频掩膜。对每个时频点，计算它属于每个源的概率（基于到各簇中心的距离），然后做软掩膜或者硬掩膜。

def tf_mask_recover(Zxx, A_est): """基于时频掩膜的源恢复""" M, N = A_est.shape F, T = Zxx.shape[1], Zxx.shape[2] S_masked = np.zeros((N, F, T), dtype=complex) for fi in range(F): for ti in range(T): x_point = Zxx[:, fi, ti] # 计算到每个源方向的投影能量 energies = np.zeros(N) for n in range(N): a_n = A_est[:, n:n+1] proj = a_n @ np.linalg.pinv(a_n) @ x_point energies[n] = np.linalg.norm(proj)**2 # 软掩膜：按能量比例分配 total = energies.sum() if total > 1e-10: mask = energies / total for n in range(N): S_masked[n, fi, ti] = mask[n] * x_point @ np.linalg.pinv(A_est[:, n:n+1]) @ A_est[:, n:n+1] return S_masked

时频掩膜法的恢复质量高度依赖混叠矩阵的精度。如果 A 估计有 5% 以上的误差，掩膜就会把源之间的能量搞混，恢复信号里出现明显的串扰。它的优势是速度快，没有迭代求解，适合做实时处理的预处理。

4.3 恢复质量评估：相关系数、信干比与主观听感

评估恢复质量不能只看一个指标。相关系数反映波形形状的相似度，但对幅度不敏感。信干比（SIR）反映源之间的串扰程度。如果是语音信号，还得靠人耳听。

def evaluate_separation(S_true, S_est): """计算各源的相关系数和信干比""" N = S_true.shape[0] results = [] for i in range(N): # 相关系数 corr = np.corrcoef(S_true[i], S_est[i])[0, 1] # 信干比：目标源能量 / 其他源泄漏能量 # 用最小二乘拟合找泄漏系数 coeffs = np.linalg.lstsq(S_true.T, S_est[i], rcond=None)[0] target_energy = (coeffs[i] * S_true[i])**2 interference = sum((coeffs[j] * S_true[j])**2 for j in range(N) if j != i) sir = 10 * np.log10(target_energy.sum() / (interference.sum() + 1e-10)) results.append({'source': i+1, 'corr': corr, 'SIR_dB': sir}) return results for r in evaluate_separation(S, S_est): print(f"源 {r['source']}: 相关系数={r['corr']:.4f}, SIR={r['SIR_dB']:.2f} dB")

相关系数高于 0.9 且 SIR 高于 15 dB 通常算合格。如果 SIR 低于 10 dB，说明串扰严重，需要回头检查混叠矩阵估计或者增加稀疏性约束。

5. 避坑与排查：欠定盲源分离翻车实录

5.1 聚类簇数设错导致混叠矩阵完全跑偏

现象：恢复出来的信号和真实源完全不相关，相关系数在 0.2 以下，但代码没有报错。

原因：K-means 的n_clusters设成了错误的源数。如果实际有 4 个源但设成 3，K-means 会把两个方向相近的源合并成一个簇，导致混叠矩阵少一列，后续恢复时这两个源被当成一个处理。设成 5 则多出一个虚假簇，恢复时引入不存在的源。

解决：源数未知时用 DBSCAN，或者用特征值分解方法估计源数——对观测协方差矩阵做特征值分解，看特征值衰减的拐点。另一个实用技巧是跑多个不同簇数的 K-means，比较聚类轮廓系数（silhouette score），选最高的那个。

5.2 时频域稀疏性不足导致方向向量弥散

现象：方向向量在单位球面上不是聚成几个清晰的簇，而是弥散成一片，聚类轮廓系数低于 0.3。

原因：源信号在所选变换域里不够稀疏。比如两个源在大部分时间频率上同时活跃，每个时频点都是多源混合，方向向量就失去了“单源主导”的几何意义。常见于源之间频率重叠严重或者源本身是宽带非稀疏信号。

解决：换变换域。STFT 不稀疏就试小波变换或者同步挤压变换。如果源在时域有稀疏性（比如脉冲信号），直接在时域做聚类。另一个方向是增加观测通道数——哪怕从 2 个增加到 3 个，稀疏性要求都会大幅降低。

5.3 噪声能量估计偏差让阈值筛选失效

现象：能量阈值筛选后有效时频点数量波动很大，换一段数据聚类结果就完全不同。

原因：用了固定的百分位阈值（比如 70%），但不同数据段的噪声水平不一样。噪声大的段，前 30% 的点里也混着大量噪声；噪声小的段，可能把有效信号也筛掉了。

解决：改用自适应阈值。先估计噪声底——取时频能量分布的最低 10% 分位作为噪声水平，然后阈值设为噪声水平的 3 到 5 倍。这样阈值随数据自适应，鲁棒性好很多。

5.4 线性规划恢复时源信号被过度稀疏化

现象：恢复出来的源信号波形被“压扁”了，峰值被削平，听起来有失真。

原因：L1 范数最小化倾向于把小的系数压到零，导致恢复信号过于稀疏。这在源信号本身不是严格稀疏时尤其明显。

解决：改用 L1-L2 混合范数，或者加一个小的 L2 正则项。另一个办法是减小sparsity参数，允许更多源同时活跃。如果对波形保真度要求高，考虑用迭代硬阈值（IHT）替代线性规划，它在稀疏度和保真度之间平衡得更好。

5.5 逆 STFT 时相位不一致导致重建波形畸变

现象：恢复的时域波形在拼接处有跳变，频谱上出现不该有的高频分量。

原因：对每个时频点独立做 L1 恢复时，相邻时频点的相位没有连续性约束。逆 STFT 时相位不连续会导致重建波形出现伪影。

解决：在恢复阶段加相位平滑约束，或者用重叠相加法时确保窗函数满足 COLA（Constant Overlap-Add）条件。更简单的做法是增大 STFT 重叠率到 87.5%（noverlap = 7/8 * nperseg），让相位不连续的影响被平滑掉。

6. 从仿真到实测：SCAN 在机械振动信号上的调参心得

仿真跑通只是第一步，实测信号才是真正的考验。我拿一个双通道加速度传感器采的齿轮箱振动信号试过，里面混着轴承故障冲击、齿轮啮合振动和电机电磁噪声三种源。仿真里表现很好的参数，实测时几乎全要重调。

第一个要改的是 STFT 窗长。仿真里 128 点窗对应 7.8 Hz 分辨率够用，但实测轴承故障冲击的重复频率只有十几 Hz，7.8 Hz 的分辨率根本分不开。我把窗长加到 512 点，频率分辨率降到 2 Hz，冲击的边频带才在时频面上清晰起来。代价是时域分辨率变差，但振动信号本来就是频域特征更关键。

第二个要改的是能量阈值。实测信号的噪声底比仿真高得多，固定百分位阈值筛出来的点里混了大量噪声。我改成先对时频能量做对数变换，然后用 Otsu 方法自动找双峰之间的谷底作为阈值。这个自适应方法在不同工况下都稳定，不用每次手动调。

第三个要改的是聚类距离度量。仿真里用欧氏距离就够了，但实测信号的方向向量在球面上分布不均匀，欧氏距离会把密集区域的点过度合并。换成余弦距离后，簇的分离度明显改善。scikit-learn 的 KMeans 不支持自定义距离，得用sklearn.cluster.SpectralClustering或者自己实现。

最后一个血泪教训：混叠矩阵的估计精度比源恢复算法重要得多。我花了大量时间调 L1 求解器的参数，效果提升有限。后来回头把聚类环节的时频点筛选和距离度量优化了一下，混叠矩阵的列方向误差从 8 度降到 2 度，恢复信号的 SIR 直接涨了 6 dB。所以如果你的 SCAN 效果不好，先检查混叠矩阵估计，别急着换恢复算法。

实测中还有一个仿真里遇不到的问题：传感器通道之间的时间延迟。两个加速度计安装位置不同，同一个冲击到达两个传感器有时间差，这在频域表现为相位随频率线性变化。如果不补偿，混叠矩阵的估计会有系统性偏差。补偿方法是在 STFT 后对每个通道乘以一个相位校正因子，校正量通过互相关估计。

这套流程跑通之后，我把代码封装成了一个可配置的管道，STFT 参数、聚类参数、恢复算法都做成配置文件。换一个工况只需要改配置，不用动代码。现在处理一段 10 秒的双通道振动信号，从 STFT 到恢复出三个源波形，全程大约 15 秒，其中聚类占 2 秒，OMP 恢复占 12 秒。如果换成线性规划，同样数据要跑 8 分钟以上。所以除非精度实在不够，我一般默认用 OMP。

希望帮到你。

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