最近做一组长时振动数据的时候，我碰到了一个很实际的问题：信号里只有少数几个谐波成分，剩下全是噪声，可数据量是百万点起步，传统SVD那套重型算法根本扛不住。后来我把随机奇异值分解（Randomized SVD）和软阈值（Soft Thresholding）组合起来，才把谐波去噪这件事真正理顺。这篇文章就围绕这套方案展开，包含原理拆解、Matlab完整实现、参数调整经验和大数据场景下的优化思路，适合做信号处理、振动分析、电力谐波治理或者任何需要从长序列里提纯谐波成分的工程师参考。
1. 谐波去噪到底在去什么：Hankel矩阵与奇异值谱
1.1 为什么要把一维信号写成矩阵
谐波去噪，目标是把周期性成分从噪声里分离出来。一个典型的谐波信号可以写成：
x(t) = Σ A_k·sin(2π·f_k·t + φ_k) + noise
这里的难点在于，谐波个数通常不多，但噪声可能很强，时域波形完全看不出来。频域里虽然能看到峰值，但FFT本身对宽带噪声比较敏感，尤其是谐波频率邻近时，旁瓣会互相干扰。
更稳健的做法是延迟嵌入，也就是把一维信号构造成Hankel矩阵：
H(i, j) = x(i + j - 1)，i = 1...m，j = 1...n
如果原始信号长度是N，取m ≈ N/2，那么n = N - m + 1，矩阵H的尺寸大约是m×n，接近方阵。这个矩阵的价值在于，无噪条件下，K个不同频率的正弦成分对应的Hankel矩阵秩不超过2K。意思是，信号的结构信息被压缩在极少数大奇异值里，而噪声能量均匀散布在所有奇异值上。
所以对H做奇异值分解，把前几个大奇异值对应的成分保留、把后面小奇异值对应的成分压掉，本质上就是在做谐波去噪。
1.2 精确SVD为什么在大数据集上吃不消
这个方法的问题在于计算量。对m×n的稠密矩阵做完全SVD，复杂度通常达到O(m·n·min(m, n))。N = 16384点数据，Hankel矩阵约8192×8193，double精度下矩阵本身就要500多MB内存，MATLAB里svd跑一次几十秒属于正常现象。一旦N到了百万级，矩阵尺寸就变成500000×500001，直接构造这个矩阵需要约2TB内存，这已经不是"慢"的问题，而是根本做不了。
MATLAB自带的svds可以算少数几个奇异值，但它本质是迭代法，对Hankel这种稠密且结构特殊的矩阵，收敛速度不稳定，参数换一次就要重新迭代，批量处理完全不行。所以需要一种把计算复杂度压到线性或近线性级别、同时对奇异值估计精度又有保障的方案。
2. 随机SVD与软阈值：提速和稳噪的两个关键
2.1 随机SVD：先压缩列空间，再做小规模分解
随机SVD的核心思想非常直白：与其对整个大矩阵做完整分解，不如先用随机投影把矩阵的列空间压缩到一个低维子空间，然后在这个小矩阵上做精确SVD，最后把结果映射回原空间。
具体步骤是这样：
- 生成一个n×K的随机高斯矩阵Ω，K取目标秩r加上过采样量p，通常p = 10。
- 计算Y = H·Ω，Y的尺寸是m×K，这相当于用随机组合采样的方式探测H的主要列方向。
- 对Y做QR分解，得到正交基Q，此时H ≈ Q·Qᵀ·H，误差由被丢弃的尾部奇异值决定。
- 计算小矩阵B = Qᵀ·H，尺寸K×n。
- 对B做精确SVD：B = Û·Σ·Vᵀ。
- 恢复原始左奇异向量：U = Q·Û。
整个过程的计算量集中在H·Ω和Qᵀ·H这两个矩阵乘积上，乘积次数是K量级，而不是min(m,n)量级，所以复杂度大幅下降。
在实际实现里，为了提升随机子空间的精度，会在第3步之前加幂迭代。也就是把Y = H·Ω替换成Y = H·(Hᵀ·H)ᵠ·Ω。每多做一次q，随机扰动就会被压低，子空间估计就越接近精确主空间。代价是每次幂迭代要额外做两次H与向量的乘积，所以q一般取0~2。谐波Hankel矩阵的奇异值谱衰减极快，前面几个信号奇异值远远大于噪声平台，这种谱结构正是随机SVD误差最小的情况，q取1已经足够稳。
2.2 软阈值优于硬阈值的三个理由
得到奇异值之后，面临一个选择：是直接截断，还是收缩？
硬阈值的做法是：s_i大于τ就保留，否则直接置零。软阈值则是让所有奇异值向零收缩，也就是s_i' = max(s_i - τ, 0)，它并不做"留还是不留"的二值决定。
从数学上看，软阈值操作对应的是核范数正则化问题：
min ½‖H - X‖F² + τ‖X‖*
这个问题的闭式解就是H的奇异值全部做软阈值收缩。硬阈值则对应秩约束问题，解是保留前k个奇异值，但需要精确知道秩k。
在实际工程中，软阈值至少有三个实打实的优势。第一，硬阈值对谱边界极度敏感，奇异值只要稍微碰到阈值线，结果就天差地别；而随机SVD估计出来的奇异值本身带有随机波动，硬阈值会把这个波动放大成完全不同的分解结果。第二，软阈值对弱谐波更友好，当一个谐波分量的奇异值刚好在噪声平台附近时，软阈值不会把它一票否决，而是保留一部分能量，这对后续重构的连续性很有帮助。第三，软阈值不需要精确的秩估计，只要目标秩给得够宽，阈值会把多余的小奇异值自然压掉，这就大大降低了参数敏感性。
2.3 阈值tau的确定：从谱平台到自适应估计
tau是整个算法里最需要经验的参数。我的做法是先用一次快速rSVD拿到奇异值谱，然后画出来看。
理想情况下，奇异值谱会呈现明显的"陡降+平台"形状：前面几个大奇异值是谐波信号，后面一大堆小奇异值构成噪声平台，平台高度基本反映噪声能量。tau取平台高度的0.8~1.2倍比较稳。如果噪声白且平稳，这个区间内结果差别不大；如果信号本身比较干净，可以把tau取到平台值附近，多压一点噪声。
工程上还可以用一个粗略的自适应估计：
tau = 1.2 × median(奇异值谱末尾20个值)
但这个公式只对纯白噪声场景靠谱。如果信号里有调幅或瞬态成分，尾部奇异值不一定来自噪声，得先看谱再下结论。千万别上来就设一个固定tau，不同量纲、不同采样率下奇异值绝对数值可以差很多。
3. Matlab代码实现与关键细节
3.1 核心函数：rSVD、软阈值、对角平均一次搞定
下面是我在实际项目里用的核心函数，把随机SVD、软阈值收缩和信号重构整个流程放在一起。为了让大家直接跑通，这里采用显式构造Hankel矩阵的版本，适合中等规模数据。
function [xr, s_orig, info] = rsvd_soft_denoise(x, r, q, tau, m) % rsvd_soft_denoise 基于随机SVD和软阈值的谐波去噪 % 输入： % x : 一维信号，列向量 % r : 随机SVD目标秩（建议过估计，例如2*谐波数+10） % q : 幂迭代次数，0-2之间常用 % tau : 软阈值，取大于0的值 % m : Hankel矩阵行数，默认floor(length(x)/2) % 输出： % xr : 去噪后的信号 % s_orig : 收缩前的奇异值序列，便于画谱 % info : 附带保留秩、参数等信息的结构体 N = numel(x); x = x(:); if nargin < 5 || isempty(m) m = floor(N / 2); end n = N - m + 1; K = r + 10; % 过采样 Omega = randn(n, K); % 显式构造Hankel矩阵（中等规模直接用） H = zeros(m, n); for j = 1:n H(:, j) = x(j:j+m-1); end Y = H * Omega; % 幂迭代，提升随机子空间精度 if q > 0 for it = 1:q [Y, ~] = qr(Y, 0); Z = H' * Y; [Z, ~] = qr(Z, 0); Y = H * Z; end end [Q, ~] = qr(Y, 0); B = Q' * H; [U, S, V] = svd(B, 'econ'); U = Q * U; s = diag(S); s_orig = s; % 软阈值收缩 s_soft = max(s - tau, 0); keep = s_soft > 0; % 用卷积做对角平均，避免构造完整重构矩阵 xr = zeros(N, 1); for l = find(keep)' seg = conv(U(:, l), V(:, l)); xr = xr + s_soft(l) * seg(1:N); end % 对角平均归一化：每条反对角线的参与元素个数 cnt = min([(1:N)', repmat(m, N, 1), repmat(n, N, 1), (N:-1:1)'], [], 2); xr = xr ./ cnt; info.r = r; info.q = q; info.tau = tau; info.kept = nnz(keep); info.s_orig = s_orig; info.s_soft = s_soft; end调用方式很简单：
load 某振动信号; % 替换成自己的数据即可 N = length(x); m = floor(N / 2); [denoised, s, info] = rsvd_soft_denoise(x, 60, 1, tau_estimate, m);tau_estimate建议先用一次tau = 0的调用拿奇异值谱，肉眼定平台后再填。
3.2 重构阶段用卷积实现对角平均，绕开O(N²)内存
这个可能是整篇文章里最值得借鉴的细节。按照常规思路，重构去噪后的Hankel矩阵需要计算Hr = U·diag(s_soft)·Vᵀ，然后沿每条反对角线取平均。但Hr的尺寸和原Hankel矩阵完全一样，大N下又是一个天文数字的内存占用。
我的做法是直接观察反对角平均的数学结构。对第l个保留的奇异三元组s_l·u_l·v_lᵀ，它的第k条反对角线元素和恰好等于：
Σ_i u_l(i) · v_l(k+1-i) = conv(u_l, v_l)(k)
卷积！Matlab里conv两条长度分别为m和n的向量，得到长度正好N = m+n-1的结果。也就是说，每个保留分量对去噪信号的贡献可以一次性算出来，完全不需要构造Hr矩阵。整个重构阶段只需要做r次卷积，计算量降到O(N)级别，内存占用只剩几个向量。
当时我在实际数据上验证过这个思路和显式构造矩阵逐元素对角平均的结果，差别在浮点误差范围内，但内存占用差了三个数量级。
3.3 超大数据集：把矩阵乘法改成算子
显式构造Hankel矩阵在N超过50万点之后基本不现实，但随机SVD本身并不需要真正拥有整个矩阵，它只需要能够计算H·v和Hᵀ·u这两类矩阵向量乘积。所以可以把H理解成一个算子，而不是一个实体矩阵。
对Hankel矩阵来说，H·v本质上是对信号x做一次滑动窗加权求和，这是一次互相关运算，可以用FFT加速到O(NlogN)，不需要构造任何中间矩阵。具体实现思路是补零后用ifft(fft(x补零) × fft(v翻转补零))，再截取有效部分，Hᵀ·u同理。这套算子化方案在百万点数据上依然流畅，耗时可从"不可接受"降到"几秒级别"。
给一个忠告：在使用算子版本之前，务必先用小规模矩阵同时跑一次显式版和算子版，对比H·Ω的结果是否一致。FFT实现里边界对齐和补零位置很容易出错，这种验证能帮你省下大量调bug的时间。
4. 实验对比：从仿真信号到百万点数据
4.1 仿真设计：三个谐波加白噪声
为了验证这套方案，我构造了一组仿真信号：
- 采样率Fs = 1024 Hz
- 信号长度N = 16384
- 三个谐波：50Hz幅度1.0，120Hz幅度0.5，230Hz幅度0.3
- 叠加高斯白噪声，使信噪比约0dB
这个信噪比下，时域波形已经完全看不出周期结构，频域也只能勉强看到50Hz的峰，120Hz和230Hz基本淹没在背景里。
运行上面的核心函数，设r = 30，q = 1，tau用谱平台法确定。去噪后的信号与真实谐波成分做对比，输出信噪比从0dB提升到17.8dB左右，三个频率的幅度估计误差都在5%以内。
4.2 方法对比：svd、svds、rSVD+软阈值
同样数据和同样机器上，我做了几组对比，结果如下表。注意具体耗时随硬件和MATLAB版本会浮动，看相对关系就行。
| 方案 | 运行时间 | 输出信噪比 | 主要问题 |
|---|---|---|---|
| 完全SVD + 硬阈值截断 | 约30秒 | 15.9dB | 需要精确知道保留几个奇异值 |
| 完全SVD + 软阈值 | 约30秒 | 18.1dB | 慢，内存占用大 |
| svds(40) + 软阈值 | 约8秒 | 17.6dB | 收敛不稳定，参数敏感 |
| rSVD + 硬阈值 | 约1秒 | 14.2dB | 阈值附近极不稳定 |
| rSVD + 软阈值 | 约1秒 | 17.8dB | 稳定，快，推荐 |
硬阈值方案在rSVD下明显吃亏，就是因为随机估计带来的奇异值波动被二值判断放大了。而rSVD加软阈值无论在速度还是稳定性上都和完全SVD加软阈值基本持平，但耗时只有后者的三十分之一。
4.3 百万点规模：内存与耗时的实测感受
为了模拟真实大数据场景，我把信号拉到N = 1000000点，谐波参数不变，噪声水平也保持0dB左右。
显式构造Hankel矩阵这条路首先就断了，矩阵内存接近2TB。所以这一轮测试我把核心函数改成算子版本，只保留H·v和Hᵀ·u的FFT实现。目标秩r = 64，q = 1，总耗时大约在4~6秒量级，保留的奇异值数量自动收敛到6个左右，正好对应三个谐波。
这个结果已经完全可以支撑"先整体去噪再分析"的工作流。如果换成完全SVD或svds，时间单位要变成分钟甚至小时，而且每次调参都要重新算一遍，交互体验完全是两个世界。
5. 工程调参与踩坑手记
5.1 秩r怎么给：高估不亏，低估全丢
随机SVD和软阈值这套组合的优点在于，r可以不那么准。谐波数未知时，把r设成谐波数可能上界的2倍再加10，也就是哪怕你怀疑有8个谐波，r设个40都行。软阈值会把多余的伪奇异值压掉，不会对结果产生实质影响。
但r千万别设小。一旦r小于真实信号子空间的维度，有用成分还没来得及被软阈值判断就被随机投影丢掉了，这个损失是不可逆的。我踩过这个坑，一开始按谐波数严格设r，结果某个弱谐波奇异值排在偏后位置，直接被截掉，事后排查很久才发现是秩低估。
5.2 幂迭代次数q怎么选
q = 0时速度最快，但对随机矩阵的采样质量要求高，偶尔会出现谱前端出现混叠的情况，也就是信号奇异值和噪声平台边界模糊。谐波场景我推荐q = 1，这一个额外代价换来的子空间稳定性非常值得。q = 2的边际收益很小，但计算量翻倍，除非你处理的是奇异值谱衰减很慢的宽带信号，否则没必要。
5.3 噪声不是白噪声怎么办
软阈值对高斯白噪声效果最好，因为噪声在Hankel矩阵里产生的奇异值会形成一个均匀平台。实际情况往往没这么干净：脉冲噪声、偶发冲击、数据缺失都会在奇异值谱上制造孤立的"假高峰"，这些高峰的幅度可能和真实弱谐波相当，软阈值单独处理不了。
我的经验是，在进Hankel矩阵之前先做一步鲁棒预处理。比如幅度超过5倍中位数绝对偏差的采样点先做限幅，或者用中值滤波把脉冲毛刺压一遍。这一步不追求净化信号，只为防止脉冲噪声污染奇异值谱的前端结构。等主成分提取完，再去处理局部细节。
5.4 边界效应和重构误差
对角平均在信号两端参与平均的元素数量很少，所以边界附近的重构误差天然比中间大。如果后续分析依赖整段信号，建议留出首尾各5%~10%的余量，分析时直接裁剪掉。另外，卷积对角平均实现里，如果u或v没有按单位向量归一化，重构幅度会整体偏移。svd输出的奇异向量本身是单位范数，但经过U = Q·Û这一步后，由于CGS或QR数值误差，极少数情况下范数会轻微偏离1。如果发现重构信号整体偏大或偏小，检查一下保留奇异向量的范数。
6. 还能往哪个方向扩展
6.1 从批量去噪到分块流式处理
百万点还能整段处理，再往上比如连续监测数据一天几千万点，整段做反而不划算。工程上的做法是分块：把信号切成若干段，段与段之间重叠50%，每段独立做rSVD加软阈值去噪，最后在重叠区做线性加权拼接。这个方案的好处是每段局部谱结构变化时能自适应，非平稳谐波也能跟上变化，缺点是拼接处可能出现幅度调制痕迹，重叠区的权重设计要平滑过渡。
6.2 把核范数正则与ADMM结合
如果去噪对象不只是简单谐波，还要同时分离稀疏冲击成分，可以把模型扩展成低秩加稀疏分解，也就是H = L + S。低秩部分代表谐波，稀疏部分代表冲击，用软阈值交替更新两个变量。核心的低秩子问题还是用随机SVD加软阈值来解，稀疏子问题用逐元素软阈值。这套思路在处理机械故障信号时非常实用，谐波和故障冲击往往叠加在一起，单靠一路去噪会互相污染。
6.3 可复现性和进一步优化
随机算法绕不开随机数，建议在调用前设置rng固定种子，否则每次分解结果虽然有微小差异，但做严格的指标对比时会引入额外方差。进一步优化可以考虑把函数改成完全算子化后用MATLAB的gpuArray直接跑，FFT和矩阵乘法在GPU上都能获得明显加速。我自己在NVIDIA入门级显卡上实验，算子版rSVD还能再快两三倍，对于实时性要求不高的批处理任务，这点加速已经非常可观。
我个人的整体体会是：谐波去噪这套组合拳，真正把"大数据"和"高鲁棒"这两件事从矛盾变成了互相成就。随机SVD负责把算力门槛降下来，软阈值负责把参数敏感度降下来，两者缺一都会让工程落地变得很别扭。如果你手头也有长序列谐波提取的需求，建议先拿仿真数据跑通，再逐步替换成自己的实测数据，过程中的谱图变化会告诉你下一步该调什么。
