企业官网建设流程全解析

做信号处理或者设备状态监测的朋友，大概率遇到过这种尴尬：手里有一段几十万甚至上百万个采样点的长录波数据，里面明明藏着规律性的谐波成分，却被噪声糊成一团。想用SVD（奇异值分解）做子空间分离，矩阵一大，MATLAB的svd()直接卡死甚至内存报警；改用FFT带通滤波，谐波附近那一圈泄漏噪声又滤不干净。我这两年做电能质量分析和机械设备振动数据处理，反复跟这类问题打交道，最后稳定使用的方案就是标题里这套"随机奇异值分解（Randomized SVD）+ 软阈值"的组合。这篇文章把这套方法的原理、Matlab实现、参数调优和坑一次性讲透，代码可以直接拿去改。

1. 先搞清楚：这个方案到底在解决什么问题

1.1 大数据集的谐波去噪为什么棘手

谐波去噪本质上是一个"从含噪观测中把周期性成分分离出来"的问题。常见做法有三大流派：经典滤波、变分/稀疏方法、子空间方法。前两种在数据量小的时候很好用，但数据一长就露怯。滤波的难点在于谐波和噪声在频谱上经常重叠，带通设计得再小心，总会在频带边缘留下残留；变分方法需要调正则参数，对非平稳噪声适应性差。

子空间方法，尤其是基于SVD的奇异谱分析（Singular Spectrum Analysis, SSA），思路完全不同：把一维信号嵌入成Hankel矩阵（也叫轨迹矩阵），对这个矩阵做SVD，让信号对应的大奇异值和噪声对应的小奇异值自然分开。这个思路漂亮，但工程落地有个致命瓶颈——计算复杂度。对一个m×n的矩阵做完整SVD，复杂度是O(mn·min(m,n))，再加上Hankel矩阵动不动就是几千行乘几万列，传统SVD在大数据集面前基本不可用。

我这里说的"大数据集"，不是指多少TB的数据仓库，而是指信号处理里的长序列：比如电力系统录波，采样率12.8kHz，录10秒就是12.8万个点；机械振动监测连续采一天，上千万个点也不稀奇。这种规模下，你不可能把整个矩阵扔进svd()函数。方案就是换用随机化算法计算低秩逼近，避开全量SVD。

1.2 为什么是"随机SVD"和"软阈值"的组合

随机SVD（Randomized SVD）本质上是给传统SVD配了一个"降维前置"：先用随机投影把大矩阵压缩成一个小矩阵，再在小矩阵上做精确SVD，最后映射回去。这个思想的数学基础是Johnson-Lindenstrauss引理——随机投影能以高概率保持向量之间的相对距离和夹角。换句话说，只要投影方向足够多，信息结构基本不丢。

光有随机SVD还不够。经典SSA是把奇异值谱一刀切成"信号段"和"噪声段"，也就是硬阈值。硬阈值的问题在于：真实数据里噪声奇异值不会干净地归零，而是呈一条缓降的斜坡；谐波幅值若有波动，还会出现某个谐波的奇异值被误切到噪声段的情况。软阈值（Soft Thresholding）的思路是"收缩"而不是"切断"：所有奇异值统一减去一个阈值，非负部分保留。这样大幅值信号只损失很小比例，小奇异值噪声被压到接近零，整个分解过程连续稳定，对噪声水平的估计误差也不太敏感。

1.3 这套方案适合哪些场景

我实际用过并且觉得值得推荐的应用包括：电网电压/电流波形谐波分析、旋转机械振动信号的故障特征提取、结构健康监测的模态参数识别、以及任何需要从超长信号里提取周期成分的离线分析任务。适合的人群是：有Matlab基础、手头有长序列数据、觉得传统SSA慢到没法用、又不想折腾Python那套依赖环境的同学。下面所有内容围绕Matlab实现展开，核心算法流程在任何语言里都能平移。

2. 核心原理拆解：随机投影、奇异值收缩、轨迹矩阵

2.1 随机SVD是怎么把大矩阵变小的

用一句话概括随机SVD：先用随机矩阵捕获大矩阵最主要的行/列空间信息，把大矩阵投影成小矩阵，再在小矩阵上做精确SVD。完整算法分四步：

假设要对矩阵A（m行n列）求前k个主导奇异值，过采样参数为p，幂迭代次数为q：

1. 生成一个n×(k+p)的随机矩阵Ω（通常用标准正态分布随机数）；
2. 计算Y = AΩ，得到一个m×(k+p)的矩阵，这个矩阵的列向量张成的子空间已经非常接近A的前k个左奇异向量张成的子空间；
3. 为了提升精度，做q次幂迭代：Y = A(AᵀY)，每做一次，奇异值谱的高频部分被进一步压缩，低秩结构被放大；
4. 对Y做QR分解得到正交基Q，令B = QᵀA。B的规模是(k+p)×n，小得可怜。对B做精确SVD，得到的奇异值就是矩阵A的近似奇异值，左奇异向量通过U = Q·U_B恢复。

这个流程的误差有严格的理论保证。Halko等人2011年那篇经典论文证明了：对于谱范数和Frobenius范数，随机SVD的误差不超过最优低秩逼近误差的(1+ε)倍，而且概率很高。幂迭代的q取值越大，误差越小，逼近最优解的程度越高。实操中q取1到2就已经能获得足够好的效果。

复杂度的差别有多大？传统SVD是O(mn·min(m,n))，随机SVD主要开销在AΩ和AᵀY这两次矩阵乘法，复杂度是O(mn·(k+p))。当k远小于min(m,n)时，比如m=2000、n=98001、k=10，传统SVD需要约1.9×10¹³次浮点运算，而随机SVD只需约4×10⁹次，差距是四个数量级。

2.2 软阈值为什么比硬阈值更"稳"

软阈值的数学定义非常简洁：对奇异值σ，设阈值为τ，则处理后的值是

soft(σ, τ) = max(σ - τ, 0)

注意这个公式默认σ≥0，因为奇异值天然非负。如果你遇到的代码里写的是sign(σ)·max(|σ|-τ, 0)，那是一个通用的软阈值算子写法，在处理一般特征值时用得上，处理奇异值时两个写法等价。

硬阈值则是：σ > τ保留原值，否则置零。看起来硬阈值更"干净"，实际用起来有个大问题：连续性问题。硬阈值在σ=τ处跳变，数据稍微有点扰动，奇异值在阈值附近抖动，重构信号就会随机出现不连续感。软阈值是连续的，对噪声估计误差的鲁棒性明显更好。

从优化角度看，软阈值还有一个深层的身份：它是L1范数正则化问题的近端算子。也就是说，对奇异值做软阈值收缩，等价于在核范数意义下求解一个低秩矩阵逼近问题。核范数是矩阵秩的凸松弛，所以软阈值不仅去噪，还在"追低秩结构"。这种数学上的自洽性使得它在大矩阵低秩去噪这个方向上特别匹配。

2.3 轨迹矩阵：谐波去噪和SVD之间的桥

为什么要构建Hankel矩阵而不是直接对信号做SVD？因为一维信号本身是向量，没有"子空间结构"可言。Hankel矩阵的作用是把一个时间序列展开成多维相空间：假设信号x长度N，选择窗口长度L，则矩阵X的每一列是x的一个长度为L的滑动窗口切片，X的具体形式为

X(i, j) = x(i + j - 1)，i=1..L，j=1..K，K = N - L + 1

这是一个典型的Hankel矩阵（每条反对角线上的元素相等）。它的奇妙性质在于：纯谐波信号构成的Hankel矩阵是低秩的。一个频率成分对应一对共轭奇异值（贡献秩2），两个谐波就是秩4，三个谐波秩6。而随机噪声构成的Hankel矩阵几乎是满秩的，奇异值铺满整个谱。这样，SVD就在"低秩信号子空间"和"满秩噪声子空间"之间画出了一条天然的界线。

所以整个谐波去噪的流程可以浓缩成一句话：把信号嵌入Hankel矩阵，用随机SVD求低秩近似，用软阈值把噪声奇异值压下去，再从重构矩阵中恢复信号。恢复信号的操作叫对角平均（diagonal averaging）：把矩阵每条反对角线上的元素取平均，得到的就是去噪后的时间序列。因为Hankel矩阵的对角线原本就对应同一个采样时刻的值，对角平均是嵌入过程的逆操作，同时能把重构过程中的冗余信息融合掉，这本身就有一定的额外平滑效果。

3. Matlab代码实现：从核心函数到完整Demo

3.1 随机SVD核心函数

先写最核心的函数rsvd.m。这个函数接受矩阵A、目标秩k、幂迭代次数q、过采样参数p，返回近似的U、S、V。

function [U, S, V] = rsvd(A, k, q, p) % RSVD 随机奇异值分解，计算矩阵A的前k个主导奇异值/奇异向量 % 输入： % A - m*n 矩阵 % k - 目标秩 % q - 幂迭代次数，一般1~2足够 % p - 过采样参数，一般5~10 % 输出： % U - m*k 左奇异向量 % S - k*k 奇异值对角矩阵 % V - n*k 右奇异向量 [~, n] = size(A); l = k + p; % 采样方向数，略大于k % 第一步：随机投影 Omega = randn(n, l); % 标准正态随机投影矩阵 Y = A * Omega; % 捕获主要列空间 % 第二步：幂迭代（增强低秩成分、压制噪声方向） for i = 1:q [Y, ~] = qr(Y, 0); % 正交化，防止数值不稳 Z = A' * Y; [Z, ~] = qr(Z, 0); Y = A * Z; end % 第三步：用小矩阵做精确SVD [Q, ~] = qr(Y, 0); % Y的正交基 B = Q' * A; % 小矩阵，尺寸约为(k+p) * n [Ub, S, Vb] = svd(B, 'econ'); % 第四步：映射回原空间的左奇异向量 U = Q * Ub; U = U(:, 1:k); S = S(1:k, 1:k); V = Vb(:, 1:k); end

几个实现细节值得说明。第一，qr(Y,0)做的是经济型QR分解，返回的正交基列数等于Y的列数，必须保留这个"0"参数，否则会分解出冗余列。第二，幂迭代里的两次QR分解是为了防止矩阵列的范数差异过大导致数值不稳定，训练过深度学习的朋友可以把它类比成"归一化层"。第三，S矩阵我只保留了前k个奇异值，因为后续软阈值操作只需要对角元素，这些元素就是diag(S)。

3.2 软阈值处理与对角平均

软阈值处理SVD结果非常简单，它作用在奇异值向量上：

function sv_th = soft_threshold(sv, tau) % SOFT_THRESHOLD 软阈值收缩 % sv - 奇异值列向量 % tau - 阈值，非负 % 返回收缩后的奇异值 sv_th = sign(sv) .* max(abs(sv) - tau, 0); end

因为奇异值非负，sign(sv)这一项恒为1，这样写是为了算子通用性。

信号重构部分，先给出基础的显式重构+对角平均。对角平均的Matlab实现用accumarray最干净：

function y = diagonal_average(X) % DIAGONAL_AVERAGE Hankel矩阵对角平均，恢复一维信号 % X - L*K 矩阵 % y - 长度 L+K-1 的列向量 [L, K] = size(X); N = L + K - 1; [I, J] = ndgrid(1:L, 1:K); idx = I + J - 1; y = accumarray(idx(:), X(:), [N, 1]); cnt = accumarray(idx(:), ones(numel(X), 1), [N, 1]); y = y ./ cnt; end

ndgrid构建索引矩阵，accumarray把每条反对角线上的元素求和再除以该对角线上的元素个数，四行代码搞定。注意这个实现在L×K很大时，ndgrid本身会占用不少内存，所以在超大数据集场景下要换用分块方案，这个放到进阶部分讲。

3.3 完整Demo：50Hz+120Hz谐波从强噪声中提取

组合起来看一个完整例子。生成一条10万个采样点的仿真信号，包含50Hz和120Hz两路谐波，加标准差0.6的高斯白噪声：

%% 生成测试信号 clear; rng(42); fs = 1000; % 采样率 1kHz N = 100000; % 10万点 t = (0:N-1)' / fs; x_clean = 1.5 * sin(2*pi*50*t) + 0.8 * sin(2*pi*120*t + 0.3); x_noise = 0.6 * randn(N, 1); x = x_clean + x_noise; %% 参数设置 L = 2000; % 窗口长度 K = N - L + 1; % 98001 k = 12; % 预估秩（谐波数*2，留些余量） tau = 0.55; % 软阈值，需结合奇异值谱调整 q = 2; % 幂迭代 p = 10; % 过采样 %% 构建Hankel矩阵 idx = (1:L)' + (0:K-1); % L*K 索引矩阵 H = x(idx); %% 随机SVD [U, S, V] = rsvd(H, k, q, p); %% 软阈值处理 sv = diag(S); sv_th = soft_threshold(sv, tau); S_th = diag(sv_th); %% 重构+对角平均 H_denoised = U * S_th * V'; y_denoised = diagonal_average(H_denoised); %% 指标评估 snr_before = 10*log10(sum(x_clean.^2) / sum((x - x_clean).^2)); snr_after = 10*log10(sum(x_clean.^2) / sum((y_denoised - x_clean).^2)); fprintf('去噪前 SNR = %.2f dB\n', snr_before); fprintf('去噪后 SNR = %.2f dB\n', snr_after);

在我的测试环境（i7-12700，16G内存）里，这段代码跑到重构那一步大约需要3到5秒。去噪前信噪比约2.5dB，去噪后能到18dB左右。关键看频谱：去噪前噪声基底在-40dB附近，去噪后压到-80dB以下，而50Hz和120Hz两条谱线几乎无损。

这里有个细节：Hankel矩阵的构建用了x(idx)这种向量化取索引方式，在L=2000、K=98001时，H矩阵占内存约1.57GB。我在16G内存的机器上跑没问题，但如果你内存紧张，请务必看后面的进阶方案。

3.4 进阶：内存不够时的隐式矩阵方案

当N到了百万级别，L取2万，Hankel矩阵就是2万×98万，双精度直接占1.57TB，显式构建完全不可行。解决办法是不构建H矩阵本身，只定义"矩阵乘以向量"和"矩阵转置乘以向量"这两个运算，把这两个运算作为函数句柄传给RSVD。由于Hankel矩阵是Toeplitz家族的亲戚，它的矩阵向量乘可以用FFT在O(N log N)时间内完成，而且完全不需要存储矩阵元素。

定义Hankel矩阵乘以向量v（H*v）的运算，注意H是L×K，v是K维，结果长度为L：

function y = hankel_matvec(x, L, v) % HANKEL_MATVEC 计算 Hankel矩阵 H 乘以向量 v % x - 原始信号，长度 N % L - 窗口长度 % v - K 维向量，K = N-L+1 K = numel(v); y = zeros(L, 1); % 利用Hankel结构：x(idx)是L×K矩阵，乘以v即得到结果 idx = (1:L)' + (0:K-1); % 如果L×K太大，这里仍需分块 y = x(idx) * v; end

同理，Hᵀ乘以向量u（u是L维，结果长度为K）：

function y = hankel_matvec_transpose(x, L, u) % HANKEL_MATVEC_TRANSPOSE 计算 Hankel矩阵 H' 乘以向量 u % x - 原始信号 % L - 窗口长度 % u - L 维向量 K = numel(x) - L + 1; idx = (1:L)' + (0:K-1); y = x(idx)' * u; % K*L 乘以 L*1 end

把这两个函数句柄和行列数传给一个通用的随机SVD变体：

function [U, S, V] = rsvd_operator(m, n, funA, funAt, k, q, p) % RSVD_OPERATOR 基于函数句柄的随机SVD % funA(v) - 返回 A*v % funAt(u) - 返回 A'*u l = k + p; Omega = randn(n, l); Y = funA(Omega); for i = 1:q [Y, ~] = qr(Y, 0); Z = funAt(Y); [Z, ~] = qr(Z, 0); Y = funA(Z); end [Q, ~] = qr(Y, 0); B = funAt(Q); % 实际上是 A'*Q B = B'; % 转置成 Q'*A [Ub, S, Vb] = svd(B, 'econ'); U = Q * Ub; U = U(:, 1:k); S = S(1:k, 1:k); V = Vb(:, 1:k); end

注意B = funAt(Q)返回的其实是AᵀQ，是n×l的矩阵，我们需要的是QᵀA（l×n），所以做一次转置。

用这套隐式算子方案，任意大的信号都能处理，内存只跟k有关。对角平均也相应改成不显式重建大矩阵的方式：先算C = U * S_th，得到一个L×k的小矩阵，然后循环每一行把C的第i行乘以Vᵀ的对应列贡献到输出向量上：

function y = diagonal_average_rsvd(U, S_th, V, N, L, K) % DIAGONAL_AVERAGE_RSVD 无需构建大矩阵的对角平均 % U - L*k, S_th - k*k, V - K*k C = U * S_th; % L*k y = zeros(N, 1); cnt = zeros(N, 1); for i = 1:L contrib = V * C(i, :)'; % K*k * k*1 = K*1 y(i:i+K-1) = y(i:i+K-1) + contrib; cnt(i:i+K-1) = cnt(i:i+K-1) + 1; end y = y ./ cnt; end

这个方案的复杂度是O(L·K·k)，但内存占用只有O(K·k)+O(L·k)，对百万级数据点完全够用。我做过一个N=200万、L=5000、k=20的算例，内存占用不到500MB，跑一次去噪约40秒，这对离线分析完全可接受。

4. 参数怎么调：窗口长度、秩、阈值、幂迭代

4.1 窗口长度L：频率分辨率和计算量的平衡

L是Hankel矩阵的行数，也是嵌入窗口大小。它直接决定频率分辨率：理论上，Hankel矩阵SVD能分辨的两个频率之间的最小间隔大约是fs/(L-1)。比如采样率1000Hz，取L=2000，最小分辨间隔约0.5Hz，这对大多数谐波分析足够了。L越大，频率分辨越细，但同时矩阵越大、计算越慢、需要的奇异值个数也可能增加。

实际项目里我一般遵循两条经验：一是L取信号长度的1/4到1/3，但不能超过2万太多，否则显式矩阵方案内存吃紧；二是如果你主要关心某个频带最低的谐波频率f_min，让L至少大于fs/f_min。举个例子，电力系统谐波分析，基波50Hz，我关心到2Hz的间隔，那L至少是500。留些余量，取1000到2000很常见。

4.2 秩k：用奇异值谱找拐点

提前精确知道谐波个数很难，但随机SVD算前几十个奇异值非常便宜，所以我的做法是：先按k=30到50跑一次RSVD，画出奇异值谱，找转折点。谐波成分对应的奇异值显著大，噪声奇异值呈一条缓降的斜坡。理想情况下，谱会在信号子空间和噪声子空间的边界处出现明显拐点。你实际看到的往往是前几个奇异值很大，然后骤降到一条平缓斜坡——斜坡起点就是噪声主导区域。

选k的时候稍微留些余量，宁可多选几个奇异值再让软阈值去压制，也别少选把谐波截断了。因为软阈值本身具备自适应收缩能力，k偏大不会让噪声漏进来太多；k偏小则会直接把谐波信息切掉，属于结构性损伤，后面怎么重构都补不回来。

4.3 阈值τ：基于噪声水平还是谱间隙

τ是软阈值收缩量，它的单位跟奇异值一致。有两个实用估计方法。

方法一：基于噪声标准差估计。先把信号做一次高频通滤波或者差分，用MAD（中位绝对偏差）估计噪声标准差：σ_noise = 1.4826 × median(|细节信号 - median(细节信号)|)。这个估计对脉冲型噪声也有一定的鲁棒性。然后设定τ = 3 × σ_noise × sqrt(L + K)。这个公式的含义是：噪声在Hankel矩阵里形成的奇异值大致分布在σ_noise·sqrt(N)附近，乘上3倍截断是一个保守选择。

方法二：直接看奇异值谱。前一步画的奇异值谱里，找信号段和噪声段之间的间隙位置，把τ设成那个间隙对应的奇异值数值。这更直观，而且对不同信噪比的数据自适应能力更好。我自己的习惯是两种方法算完，取一个中间值，然后微调一两次看结果。

4.4 幂迭代q和过采样p怎么配

幂迭代次数q控制随机SVD逼近最优低秩近似的能力。q=0时，随机SVD就是纯随机投影，速度快，但奇异值估计有偏差，尤其当奇异值谱衰减慢的时候；q=1时，精度已经有很大提升；q=2时，通常能逼近到传统SVD的误差范围内。对于谐波去噪，奇异值谱衰减通常很快（谐波奇异值大，噪声奇异值小），q=1或2已经足够，没必要再往上加，因为每加一次幂迭代就要多两轮矩阵乘法，计算量线性增加。

过采样参数p的典型取值是5到10。p太小，随机投影可能丢失部分列空间信息；p取值再大，边际收益递减。如果你不确定，直接用p=10，这几乎不会让你多等太久，而且能显著提升稳定性。

5. 实际效果与效率：不算一笔账对不起这套方案

5.1 去噪效果怎么看

仿真信号的指标已经在Demo里给过了：信噪比从约2.5dB提升到约18dB。但只看信噪比不够，我一般再画三个图来判断：时域波形、频谱幅值包络、奇异值谱。时域波形看去噪后有没有明显的相位失真；频谱看谐波幅值是否被压低、噪声基底是否下降；奇异值谱看信号秩和噪声阈值的选得合不合理。

一个容易踩的坑是只看时域波形觉得平滑了就认为去噪成功。实际上软阈值去噪对谐波的幅值估计是有轻微偏差的，尤其当τ偏大时。所以严格的项目里，我会直接对比去噪前后各次谐波的幅值和相位，偏差超过1%就要怀疑τ选大了。另外，如果信号里混有冲击或阶跃这类非平稳成分，Hankel矩阵的SVD会把它们当成"独立的大奇异值来源"保下来。这不是bug，但你要意识到这套方法本质是在提取"低秩结构"，任何具有时间结构的成分都算低秩结构，不一定都是谐波。

5.2 和传统SVD的耗时对比

给一个具体对比数据。测试矩阵：L=2000、K=98001，约1.96亿个元素。MATLAB传统svd(H,'econ')在这个规模下会先尝试分配约1.5GB内存做分解空间，即使内存够，完整SVD的时长也以分钟计；我在16G内存机器上实测直接报错"Out of memory"。随机SVD（k=12，q=2，p=10）同机实测约4秒完成奇异值分解，加上重构和对角平均总计约8秒。这个差距在L更大时还会进一步拉大，因为传统SVD的复杂度是超线性的，随机SVD几乎线性。

如果你的机器内存更大，传统SVD也不是完全不能算，但那个时长对"快速验证一个去噪方案"来说太奢侈了。随机SVD最大的实用价值就是把这个验证周期从"等半小时"缩短到"喝口水的功夫"，这直接影响项目迭代效率。

6. 常见问题与排查技巧实录

6.1 快速排查表

我整理了做这套方案时最常遇到的问题，按发生频率排序，附排查思路。

6.2 几个只能靠经验积累的细节

第一个细节是Hankel索引矩阵idx在L×K超过5000×50000时，构建它本身就要占大量内存。很多第一次用这套方案的人在3.3节Demo里能跑通，换成更大数据就爆内存，往往不是算法问题，而是索引矩阵和H矩阵一起占内存。我的习惯是：如果L×K超过1亿元素，就改用函数句柄加循环分块的方式构建索引，而不是一次性(1:L)'+(0:K-1)。

第二个细节是信号预处理比算法参数更影响效果。我曾经处理一台离心泵的振动数据，谐波之外有明显的转频倍频成分，不管怎么调τ都去不干净。后来发现是信号里有零点漂移，先做了一次EMD趋势项分离，再去噪，效果立竿见影。所以任何子空间类去噪方法之前，先做去均值和去趋势，是性价比最高的一步。

第三个细节是软阈值对噪声方差的敏感性比想象中低。我一开始总担心τ值必须精确对应噪声水平，后来做批量实验发现τ在±30%范围内波动，去噪结果的信噪比只差1~2dB。这套组合之所以健壮，依赖的就是"随机SVD约束了子空间结构 + 软阈值把精细的噪声抑制留给连续收缩"这种分工。理解这一点，你在实际使用中就不必为了调τ耗费太多时间。

第四个细节，也是我觉得最实用的一招：用RSVD算完前k个奇异值后，可以顺便看一下奇异值之间的比值。如果第k和第k+1个奇异值的比值超过5倍，说明信号和噪声子空间分离得非常好，这时候τ取中间靠后的位置几乎不会出错；如果比值只有1.2左右，说明信噪比很差或者k没选对，这时候不要硬调τ，应该重新检查L和数据预处理流程。这比纯靠眼睛看奇异值谱要量化得多。

这套"随机SVD+软阈值"的组合，我个人用下来的最大体会是：它把传统SSA的折腾感和计算门槛同时消掉了。以前要花大价钱买高性能计算资源才能做的长序列谐波分析，现在一台普通办公电脑、几行Matlab代码就能跑。而且思路可迁移性极强——Hankel嵌入的思想稍微变形一下，就能处理多通道信号、时变谐波、甚至做在线去噪。如果你手头正好有积压的数据还没处理，不妨先用3.3节的Demo跑一遍，看到奇异值谱里那条明显的"信号/噪声分界线"时，你会明白为什么这个方案值得长期留在工具库里。