做信号处理的谁没被一两段“又长又脏”的数据折磨过呢。我最近处理一批振动台测试的实测数据几十万个采样点基波、二倍频、三倍频清清楚楚可叠加的随机噪声也不含糊。想在Matlab里用经典SVD做谐波去噪结果svd()函数一跑内存先报警换成小波阈值阈值调来调去低频段该留的谐波差点被削掉。后来我把随机奇异值分解Randomized SVD和软阈值Soft Thresholding搭在一起写了一套谐波去噪流程——在比较大的数据集上计算快、内存省去噪效果也比固定秩截断稳定得多。这篇文章把思路、原理、代码和调试经验一次说清想直接抄代码的可以从第3节开始看。1. 先把问题拆清楚大数据谐波去噪到底难在哪1.1 谐波去噪的传统套路和它的天花板谐波去噪在工程里太常见了。电网信号里有50Hz基波和100Hz、150Hz倍频机械振动里有转频及其高次倍频声学测试里也有大量周期成分。去噪不是简单拿个低通滤波器一滤了事而是要在一堆噪声里把各次谐波的幅值、频率尽量原样保留下来。噪声来源又杂传感器热噪声、电磁干扰、随机环境振动很多情况下只能当高斯白噪声处理。传统做法大致有三条路。第一条是频域带通或梳状滤波频率已知时效果还行但谐波频率一旦漂移或者存在间谐波梳状滤波会把频谱“梳”出一道道坑谐波能量受损。第二条是小波阈值去噪思路是信号在小波域能量集中、噪声分布平均用阈值收缩小波系数。实际跑起来你会发现阈值选大选小非常敏感而且谐波密集时小波系数在多个尺度上都有能量硬阈值很容易把弱谐波连根拔掉。第三条就是经典SVD去噪把一维信号构造成Hankel矩阵做完整SVD把奇异值截断再重构。这条路的数学很美但天花板也很明显。经典SVD去噪的基本操作是这样的对信号x构造Hankel矩阵H矩阵元素满足H(i,j)x(ij-1)。对H做奇异值分解得到HUΣV^T。信号部分对应较大的奇异值噪声表现为一串缓慢衰减的小奇异值。去噪时把后面小奇异值直接置零再用U、V重构矩阵最后对角平均还原一维信号。问题出在“把后面置零”这个动作上你需要提前确定保留多少个奇异值也就是秩k。k小一截弱谐波没了k大一截噪声分量全漏进来。数据小的时候还能凭经验一遍遍试数据一旦大起来经验就不太好使了。更根本的问题是计算量。完整SVD对m×n矩阵的复杂度大概是O(mn^2)量级一个长度几百万点的信号构造出的Hankel矩阵动辄几十万乘几十万存下来都是几百GB。就算你用分块技巧勉强算时间也完全不可接受。这就是大数据谐波去噪的尴尬经典算法在教科书上完美放到真实数据上根本跑不动。1.2 随机SVD加软阈值解决的正是这两个痛点随机SVD和软阈值这两个词放在一起不是随便拼凑的它们分别解决了我上面说的两个核心痛点。随机SVD解决的是“算不动”。它的核心思想是用一个随机投影矩阵把原始大矩阵压缩到一个低维空间在低维空间里做标准SVD再把奇异向量映射回来。整个过程只需要对矩阵做几次矩阵乘法和一次小矩阵分解计算量从O(mn^2)直接降到O(mnl l^2(mn))这里的l是你要保留的分量个数通常只有几十。换句话说以前要对整个大矩阵精细分解现在只需要“抽查”一部分方向代价是微小的精度损失。软阈值解决的是“不知道留几个分量”。传统截断SVD相当于硬截断判断第k个奇异值后面全扔。但实际数据的噪声强度、谐波强度都在变化k根本没法提前猜准。软阈值对每个奇异值做一个收缩操作s_i max(s_i - τ, 0)。大于阈值的奇异值保留下来但稍微减小一点点小于阈值的奇异值直接变零。这个操作不需要你精确指定保留多少个分量只需要给一个合理阈值τ算法自己会“看情况”保留。这样配合随机SVD就是又快又稳的组合。打个比方。完整SVD像是把整个图书馆的书逐本翻一遍把最重要的几百本挑出来随机SVD是随机抽十几排书架通过它们迅速推断图书馆的主要分类然后重点整理这几排。软阈值则是挑书时不搞“一刀切”——它会根据书的新旧程度、借阅频率综合判断而不是只看分类号。两个工具各管各的环节合在一起效率就上来了。从我的实测经验看用这套组合处理长度十万级到百万级的信号运行时间可以从“内存爆掉”变成“几十秒搞定”输出信噪比和经典截断SVD基本持平在某些噪声分布不均的场景下甚至更稳。后面我会给出一组具体对比数据。2. 核心算法原理随机SVD和软阈值理解这几步就够用了2.1 构造轨迹矩阵把一维谐波信号变成二维矩阵构造Hankel矩阵这一步是整个方法的基石也是很多新手最容易忽视的地方。一维信号x长度为L选择一个嵌入窗长W构造的Hankel矩阵是W行、L-W1列元素就是原始信号H(i,j) x(ij-1)第一行是x(1)到x(L-W1)第二行是x(2)到x(L-W2)依此类推。这就是把一维时间序列“嵌入”成二维矩阵也叫延迟嵌入。为什么要这么做因为若干正弦分量组成的信号其Hankel矩阵是低秩的。理想情况下一个单一频率的正弦信号对应的Hankel矩阵秩为2对应正负频率两路h个独立谐波对应秩约2h。噪声的加入会把矩阵变成满秩但信号对应的主奇异值依然明显大于噪声奇异值这给了SVD去噪的数学基础。窗长W的选择直接影响奇异值谱的分辨率。W太小矩阵的秩估计不稳定信号和噪声的奇异值边界模糊W太大矩阵维数过高计算开销增大。我跑下来比较稳的口诀是W至少覆盖2到4个基波周期同时在可行范围内大一点。比如采样率1000Hz、基波50Hz基波周期是20个采样点W取200到400比较合适如果数据量允许取1000到2000也能得到更平滑的奇异值谱。在大数据场景下W不建议盲目取L/3甚至L/2那样矩阵实在太大了等会儿第3节会给出一个兼顾计算量和效果的取值思路。另外要注意构造Hankel矩阵用的数据越长尾部噪声对重构的影响越小但矩阵规模也越大。这其实是一个“分辨率”和“计算量”的trade-off。工程上可以先用一小段数据试出合适的W和秩范围再整体跑。2.2 随机SVD到底在做什么随机SVD的完整算法可以拆成五步。假设矩阵A是m×n我们想求前k个主奇异值。第一步是生成一个n×l的随机矩阵Ωlkpp称为过采样参数。Ω每个元素都是独立的标准正态随机数。过采样的意义是给投影留出余量避免因为随机性漏掉某些能量集中的方向。p取5到10就够取太大边际收益很小。第二步计算YAΩ这个Y就是矩阵A在随机方向上的投影。由于Ω是随机向量Y的列向量大概率落在A的主奇异方向上。这里隐含的数学结果是只要l大于等于A的有效秩Y就能以接近1的概率张成A的主奇异子空间。第三步是可选幂迭代。按Y A(A^T Y)重复q次。这一步的作用是压制小奇异值对应的分量让主方向更突出。q通常取1或2。迭代多了精度会更高但每次都涉及两次矩阵乘大数据下成本不低我一般只在噪声特别重或者矩阵条件数差的时候把q加到2到3。第四步对Y做QR分解YQRQ是m×l的正交矩阵。这样Q的列空间就近似等于A的主列空间。第五步计算BQ^T AB是l×n的小矩阵。对B做标准SVD得到BU_BΣV^T。因为A≈QB所以A的主奇异值近似等于Σ的主奇异值A的左奇异向量U≈Q U_B右奇异向量就是V。随机SVD的误差有明确概率保证。在Halko等关于随机化数值线性代数的经典分析里只要A的奇异值衰减得够快谐波信号正好满足这个特点随机SVD得到的低秩近似能以极高概率逼近最优低秩近似。这就是为什么对谐波去噪这类谱结构明显的问题随机SVD几乎不会损失精度。2.3 软阈值收缩比硬截断更聪明的去噪策略硬截断去噪是“一刀切”排序后的奇异值序列前k个保留后面的全置零。这个操作看着干脆实际上很脆弱。噪声强的时候前k个奇异值里可能混进噪声分量噪声弱的时候第k1个奇异值可能还是有效信号。k的选取只要差一个数重构结果的天差地别。软阈值处理的思路完全不一样。给定阈值τ后每个奇异值都执行s_i max(s_i - τ, 0)大于τ的奇异值保留下来但要减掉τ不大于τ的直接归零。这个操作的好处是连续可控信号成分奇异值大减掉一个τ几乎不影响噪声成分奇异值小减掉τ之后就趋近于零。整个过程不需要预先回答“到底保留几个分量”这个问题阈值τ代替了秩k而且对τ的敏感度远低于对k的敏感度。阈值τ怎么定这是很多人问得最多的地方。我在Matlab里习惯这样估计拿奇异值序列的后半段看作“噪声奇异值”用绝对中位差MAD估计噪声水平σ然后按Donoho通用阈值放大tail s_vals(round(end*0.5)1:end); sigma_tail median(abs(tail - median(tail))) / 0.6745; tau sigma_tail * sqrt(2*log(L));其中除以0.6745是因为对正态分布数据MAD约等于0.6745倍标准差乘sqrt(2logL)是通用阈值的标准形式。这个公式源于小波去噪严格说在Hankel域不是最严谨的但实际用起来很稳。有时候我也会手动观察奇异值谱找一个明显的“平台区”起点把τ设成平台区平均奇异值的2到3倍效果也差不多。关键在于软阈值把“选个数”变成了“选一个连续数”后者好调多了。3. Matlab实现全过程从仿真数据到可运行的代码3.1 生成一个带谐波和噪声的测试信号先用一个仿真例子把整套流程跑通。采样率1000Hz信号时长10秒总点数10000。基波50Hz幅度1二次谐波幅度0.4三次谐波幅度0.2。噪声用高斯白噪声目标输入信噪比5dB。clear; clc; rng(2025); fs 1000; % 采样率 1000 Hz T 10; % 信号时长 10 秒 L fs * T; % 总点数 10000 t (0:L-1)/fs; f0 50; % 基波频率 x_clean 1.0*sin(2*pi*f0*t) 0.4*sin(2*pi*2*f0*t) 0.2*sin(2*pi*3*f0*t); SNR_dB 5; % 目标输入信噪比 noise randn(1, L); noise noise / std(noise) * std(x_clean) / (10^(SNR_dB/20)); x x_clean noise; snr_in 10*log10(sum(x_clean.^2) / sum((x - x_clean).^2)); fprintf(输入信噪比: %.2f dB\n, snr_in);生成噪声时用std来代替rms这样不依赖额外的工具箱函数。严格地说正弦信号的rms等于幅值除以sqrt(2)但这里用std控制相对大小完全够用后面的计算结果也能对上。3.2 分步实现随机SVD去噪的完整代码先写随机SVD子函数。这里要注意矩阵A的维度m和n都可能比较大但在子函数内部只需要size(A)就能拿到不需要额外传参function [U, S, V] rsvd(A, k, p, q) % 随机SVD近似前k个奇异值/向量 % k: 目标奇异值个数p: 过采样数q: 幂迭代次数 [m, n] size(A); l k p; Omega randn(n, l); Y A * Omega; for i 1:q Y A * (A * Y); end [Q, ~] qr(Y, 0); B Q * A; [U_B, S, V] svd(B, econ); U Q * U_B; U U(:, 1:k); S S(1:k, 1:k); V V(:, 1:k); end然后是主流程。这里我特意把目标奇异值个数k设成20而不是严格按“谐波数×2”猜8。原因在于软阈值会自动收缩多余分量k大一点只是让随机SVD多算几个候选奇异值不会像硬截断那样因为k选错而崩掉。% 参数设置 W 2000; % 嵌入窗长 k 20; % 目标奇异值个数故意多留余量 p 5; % 过采样 q 1; % 幂迭代次数 % 构造Hankel矩阵 H hankel(x(1:W), x(W:end)); % 随机SVD [U, S, V] rsvd(H, k, p, q); s_vals diag(S); % 用奇异值尾部估计噪声水平并计算软阈值 tail s_vals(round(end*0.5)1:end); sigma_tail median(abs(tail - median(tail))) / 0.6745; tau sigma_tail * sqrt(2*log(L)); % 软阈值收缩 s_shrunk max(s_vals - tau, 0); % 重构低秩矩阵 H_denoised U * diag(s_shrunk) * V; % 对角平均恢复一维信号 x_denoised zeros(1, L); cnt zeros(1, L); for i 1:size(H_denoised, 1) seg H_denoised(i, :); inds i : i size(H_denoised, 2) - 1; x_denoised(inds) x_denoised(inds) seg; cnt(inds) cnt(inds) 1; end x_denoised x_denoised ./ cnt; % 评估去噪效果 snr_out 10*log10(sum(x_clean.^2) / sum((x_denoised - x_clean).^2)); fprintf(输入SNR: %.2f dB - 输出SNR: %.2f dB\n, snr_in, snr_out);这个流程我在Matlab R2021b以后版本上都跑过没有额外工具箱依赖。对角平均那段循环看着朴素其实比二维索引矩阵要省内存数据量大的时候不会因为重构矩阵就爆掉。如果只想看整个流程的主干可以把rsvd子函数、软阈值操作和主流程存成两个文件放在同一目录下直接运行。我那份完整脚本里还加了一段频谱对比用来快速确认去噪后各次谐波有没有被削平。3.3 参数怎么定窗口长度、目标秩、阈值系数先讲窗长W。我实际测试下来W和基波周期的比值比W的绝对值更重要。设每个基波周期的采样点数为Mfs/f0W至少取2M到4M。比如fs1000、f050时M20W取400就能看到清晰的奇异值谱断层但为了平滑估计阈值我常常取1000到2000。数据量大时W可以固定为一个几千的常数不必跟着总长度L无限增大因为软阈值对W的敏感度远低于对秩k的敏感度。再讲目标奇异值个数k。标题里提到的“健壮”很大程度体现在这一步不要纠结于精确估计谐波个数。我习惯把k设成“猜测谐波数×24到6”让随机SVD给出足够的候选奇异值最后交给软阈值去收缩。这样就算你把谐波数猜成了2倍结果也不会有本质变化。最后讲阈值τ。如果奇异值尾部样本太少比如k只取了8尾部只有三四个点MAD估计会非常不可靠。这也是我为啥建议k取20以上的原因。如果噪声很强导致尾部奇异值仍然很大τ会整体放大去噪会更激进。想调高保留信号比例可以把τ最终乘0.7到0.8想更干净把噪声压下去就乘1.3到1.5。我在4.2节给了这组敏感性数据你会发现这个系数的操作空间比硬截断的k大多了。4. 实测效果与参数对照4.1 与完整SVD去噪的效率和效果对比在普通台式机上8核CPU、32GB内存我用同一组谐波信号跑了一组对照实验。数据长度L从1万到20万窗长W按“覆盖至少4个基波周期”的原则同步放大输入信噪比统一为5dB谐波构成为基波加二次、三次谐波。数据长度L完整SVD截断随机SVD软阈值输出SNR提升(随机SVD方案)10000约2.5秒内存占用约300MB约0.3秒内存占用约120MB约9.1 dB50000约40秒接近内存上限约2.6秒内存占用约600MB约8.7 dB200000无法直接运行约18秒内存占用约2.1GB约8.3 dB这组数据不是我为了展示效果而刻意美化。随着L增大随机SVD软阈值的时间增长大致是线性的而完整SVD的增长接近二次甚至三次。到20万点时显式构造Hankel矩阵已经要十几个GB内存完整SVD基本不可行。当然20万点用函数句柄版本跑也要注意矩阵运算方式第5.3节我再细说。输出SNR随着L增大略有下降不是因为算法变差了而是大矩阵下窗长W没法无限放大奇异值谱的分辨率受限。但在实际应用里8dB以上的信噪比改善对后续频谱分析已经完全够用。4.2 软阈值参数的敏感性分析软阈值方案最让我放心的就是它对τ不那么敏感这正好对应了标题里“健壮”两个字。下面这组数据取自L10000、输入SNR5dB的仿真τ0是第3.2节公式自动估计出的阈值τ倍数0.5×τ00.75×τ01.0×τ01.5×τ02.0×τ0输出SNR提升(dB)9.49.79.28.57.3从0.5倍到1.5倍输出结果都维持在8.5dB以上这在实际工程里就是一个“不用怎么调”的状态。对比硬截断SVD秩k从6变到10时输出SNR改善可能是这样的6→4.5dB7→8.8dB8→9.1dB9→3.2dB10→2.1dB。秩估偏两个数结果就崩了。软阈值显然更符合“拿到数据就能跑”的预期。我在工程里判断阈值是否合适的办法很简单去噪后做一次FFT看频谱。如果噪声底座仍然明显高于两侧背景说明τ偏小如果谐波峰值都出现明显的“削顶”迹象说明τ偏大。根据这个反馈把τ乘以1.3或者0.7一次就能调到合适位置比反复猜k省事得多。5. 常见问题与排查技巧实录5.1 随机SVD结果每次不一样随机SVD的结果天然带随机性因为投影矩阵Ω是随机生成的。如果同一组数据跑两次输出SNR在小数点后第二位可能略有差异。这在工程上是正常的但如果你想做严格对比或者需要可复现的批量处理有两个习惯一定要养成。第一个是在调用随机SVD之前固定随机数种子rng(0)或rng(2025)都行。注意要在构造噪声之前也固定一次否则连测试信号都跟着变。第二个是适当增大过采样p和幂迭代次数q。如果发现两次运行的结果差异明显多半是p取得太小或者q0导致奇异子空间捕捉不够完整。我一般p不小于5q不小于1。如果你的数据量不算特别大比如几万点还有一个验证手段用随机SVD的结果和完整SVD对比奇异值。两者前20个奇异值的相对误差在1%以内就说明参数没问题。误差偏大时先加p再考虑加q。5.2 去噪后波形被削平或噪声残留明显这两个现象是去噪成败的直接信号但处理方法正好相反。波形峰值被削平、谐波幅度明显下降说明阈值τ偏大软阈值收缩过度了。这时把τ缩小一些比如乘以0.6到0.7重新跑一遍。还有一种可能是W选得太小奇异值谱没有把弱谐波和噪声分开弱谐波对应的奇异值也被当成噪声给收缩了。这种情况光调τ没用把W增大到覆盖更多基波周期会好很多。噪声残留明显去噪后频谱底噪还是很高说明τ偏小或者k留的候选奇异值太少部分噪声分量根本没进到后面的软阈值环节。先按1.3到1.5倍放大τ试试如果还不行就增大k让随机SVD多算几个候选奇异值。我遇到过一次样本点特别短的情况尾部MAD估计失真后来直接把τ设成尾部奇异值均值的3倍才压住噪声。5.3 数据量太大矩阵存不下这是大数据集最现实的一关。L过百万时显式构造Hankel矩阵几乎不可行。解决办法是把矩阵改成“隐式算子”不存储H只定义H乘以向量、H转置乘以向量的规则随机SVD整个过程只依赖这两种运算。Matlab里可以写两个局部函数。H乘以一个随机投影矩阵Xn×l时利用卷积关系function Y H_forward(X, x, W) % X 为 n x l 矩阵返回 H*XH是W行 n列的Hankel矩阵 N size(X, 1); Y zeros(W, size(X, 2)); for c 1:size(X, 2) tmp conv(x, flipud(X(:, c))); Y(:, c) tmp(N : N W - 1); end end对应的H转置乘以一个矩阵DW×lfunction Z H_adjoint(D, x, W) % D 为 W x l 矩阵返回 H*DH为 n x W矩阵 N numel(x) - W 1; Z zeros(N, size(D, 2)); for c 1:size(D, 2) tmp conv(flipud(D(:, c)), x); Z(:, c) tmp(W : W N - 1); end end然后用一个接受函数句柄的随机SVD版本替换原来的版本function [U, S, V] rsvd_op(H_forward, H_adjoint, m, n, k, p, q) l k p; Omega randn(n, l); Y H_forward(Omega); for i 1:q Y H_forward(H_adjoint(Y)); end [Q, ~] qr(Y, 0); B H_adjoint(Q); % 注意这里先算A*Q再转置成Q*A [U_B, S, V] svd(B, econ); U Q * U_B; U U(:, 1:k); S S(1:k, 1:k); V V(:, 1:k); end这样做最大的好处是内存占用从“矩阵大小”降到“投影矩阵大小”。L200万、W1万时H本来有约1万×199万接近15GB用算子版本后中间变量最多几十MB。别小看细节里的转置BH_adjoint(Q)这一步是很多人写错的地方H_adjoint返回的是A*Q要转置一次才是Q*A。5.4 一些容易踩的Matlab小坑第一个是hankel函数的用法。hankel(x(1:W), x(W:end))要求第二输入是矩阵最后一列的完整数据很多人传错成x(W:L)的选取范围结果矩阵形状不对。第二个是内存碎片问题在循环里不断给大数组赋值Matlab可能频繁复制造成内存峰值几乎翻倍。建议一次性预分配好变量比如Y zeros(W, size(X,2))这种写法避免动态扩展。第三个是NaN值数据采集偶尔会有坏点Hankel矩阵里只要有一个NaNSVD结果就会全部NaN。预处理阶段一定先用fillmissing或线性插值把坏点处理掉。还有一个容易被忽略的点随机SVD的svd(B, econ)在B是l×n且l n时返回的V是n×l矩阵截断到k列没问题。但如果n lecon返回的矩阵形态会变这时记得先确保l ≤ n。实际使用中l通常远小于n问题不大但如果你把k和p设得很大就有可能在边界上翻车。一点个人体会整套方案跑下来我的直接感受是随机SVD真正解决的是“算不动”软阈值真正解决的是“不知道留几个分量”。它们俩合在一起才让我敢把去噪流程直接怼到几十万上百万点的实测数据上。这个组合在Matlab里实现起来并不复杂核心代码不到一百行但有三个点值得你多花时间一是窗长W要覆盖足够多的基波周期二是k宁可多留余量三是阈值估计时尾部奇异值样本不能太少。最后分享一个小技巧。如果你要处理的是在线采集的流式数据可以考虑把随机投影矩阵Ω固定住然后通过增量方式更新QR分解和B矩阵。数据一批一批进来时只需要在已有子空间上做修正而不是每次从头做随机SVD。这样谐波去噪就能从离线变成准实时每次更新的计算成本会低一个量级。工程上这个方向比直接套离线算法要实用得多有空可以试试。