做信号去噪这些年，奇异值分解一直是我工具箱里优先级很高的方法。尤其是处理谐波类信号——电网电压电流波形、旋转机械的振动信号、结构响应里的周期性分量——把一段信号排成Hankel矩阵，再对奇异值动点手脚，重建出来的波形干净程度远超高通低通那套时域滤波。但最近两年数据规模涨得实在太快，几十万采样点已经是家常便饭，动不动上百万点。早期习惯用的svd(A, 'econ')在这种规模下基本就是卡死，内存和时间双双爆炸，逼着我去翻了随机矩阵投影那套理论。标题里的完整套路其实就一句话：基于随机奇异值分解和软阈值，在大数据集上做谐波去噪，用Matlab实现。这里说的“谐波去噪”不是把谐波去掉，而是从带噪观测里把谐波成分干净地捞出来，谐波是我们要保的信号，噪声才是要杀的。这篇就把原理、可复现的Matlab实现、参数调试方法，以及我踩过的几个坑一次性写清楚。
1. 为什么谐波去噪需要“随机SVD+软阈值”
1.1 谐波信号在Hankel矩阵里天然是低秩的
先解释一个基础问题：谐波信号凭什么能用奇异值分解去噪？把一段离散信号x(1), x(2), ..., x(n)按滑动窗口排成一个矩阵，行数和列数分别是m和L，满足m + L - 1 = n，矩阵的第(i, j)个元素是x(i + j - 1)。这个矩阵就是Hankel矩阵，也叫轨迹矩阵。
一个由多个正弦波叠加而成的信号，对应Hankel矩阵的秩非常低。理论上，包含H个不同频率谐波的纯净信号，其无穷维Hankel矩阵的秩不超过2H。有限长度下秩会略高一点，但远小于矩阵维度。加入白噪声以后，噪声分量会让矩阵“变得满秩”，奇异值谱上表现为：前面几个奇异值很大，代表谐波信号；后面一大堆小奇异值平缓拖尾，代表噪声。这个特性是SVD去噪的核心依据——把大的奇异值留下，把小的奇异值压掉，再反变换回时域，噪声就被滤掉了。
为什么这个方法比FFT滤波更吸引人？因为FFT需要知道谐波频率在哪里，而SVD方法不需要预先知道任何频率信息，完全由数据自适应地决定哪些成分保留。遇到频率轻微漂移、间谐波混入、多个相近频率分量叠加的情况，SVD类方法仍然能稳定工作。
1.2 经典SVD在大数据下会卡死在内存和时间上
经典SVD的问题是复杂度太高。对一个m × L的矩阵，经典SVD的时间复杂度大致是O(mL·min(m,L))，空间复杂度是O(mL)。svd(A, 'econ')在矩阵只有几千乘几千时没有任何问题，但数据量一旦上来就完全失控。
举个具体的数：我们处理10万个采样点，Hankel矩阵取m = 50000、L = 50001，显式构造这个矩阵需要存50000 × 50001 ≈ 2.5 × 10^9个元素，用double存储就是大约20GB内存，普通工作站直接爆掉。就算你能忍住不显式构造，经典SVD需要迭代的复杂度也让人绝望。这就是为什么“大数据集”这个关键词出现以后，必然要引入随机SVD。
随机SVD的思路非常直接：我们不求完整的奇异值分解，只求出矩阵最前面的几十个奇异值和奇异向量，把这些当成信号子空间。谐波去噪本来就只需要保留前若干个奇异值，随机SVD正好精确命中这个需求。
1.3 软阈值比硬截断稳得多
传统SVD去噪最粗糙的做法是把第k个之后的奇异值直接置零，这叫硬截断。硬截断的问题是：被保留的奇异值里依然藏着噪声能量，而且阈值附近一点点变化会导致重构波形产生很明显的不连续感。
软阈值的操作则是对每个奇异值做收缩：
soft(σ) = max(σ - τ, 0)
其中σ是某个奇异值，τ是阈值。大于τ的奇异值被削减了τ那么多，小于τ的直接变成0。这么做的效果是：大奇异值虽然被保留，但也被“剥掉一层皮”，这一层皮通常就对应着噪声污染；小奇异值则干脆清零。从数学上看，软阈值对应的是奇异值空间里的凸近邻算子，与矩阵核范数最小化有直接理论联系，所以它比硬截断稳定得多。
实际经验里，硬截断在谱图上如果切的位置稍微偏一点，重构波形就会出现毛刺；软阈值对这种偏差的容忍度高不少，即使τ定得不太准，波形整体也不会剧烈恶化。这也是我最终整套方案采用“随机SVD + 软阈值”而不是“随机SVD + 硬截断”的直接原因。
2. 随机SVD的核心部分：过采样与幂迭代
2.1 五步说清随机SVD在干什么
随机SVD算法是Halko等人在2011年左右系统整理出来的，思路可以压缩成五步：
- 生成一个
L × r的高斯随机矩阵Ω，其中r = k + p，k是目标秩，p是过采样参数。 - 计算
Y = AΩ。这一步把矩阵A“压”到一个低维空间里，Y的列就是随机方向上的投影。 - 对
Y做QR分解，Y = QR，得到一组正交基Q。这一步是把投影结果“掰直”，让列向量之间互不相关。 - 计算小矩阵
C = AᵀQ，对C做经典SVD，得到C = Uc·S·Vcᵀ。 - 回代：
U = Q·Vc，V = Uc，于是A ≈ U·S·Vᵀ。
为什么要多此一举先投影再分解？因为C的尺寸是L × r，远小于原始m × L矩阵。经典SVD的复杂度是O(L·r²)，相比原始矩阵的O(mL·min(m,L))，复杂度直接降了一个数量级。
在Matlab里，第4步甚至可以直接用svd(C, "econ")，第5步做两次矩阵乘法即可。整个流程对于只关心前几十个奇异值的场景，精度足够，速度却快了几十倍。
2.2 过采样参数p和幂迭代次数q为什么必不可少
随机SVD有个明显弱点：如果原始矩阵的奇异值衰减不够快，直接投影得到的Q就不够逼近真正的信号子空间。谐波信号的主奇异值通常很大，但噪声带来的小奇异值数量众多，谱尾衰减不快，这时候就需要两个补救措施。
过采样p：目标秩k之外多算几个随机方向。我常用p = 5或p = 10。多加的几个方向不是浪费，它们能把那些处于信号与噪声边界上的奇异值“兜住”，避免漏掉重要的信号成分。
幂迭代q：在随机投影之后，把Y替换成A(AᵀY)，重复若干次。数学上，(AAᵀ)ᵠ能加速矩阵奇异值谱的衰减，让前面的奇异值相对后面更大。这样再做QR分解，得到的正交基Q就更贴近前k个奇异向量张成的空间。经验上q = 1或q = 2就够了，再往上算力翻倍但收益很小，后面我会单独讲这个坑。
2.3 大数据场景下不能显式构造Hankel矩阵
这可能是整篇文章最容易被忽略、却最影响成败的一点：永远不要用zeros(m, L)去显式构造Hankel矩阵。大数据集下这一步就是自杀。
正确的做法是把Hankel矩阵当成一个线性算子，只定义两个乘法函数：A·X和Aᵀ·Y。这两个乘法在Matlab里可以用conv（卷积）高效实现，不需要存储任何完整矩阵。
推导其实很简单：因为A(i,j) = x(i+j-1)，所以A·X的每一列本质上就是信号x与X的那一列做一次相关操作；Aᵀ·Y同理。用FFT卷积一次完成的复杂度大约是O(n·log n)，比显式矩阵的O(mL)低几个数量级。下面第三节给出完整实现时，我会把这两个函数直接放出来，读者不用再推导，拿去就能用。
3. 完整Matlab实现：从数据生成到信号重构
3.1 构造一份带谐波和噪声的测试信号
先做一份可以反复复现的测试数据。采样率设成10kHz，时长10秒，即10万个采样点。信号包含50Hz基波、150Hz三次谐波、250Hz五次谐波，加高斯白噪声让信噪比大约18dB。
%% 生成测试信号：基波 + 3次谐波 + 5次谐波 + 高斯白噪声 fs = 10000; t_end = 10; t = (0 : 1/fs : t_end-1/fs).'; N = length(t); % 10万点 f0 = 50; x_true = 1.0*sin(2*pi*f0*t) ... + 0.35*sin(2*pi*3*f0*t + pi/6) ... + 0.18*sin(2*pi*5*f0*t + pi/3); SNR_before = 18; % 去噪前信噪比目标 noise_power = mean(x_true.^2) / (10^(SNR_before/10)); noise = sqrt(noise_power) * randn(N, 1); x_obs = x_true + noise;这里x_true是准确的纯净信号，仿真时可以用来计算去噪前后的信噪比。实际工程中没有真值，那就画频谱图对比去噪前后的谐波谱线清晰度。
3.2 两个Hankel矩阵快速乘法函数
这是整套实现的地基。第一个函数计算A·X，第二个计算Aᵀ·Y。它们只用conv，不显式构造矩阵。
function Y = hankel_mul(X, x, m) % 计算 Y = A*X % A 是 m 行、L 列的 Hankel 矩阵，由信号 x 生成，不显式构造 % X 是 L*r 矩阵 n = length(x); L = n - m + 1; fx = flip(x); Y = zeros(m, size(X, 2)); for c = 1:size(X, 2) tmp = conv(fx, X(:,c)); Y(:,c) = tmp(n - (1:m)' + 1); end endfunction X = hankel_adj_mul(Y, x, m) % 计算 X = A'*Y % A 同上，Y 是 m*r 矩阵 n = length(x); L = n - m + 1; fx = flip(x); X = zeros(L, size(Y, 2)); for c = 1:size(Y, 2) tmp = conv(fx, Y(:,c)); X(:,c) = tmp(n - (1:L)' + 1); end end这两个函数的核心目录是conv(fx, ...)之后取一段索引。为什么这么取，我当初推导验证过好几次，结论是：A·X的第i行对应卷积结果的第n-i+1个位置；Aᵀ·Y同理。读者如果不放心，可以先用一个3×3的Hankel小矩阵手工验证一遍，确认无误再放心跑大数据。我每次把这套代码换到新工程时都会先用小规模数据验证这俩函数，避免索引偏移这种低级错误。
3.3 封装随机SVD函数
有了乘法算子，随机SVD函数就顺理成章了。输入是信号x、Hankel行数m、目标秩k、过采样p、幂迭代次数q，输出是该Hankel矩阵的前r = k+p个奇异值及对应奇异向量。
function [U, S, V] = rrsvd_hankel(x, m, k, p, q, seed) % 随机SVD：A ≈ U*S*V' % A 是由信号 x 隐式定义的 Hankel 矩阵 if nargin < 6, seed = 1; end rng(seed); n = length(x); L = n - m + 1; r = min(k + p, m, L); % 确保r不超过矩阵维度 % 随机投影 Omega = randn(L, r); Y = hankel_mul(Omega, x, m); % 幂迭代 for it = 1:q Z = hankel_adj_mul(Y, x, m); Y = hankel_mul(Z, x, m); end % QR正交化 [Q, ~] = qr(Y, 0); % 小矩阵SVD C = hankel_adj_mul(Q, x, m); % C = A'*Q，尺寸 L*r [Uc, S, Vc] = svd(C, "econ"); V = Uc; % V: L*r U = Q * Vc; % U: m*r end这里把r限制为min(k+p, m, L)是为了防止随机投影参数超过矩阵本身维度。实际使用中k通常是个位到几十，p不超过10，基本不会触到这个下限，但写上这个保护逻辑更稳。
3.4 软阈值收缩与对角平均重构
拿到奇异值和奇异向量之后，下一步是估计噪声标准差并计算软阈值。噪声标准差我用中值绝对偏差估计，这个估计对离群点更鲁棒，比直接用标准差更稳。
% 中值绝对偏差估计噪声标准差 s_hat = median(abs(x_obs - median(x_obs))) / 0.6745; L = N - m + 1; tau = 2.0 * s_hat * sqrt(L); % 阈值系数推荐1.5~3 % 软阈值收缩 sigma = diag(S); sigma_shrunk = max(sigma - tau, 0); % 去掉被彻底压成0的分量 keep = sigma_shrunk > 0; U = U(:, keep); V = V(:, keep); ss = sigma_shrunk(keep);重构信号时我不用显式重建Hankel矩阵，而是利用一个非常省事的性质：Hankel矩阵反对角线平均等价于对U(:,r)和V(:,r)做卷积。这段代码是整套方案里速度优势最明显的地方之一：
% 权重：每个对角线上有效元素个数 w = min((1:N)', min(m, L)); w = min(w, N - (1:N)' + 1); % 对角平均重构 x_rec = zeros(N, 1); for r = 1:length(ss) x_rec = x_rec + ss(r) * conv(U(:,r), V(:,r)); end x_rec = x_rec ./ w;conv(U(:,r), V(:,r))的长度恰好是m + L - 1 = N，这正好对应信号的有效索引范围。w是对角平均的权重，也就是Hankel矩阵里每条反对角线上的元素个数。这样做重构，复杂度是O(r·N·logN)，对十万点数据完全无压力。
3.5 主脚本跑通与结果评估
把上面几段整合成完整的主脚本。为了可复现，初始化随机种子。去噪后信噪比用x_rec与x_true计算，同时输出去噪前后信噪比。
%% 主脚本：随机SVD + 软阈值谐波去噪 clear; clc; close all; rng(2025); % 生成测试信号 fs = 10000; t_end = 10; t = (0 : 1/fs : t_end-1/fs).'; N = length(t); x_true = 1.0*sin(2*pi*50*t) ... + 0.35*sin(2*pi*150*t + pi/6) ... + 0.18*sin(2*pi*250*t + pi/3); noise_power = mean(x_true.^2) / (10^(18/10)); x_obs = x_true + sqrt(noise_power) * randn(N, 1); % 参数 m = round(N / 2); % Hankel行数 k = 40; % 目标秩 p = 5; % 过采样 q = 1; % 幂迭代次数 % 随机SVD [U, S, V] = rrsvd_hankel(x_obs, m, k, p, q, 2025); % 软阈值 s_hat = median(abs(x_obs - median(x_obs))) / 0.6745; L = N - m + 1; tau = 2.0 * s_hat * sqrt(L); sigma = diag(S); ss = max(sigma - tau, 0); % 重构 keep = ss > 0; U = U(:, keep); V = V(:, keep); ss = ss(keep); w = min((1:N)', min(m, L)); w = min(w, N - (1:N)' + 1); x_rec = zeros(N, 1); for r = 1:length(ss) x_rec = x_rec + ss(r) * conv(U(:,r), V(:,r)); end x_rec = x_rec ./ w; % 评估信噪比 SNR_in = 10*log10(sum(x_true.^2) / sum((x_obs-x_true).^2)); SNR_out = 10*log10(sum(x_true.^2) / sum((x_rec-x_true).^2)); fprintf('输入信噪比: %.2f dB\n', SNR_in); fprintf('输出信噪比: %.2f dB\n', SNR_out);我实测这组参数跑10万点数据，普通笔记本电脑上几十秒内能出结果，输出信噪比相比输入能提升10dB以上。不同机器和MATLAB版本会有一些差异，但整体数量级不会有问题。
4. 参数调优经验：m、k、τ的取舍
4.1 嵌入维度m怎么选
Hankel矩阵的行数m直接决定矩阵的长宽比。奇异谱分析里最常用的选择是m = round(n/2)，这时候m和L大概相等，矩阵接近方形，低秩近似的效果最好。我也确实测试过不同m值，结论是n/2附近去噪效果最稳。
但这里有个工程上的权衡。m越大，Hankel矩阵的行数越多，即便我们用卷积实现乘法，m太大也会影响每次卷积后截取的索引范围和整体耗时。当信号长度达到几百万点，直接m = n/2会让L也接近n/2，虽然内存不会爆，但计算时间线性上涨。
我的习惯是：全局处理时m取n/2，但如果数据超过100万点，优先做分块处理，每块长度1万到5万点，块内m取块长的一半。分块之间留少量重叠，最后拼接处做交叉淡化，效果也很不错。
4.2 目标秩k怎么定
k代表我们猜测信号子空间的维度。谐波信号的理论秩上限是2H，H是谐波个数。对一个含基波、3次、5次、7次谐波的信号，2H ≈ 8，保守取k = 20已经非常宽裕。
判断方法有两条经验路径。第一，看奇异值谱的“断崖”：把奇异值从大到小画出来，如果前几十个很大，后面突然平缓，断崖点附近就是信号与噪声的分界。第二，按能量占比选：前k个奇异值平方和占总能量的99%以上，这个k通常也够用。但注意不要机械追求99.9%，噪声多的时候能量法会把秩取得很大，反而把噪声也保留了。
我还习惯在做软阈值之前故意把k调大一点，比如实际估计只需要20个分量，我给定k = 40。多余的分量不是交给秩来砍，而是交给软阈值去压，这样相当于两重保险。阈值会把多余分量压成接近0，并不会污染重构结果。
4.3 软阈值τ怎么调
软阈值τ是整个方案里最需要手感的参数。太小，噪声残留多；太大，谐波幅值被过度压缩。理论上有经典结论：对零均值噪声矩阵，最优软阈值与噪声标准差σ和矩阵维度有关，量级在σ·sqrt(L)到σ·sqrt(2·max(m,L))之间。我用的经验公式是：
τ = c · s_hat · sqrt(L)
其中s_hat是MAD估计的噪声标准差，c控制在1.5到3之间。信号纯净、谐波幅值大时用c = 1.5，噪声强、需要激进滤波时用c = 3。
一个很有效的调试技巧：先用小规模子集跑一次，把c从0.5扫到3.5，画出去噪后信噪比或波形残留随c变化的曲线，找到拐点。不同的现场数据噪声谱差异很大，这个拐点位置并不固定，但扫一遍以后基本能确定合适的区间。这个脚本解决了我不少疑难案例。
4.4 不同数据规模下的参考配置
我整理了几组自己常用的配置，供读者起步参考。注意这些不是绝对最优，但能保证快速看到可行性结果。
| 数据长度n | 块长 | m | k | p | q | 阈值系数c |
|---|---|---|---|---|---|---|
| 5,000 | 全量 | 2,500 | 30 | 5 | 1 | 2.0 |
| 50,000 | 全量 | 25,000 | 40 | 5 | 1 | 2.0 |
| 200,000 | 50,000 | 25,000 | 40 | 5 | 2 | 2.0 |
| 1,000,000+ | 50,000分块 | 25,000 | 40 | 5 | 1 | 2.0 |
幂迭代次数我很少设到2以上。理论上q越大，随机SVD越接近经典SVD，但大数据集下每多一次幂迭代就要多做两次卷积乘法，时间成本线性增加，而且迭代次数过高会让噪声方差也被放大，反而把随机性带来的稳定性抵消掉。
5. 大数据实战中的常见坑与解法
5.1 结果每次跑都不一样，随机性从哪来
随机SVD天然依赖随机投影矩阵Ω，所以两次运行结果会有轻微差异，这是正常的，不是bug。要复现结果，在调用rrsvd_hankel之前固定随机种子，甚至在主脚本开头就rng(seed)。我一般在函数内部也放一个rng(seed)，保证不管外部怎么设置，单次调用结果稳定。
如果同一个固定种子下结果仍然每次不同，检查是不是有其他函数在内部改了随机流，比如randn被并行工具箱调用时会引入额外的随机流状态。这种情况我建议把种子放到函数最前面，并且避免在循环里对randn做多次重复调用。
5.2 重构波形整体变小，幅值被压了
软阈值是收缩算子，不是硬截断，所以重构出的信号幅值天然偏小。这是软阈值的固有特性，不是代码写错了。如果后续要做谐波幅值精确测量，比如电能质量分析里求各次谐波的幅值相位，直接拿软阈值重构结果去算会偏低。
我的解决办法是用“先识别、后校正”的策略：第一步用这套随机SVD+软阈值方案得到干净频谱，确定基波和各次谐波的具体频率；第二步回到原始带噪信号，用最小二乘拟合这些频率分量的幅值和相位，相当于用原始观测数据重新“校准”参数。这样既享受了软阈值的稳定性，又不牺牲幅值精度。
5.3 内存还是爆了，问题出在哪
如果你按上面的代码实现，不会出现显式Hankel矩阵占内存的问题。但如果把svd(full(A))这种习惯带进来，或者不小心在调试时用了A = zeros(m, L)然后逐项赋值，任何优化都救不了你。
另一种隐蔽的内存问题是：随机SVD返回的U和V如果按r = k+p完整保留，在某些超大矩阵下也会占用不少内存。如果内存依旧吃紧，只保留前k列即可，多余分量本来就会被软阈值压掉。我通常会在rrsvd_hankel里加上一个可选参数，是否只返回前k列，工程上很实用。
5.4 幂迭代次数设高一点不是更精准吗
不一定。幂迭代的本意是加快奇异值谱衰减，让QR分解得到的基更贴近主奇异向量。但当噪声不轻时，(AAᵀ)ᵠ会同时把噪声的奇异值放大，迭代次数高了，噪声反而可能被“撑大”，背离去噪初衷。
我自己测试过的经验是：q = 1在绝大多数谐波场景足够，q = 2能拿到接近经典SVD的精度，q = 3以上基本就是性能倒挂。如果你追求极致精度且数据量不大，直接改用经典svd就好，没必要再挣扎参数。
5.5 超大数据集还想再快，怎么做
如果数据长度达到几百万甚至上千万点，即使随机SVD也会吃力。我有两条经验路线。
第一是分块处理。把长信号切成若干有重叠的块，每块长度几万点，块内做随机SVD+软阈值，块间重叠区长约几百点，拼接时用线性交叉淡化。这样总计算量近似与块数线性增长，内存占用稳定。
第二是降采样粗检后再细检。对超长信号先降采样到几万点做一次SVD，提取出主要谐波频率和大概幅值，再回到原始采样率只做窄带滤波或正弦拟合。这套流程特别适合电力系统那种动辄几分钟、几十分钟的连续波形记录。
分块处理时要注意：块边缘会引入不连续，软阈值重构后要避开头尾各m/4个点作为过渡区。这是我踩过多次之后总结出来最实用的经验，不加这个处理，直接拼接的波形在块边界能看到明显的台阶。
最后分享一点实际操作体会
这套代码我已经在不止一个项目里跑过，从实验室仿真的合成信号到现场采集的电能质量波形都用过。最深的感受有两个。第一，别急着调参数，先把Hankel乘法的两个函数验证对，这是所有上层逻辑的最底层地基。当初我第一次把索引写错，去噪效果看起来很合理，但重构幅值系统性偏低，排查了半天才发现是hankel_mul里取反后的索引位置偏了一位。第二，参数组合建议先用小型测试信号快速跑通全流程，再逐步放大数据规模。放大的过程中如果发现信噪比下降，多数不是算法问题，而是m和块长比例没有按n/2重新调整。这套“随机SVD加速+软阈值稳压”的组合，在当下数据规模越来越大的信号处理场景里，应该还能再用好几年。
