MATLAB从零实现经验小波变换EWT：频谱分割与Meyer滤波器组

做非平稳信号分解的同行，大概都被同一类问题磨过：一段振动信号里有几个频率成分叠在一起，想按频率把它们拆开，EMD 跑出来模态混叠，VMD 得先猜分量数，猜错了要么过分解要么欠分解。经验小波变换（EWT）走的是另一条路——先看信号的傅里叶频谱，在频域里自适应地切出一组连续区间，然后在每个区间上构造一对 Meyer 型小波滤波器，本质就是一个"照着信号频谱长出来的"自适应滤波器组。它不迭代、不搜索、一次 FFT 加一组 IFFT 就能出结果，分解层数由频谱的极大值个数直接决定，物理含义比很多黑箱方法清楚得多。这篇内容就是我在 MATLAB 里从零实现 EWT、并且拿它处理实际信号的全过程整理：频谱怎么切、Meyer 滤波器组的代码怎么写、γ 这个参数怎么取、分量数定不下来怎么办、跑一万点数据卡住又该怎么改。适合已经会用 MATLAB 基本语法、手上有一堆非平稳信号等着分解的人，也适合刚接触 EWT、被网上各种版本的函数名绕晕的人——我会把每一步的意图都写清楚，而不是甩一段代码让你自己猜。

1. 先想清楚：EWT 到底在做什么

1.1 从"切频谱"这个朴素想法说起

把一个信号的傅里叶频谱画出来，你会看到若干个鼓包，每个鼓包对应一个主导频率成分。人眼一看就知道这里有三个成分、那里有五个。EWT 干的第一件事，就是把这个"人眼看鼓包"的过程写成算法：找出频谱上的局部极大值，把它们当作各个成分的中心，相邻中心的中点作为区间边界。切完之后，整个[0, π]频段被分成 M 个连续区间，加上最靠近零频的那一段，一共 M+1 个区间。

第二件事是在每个区间上定义小波。这一点是 EWT 和固定基小波（比如 db4、sym8）最本质的区别：固定小波的频带划分是设计时定死的，跟你的信号没关系；EWT 的频带划分是运行时算出来的，信号频谱长什么样，基函数就长什么样。所以在频率成分间隔不均匀的信号上，EWT 的优势非常明显——低频成分贴得很近时它能切细，高频成分稀疏时它不会浪费基函数。

用生活化的说法打个比方：固定小波像是买一套标准尺寸的收纳盒，无论你装什么都要往这 8 个格子里塞；EWT 像是先看一眼你要装的东西体积分布，再现场定制隔板的位置。代价是定制本身要花时间——也就是那次频谱分析和极值检测。

这里有个概念必须先说透，不然后面代码看不懂。EWT 中的"经验"二字，指的是基函数由数据驱动生成，而不是指它有学习过程。它没有任何训练、优化、迭代收敛的环节，从头到尾是确定性的：同样的输入，永远得到同样的输出。这一点对工程复现极其友好——你不用设随机种子，不用担心两次跑出来结果不一样。

1.2 EWT、EMD、VMD 的取舍对比

我前后在几类信号上把三种方法都跑过一遍，下面这张表是我自己的实际感受，不是抄手册：

我的结论很直接：信号频率成分相对稳定、成分间隔比较开的时候，优先用 EWT，因为它快而且分量数好解释；成分频率随时间快速漂移的时候，EWT 的全局边界会失效，这时候要么做分段 EWT，要么老老实实用 EMD 系列。VMD 介于两者之间，参数调好了确实漂亮，但 K 和 α 的调节成本比 EWT 高不少。

有一个细节值得展开：EWT 的边界是一次性确定的，整段信号共用一套边界。这意味着如果你的信号前 0.5 秒是 50 Hz 主导、后 0.5 秒变成 200 Hz 主导，而频谱上两个峰都在，EWT 会给出两套滤波器，前 0.5 秒那个 200 Hz 滤波器输出接近于零——它不会报错，只是那个分量是空的。这个行为在故障诊断里经常被误判成"这个频段没有能量"，实际上是分割策略的问题。解决办法是滑窗分段做 EWT，或者先做时频分析看清楚成分随时间的分布再决定窗长。

1.3 我把整个实现拆成四块

写代码之前先把结构定下来，比上来就敲省事得多。我自己的实现分成四个函数，职责单一，方便单独调试：

• ewt_boundaries.m：输入信号和分量数，输出边界向量omega（从 0 到 π 单调递增），内部包含频谱计算、平滑、极值检测、边界生成。
• ewt_meyer_filterbank.m：输入信号长度 N 和边界向量，输出 N×(M+1) 的滤波器组矩阵，每一列是一个滤波器的频域响应，实对称。
• meyer_beta.m：Meyer 的过渡带函数 β(x) = x⁴(35 - 36x + 25x² - 10x³)，独立成函数是因为它被低通和带通同时调用。
• ewt1d_simple.m：主流程，串起上面三个，输出各分量时域波形和滤波器组，方便后处理画图。

注意：这四个函数名是我自己起的，跟网上流传的开源工具箱命名不一定一致。你如果混用不同来源的函数，最容易踩的坑就是边界向量的定义域不统一——有的实现用归一化到 [0, π] 的角频率，有的用 [0, 0.5] 的归一化频率（也就是除以采样率），两者差一个 2π 因子，混用之后滤波器全部错位，输出看起来"有分量"但频率完全不对。我的建议是自己写一套，从头到尾只用一种约定。

2. 核心机制拆解：频谱分割与 Meyer 滤波器组

2.1 边界是怎么"找"出来的

频谱分割是 EWT 里唯一带"经验"色彩的一步，也是最容易出问题的一步。基本流程是：算单边幅值谱 → 平滑 → 找局部极大值 → 按幅值排序取前 M-1 个 → 按频率重新排序 → 相邻取中点作为边界 → 两端补 0 和 π。

每一步都有讲究。平滑这一步很多人直接跳过，结果原始频谱上全是噪声引起的毛刺，极值检测找出几十个假峰，边界挤成一团，分出来的分量全是噪声。平滑窗口长度的经验值是 N/20 到 N/50 之间，信号短、噪声大就取大一点。我一般先用 N/30 试，看分割结果再调。

极值检测之后要丢掉直流点。频谱第一个点（k=1，对应 0 Hz）几乎总是局部极大值，因为信号有直流分量或趋势项，不丢掉它的话，第一个边界会被压到极低频，第一个分量变成一个近乎直流的"残渣带"。

选前 M-1 个这一步，排序依据是幅值而不是频率位置。原因很实际：小能量的毛刺峰不该占用边界配额，幅值大的峰才代表真实成分。选完之后必须再按频率升序排一次，否则边界不单调，后面的滤波器构造会直接崩掉。

分量数 M 的三种常见确定方式，我整理成下面这张表：

第三种的思路值得说一句：对频谱做一系列不同尺度的高斯平滑，真正的峰在各个尺度下都能存活，噪声峰在大尺度下就消失了。统计每个极值点"存活"的尺度范围，取最长的那些。这个方法理论上很优雅，但平滑尺度范围的选择本身又是个新参数，实际用起来不如直接指定 M 痛快。

2.2 Meyer 滤波器组的构造与 γ 的取值

边界定好之后，接下来是在频域直接写出每个滤波器的响应。这里用的是 Meyer 小波的频域形式，因为它在频域有紧支撑、过渡带平滑、重构误差小。

先定义过渡带函数：

β(x) = x^4 * (35 - 36x + 25x^2 - 10x^3), x ∈ [0, 1] β(x) = 0, x ≤ 0 β(x) = 1, x ≥ 1

这个多项式不是随便凑的，它满足 β(0)=0、β(1)=1，且在两端的一阶、二阶导数都为零，所以拼接处足够光滑，不会在时域引入振铃。

设边界向量为 ω₀=0 < ω₁ < ω₂ < ... < ω_M = π。尺度函数（低通）的频域响应：

φ̂₁(ω) = 1, |ω| ≤ (1-γ)ω₁ φ̂₁(ω) = cos[ (π/2) β( (|ω|-(1-γ)ω₁) / (2γω₁) ) ], (1-γ)ω₁ < |ω| < (1+γ)ω₁ φ̂₁(ω) = 0, 其他

第 n 个经验小波函数的频域响应（n = 1, ..., M-1）：

ψ̂ₙ(ω) = 1, (1+γ)ωₙ ≤ |ω| ≤ (1-γ)ωₙ₊₁ ψ̂ₙ(ω) = cos[ (π/2)β( (|ω|-(1-γ)ωₙ₊₁)/(2γωₙ₊₁) ) ], (1-γ)ωₙ₊₁ < |ω| < (1+γ)ωₙ₊₁ ψ̂ₙ(ω) = sin[ (π/2)β( (|ω|-(1-γ)ωₙ)/(2γωₙ) ) ], (1-γ)ωₙ < |ω| < (1+γ)ωₙ ψ̂ₙ(ω) = 0, 其他

最后一个 ψ̂_M 覆盖从 (1+γ)ω_M 到 π 的整段，右侧不设过渡带，因为已经到奈奎斯特频率，截断即可。

γ 是过渡带的宽度系数，取值必须满足：

γ < min_n [ (ωₙ₊₁ - ωₙ) / (ωₙ₊₁ + ωₙ) ], n = 1, ..., M-1

这个约束的物理含义很清楚：相邻两个频带的过渡带不能重叠，否则同一个频率会被两个滤波器同时放大，重构就不再是"近似完美"的了。我实际写代码时会在算出的最小值上再乘 0.9，留一点余量，纯粹是为了避免浮点误差恰好落在边界上导致的数值不稳定。

举个具体的算例，假设采样率 fs = 1000 Hz，信号长度 N = 1000，检测到的三个峰在 15 Hz、80 Hz、220 Hz，那么边界是中点：

• ω₁ 对应 (0+15)/2 = 7.5 Hz
• ω₂ 对应 (15+80)/2 = 47.5 Hz
• ω₃ 对应 (80+220)/2 = 150 Hz
• ω₄ = π 对应 500 Hz（奈奎斯特）

转换成角频率（乘 2π/fs），ω₁ ≈ 0.0471 rad，ω₂ ≈ 0.2985，ω₃ ≈ 0.9425，ω₄ ≈ 3.1416。逐项算 γ 约束：

• n=1：(0.2985-0.0471)/(0.2985+0.0471) = 0.2514/0.3456 ≈ 0.727
• n=2：(0.9425-0.2985)/(0.9425+0.2985) = 0.644/1.241 ≈ 0.519
• n=3：(3.1416-0.9425)/(3.1416+0.9425) = 2.199/4.084 ≈ 0.538

最小值 0.519，乘 0.9 后 γ ≈ 0.467。这个值偏大但不违规——它意味着过渡带很宽，第一个滤波器的上过渡带到 (1+0.467)×0.0471 ≈ 0.0691 rad，也就是 11 Hz，而 15 Hz 的成分中心在 0.0942 rad，还在通带内，没问题。但如果两个成分靠得更近，γ 会被压得很小，过渡带变窄，时域的振铃会变明显。这是一个需要权衡的地方，后面讲端点效应时再展开。

2.3 三个容易忽略的数值细节

第一个是频率轴的对齐。MATLAB 的fft输出第 k 个点（k 从 1 开始）对应角频率 2π(k-1)/N，范围是 [0, 2π)。取前 floor(N/2)+1 个点就是 [0, π] 的单边谱。构造滤波器时，我习惯先把整个 [0, 2π) 的 ω 向量算出来，再把大于 π 的部分折叠回 [0, π]（w(w>pi) = 2*pi - w(w>pi)），这样滤波器矩阵天然是实对称的，做ifft之后取实部误差极小。如果只构造半边再手工镜像，很容易在 N 为偶数时把奈奎斯特那个点搞错一位，表现出来就是重构误差莫名其妙地大。

第二个是幅值归一化。abs(fft(x))/N得到的是单边幅值谱的一半（除了直流和奈奎斯特点），做极值检测时这个系数不影响位置，但如果你要拿幅值谱去和阈值比较，就得乘 2。我一般只在极值检测时用归一化谱，其他场合都保留原始abs(fft(x))，避免自己绕晕。

第三个是序列长度。EWT 本身对 N 没有特殊要求，fft会自动补零到合适的长度。但补零会让频谱插值变密，极值检测可能多找出几个"伪峰"——它们不是真实频率成分，只是补零带来的频谱插值纹波。所以我的做法是：如果信号长度是个很别扭的数（比如 997），先裁掉尾部几个点凑成 2 的幂或至少是 2×3×5 的倍数，再做 EWT。代价是丢掉几个采样点，收益是频率分辨率确定、边界位置可预期。

实操心得：我在调试阶段习惯把abs(fft(x))的前半段画出来，把检测到的极大值位置用红色竖线叠上去，再把边界用虚线标出来。这一步花不了两分钟，但能立刻看出"分量数不对"到底是极值检测错了还是排序取错了。很多新手直接看分解结果的时域波形，那里面信息量太少，出了问题根本定位不到是哪一步。

3. MATLAB 手写实现：从频谱到分量

3.1 函数清单与工程目录

我建议的目录结构是这样的，简单粗暴但够用：

ewt_demo/ ├── main_demo.m % 演示脚本 ├── ewt1d_simple.m % 主流程 ├── ewt_boundaries.m % 边界检测 ├── ewt_meyer_filterbank.m % 滤波器组构造 ├── meyer_beta.m % 过渡带函数 └── find_local_max.m % 局部极大值检测

不用 MATLAB 的包（package）机制，是因为这类小工具函数在多个项目间复用，平铺在路径里最方便。如果你要长期维护，改成+ewt/包目录也行，调用时写ewt.ewt1d_simple即可。

依赖方面，核心功能只需要基础 MATLAB，不需要任何工具箱。只有画包络谱时用到hilbert，那需要 Signal Processing Toolbox；极值检测我也手写了，不依赖findpeaks。这样做的原因很实际——很多人是在没有全套工具箱的环境里跑代码，能少一个依赖就少一个部署障碍。

3.2 第一步：频谱极值检测与边界生成

先写最底层的两个函数。局部极大值检测，注意用>和>=的组合，避免平顶峰被重复计数：

function idx = find_local_max(v) % 找出向量 v 中所有严格局部极大值的位置（首尾点不参与） v = v(:); idx = find(v(2:end-1) > v(1:end-2) & v(2:end-1) >= v(3:end)) + 1; end

过渡带函数，注意三个分段要处理干净，尤其是边界处的赋值顺序：

function b = meyer_beta(x) % Meyer 过渡带函数，输入可为任意实数 b = zeros(size(x)); b(x >= 1) = 1; % 先处理大于 1 的部分 idx = (x > 0) & (x < 1); % 中间段单独算 xx = x(idx); b(idx) = xx.^4 .* (35 - 36*xx + 25*xx.^2 - 10*xx.^3); % x <= 0 的部分保持 0 end

然后是边界检测。这个函数里我加了平滑窗口参数，默认按 N/30 取，最小 3：

function [omega, wc] = ewt_boundaries(f, M, smoothLen) f = f(:); N = numel(f); F = fft(f); % 取 [0, pi] 单边谱 half = floor(N/2) + 1; mag = abs(F(1:half)); % 频谱平滑，抑制噪声引起的伪极值 if nargin < 3 || isempty(smoothLen) smoothLen = max(3, round(N/30)); end if smoothLen > 1 mag = movmean(mag, smoothLen); end % 极值检测，去掉直流点 idx = find_local_max(mag); idx(idx == 1) = []; if numel(idx) < M - 1 error('检测到的极值只有 %d 个，不够分出 %d 个分量，请减少 M 或减小 smoothLen', ... numel(idx), M); end % 按幅值选前 M-1 个，再按频率升序 [~, ord] = sort(mag(idx), 'descend'); idx = sort(idx(ord(1:M-1))); % 索引转角频率：第 k 个点对应 2*pi*(k-1)/N wc = 2*pi*(idx - 1) / N; % 相邻中点作为边界，两端补 0 和 pi omega = [0; (wc(1:end-1) + wc(2:end)) / 2; pi]; end

这段里有两个地方值得单独说。idx(idx == 1) = []这行看着不起眼，但它决定了第一段频带的范围。去掉直流点之后，最小的边界是"最低频率成分"和"次低频率成分"的中点，第一个分量就有了一个合理的下限，而不是被压到 5 Hz 以下。

error那句是我特意加的。不检查极值个数的话，当极值不够时，ord(1:M-1)会返回一个短的索引向量，sort之后omega的长度就不对了，后面的滤波器组会生成错误的列数，报错信息要么是矩阵维度不匹配，要么更糟糕——程序跑完了，结果是错的。工程代码里这种"提前报错"比"静默出错"值钱得多。

3.3 第二步：Meyer 滤波器组的频域构造

滤波器组是整个实现的核心，也是最容易写错的地方。我的写法是先算整个 [0, 2π) 的频率向量，再折叠，然后逐个滤波器填值：

function mfb = ewt_meyer_filterbank(N, omega) % omega: 1 x (M+1)，从 0 单调递增到 pi M = numel(omega) - 1; mfb = zeros(N, M + 1); % 全频段角频率，并折叠到 [0, pi] w = 2*pi*(0:N-1)' / N; w(w > pi) = 2*pi - w(w > pi); % 计算 gamma，留 10% 余量 gamma = 1; for n = 1:M-1 gamma = min(gamma, (omega(n+1) - omega(n)) / (omega(n+1) + omega(n))); end gamma = 0.9 * gamma; % 第 1 列：经验尺度函数（低通） mfb(:,1) = meyer_lowpass(w, omega(2), gamma); % 第 2 列到第 M+1 列：经验小波函数 for n = 1:M isLast = (n == M); mfb(:, n+1) = meyer_bandpass(w, omega(n+1), omega(n+2), gamma, isLast); end end

注意这里的索引偏移：omega(1)=0，omega(2)才是第一个非零边界 ω₁，omega(M+1)=π。循环里第 n 个滤波器用omega(n+1)作为左边界、omega(n+2)作为右边界，当 n=M 时右边界是omega(M+2)——越界了。所以我用isLast标记，最后一个滤波器直接覆盖到 π，不去读omega(M+2)。这种索引偏移是新手上手 EWT 最大的障碍，建议画个表对照着写：

两个辅助滤波器函数如下。低通：

function h = meyer_lowpass(w, w1, gamma) h = zeros(size(w)); a = (1 - gamma) * w1; b = (1 + gamma) * w1; h(w <= a) = 1; idx = (w > a) & (w < b); h(idx) = cos(pi/2 * meyer_beta((w(idx) - a) / (b - a))); end

带通：

function h = meyer_bandpass(w, wn, wn1, gamma, isLast) h = zeros(size(w)); % 左过渡带（上升沿） a1 = (1 - gamma) * wn; b1 = (1 + gamma) * wn; if wn > 0 idx = (w > a1) & (w < b1); h(idx) = sin(pi/2 * meyer_beta((w(idx) - a1) / (b1 - a1))); end % 平顶段与右过渡带（下降沿） if isLast h(w >= b1) = 1; % 一直到 pi 都通 else a2 = (1 - gamma) * wn1; b2 = (1 + gamma) * wn1; h(w >= b1 & w <= a2) = 1; idx2 = (w > a2) & (w < b2); h(idx2) = cos(pi/2 * meyer_beta((w(idx2) - a2) / (b2 - a2))); end end

if wn > 0这个判断是为了第一个小波（n=1）准备的。它的左边界是 ω₁ > 0，按说不需要判断，但如果某个版本里边界向量开头多了个零，a1 就是 0，(w - a1)/(b1 - a1)在 w=0 时会得到 0，meyer_beta(0)=0，sin(0)=0，结果正确，但会覆盖掉h(w<=a)的赋值——实际上 w=0 那个点本来就该是 0。加了判断更保险。这是我调试时踩的坑：有一版代码第一个分量在零频附近有一条奇怪的凹陷，查了半天才定位到这里。

注意：构造完成后建议做一次自检，把mfb每一列平方求和后在频率轴上画出来。对于近似完美重构的滤波器组，这个和在整个频段上应该接近常数 1（在过渡带会略小于 1）。如果看到明显的凹陷，说明 γ 太大，相邻过渡带重叠不足或者边界太近；如果看到某个区间恒大于 1，说明两个滤波器在同一频段都是通带，边界算错了。

3.4 第三步：分解、重构与误差校验

主流程反而最短：

function [ewtc, mfb, omega, wc] = ewt1d_simple(f, M, smoothLen) f = f(:); N = numel(f); [omega, wc] = ewt_boundaries(f, M, smoothLen); mfb = ewt_meyer_filterbank(N, omega); F = fft(f); ewtc = zeros(N, size(mfb, 2)); for k = 1:size(mfb, 2) ewtc(:, k) = real(ifft(F .* mfb(:, k))); end end

输出的ewtc第一列是近似分量（低频段），第 2 到 M+1 列是按频率从低到高排列的经验小波分量。列的顺序天然有序，这一点比 EMD 的 IMF 方便——EMD 的模态顺序在某些实现里是乱的，做批处理时还得自己排。

重构校验一定要做，而且要用相对误差而不是绝对误差：

xr = sum(ewtc, 2); err = norm(xr - f) / norm(f); fprintf('重构相对误差：%.3e\n', err);

正常情况这个值应该在 1e-12 到 1e-8 之间，取决于 γ 和边界是否过近。如果到了 1e-2 量级，那基本就是滤波器组有重叠或者缺失，回去查 3.3 节的平方和自检。我遇到过的最离谱的一次是相对误差 0.35，最后发现是omega向量没排序——极值按幅值选完之后忘了按频率重排，边界出现回折，滤波器互相覆盖。

3.5 一段完整的可运行脚本

把上面拼起来，用合成信号验证。信号包含 15 Hz、80 Hz、220 Hz 三个成分，采样率 1000 Hz，长度 1000 点（正好 1 秒），加一点高斯白噪声：

clear; clc; close all; fs = 1000; N = 1000; t = (0:N-1)' / fs; x = 1.0*cos(2*pi*15*t) ... + 0.6*cos(2*pi*80*t + 0.4) ... + 0.35*cos(2*pi*220*t) ... + 0.15*randn(N, 1); % 期望分 3 个成分（外加 1 个近似分量） M = 3; [ewtc, mfb, omega, wc] = ewt1d_simple(x, M, round(N/30)); fprintf('检测到的中心频率：%s Hz\n', mat2str(round(wc(:)'*fs/(2*pi)), 4)); fprintf('边界频率：%s Hz\n', mat2str(round(omega(:)'*fs/(2*pi)), 4)); xr = sum(ewtc, 2); err = norm(xr - x) / norm(x); fprintf('重构相对误差：%.3e\n', err); % 画图：原信号、各分量、滤波器组 figure('Position', [100 100 900 700]); subplot(4,1,1); plot(t, x); title('原始信号'); xlabel('时间/s'); for k = 1:size(ewtc, 2) subplot(4,1,k+1); plot(t, ewtc(:,k)); title(sprintf('分量 %d', k)); xlabel('时间/s'); end figure; plot((0:N-1)*fs/N, abs(mfb), 'LineWidth', 1.2); xlim([0 fs/2]); xlabel('频率/Hz'); ylabel('幅频响应'); title('Meyer 滤波器组'); grid on;

跑出来你应该看到：检测到的中心频率约 15、80、220 Hz，边界落在 7.5、47.5、150 Hz，重构误差在 1e-13 量级，四个子图里第一个是 15 Hz 的正弦（带一点低频噪声），第二个是 80 Hz，第三个是 220 Hz，第四个（近似分量）基本只剩下低频噪声和残余。滤波器组的图上能看到三条清晰的平顶带，过渡带宽度由 γ 决定。

如果跑出来只有一个分量有能量、其他都是平线，八成是wc的索引转换错了，或者omega没有单调。如果重构误差正常但分量频率不对，那大概率是频率轴折叠那步写错了。

3.6 结果怎么读：时域分量与包络谱

分量拆出来之后，最常接的两步后处理是时域波形观察和包络谱分析。前者看瞬态冲击，后者找周期性重复的故障特征。

包络谱的做法是：对某个分量做 Hilbert 变换取解析信号，取模得到包络，再对包络去均值后做 FFT，看低频段有没有等间隔的谱线。

k = 2; % 想分析的分量序号 env = abs(hilbert(ewtc(:,k) - mean(ewtc(:,k)))); env = env - mean(env); E = abs(fft(env)) / N; f = (0:N-1)' * fs / N; plot(f(1:round(N/10)), E(1:round(N/10))); xlabel('频率/Hz'); ylabel('包络幅值'); title('包络谱');

这段代码在旋转机械的振动分析里非常常用。分量里如果藏着周期性冲击，包络谱上会出现间隔等于冲击重复频率的一簇谱线及其谐波。关键前提是分量的频带选对了——如果 EWT 的边界正好卡在共振带中间，把冲击能量劈成两半，包络谱的特征频率会变模糊。这也是我前面强调要先把频谱和边界画出来对照着看的原因。

一个提升成功率的做法：在 EWT 之前先对原始信号做一次简单的预白化或者谱峭度筛选，把确定的共振频带找出来，再让 EWT 的分量数至少覆盖这个频带。谱峭度本身也能给出"哪个频带最冲击"的线索，两者结合起来比单用 EWT 稳得多。

4. 踩坑记录与排查手册

4.1 分量数不对的四种典型症状

分量数 M 是 EWT 里唯一需要人拍板的参数，也是最容易出问题的地方。我把实际遇到的情况整理成这张速查表：

第一种情况最隐蔽。因为它不报错，输出看着也"正常"，只是某个分量里混了两个频率。我一般的判断方法是看该分量的包络谱或者时频图，如果有两条稳定谱线，那就是没分开。这时候不是盲目加 M，而是先回去看频谱——如果两个峰之间没有明显的谷，加 M 也没用，因为极大值检测本来就检不出来，这种情况得考虑先做带通滤波或者换方法。

第二种情况在含噪信号里特别常见。我的经验是平滑窗口取 N/30 起步，噪声大的时候加倍。如果加大窗口之后真实峰也被抹平了，说明信噪比已经低到 EWT 不适用的程度，与其硬调参数，不如先做降噪。

4.2 端点效应与模态混叠

严格说 EWT 没有 EMD 那种经典的端点效应，因为它全程在频域做滤波，不存在包络拟合的边界外推问题。但它有自己的"软端点效应"：如果信号存在明显的趋势项（比如缓慢漂移的基线），这部分能量集中在 0 Hz 附近，会把第一个边界往下压，导致第一个分量变成一个又宽又低能量的带，里面的内容既不像趋势也不像成分。

处理办法很土但有效：分解前先detrend或者减均值。别小看这一步，我在处理实际采集的振动数据时，传感器温漂造成的基线漂移能量往往比目标成分还大，不做预处理的话，前两个分量全被它占着。

模态混叠在 EWT 里的表现形式和 EMD 不同。EMD 是同一个频率出现在多个 IMF 里，EWT 则是同一边界内的两个成分无法分开。判断依据很直接：看这个分量的频谱是不是只有一个主瓣。如果有两个，就是混叠了。解决办法有三个方向——提高频率分辨率（加长数据）、减小平滑窗口让两个峰都被检出、或者先做一次窄带滤波把其中一个成分挪出去。

还有一个不常被提到但很实际的问题：过渡带能量丢失。Meyer 滤波器组在过渡带上的响应是部分通过，两个相邻滤波器的过渡带严格来说不重叠（γ 约束保证），所以落在过渡带里的能量会被衰减。如果某个弱成分恰好落在边界附近，它会被削掉一部分，幅值偏小。这时候不能怪 EWT，是你把边界放错地方了——边界应该在两个峰之间的谷底，而不该压在某个峰上。这也解释了为什么"中点"策略不是最优的：当两个峰的间距不对称时，中点可能偏向其中一个。更讲究的做法是在相邻峰之间的局部最小值处取边界，代价是要多写几行代码，而且谷底可能很平，取哪个点都差不多。

实操心得：我处理过一段数据，某个 40 Hz 左右的成分总是比理论幅值小 15% 左右，查了很久以为是传感器标定问题。后来把边界画出来才发现，那个成分的频率恰好落在两个极大值的中点附近，也就是过渡带上。把边界改到两个峰之间的谷底之后，幅值立刻对上了。这个教训是：别默认中点就是最优边界，尤其是当频谱谷不对称的时候。

4.3 速度与工程化：从脚本到可部署函数

说完精度说速度。EWT 的计算量主体是一次长度为 N 的 FFT、M+1 次逐点复数乘法、M+1 次 IFFT。M 一般不超过 10，所以总开销大约是 2(M+1) 次 FFT 的量级。对 10 万点的信号，在普通笔记本上大概是零点几秒，完全可以接受。

真正拖慢速度的往往不是算法本身，而是几处实现细节：

第一，mfb矩阵大小是 N×(M+1)，N 很大时内存占用不小。10 万点、10 个分量就是 8 MB 的双精度复数（如果构造时用了复数）或者 8 MB 实数。如果信号更长，建议改成逐列计算、算完一列就ifft并释放，别一次性把整个矩阵留在内存里。

第二，real(ifft(...))里那个real在理论上是不需要的（因为滤波器对称、输入实信号，结果本来就是实的），但浮点误差会产生极小的虚部。如果你要接着做 Hilbert 变换或者求包络，这些虚部会被放大，所以保留real是必要的。代价几乎为零。

第三，如果要做批量处理（几百段信号各自做 EWT），把ewt_boundaries里的movmean换成手写的滑动平均会快一些，因为movmean对短向量有固定开销。不过在绝大多数场景下这点差异可以忽略，不值得为它牺牲可读性。

工程化方面，我建议至少做三件事。一是把参数校验写进函数入口，检查M是否为正整数、信号长度是否大于2*M、omega是否严格单调。二是把中间结果（频谱、极大值索引、边界、滤波器组）作为可选输出返回，方便出问题时复现。三是固定随机相关的操作——虽然 EWT 本身是确定性的，但如果你的预处理里有加噪、抖动之类的操作，把种子固定住，否则同一份数据两次跑出来不一样，排查起来会疯。

最后提一个扩展方向。一维 EWT 的框架可以比较自然地推到二维：把二维频谱按极坐标或直角坐标划分区域，在每个区域上构造二维 Meyer 滤波器，做二维 FFT 后逐区域滤波。图像处理里做纹理分离、条纹去除时用过这套思路，效果不错，但边界划分的策略要重新设计，一维那套"找极大值取中点"在二维频谱上会变得很复杂。如果你只是想做信号分解，一维这套已经够用了。

我个人在实际操作中的体会是，EWT 最大的价值不在于它比 EMD 或 VMD 更"先进"，而在于它的每一个中间结果都能画出来看、每一处不对劲都能追溯到具体某一步。频谱画出来、极值标出来、边界画出来、滤波器组画出来，四张图一摆，问题基本就无处可藏了。上面这套代码我用了挺长时间，改过好几版，最长的一次调试花了整个下午，最后发现只是omega少了个pi端点——所以如果你跑出来结果不对劲，先别怀疑算法，把四个中间量挨个打印一遍，八成能自己找到答案。