做信号去噪这些年奇异值分解一直是我工具箱里优先级很高的方法。尤其是处理谐波类信号——电网电压电流波形、旋转机械的振动信号、结构响应里的周期性分量——把一段信号排成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 pk是目标秩p是过采样参数。计算Y AΩ。这一步把矩阵A“压”到一个低维空间里Y的列就是随机方向上的投影。对Y做QR分解Y QR得到一组正交基Q。这一步是把投影结果“掰直”让列向量之间互不相关。计算小矩阵C AᵀQ对C做经典SVD得到C Uc·S·Vcᵀ。回代U Q·VcV 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(ij-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-i1个位置Aᵀ·Y同理。读者如果不放心可以先用一个3×3的Hankel小矩阵手工验证一遍确认无误再放心跑大数据。我每次把这套代码换到新工程时都会先用小规模数据验证这俩函数避免索引偏移这种低级错误。3.3 封装随机SVD函数有了乘法算子随机SVD函数就顺理成章了。输入是信号x、Hankel行数m、目标秩k、过采样p、幂迭代次数q输出是该Hankel矩阵的前r kp个奇异值及对应奇异向量。function [U, S, V] rrsvd_hankel(x, m, k, p, q, seed) % 随机SVDA ≈ 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(kp, 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代表我们猜测信号子空间的维度。谐波信号的理论秩上限是2HH是谐波个数。对一个含基波、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块长mkpq阈值系数c5,000全量2,50030512.050,000全量25,00040512.0200,00050,00025,00040522.01,000,00050,000分块25,00040512.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 kp完整保留在某些超大矩阵下也会占用不少内存。如果内存依旧吃紧只保留前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加速软阈值稳压”的组合在当下数据规模越来越大的信号处理场景里应该还能再用好几年。