企业官网建设流程全解析

做谐波去噪做久了，你一定会遇到一个尴尬的坎：数据量一大，传统的奇异值分解（SVD）直接算不动，或者算一次要等到天荒地老。我当初处理一组 10 万点的振动谐波信号，构造出轨迹矩阵后，精确 SVD 跑了整整三分钟，而后续参数调优又要反复重算，那种感觉就像抱着个大石头走路。后来换成了随机奇异值分解（Randomized SVD）配合软阈值（Soft Thresholding）做去噪，同样规模的数据，整个流程压缩到了几秒，而且去噪效果和精确分解基本没有肉眼可见的差别。这篇文章我就把这条技术路线完整拆开，从原理到 Matlab 代码实现，再到大数据集下的工程化处理，一次性讲透。

这篇文章适合谁？如果你正在做电力谐波分析、机械振动信号处理、水声信号清洗，或者手头有大量含噪周期信号需要批量去噪，同时又对计算资源和实时性有要求，那这套"随机 SVD + 软阈值"的组合拳就是为你准备的。内容不会只停留在给代码层面，我会把每个关键选择的"为什么"也讲清楚——为什么是随机 SVD 而不是精确 SVD，为什么用软阈值而不是硬阈值，为什么轨迹矩阵的维度这样取。读完你不仅能用 Matlab 复现，还能在真实项目里根据数据特征自己调参。

1. 谐波去噪的本质：为什么问题可以归结为低秩近似

1.1 谐波信号的数学结构

谐波信号说白了就是一组固定频率的正弦波叠加。电力系统里的 50Hz 基波加 3、5、7 次谐波，旋转机械里的转频及其倍频成分，水声通信里的窄带载波，这些都是典型的谐波结构。这类信号有个非常重要的共同点：它们在一个足够长的观测窗口内，本质上被少数几个频率参数完全确定。

写成离散形式：

[ x(t) = \sum_{i=1}^{r} A_i \sin(2\pi f_i t + \phi_i) + n(t) ]

其中 ( n(t) ) 是宽带噪声，( r ) 是谐波个数。去噪的目标就是把这 ( r ) 个正弦成分从混合物中干净地拎出来。

如果只是单频正弦波，用一个二阶线性预测或者直接 FFT 就能搞定。但当谐波分量多、频率接近、噪声强的时候，频域方法就容易出现谱泄漏和旁瓣干扰。这时候矩阵方法反而更有优势——因为它把时间序列的整体结构当作一个对象来处理，而不是逐个频率点去猜。

1.2 轨迹矩阵：把一维时间序列变成二维低秩矩阵

把一维时间序列 ( x(1), x(2), ..., x(N) ) 变成一个矩阵，最经典的做法是构造 Hankel 矩阵（也叫轨迹矩阵）。选定一个窗口长度 ( K )，令 ( L = N - K + 1 )，构造：

[ H = \begin{bmatrix} x(1) & x(2) & \cdots & x(L) \ x(2) & x(3) & \cdots & x(L+1) \ \vdots & \vdots & \ddots & \vdots \ x(K) & x(K+1) & \cdots & x(N) \end{bmatrix} ]

这个矩阵的维度是 ( K \times L )。关键性质来了：如果原始信号是 ( r ) 个正弦分量的叠加（无噪声），那么理论上这个 Hankel 矩阵的秩不超过 ( 2r )。这是因为正弦函数可以表示为两个指数函数的线性组合，而这样构造的 Hankel 矩阵能够写成 ( r ) 个秩 2 矩阵的和。

换句话说，谐波信号的多维能量集中在一个极低维的子空间里。加性噪声则均匀地散布在所有奇异方向。这个差异就是低秩去噪的物理基础——把矩阵分解成"低秩成分 + 扰动成分"，低秩部分就是去噪后的信号。

1.3 精确 SVD 的复杂度困境

有了低秩这个先验，直觉上的做法就是对轨迹矩阵做 SVD，截断前 ( 2r ) 个奇异值重构。但问题出在规模上。

假设 ( N = 100000 )，窗口 ( K = 200 )，那么轨迹矩阵是 ( 200 \times 99801 )。精确 SVD 的复杂度约 ( O(KL^2) )，也就是 ( 200 \times 99801^2 )，这个是接近 ( 2 \times 10^{12} ) 量级的运算。虽然现代 Matlab 的 svd 函数用了分块算法和 LAPACK 优化，但这个量级依然需要几分钟甚至更久。而如果数据是批量的——比如 1000 组传感器数据——那这个方案就彻底不现实了。

随机 SVD 的思路非常直接：我们并不需要完整的 ( 200 \times 99801 ) 维分解，只需要前几个大奇异值和对应奇异向量。既然如此，为什么要把所有奇异方向都算一遍？用随机投影把矩阵"压扁"到一个低维子空间，在这个小矩阵上做精确 SVD，然后映射回去——这就是随机 SVD 的核心思想。

2. 随机 SVD 的三步走：投影、正交、小矩阵分解

随机 SVD 的算法流程并不复杂，我直接写出每一个步骤的线性代数含义。

2.1 算法流程拆解

给定矩阵 ( A )（维度 ( m \times n )），目标秩 ( r )，过采样参数 ( p )，幂迭代次数 ( q )：

第一步：随机投影。生成一个 ( n \times (r+p) ) 的高斯随机矩阵 ( \Omega )，计算 ( Y = A\Omega )。这个 ( Y ) 的列是矩阵 ( A ) 在随机方向上的投影，它捕获了 ( A ) 的主要列空间信息。因为随机方向几乎不会与 ( A ) 的低秩子空间正交，所以 ( Y ) 里保留了前 ( r+p ) 个主奇异方向的信息。

第二步：列空间正交基。对 ( Y ) 做 QR 分解，得到 ( Q )（维度 ( m \times (r+p) ) ）。( Q ) 的列向量张成了 ( A ) 主要列空间的近似正交基。

第三步：小矩阵精确分解。令 ( B = Q^T A )，维度是 ( (r+p) \times n )，远小于原始矩阵的规模。对 ( B ) 做精确 SVD：( B = U_B \Sigma V^T )。然后令 ( U = Q U_B )，奇异值矩阵就是 ( \Sigma )。

这里有一个值得注意的细节：常规的公式是 ( Y = (AA^T)^q A\Omega )，其中 ( q ) 是幂迭代次数。直接这样计算涉及两次大矩阵乘法。更聪明的做法是利用结合律：

[ Y = A(A^T(A(\cdots(A^T(A\Omega))\cdots))) ]

先算 ( A\Omega )，再算 ( A^T(\cdot) )，如此交替。这样可以避免显式构造 ( AA^T )，在内存上友好得多。

2.2 三个超参数怎么选

随机 SVD 有 ( r )、( p )、( q ) 三个超参数，选择一个合适参数组合非常重要。

目标秩 ( r )：对谐波去噪而言，( r = 2 \times )（谐波个数）是理论值。如果你不确定谐波个数，可以稍微取大一点。比如估计有 3~5 个谐波，就取 ( r = 12 ) 左右，给奇数个频率分量留出余量。取大了只会多算几个小奇异值，对最终结果影响不大，代价是计算时间增加。

过采样参数 ( p )：通常取 5~10，或者 ( r ) 的 10%。过采样的作用是给随机投影一点"容错空间"，避免因为随机方向与某几个奇异向量恰好正交而丢失信息。这在谱能量衰减不够快的时候尤其重要。我自己习惯 ( p = 10 )，简单又稳。

幂迭代次数 ( q )：这是个容易忽略但很关键的参数。当矩阵的奇异值谱衰减较慢（即后面的奇异值没有迅速趋近于零）时，随机投影捕获前 ( r ) 个主方向的能力会下降。幂迭代的作用是把奇异值差距"拉开"——每做一次 ( Y = (AA^T)^q A\Omega ) 的等价操作，就相当于把谱的衰减速度提升 ( 2q ) 倍。

我建议：如果信号比较干净（SNR 高于 10dB），( q = 1 ) 就够了；如果噪声很强，或者信号频带很宽导致奇异值衰减慢，用 ( q = 2 ) 或 ( q = 3 )。注意 ( q ) 每增加 1，计算量大约翻一倍，所以别贪。

2.3 误差直觉

很多人担心随机 SVD 是"近似方法"，会丢精度。实际上它的误差是有理论界的。对于矩阵 ( A )，取 ( \Omega ) 为高斯随机矩阵，那么：

[ |A - Q Q^T A|2 \approx \sigma{r+1} \cdot (1 + \text{小修正}) ]

其中 ( \sigma_{r+1} ) 是第 ( r+1 ) 个奇异值。也就是说，随机 SVD 的误差主要取决于被截断掉的奇异值有多大，而随机性带来的额外误差可以通过过采样和幂迭代压到极小。

这个误差界与时谐波去噪场景配合得很好——谐波信号的低秩部分能量占绝对主导，噪声对应的小奇异值虽然数量多但单个都小。所以误差的绝对值非常有限。噪声引起的最大奇异值通常远小于谐波对应的奇异值，截断误差自然就不会大。

3. 软阈值操作：为什么在奇异值上做"温柔的截断"

拿到了 SVD 结果，接下来就是去噪的关键动作。天真做法是把小于某个阈值的奇异值直接置零，大于阈值的保留原值，然后重构。这就是硬阈值。但它有个毛病：对于含噪的奇异值，直接置零会产生不连续的突变，重构信号容易出现振铃和毛刺，俗称"伪吉布斯现象"。谐波信号是平滑的周期函数，最怕这种人工突变。

软阈值的做法就温和得多。对每个奇异值 ( \sigma_i )，应用：

[ \sigma_i^{\text{new}} = \text{sign}(\sigma_i) \cdot \max(|\sigma_i| - \tau, 0) ]

因为奇异值都是非负的，sign 部分可以忽略，实际就是 ( \max(\sigma_i - \tau, 0) )。小于阈值的直接抹掉，大于阈值的整体收缩一个 ( \tau )。这保证了奇异值序列的连续性，重构信号更平滑。

我把两者的差异总结成表格：

软阈值唯一的副作用是把保留下来的奇异值也收缩了，相当于对信号整体能量打了折扣。这在去噪里是可接受的，因为噪声对奇异值的贡献本身就偏大，收缩相当于做了折中。如果你关心最终的幅值精度，可以在重构后对去噪信号做一次幅值校准——比如用最小二乘拟合各谐波分量的实际幅值。后文我会给具体做法。

阈值 ( \tau ) 的确定是整个流程里最需要经验的地方。我的做法是两步：

第一步，估计噪声标准差 ( \sigma_n )。用残差矩阵的奇异值，取中位数绝对值（MAD）估计：

[ \sigma_n \approx \frac{\text{median}(|\Sigma_{\text{rejected}}|)}{0.6745} ]

实际操作中更简单：对 Hankel 矩阵的每个元素减去列中位数，然后把所有元素绝对值取中位数。这个估计对异常脉冲稳健，不依赖模型假设。

第二步，设定阈值。

[ \tau = \alpha \cdot \sigma_n \cdot \sqrt{m \cdot n} ]

其中 ( m, n ) 是轨迹矩阵维度，( \alpha ) 是经验系数，一般在 0.5 到 1.5 之间。这个公式是有量纲依据的：高斯随机矩阵的最大奇异值近似在 ( \sigma_n(\sqrt{m} + \sqrt{n}) ) 量级，乘一个安全系数后，能把纯噪声的奇异值全部压掉，又不会伤及有效成分。

如果你用的是 Matlab，可以直接调用wthresh(s, 's', tau)做软阈值。但要提醒一句，这个函数需要 Wavelet Toolbox，不是所有人都有这个工具箱。我一般直接手写一行代码：

s_new = sign(s) .* max(abs(s) - tau, 0);

零依赖，执行速度快，效果完全一致。

4. Matlab 实现：完整的"轨迹矩阵 + 随机SVD + 软阈值"流程

下面给出完整的可运行代码框架。为了可读性，我拆成几步，每一步都加了注释。

4.1 主流程代码

function [x_clean, info] = rsvd_soft_harmonic_denoise(x, fs, n_harmonic, K, opts) % RSVD_SOFT_HARMONIC_DENOISE 基于随机SVD和软阈值的谐波去噪 % 输入： % x - 含噪的时间序列，N x 1 % fs - 采样率（Hz） % n_harmonic- 估计的谐波个数 % K - 轨迹矩阵窗口长度 % opts - 结构体，可选字段：p, q, tau_alpha % 输出： % x_clean - 去噪后的时间序列 % info - 附加信息（奇异值、时间等） N = length(x); % 1. 构造轨迹矩阵 L = N - K + 1; H = zeros(K, L); for i = 1:K H(i, :) = x(i : i+L-1); end % 2. 随机SVD参数 r = 2 * n_harmonic + 2; % 留一点余量 p = 10; q = 2; tau_alpha = 1.0; if nargin >= 5 && isfield(opts, 'p'), p = opts.p; end if nargin >= 5 && isfield(opts, 'q'), q = opts.q; end if nargin >= 5 && isfield(opts, 'tau_alpha'), tau_alpha = opts.tau_alpha; end % 3. 随机SVD tic; [U, S_vals, V] = rsvd(H, r, p, q); info.time_svd = toc; % 4. 噪声估计与阈值 sigma_noise = median(abs(H(:))) / 0.6745; tau = tau_alpha * sigma_noise * sqrt(K * L); % 5. 软阈值收缩奇异值（保留最大的一个，防过度收缩） S_new = sign(S_vals) .* max(abs(S_vals) - tau, 0); % 如果第一个奇异值也被削了，至少保留一定能量 if S_new(1) < tau S_new(1) = S_vals(1) * 0.9; end % 6. 重构矩阵并反变换回时间序列 H_clean = U * diag(S_new) * V'; % 7. 对角平均（Hankel矩阵的反操作） x_clean = hankel_to_series(H_clean, N); info.S_vals = S_vals; info.S_new = S_new; info.sigma_noise = sigma_noise; info.tau = tau; info.time_total = toc; end

4.2 随机 SVD 的具体实现

function [U, S_vals, V] = rsvd(A, r, p, q) % RSVD 随机奇异值分解 % A - m x n 矩阵 % r - 目标秩 % p - 过采样 % q - 幂迭代次数 [m, n] = size(A); ell = r + p; Omega = randn(n, ell); % 使用结合律做幂迭代，避免显式计算 A*A' Y = A * Omega; for i = 1:q Y = A * (A' * Y); end % QR分解获取列空间基 [Q, ~] = qr(Y, 0); % 小矩阵精确分解 B = Q' * A; [U_B, S_B, V] = svd(B, 'econ'); S_vals = diag(S_B); U = Q * U_B; % 截断到前 r 个（严格来说随机SVD给了 ell=r+p 个，我们截到r） U = U(:, 1:r); S_vals = S_vals(1:r); V = V(:, 1:r); end

这里有个容易踩的坑：qr(Y, 0)是经济型 QR，返回的 Q 列数为min(m, ell)。当m远大于ell时没问题。但如果矩阵的m < ell（比如窗口长度小于目标秩加过采样），需要先转置处理，否则 Q 的列数不够。实际问题中窗口长度 K 通常远大于 2 倍谐波个数，所以这个坑不常遇到，但我在写通用工具时会加一个判断：

if m >= ell [Q, ~] = qr(Y, 0); else [Q, ~] = qr(Y' , 0); Q = Q'; end

4.3 Hankel 矩阵反向拼接

去噪后的矩阵H_clean是一个近似 Hankel 矩阵，但因为有截断误差，矩阵的对角线元素并不完全相等。标准的做法是"对角平均"——把每条反对角线上的元素取平均，作为时间序列在该时刻的输出值。

function x = hankel_to_series(H, N) % HANKEL_TO_SERIES 对角平均还原时间序列 [K, L] = size(H); x = zeros(N, 1); count = zeros(N, 1); for i = 1:K for j = 1:L idx = i + j - 1; x(idx) = x(idx) + H(i, j); count(idx) = count(idx) + 1; end end x = x ./ count; end

这个双重循环在数据量大时会有一点慢。我测试过，对于一个 200×99801 的矩阵，单纯对角平均大概耗时 0.5~1 秒左右。如果你要追求极致性能，可以把它改成矩阵运算——用spdiags或者accumarray实现。不过考虑到整个流程已经从几分钟降到了几秒，这里的一秒代价我选择接受，代码清晰更重要。

4.4 窗口长度 K 怎么定

窗口长度 ( K ) 是轨迹矩阵方法里对结果影响最直接的参数。

从秩的角度看，( K ) 至少要大于 ( 2r )，否则矩阵的秩根本容不下所有谐波成分。从频率分辨率的角度看，( K ) 太小时，Hankel 矩阵的"视野"太短，无法区分频率接近的谐波；( K ) 太大时，矩阵规模变大，计算量上升，而且随机 SVD 的优势会被弱化。

我的经验法则是：

• 保证 ( K ) 至少包含最低频率谐波的一个完整周期。
• 如果数据里有频率 ( f_{min} )，则建议 ( K \geq 2 \cdot \text{round}(f_s / f_{min}) )。
• 对于长序列，( K ) 取 200~500 通常是合理的平衡点。

举个例子：采样率 1000Hz，最低谐波频率 50Hz，那么一个周期是 20 个点，K 至少 40；实际我会取 ( K = 200 )，这样能覆盖 10 个周期，频率分辨率和统计稳定性都够了。

5. 实测对比：随机SVD + 软阈值 vs 精确SVD + 硬阈值

光说不练假把式。我构造了一组仿真信号来验证整个流程：

采样率 1000Hz，时长 100 秒（N=100000）。三个谐波：50Hz（幅度1.0）、150Hz（幅度0.5）、250Hz（幅度0.25），相位分别是 0、π/4、π/3。叠加高斯白噪声，信噪比约 5dB。

窗口 ( K = 200 )，目标秩 ( r = 8 )，过采样 ( p = 10 )，幂迭代 ( q = 2 )。

下面是三种方案的结果对比：

可以明显看到，随机 SVD 的计算时间比精确 SVD 低两个数量级，而去噪质量的损失几乎可以忽略。软阈值确实比硬阈值好，输出 SNR 高了 1.3dB，而且重构信号的波形更光滑。

我还特意检查了重构信号末尾段和开头段，硬阈值方案在信号幅度跃变处出现了轻微振荡，软阈值方案则非常干净。这个差异在肉眼观察时不易察觉，但在后续做谐波幅值精确提取时，会直接影响测量精度。

5.1 噪声强度变化时的表现

改变噪声水平，观察随机SVD+软阈值的输出 SNR：

不同输入 SNR 下，输出 SNR 的提升稳定在 14dB 上下。这个结果说明软阈值收缩能稳定剥离约 95% 的噪声功率。当然，如果噪声模型换成非高斯（比如脉冲噪声），这个提升幅度会下降，"健壮性"就体现在这里——用 MAD 估计噪声时，个别大的异常值对中位数影响不大，所以阈值不会因为几个离群点而乱跳。

5.2 超参数敏感性测试

我做了个简单的网格搜索：r 从 6 到 20，p 从 5 到 15，q 从 1 到 3。输出 SNR 的变化范围最大只有 0.8dB。说明这套方案的性能对于参数选择不敏感，这对实际应用很重要——因为大多数场景下，你不可能提前精确知道谐波个数。

唯一需要注意的是 r 取得太小。如果真实谐波数是 3（需要秩 6），但你把 r 设为 4，那么前两个谐波会被保住，第三个谐波会丢失一部分能量。因为它的奇异值和噪声混在一起被软阈值压掉了。这是一个不可逆的信息损失，比参数取大严重得多。所以我的原则是：r 宁可大不可小，用 p 和 q 来控制计算精度。

6. 大数据集的工程化处理：批量信号去噪

实际项目里很少只处理一条时间序列。传感器阵列、多通道振动数据、批量离线文件，都是成百上千条信号堆在一起。逐条调用上面的函数虽然可行，但效率明显偏低。我总结了三个工程化技巧。

6.1 批量数据的并行计算

Matlab 的parfor可以直接套在我的主函数外面。由于每条时间序列的 Hankel 矩阵构造和 SVD 分解彼此独立，这个场景天然适合并行。

% 假设 X 是 N x C 的矩阵，C 是通道数 parfor c = 1:C X_clean(:, c) = rsvd_soft_harmonic_denoise(X(:, c), fs, n_harmonic, K, opts); end

不过要注意，并行池的启动和进程间数据传递也会带来开销。如果单条信号不大，建议换个思路：把多条信号叠成一个三维数组，一次批量构造 Hankel 块对角矩阵，用一次随机 SVD 同时分解。但这个方案实现复杂度高，除非数据量大到单条处理真的不可接受，否则我不建议这样做。

6.2 流式处理长数据

还有另一种情形：单条时间序列特别长，比如连续监测一个小时的高频振动数据，N 可能到了百万甚至千万量级。这时构造完整轨迹矩阵的内存开销已经很可观，一次全部读入不现实。

我的做法是分段处理。将原始信号切成有重叠的片段，每段长度 2~5 万点，分别去噪后再用重叠相加（OLA）拼接。重叠率取 50%，两端各加窗（推荐汉宁窗）抑制边界效应。

具体分段参数就根据你的实际数据来定。核心逻辑是：让每段内的谐波频率保持相对稳定，同时保证段长足够覆盖多个周期。这个方法简单可靠，唯一要注意的是段间拼接处可能出现相位不连续，但 50% 重叠的汉宁窗加法可以很好地解决。

6.3 内存控制

随机 SVD 虽然计算快，但如果构造出完整的 Hankel 矩阵再传进函数，内存峰值依然可能很高。一个 200×99801 的 double 矩阵大约是 1.6GB。

一个实用技巧是分块构造 Y。先初始化Y = zeros(K, ell)，然后循环计算 ( A\Omega ) 的每一部分，而不需要一次性实例化整个 H 矩阵。

Y = zeros(K, ell); Omega = randn(L, ell); for i = 1:ell % 用快速卷积或部分矩阵乘法实现 A*Omega(:,i) Y(:, i) = partial_hankel_mult(x, Omega(:, i), K); end

partial_hankel_mult可以利用 FFT 加速 Hankel 矩阵的乘法——Hankel 矩阵乘以向量本质是一个卷积。这是一个更高级的优化手段，有兴趣的读者可以自己研究。我在实际项目中就靠这一手，把 500 万点数据的去噪完整跑进了 20 秒以内。

7. 实操中容易踩的坑和我的调参经验

这部分是我最想跟读者分享的内容。理论再漂亮，落地时总有几个细节会坑到你。

7.1 软阈值过度收缩的问题

我在最初测试时发现，当噪声很强（SNR 低于 0dB），通过 MAD 估计出的噪声标准差会偏大，导致阈值 τ 过大，连第一个奇异值也被削掉了一大半。结果是去噪后的信号虽然干净，但谐波幅值严重缩水，和真实值差了 20% 以上。

解决方法是给第一奇异值一个保护机制。前文代码里那个if S_new(1) < tau的判断就是从实际经验里来的。或者更精细一点：如果前 ( 2r ) 个奇异值中有超过一半都被压到零，就适当降低 ( \alpha ) 重跑一次。这个简单策略能让输出 SNR 额外提升 1dB 左右。

7.2 随机种子和可复现性

随机 SVD 里用到了randn。如果你不固定随机种子，同样的代码跑两次结果会有细微差异。这在开发测试阶段会让人抓狂——你明明什么参数都没改，第三次运行的结果和前两次略有不同，你会怀疑是自己代码出 bug 了。

我的建议是：在调用随机 SVD 之前固定种子：

rng(42);

注意rng最好只在测试和调试时固定。正式的大规模处理中，随机种子的影响会随矩阵规模增大而迅速减小，不固定反而能避免某些极端随机矩阵带来的不利情况。

7.3 谐波幅值如何精确恢复

软阈值收缩会整体压低奇异值，导致重构信号的谐波幅值偏小。如果你要做的是谐波检测，而不仅仅是波形清洗，那么后续的幅值校准步骤不能省。

最简单有效的方法：对去噪后的信号做 FFT，提取各峰值频率处的幅值 ( A_{meas} )，再和原始含噪信号在同一频率处的 FFT 幅值做对比。因为谐波成分在窄带内的 SNR 通常远高于宽带平均 SNR，原始信号的该频率幅值反而是可信的。

也可以反过来验证去噪效果：如果去噪后的 FFT 幅值和去噪前差得太多（比如超过 10%），说明你的软阈值收缩过度了，需要调低 ( \alpha )。我一般用这个方法作为快速调参手段，几秒钟就能判断阈值设得合不合理。

7.4 与 FFT 谐波分析的衔接

有人会问：既然最后还是用 FFT 提幅值，那干嘛还要做 SVD 去噪？我的体会是，这两者解决的其实是不同层面的问题。FFT 适合在信噪比尚可的情况下快速提取频谱峰，但当噪声很大时，频谱泄漏和旁瓣抬升会让小谐波被噪声淹没。SVD 去噪后，噪声底被压低，谐波的频谱峰变得尖锐且孤立，FFT 提取的幅值精度自然更高。

另外一个实用衔接方式是：先用 SVD 去噪得到干净的时域波形，再对波形分段做 FFT 看频谱的时变性。你可以清晰地观察各次谐波幅值随时间的变化，这是直接用含噪信号做 FFT 很难做到的。

7.5 一个关于 q 的细节

幂迭代 q 有一个容易忽略的副作用：它会放大主要成分的能量差距。对于谐波去噪，这意味着小谐波（比如 5 倍频之后的高次谐波）对应的奇异值可能被过度压缩。如果你的数据里存在幅度很小但真实存在的谐波成分，建议 q 不要超过 2。在类似场景下，我给客户的推荐配置一直是 q=1 起步，根据结果再决定是否增加。低信噪比时 q=2，很少会用到 q=3，因为收益太小，计算代价却不小。

这套"随机SVD + 软阈值"的谐波去噪流程我已经在多个项目里落地用过。最初吸引我的是它把计算时间从分钟级压到秒级，真正做久了之后，反而觉得它最可贵的地方是稳定——参数不敏感、对噪声模型不敏感、对数据规模不敏感。你不用担心某一天换了一批数据就要重新调一整天参数。如果你手头也有大批量含噪谐波数据需要清洗，我建议你直接把这套代码拿去跑一遍，把第一版本的参数设成我在文中给的默认值，然后观察一下输出。大概率你会发现，那些之前被噪声盖住的细节，现在能看清了。