基于小波模极大值的信号奇异点检测与MATLAB实现

简介：面向小波模极大值方法在MATLAB中的实现需求，这份资源以单脚本形式提供了完整的算法演示，适用于信号检测、特征提取与图像分析方向的学习者，以及需要处理突变信号或噪声背景的研究人员。压缩包内仅含1个m文件，整体体积仅563B，但脚本流程清晰，依次涵盖小波分解、各尺度模极大值计算与定位、阈值去噪处理以及基于极大值点的信号重构或特征输出，可直接在MATLAB中运行修改。已有227人学习了该资源，适合希望理解“cgau”连续高斯小波特性、掌握wavedec与waverec函数配合使用的读者，尤其对想快速入门小波模极大值原理的初学者颇具参考价值。通过研读并运行该脚本，读者既能观察到信号突变点的定位效果，也能借鉴其代码框架，将小波模极大值方法迁移至地震数据、医学成像或金融时间序列等实际场景中，是一份轻量而实用的动手学习资料。

1. 小波模极大值检测奇异点：先看信号在哪一刻变了

把压缩包名称拆开，“qiyidian”是“奇异点”的拼音，“小波模极大值”才是核心算法：信号在某一瞬间发生跳变、脉冲或斜率突变时，普通差分和FFT只能告诉你“有大变化”，却回答不了变化发生在第几个采样点、突变有多尖锐。小波模极大值方法把原始信号放到多尺度空间里观察，每个突变点的小波系数模极大值会沿尺度方向形成一条极值线，其衰减速率由Lipschitz指数决定。把这个指数估算出来，就能把“突变”从定性变成定量。下面这套流程面向故障诊断、信号处理和MATLAB工程实践，从算法原理、代码实现、参数调整到验证技巧，完整走一遍从原始波形到奇异点列表的路线。新手能按步骤复现，熟手也能看到噪声条件下参数取舍的边界。

2. 小波模极大值的数学基础：Lipschitz指数与极值线传播

2.1 Lipschitz指数度量奇异强度

奇异点不是非黑即白的“跳变”：同样的电压突降，电网上相邻节点捕捉到的陡峭程度不同；同样的边缘像素，经图像退化后梯度也会变缓。Lipschitz指数把“陡峭程度”收敛成一个实数。若x(t)在t0附近满足

|x(t) - x(t0)| ≤ K·|t - t0|^α

则称x(t)在t0处具有Lipschitz指数α。α=1对应光滑可微点，α=0对应有界间断（阶跃），α=-1对应比阶跃更尖锐的冲击（理想脉冲），0和1之间的值对应斜坡状突变。指数越小，奇异越强。

在离散信号里直接算这个指数很困难，因为采样点落在突变处的相位会严重影响差分结果。小波变换的优势在于，它相当于用可伸缩的窗口对信号做带通卷积，每个尺度的输出都携带该尺度下突变轮廓的信息。Mallat证明了小波系数模极大值满足

log|Wx(s,t)| ≤ log K + α·log s

这给了工程上一个可操作的估计办法：沿极值线取多个尺度的小波系数模，在双对数坐标里拟合直线，斜率即α。实现这一步所需的全部工作，就是先把小波变换系数算对，再把极值点从二维矩阵里准确拎出来。

2.2 小波基的选择影响极值线含义

先将数学结论落到小波基选择上。用于奇异点检测的小波必须满足两个条件：具有一阶消失矩或以上，能够反映信号的局部变化；尽可能平滑，避免小波自身的高频抖动在模极大值提取阶段产生伪点。最常用的是高斯函数的一阶导数，即gauss1小波。gauss1在时域是单峰形态的疏波，奇数对称，模板本身没有直流分量，卷积结果近似于对信号做一阶微分后再平滑。

不同小波基对应不同的“极值位置含义”。表2-1列出的定位关系决定后续追踪逻辑，建议在写代码前先想清楚自己的目标是找信号本身突变（用gauss1），还是找信号拐点（用gauss2）。

表2-1 不同小波基对奇异点检测的定位差异

为什么不建议直接用MATLAB的cwt默认morse小波估Lipschitz指数：morse小波在尺度频带上的定位更偏向时频分析，其相位特性和极值衰减规律比gauss1复杂，拟合出的斜率波动大，解释起来也更困难。gauss1的另一个好处是模板解析式简单，便于在调试时手算特定尺度下的期望模值。如果需要确认当前MATLAB版本是否支持直接调用gauss1，在命令行执行waveinfo('gaus')查看即可。

2.3 极值线在尺度空间中的传播与偏移

一个小窍门：若在最小尺度上检测到模极大值，并且随着尺度增大该极大值点持续存在、位置漂移不超过各自尺度下小波支撑半径，则基本可判定这是一个真实奇异点。理论来源是，奇异点周围小波系数构成锥形影响区，极值线穿过该区域一直延伸到最大尺度。反过来，噪声形成的极值线覆盖尺度很短，通常在2到3个尺度就消失。

定位和偏移的量级需要估计。尺度s下gauss1的等效支撑半径大约为2.5s个采样间隔；若该尺度下有另一个同极性奇异点在附近，两条极值线会在大尺度方向汇合。实际检测时，时间坐标取极值线最细尺度端的位置，而不是把整条线平均，这样能避免大尺度偏移引入的定位误差。

下面这段代码是极值线追踪和Lipschitz拟合的基础，也解释了为什么会用双对数坐标：

scales = 2 .^ (1:0.5:7); % 指数增长的尺度序列 amp = 0.8 * scales .^ (-0.35); % 假设α=-0.35 loglog(scales, amp, 'o-'); xlabel('尺度 s'); ylabel('模极大值 |W|');

若模值随尺度呈幂律衰减，在双对数坐标里就是一条直线，直线斜率即α。写代码之前先造出这样一组数据做基准，实现阶段会少走很多弯路。

3. MATLAB实现：从连续小波变换到模极大值提取

3.1 构造测试信号与gauss1小波核

先搭建一个最小可运行流程。测试信号应包含至少两种尺度的突变：阶跃和冲击。阶跃突变对应α=0，冲击对应α≈-1，二者在模极大值提取中的行为差异足够明显。

Fs = 1000; t = (0:999) / Fs; x = sin(2*pi*50*t); % 干净背景 x(300:end) = x(300:end) + 0.6; % 阶跃突变，位于第300点 x(700) = x(700) + 1.0; % 冲击突变，位于第700点

信号构造本身很简单，第300点开始整体抬高形成阶跃，第700点是单点脉冲。测试信号的作用是给整个流程提供一个可对照的真值，检测完成后能直接拿检测位置与构造位置做差，定位误差一目了然。

接下来生成gauss1小波核。常见做法是先固定时间轴分辨率dt，再按尺度s决定核长度。这样生成的每个尺度核都覆盖了高斯函数的完整支撑范围，不会因为核太短把高频细节切掉：

dt = 1 / Fs; nmax = ceil(8 * s / dt); % 尺度s下核的半长度 nu = -nmax:nmax; psi = -(nu * dt / s) .* exp(-(nu * dt / s).^2 / 2); W = conv(x, psi, 'same') * (dt / sqrt(s));

这段代码的核心是把连续小波定义离散化：小波模板ψ(u)的自变量u取为(n·dt)/s，尺度越大，核越长，对应观察的波形窗口越宽。末尾乘上dt/√s是把积分写成黎曼和后的小波能量归一化，保证不同尺度下的模值处于同一量级，后面Lipschitz拟合才有意义。完整循环就是逐尺度做卷积，把每一层的输出作为一行，堆成W矩阵。

信号两端用conv的same模式会默认补零，制造边界假极值。建议先对x做镜像延拓：

xext = wextend('sym', 2, x, floor(nmax/2)+2); Wext = conv(xext, psi, 'same'); W = Wext(floor(nmax/2)+3 : end - floor(nmax/2)-2);

延拓后出现假极值的位置会被裁剪掉，这是数字实现里最容易忽略的一层，也是“同样代码别人结果好、自己结果乱”的最常见原因。完整代码中需要把延拓放在尺度循环外一次完成，避免每个尺度重复延拓造成边界不一致。

3.2 在每个尺度行上提取模极大值候选点

小波系数矩阵建好后，奇异点检测的下一步是把每个尺度上|W|的局部极大值抽出来。这里要区分两个概念：单尺度内的局部极大值，和跨尺度同一条极值线上的点。前者用一维滑动比较或islocalmax即可：

threshold = 0.1 * max(abs(W(:))); cand = cell(length(scales), 1); for k = 1:length(scales) w = abs(W(k, :)); lm = islocalmax(w, 'MinProminence', 0.05 * max(w)); lm = lm & (w > threshold); pos = find(lm); cand{k} = [pos(:), w(lm(:))]; % [位置，模值] end

islocalmax的MinProminence参数值得展开说。没有这个参数时，只要某点两侧比邻点大就会判定为极大值，含噪信号里会出现大量假峰。MinProminence设为当前尺度最大模值的5%，能滤掉明显的浅峰。threshold阈值再兜底一次，把整体模值很低的候选点全部去掉。如果你的MATLAB版本较老，没有islocalmax，可以用diff(w)完成等价的局部极大值判定，核心逻辑不变。

cand{k}里存放的是二维数组：第一列是时间位置，第二列是对应的模值。后续极值线追踪完全依赖这两个字段，所以提取阶段的过滤宁可偏严，也不要偏松，否则追踪代码会因为这些太密的候选点而连线错乱。

3.3 跨尺度极值线追踪

追踪的目标是把同一奇异点在不同尺度产生的模极大值连成线。常见做法是从最小尺度出发，向大尺度找数最近邻：当前尺度的候选点pos，在上一尺度所有候选点中找时间差最小的点，若差值小于允许范围就判定为同一条线。允许范围一般取0.2倍当前尺度对应的小波支撑宽度，或干脆固定为3到5个采样点：

dtmax = 0.2 * scales(k) * dt + 2 * dt; ids = zeros(size(cand{k}, 1), 1); for k = 2:length(scales) c1 = cand{k-1}; c2 = cand{k}; for i = 1:size(c1, 1) [mind, j] = min(abs(c2(:, 1) - c1(i, 1))); if mind <= dtmax id = assign_line(c1(i,:), c2(j,:)); % 把两点并到一条线 end end end

此处的assign_line是你要自己维护的线id分配函数；在正式工程里，可以用图论里的连通分量来做，避免多条线在赋值时互相覆盖。更简单的实现是把lines存成元胞数组，每条线保存posSeq和ampSeq，匹配成功就在线尾追加节点；如果当前候选点与任何已有线都不匹配，则新开一条线。追踪方向也可以从大尺度向小尺度反向进行，但对噪声信号，从小尺度出发更容易保持定位精度。

追踪完成后，一条有效的极值线至少要覆盖60%以上尺寸。这个比例是过滤噪声伪线的基本标准，具体取值在第4章结合参数展开。

3.4 用双对数拟合估计Lipschitz指数

每条线的posSeq对应的尺度系数已知，ampSeq由cand第二列得到。接下来只需要取对数后做线性拟合：

function alpha = est_lipschitz(scales, amps, kRange) p = polyfit(log(scales(kRange)), log(amps(kRange)), 1); alpha = p(1); end

调用例子：如果某条线在尺度索引5到12之间均有数据，就令kRange=5:12，得到该点的α估计值。为什么要舍弃两端的尺度？最小尺度的模值容易受噪声扰动，最大尺度的模值又可能被邻近奇异点通过极值线合并而干扰；中间段更接近“只受该奇异点控制”的理论区域，拟合稳定性最高。

实际操作时通常做两次过滤：第一次按覆盖率过滤线，第二次按拟合残差过滤。若拟合残差的均方根超过0.3，说明该线位置附近可能有两个互相影响的奇异点，或者小波变换数值精度有问题，需要回到第3.1节的延拓部分检查。拟合得到α后，按表3-1做物理含义映射。

表3-1 估计α与奇异类型的对应参考

4. 参数选择与噪声场景下的小波模极大值调节

4.1 尺度范围与尺度步长怎么定

尺度范围的选择同时影响计算量和检测精度。最小尺度决定定位精度：尺度1对应的gauss1核只覆盖约8个采样点，定位误差也在1-2个采样点之内。如果信号采样率是1000Hz，最小尺度取1就够；若信号本身带宽有限（例如传感器前端内置低通滤波），则需要把最小尺度抬高到2-4，避免检测出大量高频伪突变。

最大尺度随着噪声水平和信号长度的变化而变化。理论上尺度取到信号长度的1/4以上没有意义，因为此时整个信号基本都在小波支撑内，奇异点已无法在时域上区分。经验上取信号长度的1/16到1/8。表4-1给出典型的初始参数组合，实际使用时按信号能量分布微调最大尺度即可。

表4-1 不同采样率下的尺度范围建议

步长决定极值线匹配成功率。指数增长序列2^(a:b:c)中，b取0.5时相邻尺度之间小波核长度相差约1.4倍，核形状变化平缓，追踪稳定。b取1时计算量小，但相邻尺度间的模值可能突变，对噪声敏感。推荐先从b=0.5开始，若信号很短或者只需要做粗检测，再换到b=1。

4.2 阈值与显著性过滤

固定百分比threshold=0.1*max(|W|)在背景噪声很小时工作正常，一旦噪声水平升高，真实微弱奇异点的模值可能远低于背景峰值，固定百分比就直接漏检。更稳健的做法是从中位数估计噪声能量：

noiseLevel = median(abs(W(:))) / 0.6745; thr = 4 * noiseLevel;

0.6745是高斯噪声标准差与中位数绝对偏差的换算常数，4倍噪声标准差对应的虚警率已经很低。该方法不需要知道信号的真实信噪比，对工程应用更友好。如果信号中存在明显的有色噪声，中位数会被低频成分抬高，推荐改用每一尺度行单独计算中位数，得到每层不同的阈值。

噪声本身也会在尺度平面里形成极值点，但这些极值点的显著性随尺度增大迅速衰退。若某条极值线覆盖尺度数少于总数一半，其拟合出的α可信度非常低，建议直接丢弃。阈值和覆盖率两个参数要一起调节，只动一个往往压不住噪声。

4.3 伪极值线的判别：覆盖率与线长

追踪完成后按线覆盖的尺度数过滤：

valid = false(length(lines), 1); for i = 1:length(lines) valid(i) = numel(lines(i).amp) >= 0.6 * length(scales); end lines = lines(valid);

0.6这个比例需要结合尺度范围做微调。当尺度步长b=0.5、总尺度数是13时，真实极值线通常能覆盖全部13层；噪声线一般只覆盖2到4层。如果最小尺度或最大尺度设置过短，真线的覆盖率也会下降，因此判断标准以“至少覆盖总尺度数的60%”为准，不要直接写死为某个固定数量。

白噪声在尺度谱上还有一个明显特征：极值线分布均匀且密度高。过滤时可以把尺度平面的二维密度也纳入考量，利用histcounts2统计单位区域的线数，密度超标的区域整块剔除。这样处理对突发性强噪声更有效，因为突发噪声通常只集中在某几个尺度上。

4.4 低信噪比下的预平滑替代方案

当信噪比降到10dB以下时，只靠提高阈值往往压制不住噪声极值线，需要在对小波系数做极值提取之前先做时间方向的中值滤波。常见做法是对每一尺度行w做medfilt1(w, m)，m取该尺度小波主瓣宽度的一半。中值滤波在保边缘的同时能显著降低噪声峰，缺点是小尺度端的极值位置会偏移几个样本点，定位精度会略微下降。

另一种方案是提高阈值的同时放宽dtmax。噪声让极值线位置抖动，如果追踪窗口过紧，真实极值线会被切碎成若干短线。此时可以把dtmax从固定3个采样点调整为0.4*scales(k)dt+3dt。放宽追踪窗口的代价是邻近两个奇异点的极值线可能被误并为一条，解决方法是检查线内模值的平滑性：若同一条线内模值出现非单调反复，则拆线重连。

5. 验证技巧：用合成信号检验小波模极大值的整条链路

5.1 构造三类标准奇异点

写完代码后，第一件事不是跑真实数据，而是用理论α已知的三类信号验证。阶跃点的Lipschitz指数理论为0；脉冲点理论接近-1；斜坡突变的理论在0.5附近。构造这样一支信号后运行前述函数，比对检测位置与估计α。

t = (0:999) / 1000; x = zeros(1, 1000); x(400) = 1; % 脉冲 x(600:end) = x(600:end) + 1; % 阶跃 x(800:end) = x(800:end) + (0:199) / 200; % 斜坡切变

脉冲在单点注入能量，阶跃从第600点开始持久改变电平，斜坡则在第800点开始引入斜率。三者的理论α分别接近-1、0和0.5，足够覆盖多数故障信号中的突变形态。构造完成后分别调用极值提取和追踪函数，观察三条极值线是否都出现。

5.2 与理论指标做定量对照

验证时要看的三个指标：定位误差、α估计误差、极值线覆盖率。定位误差应小于最小尺度下的核半径；α误差在±0.15以内算通过；覆盖率应接近100%。任何一个不合格都优先回溯到边界延拓与阈值设置，而不是直接改拟合代码。把这三个指标打印出来作为回归基准，之后每改动一次参数，跑一遍合成信号对比结果，能快速发现是参数还是实现出了问题。若定位误差持续偏大，可以在极值线追踪后用一次一维抛物线插值对峰值位置做亚像素修正，这是不增加计算量就能提升定位精度的常用手段。

5.3 按qiyidian.rar的组织习惯拆分函数

标题里的qiyidian.rar让人联想到一个现成的MATLAB包。不管解压后原始文件怎么命名，按“数据读取-尺度变换-极值提取-线追踪-指数估计”拆成五个函数去做，以后换数据、换小波基时只需改对应层。例如换用gauss2，只需要替换尺度变换里的模板函数，其他模块完全不动。推荐的调用顺序：

W = wt_sgram(x, scales, Fs); cand = wt_findmax(W, 'MinProminence', 0.05); lines = wt_track(cand, scales, 'dtmax', 0.3); result = wt_lipschitz(lines, scales, 'MinCoverage', 0.6);

参数名在函数接口里显式暴露，比在脚本里写散落的魔数更容易维护。把合成信号验证脚本放到最前面，作为一个独立的testSmoke.m，形成每次改代码后必跑一次的回归习惯。这样即使后面接入的工程信号形态复杂，也能保证小波模极大值这层核心算法始终处于可信状态。

本文还有配套的精品资源，点击获取