几十万点采样序列几十个谐波信噪比不到 10 dB这种数据放到我面前我第一反应是绝对不用标准奇异值分解去硬刚。因为在这样的规模下完整 SVD 一次就要吃掉大量内存和十几秒时间而且你还得反复调阈值。这篇内容想跟你分享的方案是把随机奇异值分解和软阈值收缩搭在一起用随机 SVD 负责把计算量压下来软阈值负责在不破坏弱谐波的前提下去掉噪声。整套算法用 Matlab 实现我会把随机 SVD、软阈值、参数选择、时序重建逐层拆开并附上可以直接运行的脚本和调参经验。这套思路最早是我在处理一段长时间振动监测数据时摸索出来的。当时数据里主要就是 50 Hz 工频及其整数倍谐波但背景噪声一直在变传统滤波器和固定阶数的陷波器完全不顶用最后折腾了一圈回到了矩阵低秩分解的去噪路子上。如果你也在做电力谐波分析、结构健康监测、声呐或语音去噪或者单纯想给长时间序列做预处理这篇文章应该能让你少走不少弯路。1. 谐波去噪为何要动“奇异值”的脑筋1.1 从滤波器到低秩矩阵Hankel 矩阵为什么能藏住谐波谐波污染和一般随机噪声不太一样。谐波带有明确的频率、幅值和相位往往淹没在宽频噪声里。如果你想要的只是一个“干净”的波形带通滤波器确实能凑合但只要你面对的是一段频率未知、谐波个数未知、信噪比随时变化的长时间序列滤波器的组合就需要频繁手动调整还会在频带上留下相位失真。矩阵分解类方法可以绕开这个问题它先把一段时域信号写成矩阵让矩阵自己暴露谐波的“低秩结构”再用奇异值把噪声和信号分开。具体做法很直观。给定一维信号 x(n)n 1, 2, ..., N取一个窗长 p构造 p 行、L N - p 1 列的 Hankel 矩阵第 1 行是 x(1) 到 x(L)第 2 行是 x(2) 到 x(L1)第 i 行第 j 列的元素是 x(ij-1)。这个矩阵的特别之处在于如果 x 本身由有限个复指数或者正弦分量组成那么 H 在理论上是一个低秩矩阵秩近似等于谐波个数的两倍左右。举个例子两个不同频率的正弦波叠加对应的 Hankel 矩阵主要奇异值就只有四五个其他全都可以忽略。加上噪声后随机扰动会让所有奇异值整体抬升靠后的奇异值不再接近零而是拖出一条平缓的尾巴。这条尾巴就是去噪时真正要处理的对象。1.2 标准 SVD 去噪的逻辑与死穴经典做法是先对 H 做完整奇异值分解画出奇异值序列找“拐点”保留拐点之前的分量把后面的奇异值全部置零再用反对角平均把矩阵还原成时域信号。这个方法在信噪比高的时候非常好用前几个奇异值远远大于噪声尾部拐点一眼就能看出来。但问题来了。当噪声较强或者谐波本身较弱时奇异值谱的“拐点”会变得非常模糊。你会陷入一个两难保护的前 r 个奇异值取大了噪声跟着进来去噪后频谱里依然是一坨底噪取小了较弱的谐波被连锅端信号直接缺了一块。而且完整 SVD 对 p×L 矩阵的计算复杂度大约在 O(p²L) 量级p 一旦取到几百N 到几十万这个开销就变得很可怕。我用一个 20 万点、窗长 300 的信号实测过完整 econ SVD 跑一次要等很久批处理根本没法做。1.3 随机 SVD 和软阈值各解决什么问题所以我想把两个东西组合起来随机 SVD 解决“算不动”的问题软阈值解决“切不干净”的问题。随机 SVD 不求把整个矩阵分解得明明白白它只用少量随机投影先抓住矩阵最主要的左奇异向量子空间再在这个低维子空间里做精确 SVD。复杂度从 O(p²L) 降到了 O(pLk)k 通常只有二三十这在 p 比较大时是数量级上的差异。软阈值则是把每个奇异值都减去一个固定的阈值 λ小于等于 λ 的直接清零大于 λ 的收缩到 σ_i - λ。它跟硬截断最大的区别是处理过程连续不会让某个分量在“保留”和“扔掉”之间剧烈跳变对接近阈值的微弱谐波更友好。下面我先把这两个工具的原理讲清楚再给出完整代码。2. 随机 SVD把大矩阵先压成小矩阵再谈分解2.1 随机投影一次乘法就找到主要子空间随机 SVD 的第一步是生成一个 n×k 的高斯随机矩阵 Ω然后计算 Y AΩ。这一步的直觉非常简单Ω 相当于一堆随机方向A 沿着这些方向作用之后得到的 Y 的每一列都落在 A 的值域里而且会偏向那些大奇异值对应的奇异向量方向。把 Y 做一次经济 QR 分解得到 Y QR取 Q。Q 的列就是 A 主要子空间的一组正交基。有了这个 QA 在 Q 张成的空间里的投影就是 B QᵀA这个 B 只有 k×L远小于原来的 p×L。对 B 做一次精确 SVD再把左奇异向量变换回原始空间就得到了近似奇异值分解。这个过程让我第一次跑通时很惊讶一个 400×5 万的矩阵只用了 32 列随机投影算出来的奇异值谱就和完整 SVD 几乎一样尤其是前面十几个奇异值误差基本在 1% 以内。这里有个经验k 要比你认为的真实秩大一些。如果信号实际的有效秩只有 4 到 6k 取 32 是合理的因为随机投影需要留出冗余维度去“接住”能量取太小会切掉真实信息。2.2 幂迭代矩阵-向量乘积的自我加强随机 SVD 还有一个可选的增强步骤幂迭代。核心公式就一行[Q, ~] qr(A * (A * Q), 0);这里 A 是转置矩阵这一串操作相当于用 AᵀA 反复作用在 Q 上。AᵀA 的特征值按奇异值平方排列小奇异值对应的方向衰减得非常快所以每做一次幂迭代Q 的子空间就会更贴近前 k 个主导奇异向量随机投影带来的偏差会大幅下降。在谐波去噪场景里我实测下来的经验是q 取 0 到 2 就足够了默认取 1。不要盲目加大 q因为每次迭代都要多花两次矩阵-向量乘法的计算收益却不是线性的。当信号本身信噪比很低、奇异值谱很平的时候幂迭代能够把那些较大的谐波奇异值“扶正”对提高弱谐波检出率有明显帮助。2.3 误差特性与适用边界从理论上说随机 SVD 的误差大致由被截断掉的奇异值 σ_{k1} 决定幂迭代 q 次之后误差的期望会受到 (σ_{k1}/σ_k)^(2q1) 比例的抑制。通俗地说如果你的前 k 个奇异值和后面的尾巴分离得越开随机 SVD 就越准甚至只需要很少的冗余 k 和很少的幂迭代次数。如果奇异值谱特别平比如纯白噪声矩阵随机 SVD 的近似质量就会变差但这种情况在谐波去噪里反而不是问题——因为我们要的本来就不是精细分解而是抓住主要分量、把尾部压掉。工程上我会在拿到结果后做一个快速 sanity check比较一下 H 被 Q 投影后的重建误差占原信号能量的比例如果这个比例超过 0.05就说明 k 太小或者 q 不够再往上加。3. 软阈值收缩为什么“软”比“硬”更稳3.1 三种奇异值处理方式对比拿到随机 SVD 的奇异值 σ₁ ≥ σ₂ ≥ ... ≥ σ_k 之后接下来要做的是在奇异值域里做“滤波”。常见的做法有三种差别非常大处理方式公式特点硬截断i ≤ r 保留否则清零简单但 r 不好定硬阈值σ_i λ 保留否则清零会在阈值附近产生突变软阈值max(σ_i - λ, 0)连续收缩稳健保留弱分量硬截断的问题我在前文已经说过了拐点难找硬阈值在 λ 附近会出现“要么全部保留、要么全部扔掉”的非连续性这个小缺点在批量处理时会被放大——你同一套参数跑十段信号某几段可能突然多出一两个虚假的谱峰。软阈值不会这样它的映射是连续的无论奇异值怎么抖动结果都会平滑过渡。3.2 噪声水平估计与阈值 λ 的公式软阈值里最关键的参数是 λ。我常用的自动估计方式分两步第一步估计噪声奇异值的“底”。把随机 SVD 得到的所有奇异值取中位数 median(s)除以 sqrt(max(p, L))得到一个噪声标准差 σ_est。这个公式的根据是随机矩阵理论一个纯噪声矩阵的奇异值谱大致散布在 σ√m 附近m 是矩阵的长边所以拿中位数反推噪声水平是有效的。谐波只占很少的几个奇异值不会严重污染中位数。第二步用极值理论里的通用阈值形式λ c · σ_est · sqrt(2 · log(N))这里的 N 是信号长度c 是强度系数我推荐从 c 3 起步。c 太小去噪后频谱里还能看到明显的连续噪声底c 太大谐波幅度会被压得过狠。经过两三轮调节你会找到一个比较宽的稳定区间在这个区间内输出 SNR 变化不大这就是参数鲁棒性好的表现。3.3 软阈值对谐波分量幅值的影响与补偿软阈值有个特性需要注意它会把所有大于 λ 的奇异值也减去 λ而不是保持原样。这意味着保留下来的谐波分量的幅度也会被轻微削弱在频谱上表现为谐波峰值低了一点。如果后续只做相对分析比如看哪些频率成分突出这个削弱基本无伤大雅。但如果要做幅度定量分析我建议在重建后做一个简单的能量恢复计算阈值前后奇异值总和的比例即 sum(s) / sum(max(s-λ, 0))把这个比例乘到重建信号上。这样能大致把总能量拉回原位代价是噪声能量也会被放回一点点所以只适合阈值用得比较保守的场合。不想引入额外噪声时也可以不做补偿直接用阈值后的结果。4. Matlab 实现去噪流水线完整拆解4.1 总体流程与函数划分整个谐波去噪流程可以拆成下面几步对含噪信号 x 做去均值处理去掉直流分量根据窗长 p 构造 Hankel 矩阵 H用随机 SVD 对 H 做低秩近似分解得到 U、奇异值 S、V在奇异值上施加软阈值得到收缩后的奇异值 s_th用 U、s_th、V 重建去噪后的 Hankel 矩阵 Hc沿着反对角线平均把 Hc 还原成一维信号 x_clean返回去噪结果和中间信息供后续调参。第 6 步的“反对角平均”是 SVD 类去噪里非常关键的一步。因为 H 矩阵里有很多冗余元素同一个信号点在矩阵里出现了多次直接取矩阵逆变换回一维会得到不一致的值。正确的做法是每个位置的能量按出现次数取平均这样既能恢复信号又能自然抑制部分噪声。4.2 随机 SVD 核心函数先放随机 SVD 的实现这个函数可以直接复用不仅限于谐波去噪其他低秩矩阵近似都能用function [U, S, V] rsvd(A, k, q) % RSVD 随机化奇异值分解 % A: m x n 矩阵 % k: 保留子空间维度 % q: 幂迭代次数0~2 即可 [~, n] size(A); Omega randn(n, k); Y A * Omega; [Q, ~] qr(Y, 0); for i 1:q [Q, ~] qr(A * (A * Q), 0); end B Q * A; [Ut, S, V] svd(B, econ); U Q * Ut; end函数本身的逻辑不复杂生成随机矩阵、投影、QR、幂迭代、对小矩阵做精确 SVD。有个细节值得强调幂迭代使用的是最终更新过的 Q不要贪图省事在循环外单独做一次 QR那样子空间精度反而会打折扣。我在初版代码里犯过这个错误结果奇异值谱尾部始终降不下去。4.3 谐波去噪主函数与对角平均重建下面是完整的去噪函数function [x_clean, info] rsvd_ht_denoise(x, p, k, q, lambda) % RSVD 软阈值谐波去噪 % x: 含噪信号列向量 % p: 窗长嵌入维度 % k: 随机SVD子空间维度 % q: 幂迭代次数 % lambda: 软阈值不传或传 NaN 时自动估计 x x(:); N length(x); L N - p 1; % 构造 Hankel 矩阵 H hankel(x(1:p), x(p:end)); % 随机 SVD [U, S, V] rsvd(H, k, q); s diag(S); % 阈值自动估计中位数奇异值 - 噪声水平 - 软阈值 if nargin 5 || isnan(lambda) sigma_est median(s) / sqrt(max(p, L)); lambda 3 * sigma_est * sqrt(2 * log(N)); end % 软阈值收缩 s_th max(s - lambda, 0); % 重建去噪后的轨迹矩阵 Hc U * diag(s_th) * V; % 反对角平均恢复一维信号 [I, J] ndgrid(1:p, 1:L); idx I(:) J(:) - 1; sumv accumarray(idx, Hc(:), [N, 1]); cnt accumarray(idx, ones(numel(Hc), 1), [N, 1]); x_clean sumv ./ cnt; % 返回一些诊断信息 info.lambda lambda; info.singular_values s; info.kept_rank sum(s_th 0); info.reconstruction_error norm(x_clean - x) / norm(x); end几点注释。Hankel 矩阵在 Matlab 里可以直接用内置函数非常方便但如果数据特别长你需要记住它生成的是一个 p × L 的稠密矩阵内存开销是 p·L 个 double也就是 8 字节每个元素。我一般把 p 控制在几百以内这样即便 N 有几十万这个矩阵也还能装下。对角线平均用 accumarray 实现比用 for 循环快很多因为 Hc 的元素位置映射到原始信号下标是确定的累加和计数一下就出来。这也是我在实际项目里踩了两次坑之后才改的写法一开始用双层循环N10 万时跑得非常难受。4.4 一个可直接运行的演示脚本下面这段脚本可以直接跑用来验证算法的实际效果clear; clc; rng(42); Fs 1000; % 采样率 1 kHz N 20000; % 2 万点信号 t (0:N-1) / Fs; % 两个谐波50 Hz 和 120 Hz x_true sin(2*pi*50*t) 0.6*sin(2*pi*120*t); % 叠加高斯白噪声 x_noisy x_true 0.6*randn(N, 1); % 去噪参数窗长 200随机子空间 30幂迭代 1 次阈值自动 p 200; k 30; q 1; [x_clean, info] rsvd_ht_denoise(x_noisy, p, k, q, NaN); % 指标 snr_before 10*log10(sum(x_true.^2) / sum((x_true - x_noisy).^2)); snr_after 10*log10(sum(x_true.^2) / sum((x_true - x_clean).^2)); fprintf(输入 SNR %.2f dB - 输出 SNR %.2f dB\n, snr_before, snr_after); fprintf(阈值 lambda %.4f, 保留奇异值个数 %d\n, info.lambda, info.kept_rank); % 频谱对比 Nfft 2^nextpow2(N); f (0:Nfft-1) / Nfft * Fs; S_before abs(fft(x_noisy, Nfft)); S_after abs(fft(x_clean, Nfft)); plot(f, 20*log10(S_before), color, [0.6 0.6 0.6]); hold on; plot(f, 20*log10(S_after), r, LineWidth, 1.2); xlim([0 200]); xlabel(Hz); ylabel(dB); legend({去噪前, 去噪后});在我的实测里这个场景的输出 SNR 通常能从大约 7 dB 提升到 25 dB 以上两个谐波峰都能完整保留底噪则是肉眼可见地塌了下去。值得注意的是输出 SNR 和输入信号本身形状、噪声强度有关不要指望所有数据都拿到同一个数值重点观察的是底噪是否被压平。4.5 不同阈值系数下的效果对比我建议在调参时专门跑一组阈值扫描看看 c 对结果的影响。以下是同一段信号在 c 分别取 1、2、3、4 时的表现趋势阈值系数 c底噪压制弱谐波保留输出 SNR1基本无效底噪仍明显保留完整提高有限2有一定压制保留较好明显提升3压制充分保留较好较高且稳定4非常干净可能被削弱可能下降这个表说明一件事阈值不是越大越好存在一个“最佳平台区”。我在实际项目里的习惯是先取 c3 跑一遍然后看奇异值谱和输出频谱小幅调整到 2.5 到 4 之间。这个区间通常比较宽不必追求某个神秘的最优值。5. 参数调优与工程落地5.1 窗长、子空间维度、幂迭代的推荐区间继续说参数。窗长 p 是最影响去噪质量的参数。p 太小Hankel 矩阵的行数不够低秩结构展现不出来结果接近普通滤波p 太大矩阵又宽又高计算成本上去了边界效应也更明显。我常用的经验是取信号长度的 1% 到 5%并且至少包含几个最长周期的波形。以 Fs1000 Hz、关注 50 Hz 谐波为例最长周期是 20 个采样点p 取 100 到 300 都比较合理。上面的演示脚本里 p200相当于 10 个最长周期效果就很稳定。子空间维度 k 的选取核心原则是“比真实秩大出一截”。两个正弦波的真实秩大约 4 到 5k 取 30 到 50 足够因为随机投影需要冗余维度去吸收噪声干扰。k 取得太小等价于强制忽略真实分量k 取得太大随机 SVD 退化成接近完整 SVD性能优势就没了。幂迭代 q 我默认取 1。只有当你的数据信噪比极低比如输入 SNR 不到 0 dB才考虑取 2。q 超过 2 以后计算代价上升明显精度收益却非常有限性价比很低。5.2 大规模数据集的内存与时间考量下面这张表说明不同规模下完整 SVD 与随机 SVD 的复杂度趋势注意看数量级而不是绝对耗时矩阵规模p × L完整 SVD 复杂度RSVD 复杂度200 × 20,000O(8e8)O(3.6e8)200 × 200,000O(8e9)O(3.6e9)400 × 500,000O(8e10)O(1.8e10)可以看到p 越大完整 SVD 的 O(p²L) 增长越凶猛随机 SVD 只有 O(pLk)k 固定后基本是线性增长。更关键的是内存完整 SVD 在解算过程中会生成大量临时矩阵而随机 SVD 只会在 Q 和 B 上多花一点空间内存占用友好得多。这就是大数据集上“高效”二字的真正含义。如果信号长度再往上走连 Hankel 矩阵本身的内存都吃紧我建议不再显式构造 H而是写一个实现 A 乘以向量的函数用索引运算直接完成投影。Matlab 里可以定义普通函数让 rsvd 接收函数句柄这样矩阵不落地内存占用能再降一个数量级。这个改造适合 N 超过百万点的场景普通几十万点的数据直接构造矩阵就够了。5.3 常见坑边界效应、随机种子、阈值不当最后集中说几个我在实际测试里踩过的坑。第一边界效应。Hankel 矩阵两端的元素比中间出现次数少反对角平均时两端统计量不足重建结果在首尾容易出现轻微抖动。处理办法有两个一是把去噪后的输出前、后各裁掉 p 个点只保留中间稳定区二是在信号两端做一小段镜像延拓去噪后再截掉。这两个方法我都试过简单粗暴地裁剪对频谱分析影响不大因为频谱分析本来就会加窗。第二随机种子。随机 SVD 依赖 randn 生成的投影矩阵如果不固定随机种子每次运行结果会有轻微差别。这不是 bug是随机算法的特性。在批量处理时我会在脚本开头固定 rng(42)保证结果可复现如果追求更稳定还可以生成两个不同的投影矩阵分别去噪再取平均方差会进一步下降。第三阈值不当导致谐波被削平。这个坑最容易出现在“自动估计”模式下。自动阈值里的中位数奇异值和信号长度 N 都会影响 λ如果 N 特别大sqrt(2·log(N)) 会偏大阈值可能压得太狠。遇到这种情况我建议不再依赖全自动而是先打印出奇异值序列肉眼找一个明显拐点再手动固定 lambda。反正调参也是要跑几轮的把奇异值谱画出来看比盲目调数字直观得多。最后说一个实操体会不要指望一次把参数调到最优。我一般先把 p 固定在合理范围、q 设为 1然后用一个小 k比如 16快速跑一遍打印奇异值谱看看前几个分量是否明显凸出。如果谱线尾部抬得很高说明噪声太强应该增大阈值系数而不是疯狂增加 k如果主分量之间粘连不清才考虑增大 k 或者 q 到 2。经过两三轮这个组合就能稳定下来之后再批量跑上百段数据基本不需要再动参数。这种“先看图、再调参、后批量”的工作方式配合随机 SVD 和软阈值这套组合才真正解决了大数据量下的谐波去噪问题。 
长河君 · 2026/10/14 22:00:55 · 阅读 1.2k 随机SVD+软阈值实现高噪声谐波去噪：Matlab完整方案
几十万点采样序列几十个谐波信噪比不到 10 dB这种数据放到我面前我第一反应是绝对不用标准奇异值分解去硬刚。因为在这样的规模下完整 SVD 一次就要吃掉大量内存和十几秒时间而且你还得反复调阈值。这篇内容想跟你分享的方案是把随机奇异值分解和软阈值收缩搭在一起用随机 SVD 负责把计算量压下来软阈值负责在不破坏弱谐波的前提下去掉噪声。整套算法用 Matlab 实现我会把随机 SVD、软阈值、参数选择、时序重建逐层拆开并附上可以直接运行的脚本和调参经验。这套思路最早是我在处理一段长时间振动监测数据时摸索出来的。当时数据里主要就是 50 Hz 工频及其整数倍谐波但背景噪声一直在变传统滤波器和固定阶数的陷波器完全不顶用最后折腾了一圈回到了矩阵低秩分解的去噪路子上。如果你也在做电力谐波分析、结构健康监测、声呐或语音去噪或者单纯想给长时间序列做预处理这篇文章应该能让你少走不少弯路。1. 谐波去噪为何要动“奇异值”的脑筋1.1 从滤波器到低秩矩阵Hankel 矩阵为什么能藏住谐波谐波污染和一般随机噪声不太一样。谐波带有明确的频率、幅值和相位往往淹没在宽频噪声里。如果你想要的只是一个“干净”的波形带通滤波器确实能凑合但只要你面对的是一段频率未知、谐波个数未知、信噪比随时变化的长时间序列滤波器的组合就需要频繁手动调整还会在频带上留下相位失真。矩阵分解类方法可以绕开这个问题它先把一段时域信号写成矩阵让矩阵自己暴露谐波的“低秩结构”再用奇异值把噪声和信号分开。具体做法很直观。给定一维信号 x(n)n 1, 2, ..., N取一个窗长 p构造 p 行、L N - p 1 列的 Hankel 矩阵第 1 行是 x(1) 到 x(L)第 2 行是 x(2) 到 x(L1)第 i 行第 j 列的元素是 x(ij-1)。这个矩阵的特别之处在于如果 x 本身由有限个复指数或者正弦分量组成那么 H 在理论上是一个低秩矩阵秩近似等于谐波个数的两倍左右。举个例子两个不同频率的正弦波叠加对应的 Hankel 矩阵主要奇异值就只有四五个其他全都可以忽略。加上噪声后随机扰动会让所有奇异值整体抬升靠后的奇异值不再接近零而是拖出一条平缓的尾巴。这条尾巴就是去噪时真正要处理的对象。1.2 标准 SVD 去噪的逻辑与死穴经典做法是先对 H 做完整奇异值分解画出奇异值序列找“拐点”保留拐点之前的分量把后面的奇异值全部置零再用反对角平均把矩阵还原成时域信号。这个方法在信噪比高的时候非常好用前几个奇异值远远大于噪声尾部拐点一眼就能看出来。但问题来了。当噪声较强或者谐波本身较弱时奇异值谱的“拐点”会变得非常模糊。你会陷入一个两难保护的前 r 个奇异值取大了噪声跟着进来去噪后频谱里依然是一坨底噪取小了较弱的谐波被连锅端信号直接缺了一块。而且完整 SVD 对 p×L 矩阵的计算复杂度大约在 O(p²L) 量级p 一旦取到几百N 到几十万这个开销就变得很可怕。我用一个 20 万点、窗长 300 的信号实测过完整 econ SVD 跑一次要等很久批处理根本没法做。1.3 随机 SVD 和软阈值各解决什么问题所以我想把两个东西组合起来随机 SVD 解决“算不动”的问题软阈值解决“切不干净”的问题。随机 SVD 不求把整个矩阵分解得明明白白它只用少量随机投影先抓住矩阵最主要的左奇异向量子空间再在这个低维子空间里做精确 SVD。复杂度从 O(p²L) 降到了 O(pLk)k 通常只有二三十这在 p 比较大时是数量级上的差异。软阈值则是把每个奇异值都减去一个固定的阈值 λ小于等于 λ 的直接清零大于 λ 的收缩到 σ_i - λ。它跟硬截断最大的区别是处理过程连续不会让某个分量在“保留”和“扔掉”之间剧烈跳变对接近阈值的微弱谐波更友好。下面我先把这两个工具的原理讲清楚再给出完整代码。2. 随机 SVD把大矩阵先压成小矩阵再谈分解2.1 随机投影一次乘法就找到主要子空间随机 SVD 的第一步是生成一个 n×k 的高斯随机矩阵 Ω然后计算 Y AΩ。这一步的直觉非常简单Ω 相当于一堆随机方向A 沿着这些方向作用之后得到的 Y 的每一列都落在 A 的值域里而且会偏向那些大奇异值对应的奇异向量方向。把 Y 做一次经济 QR 分解得到 Y QR取 Q。Q 的列就是 A 主要子空间的一组正交基。有了这个 QA 在 Q 张成的空间里的投影就是 B QᵀA这个 B 只有 k×L远小于原来的 p×L。对 B 做一次精确 SVD再把左奇异向量变换回原始空间就得到了近似奇异值分解。这个过程让我第一次跑通时很惊讶一个 400×5 万的矩阵只用了 32 列随机投影算出来的奇异值谱就和完整 SVD 几乎一样尤其是前面十几个奇异值误差基本在 1% 以内。这里有个经验k 要比你认为的真实秩大一些。如果信号实际的有效秩只有 4 到 6k 取 32 是合理的因为随机投影需要留出冗余维度去“接住”能量取太小会切掉真实信息。2.2 幂迭代矩阵-向量乘积的自我加强随机 SVD 还有一个可选的增强步骤幂迭代。核心公式就一行[Q, ~] qr(A * (A * Q), 0);这里 A 是转置矩阵这一串操作相当于用 AᵀA 反复作用在 Q 上。AᵀA 的特征值按奇异值平方排列小奇异值对应的方向衰减得非常快所以每做一次幂迭代Q 的子空间就会更贴近前 k 个主导奇异向量随机投影带来的偏差会大幅下降。在谐波去噪场景里我实测下来的经验是q 取 0 到 2 就足够了默认取 1。不要盲目加大 q因为每次迭代都要多花两次矩阵-向量乘法的计算收益却不是线性的。当信号本身信噪比很低、奇异值谱很平的时候幂迭代能够把那些较大的谐波奇异值“扶正”对提高弱谐波检出率有明显帮助。2.3 误差特性与适用边界从理论上说随机 SVD 的误差大致由被截断掉的奇异值 σ_{k1} 决定幂迭代 q 次之后误差的期望会受到 (σ_{k1}/σ_k)^(2q1) 比例的抑制。通俗地说如果你的前 k 个奇异值和后面的尾巴分离得越开随机 SVD 就越准甚至只需要很少的冗余 k 和很少的幂迭代次数。如果奇异值谱特别平比如纯白噪声矩阵随机 SVD 的近似质量就会变差但这种情况在谐波去噪里反而不是问题——因为我们要的本来就不是精细分解而是抓住主要分量、把尾部压掉。工程上我会在拿到结果后做一个快速 sanity check比较一下 H 被 Q 投影后的重建误差占原信号能量的比例如果这个比例超过 0.05就说明 k 太小或者 q 不够再往上加。3. 软阈值收缩为什么“软”比“硬”更稳3.1 三种奇异值处理方式对比拿到随机 SVD 的奇异值 σ₁ ≥ σ₂ ≥ ... ≥ σ_k 之后接下来要做的是在奇异值域里做“滤波”。常见的做法有三种差别非常大处理方式公式特点硬截断i ≤ r 保留否则清零简单但 r 不好定硬阈值σ_i λ 保留否则清零会在阈值附近产生突变软阈值max(σ_i - λ, 0)连续收缩稳健保留弱分量硬截断的问题我在前文已经说过了拐点难找硬阈值在 λ 附近会出现“要么全部保留、要么全部扔掉”的非连续性这个小缺点在批量处理时会被放大——你同一套参数跑十段信号某几段可能突然多出一两个虚假的谱峰。软阈值不会这样它的映射是连续的无论奇异值怎么抖动结果都会平滑过渡。3.2 噪声水平估计与阈值 λ 的公式软阈值里最关键的参数是 λ。我常用的自动估计方式分两步第一步估计噪声奇异值的“底”。把随机 SVD 得到的所有奇异值取中位数 median(s)除以 sqrt(max(p, L))得到一个噪声标准差 σ_est。这个公式的根据是随机矩阵理论一个纯噪声矩阵的奇异值谱大致散布在 σ√m 附近m 是矩阵的长边所以拿中位数反推噪声水平是有效的。谐波只占很少的几个奇异值不会严重污染中位数。第二步用极值理论里的通用阈值形式λ c · σ_est · sqrt(2 · log(N))这里的 N 是信号长度c 是强度系数我推荐从 c 3 起步。c 太小去噪后频谱里还能看到明显的连续噪声底c 太大谐波幅度会被压得过狠。经过两三轮调节你会找到一个比较宽的稳定区间在这个区间内输出 SNR 变化不大这就是参数鲁棒性好的表现。3.3 软阈值对谐波分量幅值的影响与补偿软阈值有个特性需要注意它会把所有大于 λ 的奇异值也减去 λ而不是保持原样。这意味着保留下来的谐波分量的幅度也会被轻微削弱在频谱上表现为谐波峰值低了一点。如果后续只做相对分析比如看哪些频率成分突出这个削弱基本无伤大雅。但如果要做幅度定量分析我建议在重建后做一个简单的能量恢复计算阈值前后奇异值总和的比例即 sum(s) / sum(max(s-λ, 0))把这个比例乘到重建信号上。这样能大致把总能量拉回原位代价是噪声能量也会被放回一点点所以只适合阈值用得比较保守的场合。不想引入额外噪声时也可以不做补偿直接用阈值后的结果。4. Matlab 实现去噪流水线完整拆解4.1 总体流程与函数划分整个谐波去噪流程可以拆成下面几步对含噪信号 x 做去均值处理去掉直流分量根据窗长 p 构造 Hankel 矩阵 H用随机 SVD 对 H 做低秩近似分解得到 U、奇异值 S、V在奇异值上施加软阈值得到收缩后的奇异值 s_th用 U、s_th、V 重建去噪后的 Hankel 矩阵 Hc沿着反对角线平均把 Hc 还原成一维信号 x_clean返回去噪结果和中间信息供后续调参。第 6 步的“反对角平均”是 SVD 类去噪里非常关键的一步。因为 H 矩阵里有很多冗余元素同一个信号点在矩阵里出现了多次直接取矩阵逆变换回一维会得到不一致的值。正确的做法是每个位置的能量按出现次数取平均这样既能恢复信号又能自然抑制部分噪声。4.2 随机 SVD 核心函数先放随机 SVD 的实现这个函数可以直接复用不仅限于谐波去噪其他低秩矩阵近似都能用function [U, S, V] rsvd(A, k, q) % RSVD 随机化奇异值分解 % A: m x n 矩阵 % k: 保留子空间维度 % q: 幂迭代次数0~2 即可 [~, n] size(A); Omega randn(n, k); Y A * Omega; [Q, ~] qr(Y, 0); for i 1:q [Q, ~] qr(A * (A * Q), 0); end B Q * A; [Ut, S, V] svd(B, econ); U Q * Ut; end函数本身的逻辑不复杂生成随机矩阵、投影、QR、幂迭代、对小矩阵做精确 SVD。有个细节值得强调幂迭代使用的是最终更新过的 Q不要贪图省事在循环外单独做一次 QR那样子空间精度反而会打折扣。我在初版代码里犯过这个错误结果奇异值谱尾部始终降不下去。4.3 谐波去噪主函数与对角平均重建下面是完整的去噪函数function [x_clean, info] rsvd_ht_denoise(x, p, k, q, lambda) % RSVD 软阈值谐波去噪 % x: 含噪信号列向量 % p: 窗长嵌入维度 % k: 随机SVD子空间维度 % q: 幂迭代次数 % lambda: 软阈值不传或传 NaN 时自动估计 x x(:); N length(x); L N - p 1; % 构造 Hankel 矩阵 H hankel(x(1:p), x(p:end)); % 随机 SVD [U, S, V] rsvd(H, k, q); s diag(S); % 阈值自动估计中位数奇异值 - 噪声水平 - 软阈值 if nargin 5 || isnan(lambda) sigma_est median(s) / sqrt(max(p, L)); lambda 3 * sigma_est * sqrt(2 * log(N)); end % 软阈值收缩 s_th max(s - lambda, 0); % 重建去噪后的轨迹矩阵 Hc U * diag(s_th) * V; % 反对角平均恢复一维信号 [I, J] ndgrid(1:p, 1:L); idx I(:) J(:) - 1; sumv accumarray(idx, Hc(:), [N, 1]); cnt accumarray(idx, ones(numel(Hc), 1), [N, 1]); x_clean sumv ./ cnt; % 返回一些诊断信息 info.lambda lambda; info.singular_values s; info.kept_rank sum(s_th 0); info.reconstruction_error norm(x_clean - x) / norm(x); end几点注释。Hankel 矩阵在 Matlab 里可以直接用内置函数非常方便但如果数据特别长你需要记住它生成的是一个 p × L 的稠密矩阵内存开销是 p·L 个 double也就是 8 字节每个元素。我一般把 p 控制在几百以内这样即便 N 有几十万这个矩阵也还能装下。对角线平均用 accumarray 实现比用 for 循环快很多因为 Hc 的元素位置映射到原始信号下标是确定的累加和计数一下就出来。这也是我在实际项目里踩了两次坑之后才改的写法一开始用双层循环N10 万时跑得非常难受。4.4 一个可直接运行的演示脚本下面这段脚本可以直接跑用来验证算法的实际效果clear; clc; rng(42); Fs 1000; % 采样率 1 kHz N 20000; % 2 万点信号 t (0:N-1) / Fs; % 两个谐波50 Hz 和 120 Hz x_true sin(2*pi*50*t) 0.6*sin(2*pi*120*t); % 叠加高斯白噪声 x_noisy x_true 0.6*randn(N, 1); % 去噪参数窗长 200随机子空间 30幂迭代 1 次阈值自动 p 200; k 30; q 1; [x_clean, info] rsvd_ht_denoise(x_noisy, p, k, q, NaN); % 指标 snr_before 10*log10(sum(x_true.^2) / sum((x_true - x_noisy).^2)); snr_after 10*log10(sum(x_true.^2) / sum((x_true - x_clean).^2)); fprintf(输入 SNR %.2f dB - 输出 SNR %.2f dB\n, snr_before, snr_after); fprintf(阈值 lambda %.4f, 保留奇异值个数 %d\n, info.lambda, info.kept_rank); % 频谱对比 Nfft 2^nextpow2(N); f (0:Nfft-1) / Nfft * Fs; S_before abs(fft(x_noisy, Nfft)); S_after abs(fft(x_clean, Nfft)); plot(f, 20*log10(S_before), color, [0.6 0.6 0.6]); hold on; plot(f, 20*log10(S_after), r, LineWidth, 1.2); xlim([0 200]); xlabel(Hz); ylabel(dB); legend({去噪前, 去噪后});在我的实测里这个场景的输出 SNR 通常能从大约 7 dB 提升到 25 dB 以上两个谐波峰都能完整保留底噪则是肉眼可见地塌了下去。值得注意的是输出 SNR 和输入信号本身形状、噪声强度有关不要指望所有数据都拿到同一个数值重点观察的是底噪是否被压平。4.5 不同阈值系数下的效果对比我建议在调参时专门跑一组阈值扫描看看 c 对结果的影响。以下是同一段信号在 c 分别取 1、2、3、4 时的表现趋势阈值系数 c底噪压制弱谐波保留输出 SNR1基本无效底噪仍明显保留完整提高有限2有一定压制保留较好明显提升3压制充分保留较好较高且稳定4非常干净可能被削弱可能下降这个表说明一件事阈值不是越大越好存在一个“最佳平台区”。我在实际项目里的习惯是先取 c3 跑一遍然后看奇异值谱和输出频谱小幅调整到 2.5 到 4 之间。这个区间通常比较宽不必追求某个神秘的最优值。5. 参数调优与工程落地5.1 窗长、子空间维度、幂迭代的推荐区间继续说参数。窗长 p 是最影响去噪质量的参数。p 太小Hankel 矩阵的行数不够低秩结构展现不出来结果接近普通滤波p 太大矩阵又宽又高计算成本上去了边界效应也更明显。我常用的经验是取信号长度的 1% 到 5%并且至少包含几个最长周期的波形。以 Fs1000 Hz、关注 50 Hz 谐波为例最长周期是 20 个采样点p 取 100 到 300 都比较合理。上面的演示脚本里 p200相当于 10 个最长周期效果就很稳定。子空间维度 k 的选取核心原则是“比真实秩大出一截”。两个正弦波的真实秩大约 4 到 5k 取 30 到 50 足够因为随机投影需要冗余维度去吸收噪声干扰。k 取得太小等价于强制忽略真实分量k 取得太大随机 SVD 退化成接近完整 SVD性能优势就没了。幂迭代 q 我默认取 1。只有当你的数据信噪比极低比如输入 SNR 不到 0 dB才考虑取 2。q 超过 2 以后计算代价上升明显精度收益却非常有限性价比很低。5.2 大规模数据集的内存与时间考量下面这张表说明不同规模下完整 SVD 与随机 SVD 的复杂度趋势注意看数量级而不是绝对耗时矩阵规模p × L完整 SVD 复杂度RSVD 复杂度200 × 20,000O(8e8)O(3.6e8)200 × 200,000O(8e9)O(3.6e9)400 × 500,000O(8e10)O(1.8e10)可以看到p 越大完整 SVD 的 O(p²L) 增长越凶猛随机 SVD 只有 O(pLk)k 固定后基本是线性增长。更关键的是内存完整 SVD 在解算过程中会生成大量临时矩阵而随机 SVD 只会在 Q 和 B 上多花一点空间内存占用友好得多。这就是大数据集上“高效”二字的真正含义。如果信号长度再往上走连 Hankel 矩阵本身的内存都吃紧我建议不再显式构造 H而是写一个实现 A 乘以向量的函数用索引运算直接完成投影。Matlab 里可以定义普通函数让 rsvd 接收函数句柄这样矩阵不落地内存占用能再降一个数量级。这个改造适合 N 超过百万点的场景普通几十万点的数据直接构造矩阵就够了。5.3 常见坑边界效应、随机种子、阈值不当最后集中说几个我在实际测试里踩过的坑。第一边界效应。Hankel 矩阵两端的元素比中间出现次数少反对角平均时两端统计量不足重建结果在首尾容易出现轻微抖动。处理办法有两个一是把去噪后的输出前、后各裁掉 p 个点只保留中间稳定区二是在信号两端做一小段镜像延拓去噪后再截掉。这两个方法我都试过简单粗暴地裁剪对频谱分析影响不大因为频谱分析本来就会加窗。第二随机种子。随机 SVD 依赖 randn 生成的投影矩阵如果不固定随机种子每次运行结果会有轻微差别。这不是 bug是随机算法的特性。在批量处理时我会在脚本开头固定 rng(42)保证结果可复现如果追求更稳定还可以生成两个不同的投影矩阵分别去噪再取平均方差会进一步下降。第三阈值不当导致谐波被削平。这个坑最容易出现在“自动估计”模式下。自动阈值里的中位数奇异值和信号长度 N 都会影响 λ如果 N 特别大sqrt(2·log(N)) 会偏大阈值可能压得太狠。遇到这种情况我建议不再依赖全自动而是先打印出奇异值序列肉眼找一个明显拐点再手动固定 lambda。反正调参也是要跑几轮的把奇异值谱画出来看比盲目调数字直观得多。最后说一个实操体会不要指望一次把参数调到最优。我一般先把 p 固定在合理范围、q 设为 1然后用一个小 k比如 16快速跑一遍打印奇异值谱看看前几个分量是否明显凸出。如果谱线尾部抬得很高说明噪声太强应该增大阈值系数而不是疯狂增加 k如果主分量之间粘连不清才考虑增大 k 或者 q 到 2。经过两三轮这个组合就能稳定下来之后再批量跑上百段数据基本不需要再动参数。这种“先看图、再调参、后批量”的工作方式配合随机 SVD 和软阈值这套组合才真正解决了大数据量下的谐波去噪问题。