做信号处理的人应该都有过这种体验现场采集的电压、电流、振动或声学信号里谐波成分就藏在噪声底下想提取特征、做分析、诊断问题第一步往往是先去掉噪声。传统做法里奇异值分解SVD是相当经典的一招把信号构造成矩阵做一次完整SVD留下前面几个大的奇异值重构回去就能把主成分捞出来。这个方法效果确实不错但有一个很现实的问题——数据量一旦涨上去比如连续采集了几百万个采样点、矩阵动辄几千乘几千甚至上万乘上万直接调系统的svd()函数会非常吃力内存和耗时都容易失控。而我最近在Matlab里反复验证下来一套“随机奇异值分解 软阈值”的组合拳能在大数据集上把谐波去噪这件事做得又快又稳。这篇文章就把这套方法的原理、代码实现、参数设置和踩坑经验完整梳理一遍。这套方案适合谁如果你在处理谐波、振动、电力质量等时间序列数据矩阵规模大到你不敢直接做全SVD又不想牺牲太多去噪精度那这条路子很值得试试。我尽量按一个从业者的口吻写把方法背后的“为什么”也讲透不是只丢给你一段能跑的代码。1. 谐波去噪到底在做什么为什么SVD能派上用场1.1 谐波去噪的场景与难点谐波去噪不是一个抽象的概念它在很多实际场景里都在出现。电力系统里电网电流和电压中含有基波和整数倍频的谐波分量被干扰和测量噪声污染后需要把谐波参数估计出来机械设备状态监测里齿轮箱和轴承的振动信号中往往存在轴频及其倍频成分也就是谐波族但现场采集信号夹杂了大量随机冲击和环境噪声通信和水声领域同样需要对带噪信号里的周期、准周期成分做分离。这类问题的核心难点在于谐波分量通常只占据整个频带的一小部分能量但是它在时域里是持续存在的周期成分。噪声则可能是宽带随机噪声、脉冲干扰或者是同频段的非平稳成分。简单的高通或低通滤波很难在保住谐波形状的同时完全去掉噪声而自适应滤波又需要参考信号或迭代收敛在数据量大的时候并不友好。SVD去噪的切入点是把一维信号重构成二维矩阵然后利用矩阵的奇异值分布来区分“信号子空间”和“噪声子空间”。谐波对应的周期成分会使矩阵产生少数几个较大的奇异值而随机噪声经过矩阵化之后能量被摊到大量小奇异值上。用大奇异值重构信号实际上就是在做一个基于数据本身的自适应滤波不需要额外知道谐波频率和幅值因此适合频谱未知、谐波数量较多的场景。1.2 为什么大数据集让传统SVD难以为继我自己踩过这个坑。早先处理一个长度为几十万点的振动信号按经典的奇异谱分析SSA方法构建Hankel矩阵后矩阵维度大概是5000乘5000左右Matlab里直接调svd()跑完大概需要几百秒内存峰值也高得吓人。后来数据变成几百万点矩阵马上变成上万乘上万我当时想都没想就调svd()结果机器风扇狂转内存告急最后只能强行取消任务。根本原因是完整SVD的计算复杂度大致是O(mnmin(m,n))对于2000乘2000的矩阵还能忍受但对于10000乘10000矩阵运算量是前者的上百倍。而且svd()对整个奇异谱全部求解但我们做去噪时真正需要的其实只有前面几十个甚至几个奇异值和对应向量。这就好比你去超市买瓶酱油却非要把整个货架打包搬回家资源和时间都浪费了。所以问题就变成能不能用一种近似但足够精确的方式只算出信号子空间的前k个奇异值和奇异向量然后在这个低秩近似的基础上做软阈值收缩这正是随机奇异值分解Randomized SVD解决的痛点。2. 核心方法拆解随机SVD加软阈值组合起来为什么这么顺2.1 随机奇异值分解的基本原理随机SVD的核心思路是用一个低维随机投影先把原始大矩阵“压”到一个很小的子空间上再在这个小子空间里做标准SVD。通俗点说一个身高体重特征特别多的大型数据矩阵先把它的“最显著方向”捕捉到一个低维切片里然后只对这个切片做精细分解最后再还原回原空间。具体实现是一个很经典的流程。假设我们面对的是m乘n的矩阵A希望得到前k个左奇异向量U、奇异值S和右奇异向量V。做法是先构造一个n乘L的高斯随机矩阵Omega其中L k pp是过采样参数。然后计算Y A * Omega相当于把A的列空间通过随机投影采样了一遍。接着对Y做QR分解得到正交基Q再把A投影到Q的列空间上得到小矩阵B Q * A。B的维度只有L乘n已经小了很多对它做标准SVD几乎不费什么时间。最后把B的左奇异向量通过Q映射回原空间就得到了A的近似左奇异向量。这里有一个很重要的细节如果只做一次投影随机误差可能比较大。工程上通常会加入一步幂迭代power iteration把Y更新为A*A*Y的Q因子循环一两次。幂迭代能增强大奇异值对应的方向抑制小奇异值对应的随机成分相当于让低秩近似更贴近真正的信号子空间。代价只多了一两次矩阵乘法和QR分解对大数据集来说依然比全SVD划算得多。我之前看到很多讲解随机SVD的帖子会把“随机投影幂迭代”当成一个黑盒但对于谐波去噪场景幂迭代尤其有用。因为谐波对应的奇异值虽然相对较大但和某些强噪声相比并没有绝对优势直接一次投影的话被投影丢失的方向可能正好是某个谐波分量。加上一次幂迭代之后结果会稳定很多。2.2 软阈值收缩在去噪中的角色得到随机SVD分解后通常的做法是直接截断只保留最大的k个奇异值其余置零然后重构信号。这一步是硬阈值。但硬阈值有一个让人头疼的问题它对奇异值“一刀切”保留的保留删掉的删掉重构时因为截断不连续信号里容易出现局部振荡特别是谐波和噪声之间界限模糊的时候。软阈值做得更细腻一些。它把每个奇异值都向零方向收缩固定的量tau公式是s_tilde max(s - tau, 0)。大于tau的奇异值不只是保留而是被缩小了tau小于等于tau的奇异值直接置零。这样做的好处是它在横跨阈值时是连续变化的不会出现从保留到删除的突然跳变重构出来的波形更平滑对谐波基波的相位保持也更好。我们可以打个比方。硬阈值像剪树枝对着一条线“咔”一下剪断断口容易开裂软阈值像修枝剪完之后还把切口修掉一圈树干表面更平滑。对谐波信号这种对波形连续性要求很高的场景软阈值的这种“收缩”动作能减少重构信号的畸变。这里还要说清楚一个容易混淆的点软阈值通常用在稀疏表示和小波去噪里针对的是系数域。把软阈值用在SVD奇异值上道理是一样的——奇异值本身就是矩阵的“幅值系数”噪声在奇异值谱上表现出较小的数值收缩大数值、抹掉小数值天然就是一种去噪操作。2.3 为什么说这套组合“健壮且高效”标题里“健壮”和“高效”两个词不是我乱贴的。高效来源于随机SVD只关心低秩部分计算量从O(mnmin(m,n))降到O(mnL)加上几次QR分解实际测试下来对一个大矩阵而言往往有一个数量级以上的提升。更关键的是随机SVD还可以用函数句柄形式实现A本身不用显式构造只需要定义A*x和A*y两个乘法规则。这样即使Hankel矩阵大到没有内存能完整放下算法也照样可以跑这是传统svd()完全做不到的。健壮性则体现在两个方面。一方面随机SVD使用随机投影只要过采样参数p和幂迭代次数q设置合理它的近似误差有理论保证而且对奇异值分布比较“平”的情况也能稳得住另一方面软阈值比硬阈值对噪声水平的估计误差更不敏感即使tau稍微偏离最优值重构结果也不会突然崩坏。这意味着在工程数据上不需要每次反复调参只要参数在一个合理范围里结果都可用。我自己把随机SVD和传统SVD放在同一组数据上对比过当矩阵是8000乘8000、需要前30个分量时直接svd()跑了大概几十秒随机SVD加两次幂迭代只用了不到两秒重构信号和完整SVD截断结果做相关分析相关系数在0.999以上。这个精度换来的速度收益是值得的。3. Matlab代码实现从矩阵构建到随机SVD加软阈值3.1 整体代码框架整个实现可以拆成四块第一把一维信号构造成Hankel矩阵第二用随机SVD求前k个奇异值和奇异向量第三对奇异值做软阈值收缩第四从收缩后的矩阵重构回一维信号。这四块可以做成独立函数便于复用。下面的代码是我在Matlab里整理过的一个完整框架。为了不让示例太复杂我用一组包含基波、三次谐波和五次谐波的模拟信号来演示。Hankel矩阵的构建方法对一维时间序列是通用的。% 生成模拟谐波信号 fs 1000; % 采样率 1000 Hz t (0:1/fs:2-1/fs); % 2秒信号 N length(t); x 1.0*sin(2*pi*50*t) 0.4*sin(2*pi*150*t 0.5) 0.2*sin(2*pi*250*t - 0.8); x x 0.25*randn(N,1); % 添加高斯白噪声 % 构建Hankel矩阵 num 2000; % 窗口长度决定矩阵行数 H zeros(num, N - num 1); for i 1:num H(i,:) x(i:iN-num); endHankel矩阵的每一行都是原信号的滞后版本相当于把原始的时域信号嵌入到一个高维相空间里。这里窗长num选多大比较有讲究。它直接决定矩阵行数也决定算法能捕捉的最长周期成分。一般来说至少要有两到三个完整基波周期比如处理50Hz工频信号、采样率1000Hz时一个周期是20个点选2000点就是100个周期足够充分。如果窗长取得太小谐波信息会被切碎重构信号连续性差取得太大矩阵规模又上去了大数据集压力增大。我自己一般按“基波周期的50到100倍”来初选。3.2 随机SVD的Matlab实现要点下面是随机SVD的标准实现。我把它写成函数方便调换不同的k、p和q参数。在矩阵A特别大的时候可以把A替换成“函数句柄乘法”的形式但在这里先展示常规数值矩阵版本简单好懂。function [U, S, V] rsvd(A, k, p, q) % A: m x n 矩阵 % k: 保留的奇异值个数 % p: 过采样个数默认 k 数量的50%左右 % q: 幂迭代次数默认1或2 if nargin 3, p 5; end if nargin 4, q 1; end L k p; % 高斯随机投影Omega 是 n x L Omega randn(size(A,2), L); Y A * Omega; % 幂迭代增强低秩逼近的稳定性 if q 0 [Y, ~] qr(Y, 0); for i 1:q [Z, ~] qr(A * Y, 0); [Y, ~] qr(A * Z, 0); end end % 投影到低维空间 [Q, ~] qr(Y, 0); B Q * A; % 对B做标准SVD只取前k个 [Uhat, S, V] svd(B, econ); Uhat Uhat(:, 1:k); S S(1:k, 1:k); V V(:, 1:k); % 映射回原空间 U Q * Uhat; end这里有几个细节值得强调。qr(Y,0)是经济QR返回的是m乘L的正交矩阵能避免生成不必要的m乘m大矩阵。B Q * A这一步是核心B的行数只有L之后svd(B, econ)算的就是一个L乘n的小矩阵计算量比原矩阵小很多。如果A本身大到无法承载就需要用函数句柄来避免构造A后面我会单独写。幂迭代的次数不是越多越好。q越大结果越接近完整SVD但每多一次就要多算两次矩阵乘积和两次QR分解。我的实际经验是q取1到2就够了继续加到4或5精度提升有限耗时却成倍增加。尤其是当数据本身信噪比已经比较好的时候q1就能满足后续去噪要求。3.3 软阈值去噪与信号重构得到随机SVD的U、S、V之后对奇异值做软阈值收缩。% 设置软阈值参数 tau s diag(S); tau 0.1 * max(max(s)); s_soft max(s - tau, 0); % 去噪后矩阵重建 S_soft diag(s_soft); H_denoised U * S_soft * V; % 对Hankel矩阵做反对角线平均恢复一维信号 x_denoised zeros(N,1); cnt zeros(N,1); for idx 1:N val 0; total 0; for i 1:num j idx - i 1; if j 1 j size(H_denoised,2) val val H_denoised(i,j); total total 1; end end if total 0 x_denoised(idx) val / total; cnt(idx) total; end end这段代码的tau是临时用一个比值估的并不通用。因为奇异值的绝对量级和数据幅度、窗长都有关系。实际工程里更好的做法是先画出奇异值谱看分布再决定收缩量。不要以为软阈值就像小波去噪那样固定一个值奇异值尺度完全不一样。重构回一维信号的方法叫“对角平均”。因为Hankel矩阵的每个元素来自原信号的某个时间段同一个时间点会被多次用到重建后要对矩阵的每条反对角线取平均才能还原出原始时间轴上的一个值。这个步骤很朴素但也很容易写错。一个常见错误是循环里索引搞混导致去噪后的信号出现时间错位。强烈建议在写完重构后先跑一组已知信号做验证比如直接用正弦信号加噪声看重构结果是否对齐原时间轴。4. 参数到底怎么调大量实验之后的实测经验4.1 保留分量数k的选择随机SVD的第一步是确定保留的奇异值个数k。这个值决定了随机SVD要算多少个分量也决定了去噪后信号包含多少自由度。如果k太小谐波的某个分量会被丢掉如果k太大随机SVD计算量增大而且超过信号子空间的部分会给噪声留出口子。我总结出一个很实用的方法先用小窗长、小规模子集对Hankel矩阵做一次随机SVD看一下奇异值谱。理想情况下奇异值谱会呈现明显的“肘部”前几个大奇异值对应谐波和基波之后有一串平缓的小奇异值对应噪声。保留数k就取肘部之前的分量个数通常等于“基波和主要谐波数量的两倍”。为什么要两倍因为谐波在构造的复矩阵或实Hankel矩阵中往往对应一对共轭奇异值单独取一个会把半边相位丢掉。如果数据里谐波数量是3个k6就是一个很自然的起点。然后在6附近做小范围的扫描对比重构信号的频谱和实测频谱。谐波峰值都清晰可见且高次谐波不被压掉就可以确定下来。这个方法不需要依赖自动判断直观又可靠。4.2 软阈值tau的确定与过采样、幂迭代的搭配软阈值tau是这套方法里最需要手动把握的参数。很多刚接触的人喜欢直接套小波去噪里的公式但奇异值分布并非高斯分布Donoho阈值并不直接适用。我的经验是tau一般设置为噪声子空间奇异值平均值的一点五到三倍。可以先取前k个之后的奇异值也就是第k1个到第L个视为噪声子空间的代表性样本计算它们的平均值。把tau定在这个平均值附近能保持信号分量收缩不过度也能压低噪声基底。要注意tau太大和太小的结果完全不同。tau设太小噪声残余明显去噪后频谱里依然有一条杂散宽带tau设太大信号主成分的幅值会被压小谐波幅值估计会偏低后续做幅值诊断时就会出错。我这里也踩过坑曾经为了追求干净频谱把tau调得很高结果去噪后的基波幅值比真值低了快一半。所以调tau时一定要同时关注重构信号的幅值保持情况而不只看频谱干不干净。为方便复现我给一组常用的起始参数k取估计主分量数的两倍p取5到10q取1tau取0.05到0.2倍的第一个奇异值。在这个范围里大多数谐波去噪数据都能得到一个可接受的结果。如果发现重构信号有明显噪声泄漏再把p提高到10q提高到2tau稍微上调。不要一上来就盲目拉高q这样浪费时间也不一定提升多少。% 示例从一个完整SVD结果中观察奇异值分布辅助定tau [Ufull, Sfull, Vfull] svd(H, econ); sigma diag(Sfull); semilogy(1:length(sigma), sigma, o-); grid on;这种观察方法在大矩阵上可能不现实可以先在小规模子集上做一次。大数据集上直接算完整SVD并不划算但在调试阶段用小窗长构建一个小矩阵来做参数标定完全可行。4.3 用函数句柄处理超大矩阵避免内存爆炸前面提到随机SVD最大的优势之一就是可以配合函数句柄工作。这在处理真正的“大数据集”时几乎是必选项。因为Hankel矩阵的规模可能达到上万乘上百万显式构建它内存直接爆掉。解决办法是只定义两个乘法规则让随机SVD算法在需要的时候动态计算A*x和A*y而无需把矩阵元素全部存下来。拿Hankel矩阵来说H的第i行其实是从信号x的第i个位置开始截取的一段。那么H乘以一个向量v结果y的第i个元素就是对x的一段做内积而H乘以向量u结果z的第j个元素就是u和x对应窗口片段的内积。用Matlab的conv、convmtx或filter都可以实现不过最简单直接的方式还是写两个匿名函数接受信号x、窗长num和向量v返回乘积结果。Afun (v) hankel_multiply(x, num, v); Atfun (u) hankel_adjoint_multiply(x, num, u); Omega randn(N - num 1, k p); Y Afun(Omega); % 之后所有A*Z和A*Y操作都换成函数句柄调用hankel_multiply的具体写法我这里不展开了核心是理解它不需要存储矩阵。这种隐式矩阵处理方式在Matlab大数据处理里是非常好用的技巧。你可以先在小矩阵上用显式矩阵验证随机SVD结果再把乘法函数替换成隐式版本两个结果完全一致的话说明算子定义没有错。5. 实测效果随机SVD去噪和大规模数据到底合不合拍5.1 一组模拟信号的完整对比我构造了一组仿真测试来检查整套流程信号为50Hz基波加150Hz三次谐波加白噪声数据长度20万点窗长2000因此Hankel矩阵是2000乘198001规模不算特别大但已经超出直接svd()的舒适区。用随机SVD取k6、p10、q1奇异值软阈值收缩后重构。对比去噪前后的信噪比含噪信号SNR约12dB去噪后SNR提高到24dB左右。频谱上50Hz、150Hz处的谐波峰值清晰保留整个频带底噪明显降低。更重要的是重构的时域波形和原始无噪信号的对齐程度很好没有出现硬阈值常见的那种峰值附近振铃。这个结果印证了软阈值在这里确实比硬阈值更合适。作为对照我用同样的数据直接截断随机SVD结果也就是用硬阈值忽略小奇异值重构后SNR只到19dB左右而且波形在谐波突变处有一点点过冲。虽然硬阈值也不是不能用但对谐波去噪这种对幅值和相位都有要求的场景软阈值明显更好。5.2 大数据集上的耗时与内存变化我还测试了一个更大规模的数据100万点信号窗长5000Hankel矩阵是5000乘995001。显式矩阵的话单是存储就要接近40GB根本不可能。但用函数句柄乘以随机投影配合k10、p10、q1整个过程内存占用不超过几百MB耗时约十几秒。如果把同样场景换成传统svd()不用说完整SVD了光是把矩阵写出来内存就已经爆了。即便强行分块计算量也比随机SVD高好几个数量级。这组对比也解释了为什么在数据规模上来后随机SVD几乎是唯一实际可用的选择。省下来的内存和计算时间可以用来调参数、跑多次验证这对工程调试来说价值巨大。5.3 不同信噪比下的稳定表现我又把噪声强度从0.1逐步加到1.0观察软阈值方案的去噪表现。在低噪声下随机SVD和传统SVD的结果几乎一致在高噪声下传统截断方法会出现奇异值排序不稳定、噪声分量混入前几个分量的问题但随机SVD加一次幂迭代依然能稳定地把前几个谐波分量找出来软阈值收缩也一直保持连续没有出现结果反复跳变的情况。这一点的解释是随机SVD的随机投影在构造低维子空间时对奇异值间隔有一种“平均化”的稳定作用。配合幂迭代后大奇异值的方向被强化小奇异值的随机干扰被抑制所以算法对噪声的容忍度比直接截断高。这个特性很关键实际现场数据的噪声特性往往不理想不是教科书里那种漂亮的高斯噪声能有这种健壮性我会更放心用它。6. 常见问题与排查技巧实录6.1 随机SVD结果每次跑出来都不一样怎么办这是绕不过去的第一个问题。因为随机SVD里面用到了randn必然会有随机性。如果你每次跑结果都有差异可以先确定是不是没有固定随机种子。在函数调用前加一行rng(2024)就能保证实验可复现。但要注意即使固定了随机种子矩阵尺寸、过采样参数和幂迭代次数也会影响误差水平。如果固定种子后结果还不太稳定优先做法是把p从5提高到10把q从1提高到2误差会显著下降。有一个更隐蔽的情况你的矩阵A奇异值分布很“平”也就是信号和噪声在奇异值谱上混在一起随机投影容易丢失部分方向。这时候幂迭代的作用特别明显。我建议遇到这种情况时先试q2或者q3然后再看奇异值谱是否稳定。如果依然不稳定很可能是k取得太小信号子空间没有完全覆盖。6.2 重构出来的信号有边缘畸变或首尾失真Hankel矩阵和反对角线平均天然对信号两端处理不充分。因为边界处参与平均的元素少所以去噪后信号首尾段可能失真。这不是随机SVD的问题而是基于轨迹矩阵的去噪方法本身固有的边界效应。解决办法可以在构建Hankel矩阵前把信号两端各扩展一段去噪后再裁掉或者把重复计算的重构部分只取中间稳定区段来分析和诊断。我以前也天真地以为去噪后时间轴所有点都一样可靠直到一次做振动特征提取时发现首尾段幅值和理论上对不上排查半天才发现是边界效应。从那以后我所有去噪结果的下游分析都会刻意跳过前200个点和最后200个点。这个经验虽然小但在工程里很能避免结论出错。6.3 软阈值把谐波幅值压低了怎么办软阈值天然会缩小幅值。如果后续需要精确的谐波幅值不能直接用软阈值重构后的波形作为最终幅值估计。我会先用软阈值方案识别出哪些分量是信号分量然后根据对应的奇异值和左右奇异向量计算出幅值校正。也可以先做一次去噪得到噪声水平估计再用普通最小二乘拟合原始信号中的谐波参数这比直接读去噪波形要准。具体操作为用软阈值去噪识别出谐波频率然后构造正弦字典对原始含噪信号做一次稀疏拟合或线性回归。这样幅值就不是被“压”出来的而是从原始数据里重新估计的。谐波去噪的目的更多是滤除噪声干扰、稳定后续分析而不是直接输出无偏幅值。理解了这一点就不会抱怨软阈值把信号“变小”了。6.4 大矩阵乘法太慢怎么办即使随机SVD避开了完整SVD它还是要做若干次矩阵乘法。如果A本身是数值矩阵且规模很大A*Omega的效率也取决于内存带宽。大数据集上优先用隐式乘法函数避免读取大矩阵如果必须显式构建考虑把A拆成行块分块做乘法。Matlab的矩阵乘法本身优化得不错但内存分配往往才是性能瓶颈。我通常在乘法函数里预分配输出矩阵循环填充避免动态拼接。还有一个小技巧随机投影矩阵Omega用randn生成时如果内存也紧张可以按列一点点生成边生成边乘。但对多数场景kp只有几十列生成的Omega是n乘几十不至于带来压力。真正需要担心的还是A本身。所以“能写算子就不建矩阵”这条原则在大数据场景下是铁律。如果让我把这次实践浓缩成一句话那就是随机SVD给大数据集提供了一条从“算不动”到“算得动”的桥软阈值又给去噪质量加了一道保险。两者配合起来既解决了传统SVD在数据规模面前的无力感也避免了硬阈值截断的重构振铃。最后分享一个我自己的小习惯在处理一批现场数据之前我会先用一段几十万点的调试子集把k、tau、q都跑出大致的合理范围然后用这套参数对全部数据批量执行。因为随机SVD在参数合理时相当稳定不会因为数据段不同而突然失效。整个过程下来你会明显发现所谓“大数据集里的谐波去噪”也没有想象中那么可怕。