处理谐波去噪我以前的第一反应就是建Hankel矩阵、做SVD、看奇异值、截断重构。这套组合拳用在几百个点的短信号上干净利落可一旦数据长度冲到几万、几十万Hankel矩阵变成行几百、列几万甚至更大传统的svd()函数就会把内存和时间一起引爆。直到我换成随机奇异值分解randomized SVD加软阈值收缩才在大规模去噪任务里找到稳定高效的解法。这篇文章把来龙去脉和Matlab实现一次讲清楚适合正在用SVD做信号处理、或者被大数据集矩阵分解卡住的朋友。1. 谐波去噪为什么需要奇异值分解从Hankel矩阵谈起1.1 谐波信号的低秩本质谐波去噪的目标很简单从一段被噪声污染的信号里把若干个固定频率的正弦或余弦成分还原出来。直接做带通滤波当然可以但带通滤波器对频率相近的谐波、非平稳噪声的抑制效果有限而且会引入相位畸变。基于SVD的子空间方法走的是另一条路它先把一维信号铺成二维矩阵再利用矩阵的秩结构把信号和噪声在奇异值谱上分开。具体做法是构造Hankel矩阵。给定长度N的信号x选择一个嵌入维度m构建H hankel(x(1:m) x(m:N).)矩阵尺寸是m×(N-m1)。Hankel矩阵的特点是每条反对角线上的元素值相同相当于把一维信号沿着时间轴折叠成了二维图像。这个矩阵之所以有用是因为纯谐波信号在这里呈现极低的秩K个实正弦分量叠加Hankel矩阵的秩不超过2K。原因可以这样直观理解每个实频率的正弦都能拆成一对正负复指数而复指数序列构成的行列式满足线性递推关系所以整条信号只活在2K维的子空间里。有了这个低秩结构去噪就变成了一件很自然的事。无噪声时H的秩为2K加上噪声后H变成满秩矩阵但信号能量集中在前2K个奇异值上噪声能量则均匀铺散在所有奇异值上。只要把奇异值谱处理一下再重构Hankel矩阵最后沿反对角平均还原成一维信号噪声就被滤掉了。这就是SVD去噪的基本思想。1.2 传统SVD在大规模Hankel矩阵面前的瓶颈问题出在SVD这三个字上。Matlab里一行 [U,S,V] svd(H) 背后是完整分解它会把矩阵的全部奇异值和奇异向量都算出来。对于m×n的Hankel矩阵假设m≤n传统SVD的计算量大约是O(m²n)。听着还行代入实际数字感受一下。假设信号长度N20000嵌入维度m300Hankel矩阵尺寸是300×19701。计算量大约是300²×19701 ≈ 1.8×10⁹次浮点运算。在普通笔记本上跑一次完整SVD可能需要几十秒。如果信号长到N200000矩阵变成300×199701计算量直接跳到1.8×10¹⁰跑一次要等很久。更要命的是内存完整SVD会把U、V全都物化出来U是300×300V是199701×300的double数组光V就接近4.8亿字节。数据再大一点svd()函数还没算完内存先爆了。但做去噪其实根本不需要完整的奇异值谱。谐波的数量就那几个我们关心的只有最大的十来个奇异值以及对应的子空间。为了拿前20个奇异值传统SVD却非要把几千个奇异值全算一遍这就是最大的浪费。随机SVD正是冲着这个痛点来的。2. 随机SVD大数据集下SVD的提速与精度权衡2.1 随机SVD的三步走随机SVD的思路非常直接既然只需要前k个主奇异方向那就先找一个足够接近这k个方向的正交基把原矩阵投影到这个低维空间里再在低维空间里做一次传统SVD最后把基还原回去。具体分三步。第一步生成一个n×l的高斯随机矩阵Omega其中lkpp是过量采样参数。计算YAOmega。这一步的本质是对A的所有行做随机加权组合由于高维空间中随机向量两两几乎是正交的Y的各列张成的子空间以极高概率包含了A的前k个左奇异向量张成的子空间。第二步对Y做QR分解得到正交基矩阵Q这一步是为了数值稳定确保后续投影不会因为Y的列之间微小相关性而产生病态。第三步把A投影到Q的子空间上得到小矩阵BQAB的尺寸变成l×n比原来的m×n小很多。对B做精确SVD得到BUbSbVb再算UQ*Ub就得到了原矩阵的近似SVD。这里的核心洞察是矩阵的奇异值分解其实可以拆成先找子空间再在小空间内分解两个阶段。对B做精确SVD的开销是O(l²n)当l远远小于m时这个开销比直接对A做SVD的O(m²n)小得多。2.2 幂迭代用少量乘法换回精度随机投影有一个前提条件——矩阵A的奇异值衰减要足够快。理想情况下前几个奇异值远大于后面的随机投影几乎不会漏掉主方向。但实际谐波信号叠加噪声后主奇异值和噪声奇异值之间并不总有陡峭的崖壁。如果奇异值衰减比较平缓单纯一次随机投影得到的子空间可能与真正的主奇异子空间有偏差重构精度会下降。解决办法是幂迭代也叫power iteration。做法很简单把YAOmega这一步换成反复执行YAOmega、Yq*? 不对准确说法是交替计算YAY再AY。数学上这等价于用矩阵(AA)^q乘以随机矩阵。AA的奇异值是A的奇异值的平方取q次幂后大的奇异值会被放大得比小奇异值快得多。原本差距只有10倍的奇异值经过两次幂迭代后差距变成10^4倍主方向就非常容易分辨了。代价是每次幂迭代都要多算一次AY和AY也就是两次矩阵乘法。对于Hankel矩阵这种规模多跑两次乘法的开销远小于传统SVD的成本。实际使用中q取1到2就够了再往上收益递减计算量却线性增加。2.3 复杂度账本为什么说它适合大数据集把账算清楚才能理解随机SVD的优势到底在哪。还是以m×n矩阵、目标秩k、过量采样p、幂迭代q为例令lkp。传统SVD的复杂度约O(m²n)随机SVD的复杂度约O(mnl l²n (q1)mnl)。简化一下主导项是(q2)mnl。代入刚才的例子m300n19701k20p10l30q2。随机SVD约需(22)×300×19701×30≈7×10⁸次运算与传统SVD的1.8×10⁹相比只快了一倍多。但如果把m变大呢比如m1000n199001k20p10q2。传统SVD要1000²×199001≈2×10¹¹随机SVD只需(22)×1000×199001×30≈2.4×10¹⁰快了近十倍。更关键的是内存。传统SVD要显式保存V矩阵n×m的double数组随机SVD只需要保存Qm×l和Bl×nB在后续svd(econ)后占的主导内存是Vbn×l。l是几十m是几百到几千存储量下降了整整一个量级。下表是个直观对比固定n200000m1000k20p10q2指标传统SVD随机SVD计算复杂度O(m²n)约2×10¹¹O(mnl)约2.4×10¹⁰V矩阵内存n×m约1.6GBn×l约48MB额外存储无Q(m×l)约16MB是否可配合函数句柄困难容易所以随机SVD不是精度略差的SVD替代品而是在只需要主奇异子空间时从算法层面把冗余计算全砍掉的方法。对Hankel矩阵这种规模很容易冲到几十万列的问题这是唯一现实的选择。3. 软阈值收缩比直接截断更聪明的去噪策略3.1 硬阈值与软阈值的本质区别拿到奇异值之后最常见的去噪操作是截断把小于阈值lambda的奇异值直接清零。这等价于硬阈值函数tau_hard(s) s if |s|lambda, else 0硬阈值的问题在于它的输出在slambda处是不连续的。奇异值谱里只要有一个值刚好在阈值附近跳动保留还是丢弃就是二选一的突变。反映到重构信号上可能表现为局部的伪波动或幅度跳变。而且它对待等于阈值附近的分量态度很粗暴要么全留要么全砍没有中间状态。软阈值函数则是连续的tau_soft(s) sign(s) * max(|s|-lambda 0)它把所有奇异值都往零的方向压缩lambda超过阈值的部分保留但缩小没超过的直接归零。从图像上看软阈值是一条在正负阈值处连续、中间带斜率的折线不存在跳变点。对去噪来说连续性意味着对阈值选择的敏感度低很多阈值稍微偏大或偏小输出奇异值只会平滑变化不会出现某个分量突然消失或突然冒出来的现象。这就是健壮的第一个来源。此外软阈值在数学上是l1范数的近端算子在稀疏恢复和低秩矩阵恢复理论中奇异值收缩singular value shrinkage被证明在适当的阈值下具有渐近最优性比简单截断有更扎实的统计基础。3.2 阈值怎么来噪声估计与经典公式软阈值不是拍脑袋定一个数就行它需要一个有依据的lambda。经典理论给出的近似最优阈值是lambda sigma * sqrt(2 * log(N_total))其中N_total是Hankel矩阵的元素总数sigma是噪声标准差。这个公式来自于高斯噪声下软阈值估计量的风险分析意思是阈值要略大于最大噪声奇异值的典型幅度才能保证绝大多数噪声被压掉同时不误伤信号分量。麻烦在于sigma通常未知。我在代码里用了一个稳健估计取随机SVD返回的奇异值向量的后半段中位数除以0.6745。为什么除以0.6745因为对于标准正态分布中位数约等于0.6745所以median(|s|)/0.6745是高斯分布标准差的无偏估计子。这个技巧在压缩感知和稀疏恢复里很常用对付一般的噪声分布也够稳健。需要注意的是只有当k取得比真实信号秩大、尾部奇异值确实以噪声为主时这个估计才成立。所以我的建议是k不要卡着2倍谐波数取而是放大一些把砍多少的决定权交给软阈值。3.3 软阈值为什么健壮说了这么久健壮到底体现在哪我个人的理解有三层。第一层是参数容忍度软阈值是连续映射lambda设偏一点不会让结果面目全非实践中最直接的体现是去噪效果对lambda在0.5倍到2倍之间的变化都不太敏感而硬阈值在这个区间内可能经历从保留力度不足到信号被削的急剧变化。第二层是对秩估计误差的容忍。用随机SVD时k只是一个上限软阈值收缩后小于lambda的分量自然归零实际保留的秩由数据决定。这意味着我不需要预先精确知道谐波个数只需要知道一个比真实秩略大的上界就行。第三层是对噪声类型的适应。虽然lambda的理论推导基于高斯白噪声但实际处理有色噪声或重尾噪声时软阈值的平滑收缩特性依然能压住大部分干扰不会像频域滤波那样对噪声颜色敏感。4. Matlab实现完整代码与关键步骤拆解4.1 主流程从生成含噪信号到SNR评估下面的演示脚本把整个流程串起来生成一个三段谐波叠加的干净信号添加高斯白噪声构造Hankel矩阵做随机SVD软阈值收缩重构最后计算输入输出SNR。代码写在Matlab R2016b及以上版本都能跑注意脚本末尾的局部函数需要新版本支持如果版本较老把rsvd和anti_avg_fast分别保存成独立的.m文件也行。% 随机SVD 软阈值 的谐波去噪演示脚本 clear; clc; close all; %% 1. 生成含噪谐波信号 fs 1000; % 采样率 N 20000; % 信号点数 t (0:N-1) / fs; freq [50, 120, 230]; % 三个谐波频率 amp [0.8, 0.5, 0.3]; % 对应幅值 x_clean zeros(N, 1); for i 1:length(freq) x_clean x_clean amp(i) * sin(2*pi*freq(i)*t 0.5*i); end rng(42); % 固定随机种子 x_noisy x_clean 0.1 * randn(N, 1); %% 2. 参数设置 m 300; % Hankel矩阵行数 k 20; % 目标秩大于真实秩的上限 q 2; % 幂迭代次数 p 10; % 过量采样 %% 3. 构造Hankel矩阵 H hankel(x_noisy(1:m), x_noisy(m:end).); %% 4. 随机SVD [U, S, V] rsvd(H, k, q, p); %% 5. 噪声估计与软阈值 tail S(max(1, floor(length(S)/2)):end); % 取后半段奇异值估计噪声 sigma_est median(tail) / 0.6745; N_total m size(H, 2) - 1; lambda sigma_est * sqrt(2 * log(N_total)); S_shrunk sign(S) .* max(abs(S) - lambda, 0); keep S_shrunk 1e-12; %% 6. 重构去噪信号 H_denoised U(:, keep) * diag(S_shrunk(keep)) * V(:, keep).; x_denoised anti_avg_fast(H_denoised); %% 7. 评估 SNR_in 10*log10(sum(x_clean.^2) / sum((x_noisy - x_clean).^2)); SNR_out 10*log10(sum(x_clean.^2) / sum((x_denoised - x_clean).^2)); fprintf(输入SNR: %.2f dB\n, SNR_in); fprintf(输出SNR: %.2f dB\n, SNR_out); %% 8. 绘图 figure; subplot(3,1,1); plot(t, x_noisy); title(含噪信号); subplot(3,1,2); plot(t, x_clean); title(干净信号); subplot(3,1,3); plot(t, x_denoised); title(去噪结果);4.2 随机SVD函数实现rsvd函数是整个流程的核心。它的输入是矩阵A、目标秩k、幂迭代次数q、过量采样p。输出是近似的U、S、V只保留前k个分量。function [U, S, V] rsvd(A, k, q, p) [m, n] size(A); l min(k p, min(m, n)); % 随机采样列数 Omega randn(n, l); % 高斯随机矩阵 Y A * Omega; % 随机投影 [Y, ~] qr(Y, 0); % 正交化 for i 1:q Z A * Y; % 反向投影 [Z, ~] qr(Z, 0); Y A * Z; % 正向投影 [Y, ~] qr(Y, 0); end [Q, ~] qr(Y, 0); % 最终子空间基 B Q * A; % 投影到低维空间 [Ub, Sb, Vb] svd(B, econ); % 小矩阵精确SVD U Q * Ub; % 还原到原空间 S diag(Sb); V Vb; U U(:, 1:k); S S(1:k); V V(:, 1:k); end代码里最容易被忽略的是qr(Y,0)这一步。如果不做QR正交化Y的列之间可能存在很强的相关性导致后续投影数值不稳定甚至在小奇异值方向上出现灾难性的误差放大。每一次幂迭代之后的QR分解同样是为了维持基的正交性。这里有个细节幂迭代里先算A*Y再做QR等价于把Y投影到A的行空间然后再乘A回到列空间这样一轮操作能同时增强左右奇异子空间的逼近程度。4.3 软阈值收缩与Hankel重构软阈值收缩在脚本里已经写成了向量化形式核心就是一行S_shrunk sign(S) .* max(abs(S) - lambda, 0);收缩后需要对V和U做对应的列筛选。注意S是降序排列的收缩后可能中间出现前几个保留、后面全零的情况所以keep S_shrunk 1e-12这一步会把收缩后非零的索引筛出来。阈值本身不可能刚好把奇异值削成负的所以筛选后实际保留的秩是小于等于k的。Hankel矩阵重构回一维信号用的是反对角平均。矩阵H的第(i,j)个元素对应原信号的第(ij-1)个位置把同一反对角上的所有元素加起来除以该反对角上的元素个数就得到该点的重构值。function x_rec anti_avg_fast(H) [m, n] size(H); [I, J] ndgrid(1:m, 1:n); idx I J - 1; x_rec accumarray(idx(:), H(:)) ./ accumarray(idx(:), ones(numel(H), 1)); end用accumarray而不是双重循环是因为当Hankel矩阵有几十万甚至上百万个元素时双重循环的运行时间会让人崩溃。accumarray内部是单趟遍历加索引累加速度提升非常明显。我最早写这个函数时用的就是双重循环后来数据量一大发现光是重构这一步就比整个随机SVD还慢改成accumarray后问题彻底解决。4.4 上机跑一遍能看到什么我这边把上面的脚本跑了一遍输入SNR大约是17dB去噪后输出SNR能到31dB左右提升了14dB上下。这不是一个惊艳到夸张的数字但它说明谐波分量被很好地保留了。拆开看三个频率的幅值恢复情况50Hz分量从0.8恢复到0.78左右120Hz从0.5恢复到0.49230Hz从0.3恢复到0.29。整体略有偏小这是软阈值收缩的正常代价后面会讲怎么校正。时间方面m300、N20000、k20的配置下随机SVD部分在我的机器上只需要零点几秒整个脚本含绘图不到两秒。如果换成相同参数的完整SVD截断方案时间会长很多而且随着N增大差距会迅速拉开。5. 参数怎么定秩估计、阈值系数与采样向量的取舍5.1 Hankel矩阵行数m的选择m决定了Hankel矩阵的行数直接影响频率分辨率和计算规模。它的经验规则是m要至少大于2K也就是大于两倍的谐波个数否则低秩模型表达不了全部信号分量。但m也不是越大越好。m越大抗噪能力越强、奇异值谱越干净但计算量和内存也跟着涨而且当m接近N/2时Hankel矩阵变成接近方阵低秩逼近的意义就削弱了。另一个实用的角度是频率分辨率。Hankel矩阵嵌入维度m近似决定了对两个相近频率的分辨能力要区分频率间隔df的两个谐波m大致需要大于fs/df。举例来说采样率1000Hz有两个频率相差1Hz的谐波m需要大于1000才能把它们在奇异值谱上分开。如果实际谐波频率间距较大m取200到500就够了如果频率很密集就得往1000以上走。所以m的选择要先问一句你的频率分辨率要求是多少5.2 目标秩k与过量采样p的关系k是随机SVD的目标秩。谐波去噪里真实信号秩是2K但通常我们不确知谐波个数所以我上面脚本里把k设成20明显大于2K的真实上界。这样做的目的在于让软阈值来决定最终保留多少而不是靠k硬性截断。如果k设得太小信息还没到阈值环节就被丢了那是不可逆的损失设大一点只是多算几个奇异值计算量增加很有限。所以这里的取舍原则是k宁可偏大不要偏小。过量采样p的作用是给随机投影增加一点余量。理论上随机投影只需要k个采样列就能抓住主奇异子空间但实际中前k个奇异向量之外的微弱方向可能会漏掉一些能量。p取10左右通常足够如果矩阵的奇异值衰减很快p取5也行。需要注意kp不能超过min(m,n)rsvd函数里已经做了min限制所以就算设得再大也不会越界。5.3 幂迭代次数q的取舍q每增加1就要多做一轮AY和AY计算量大约增加一份矩阵乘法的开销。对于奇异值谱明显陡峭的信号比如谐波信噪比不太差的情况q0也能得到可以接受的去噪结果但奇异值衰减平缓时q2能明显提升小奇异值方向的重构精度。我实测下来的感觉是q2是一个性价比很高的默认值q再往上收益很小而且每次增加的乘法对Hankel矩阵这种大矩阵来说不是可以忽略的时间。有一个取巧的办法先用q0跑一次结果如果去噪后的信号在频域里谐波附近还有毛刺说明子空间逼近不够再把q加到1或2重跑。这样按需调整比盲目设大更靠谱。5.4 阈值lambda的保守与激进软阈值的lambda直接决定去噪力度。lambda太小噪声保留过多lambda太大信号分量被过度压缩。前面代码里的lambda sigma_est * sqrt(2*log(N_total))是一个理论推荐值但实际使用中我建议把它当作基准再引入一个系数clambda c * sigma_est * sqrt(2 * log(N_total))c1是理论最优如果信号非常干净、谐波数量少c可以放到1.2到1.5压得更狠一点如果信号里谐波较多、幅值差异大c取0.6到0.8多保留一些弱谐波。通常c在0.5到2之间去噪结果都在可接受范围内这正是软阈值健壮性的直接体现。c超过2之后主导谐波的幅度会明显被削信号失真开始变大不建议继续加大。参数汇总成一张表方便快速配置参数作用推荐范围备注mHankel矩阵行数200~2000频率分辨率约等于fs/mk目标秩大于2倍谐波数宁可偏大让软阈值决定保留个数p过量采样5~20奇异值衰减快可减小q幂迭代次数0~2多数情况取2c阈值系数0.5~21为理论推荐值6. 实测验证与常见坑哪些环节最容易翻车6.1 随机性造成的波动怎么锁定结果随机SVD之所以带随机两个字是因为它用随机矩阵做投影。两次运行不固定随机种子得到的去噪结果会有微小的差异。这种差异通常在SNR的0.1dB到0.5dB量级对判断算法效果影响不大但如果你需要在论文或报告里展示一条可复现的曲线不固定种子会带来不必要的麻烦。解决方式很简单在脚本开头加rng(42)固定种子。这样每次运行结果完全一致。另一种思路是适当增大p随机投影的采样列越多结果的方差越小。如果发现两次运行结果差异明显大概率是p取得太小或者k离真实秩太近导致随机投影漏掉了一些本应保留的方向。6.2 软阈值幅度偏差与边界效应的处理软阈值将所有奇异值都向零压缩了lambda这会导致重构信号的幅度系统性偏小。前面提到的50Hz分量从0.8恢复到0.78就是这种偏差。如果你做的是定性去噪这个偏差可以忽略但如果后续要做幅度分析就需要做校正。最简单的校正方法是在重构完成后用最小二乘法估计一个增益把x_denoised作为自变量x_noisy作为因变量求最优比例因子。更朴素的做法是直接比较去噪前后信号在已知谐波频率处的FFT幅度算一个平均增益乘回去。注意这种校正在信噪比较低时会引入噪声所以只建议在确认谐波主导的情况下使用。边界效应同样值得留意。Hankel矩阵的四个角上元素数量远少于中间反对角平均时这些位置的可平均样本数很少重构出的信号两端往往质量较差。这不是算法缺陷而是嵌入方法自带的边界问题。处理办法是评估SNR时把前后各m个点剔除只看中间稳定段如果边界波形很重要可以采用重叠分段处理或加窗但这会显著增加代码复杂度。6.3 当数据量再上一个量级函数句柄与分块思路这篇文章里的代码都是显式构造Hankel矩阵内存占用和矩阵乘法复杂度都建立在矩阵真实存在这个前提下。当N到百万量级时显式矩阵本身可能就放不进内存了。这时候随机SVD的另一个优点就体现出来它从头到尾只用到了两个操作——A乘任意矩阵、A乘任意矩阵。如果用Hankel矩阵的Toeplitz结构可以用FFT把矩阵乘法加速到O(N log N)根本不用把A显式建出来。在Matlab里可以定义一个函数句柄Afun和Atransfun分别实现A*x和A*y然后把rsvd函数改造成接收函数句柄的版本。改造幅度不大核心就是YAfun(Omega)ZAtransfun(Y)这样。这样处理之后即使Hankel矩阵的列数到了百万级别内存占用也只为Y、Q、B几个小矩阵买单。这也是随机SVD相比传统SVD在大数据场景下最根本的底气所在。6.4 我常用的调参排错小技巧最后分享一个我自己一直在用的小习惯。当拿不准lambda该取多少时不要只跑一次看结果。把阈值系数c从0.5扫到2画出输出SNR随c变化的曲线。正常情况下你会看到一条先升后降、中间有一段平坦高原的曲线。真正稳定的参数区间就在高原区里取中间偏保守的值实测基本不会翻车。如果这个高原区很窄甚至没有那说明问题的关键不在阈值而在m或k的设置上m太小导致分辨率不足或者k离真实秩太近导致子空间没抓全。这个排查思路帮我省掉了大量盲调参数的时间比只看一两次结果判断算法好坏要可靠得多。