简介一套面向信号处理与工程应用的小波模极大值程序聚焦小波降噪与模态参数识别适合需要借助MATLAB开展小波分析的研究人员、工程师及相关专业学生。程序包内含6个文件以5个.m脚本和1个txt说明文档为主覆盖小波系数计算、极大值点提取、信号重构与模态参数识别等核心环节便于直接运行与二次修改。文件体积仅4KB结构精简适合快速理解小波分析中阈值去噪和特征提取的实现思路。小波模极大值对应各尺度上的关键系数通过该算法可有效定位信号突变与特征成分是处理非平稳信号的有力工具。目前已有495人学习下载资源附带的说明文档和脚本注释可辅助初学者掌握Haar、Daubechies等小波基的选择与阈值策略并进一步用于科研实验或工程项目中的数据清洗与动态特性分析。1. 小波模极大值把信号突变点从噪声里捞出来的程序化方法“小波模极大值”这个名字听起来像某个数学库里的冷门函数实际上它解决的是很具体的问题在一长段信号里定位突变时刻。工程场景中突变往往是最有价值的信息——旋转机械的故障冲击、电力系统暂态、心电信号里的早搏。常规做法是把信号和小波做卷积后直接看系数但系数会被噪声淹没阈值也分不清“边缘”和“毛刺”。模极大值提供了一条跨尺度传播的线索真实突变的大系数在多个尺度上连续重现噪声则快速断链。下面先把数学原理压到最小可用程度再直接给出一个能跑的小波模极大值提取程序把尺度、阈值、链长这三个关键旋钮调明白。适合正在做故障诊断、信号检测或边缘提取的工程师。2. 小波模极大值的判定原理为什么看“模的极大值”比看系数大小可靠2.1 连续小波变换与高斯一阶导小波对连续信号 f(t)取尺度参数 s 0、位移参数 b小波 ψ(t)连续小波变换写作Wf(s,b) (1/√s) ∫ f(t) ψ*((t-b)/s) dt当母小波选用高斯函数 θ(t)e^{-t²/2} 的一阶导数 ψ(t)-t·e^{-t²/2} 时变换结果等价于先对信号做高斯平滑、再求一阶导数。原因是卷积与微分可交换次序∂/∂b f * θ_s f * (∂θ_s/∂b) 而 ∂θ_s/∂b 正比于某个尺度下的 ψ。这个等价关系正是模极大值方法的基础小波系数处处反映信号在尺度 s 下的局部斜率斜率最陡的地方|Wf| 取极大就是疑似突变点。pywt 里把这一族小波命名为 gaus1、gaus2、gaus3…gaus1 就是一阶高斯导数gaus2 是二阶。做模极大值检测用 gaus1 最常见它把阶跃或脉冲突变表现为单个正峰或负峰取绝对值后是一个单峰容易定位。如果用 morl 这类振荡型小波突变位置会留下多个正负交替的极值跨尺度连接时会出现成对的假链定位前还要额外做极值配对得不偿失。2.2 Lipschitz 指数模极大值随尺度变化的斜率模极大值之所以能区分信号与噪声依据是信号局部奇异性的衰减规律。若 f 在 b0 处是 Lipschitz α 的则存在常数 K使得|Wf(s, b0)| ≤ K · s^α两边取以 2 为底的对数log2|Wf(s, b0)| ≤ log2K α·log2s这意味着把同一个极值点在不同尺度上的模值画在 log-log 坐标里拟合斜率就是 α。阶跃信号 α≈0光滑可导信号 0α1白噪声 α0。这个定量关系直接变成程序里的筛选判据一条“真链”的模值跨尺度应基本持平或缓慢上升噪声链则快速衰减。信号局部形态Lipschitz α模极大值随尺度的走向阶跃、边缘α≈0基本持平链延伸到最大尺度光滑信号中缓慢变化0 α 1缓慢上升链长白噪声、高频毛刺α 0快速下降链短且位置随机白噪声在最小两个尺度上的模值可能比真实突变还大这就是单看某一尺度系数幅值不可靠的原因而一旦把尺度轴拉起来噪声链的斜率立刻露出马脚。2.3 为什么不能直接软阈值替代软阈值去噪在系数层面操作前提是真实系数比噪声系数幅值大。但阶跃点的系数随尺度增大而增大在最小尺度上它和小尺度噪声重叠固定阈值要么放过噪声、要么削掉阶跃的位置精度。模极大值把检测从幅度问题变成几何问题候选点在连续多个尺度上能否成为一条极值线。幅度只是第一道门槛跨尺度存活才是决定性判据。检测任务根本不需要把信号重构出来因此这里走的是“定位优先”的路子。需要完整恢复波形时再考虑 Mallat 交替投影重构那是另一套程序流程本文只解决“突变在哪儿”以及“几个突变”的问题。3. 用 Python 实现一个小波模极大值最小程序3.1 用 pywt.cwt 算连续小波变换并扫描模极大值先给一个可以直接跑的 demo 程序。它不做链管理先把每层尺度的候选极值点扫出来。import numpy as np import pywt def build_scales(s_min, s_max, num_scales): # 按2的幂对数均匀取尺度小尺度分辨率密大尺度稀疏 return 2 ** np.linspace(np.log2(s_min), np.log2(s_max), num_scales) def modmax_scan(signal, fs, s_min2, s_max128, num_scales24): scales build_scales(s_min, s_max, num_scales) coef, freqs pywt.cwt(signal, scales, gaus1, sampling_period1.0 / fs) mod np.abs(coef) # 取模正负极值统一成峰 cand [] # 每个尺度一条候选位置列表 for j, row in enumerate(mod): pos [] for i in range(1, len(row) - 1): if row[i] row[i - 1] and row[i] row[i 1]: pos.append(i) cand.append(pos) return scales, freqs, mod, cand逻辑说明pywt.cwt 对每个尺度输出一条与信号等长的系数序列row[i] 同时大于左右相邻点即为局部极大。这里只做相邻点比较采样率很低时可能把平台误判成峰后续应叠加幅值阈值。参数上s_min2 表示小波时窗约覆盖信号最短周期的两倍num_scales24 在默认配置下对几千点长度的信号运行时间可以接受。freqs 由 pywt 根据采样周期自动换算画图时标定纵轴用。3.2 跨尺度链连接把零散候选点连成“根系”单层候选点毫无意义真正的模极大值程序要把相邻尺度的候选点连接成链。连接规则基于一个物理事实尺度越大小波时窗越宽极值点在小波有效支撑内的位置漂移越明显所以容差半径要随尺度放大。def link_across_scales(cand, scales, base_shift3): links [] for j in range(1, len(scales)): if not cand[j] or not cand[j - 1]: continue # 容差随尺度线性放宽尺度越大位置漂移越大 radius max(1, int(base_shift * (scales[j] / scales[0]))) for p in cand[j]: dist [abs(q - p) for q in cand[j - 1]] if min(dist) radius: q cand[j - 1][int(np.argmin(dist))] links.append((j - 1, q, scales[j - 1], j, p, scales[j])) return links逻辑说明每次只把当前尺度与上一个较小尺度“配对”返回的元组包含上一层位置和尺度、当前层位置和尺度。这个简化版不做链 ID 管理适合先观察连接关系。逻辑上base_shift3 表示最小尺度处的容差为 3 个采样点最大尺度处按比例放大。完整工程上要给候选点挂 chain_id统计每个 ID 出现的次数出现次数达到阈值才保留下一步筛选这个链长统计放在后续流程里做。3.3 可视化让模极大值线直接可检视没有图参数调整就是瞎猜。下面的函数把系数矩阵画成热力图白色折线即跨尺度连接线。import matplotlib.pyplot as plt def plot_modmax(t, mod, scales, linksNone): fig, ax plt.subplots(figsize(10, 5)) im ax.imshow(mod, aspectauto, originlower, extent[t[0], t[-1], scales[0], scales[-1]], cmapjet) ax.set_yscale(log) ax.set_xlabel(time (s)) ax.set_ylabel(scale) fig.colorbar(im, axax) if links: for _, q, s1, _, p, s2 in links: ax.plot([t[q], t[p]], [s1, s2], w-, lw0.5, alpha0.7) return fig看图的重点不是颜色深浅而是白色折线是否从最小尺度一路延伸到最大尺度。噪声链通常只有两三条短线信号链则像树根一样扎进大尺度区域。如果某条链在中途断开先怀疑 num_scales 太少或 base_shift 太小其次才怀疑是噪声。关于小波族的选型这里整理一个对比方便以后换场景小波突变显示定位精度适用gaus1单峰好阶跃、脉冲冲击首选gaus2正负双峰中平滑起伏中的拐点db4多峰中需要近似正交基的场景morl振荡多峰差瞬时频率分析不做模极大值定位4. 小波模极大值程序的四个必调参数4.1 尺度范围 s_min 与 s_max先圈定目标频带s_min 越小越像一阶差分算子高频量化噪声和测量噪声会被放大。若目标故障脉冲的脉宽约为 20ms采样率 1kHzs_min 取 816 比较合理s_min2 更适合检测极窄的瞬态尖峰。s_max 决定了能确认的最慢变化。s_max 过大时小波时窗覆盖的目标早超过信号中的最低频成分而且边界效应从两端大范围侵入。常见做法是 s_max ≤ N/4N 为信号长度信号只有 2000 点而 s_max512 时几乎整条纵轴都被边界污染链长统计失去意义。先定 s_max再定 s_min中间隔两三个倍频程即可。4.2 尺度分辨率与链接半径链是否连续的真正决定因素num_scales 太少相邻两尺度的小波时窗差异过大极值点漂移超出 radius链断成一截截短线。经验上 num_scales 取 1632信号很长可到 64。判断依据是可视化图中纵向是否有明显断裂。radius 需要与尺度序列匹配。base_shift1 时链往往断在较深处因为大尺度极值位置偏移容易超过 1 个采样点base_shift3 通常能把信号链接上。但 base_shift 过大也会有副作用两条空间接近的噪声链被错误并成一条图上表现为白色连线在一个时刻附近横跨一大段空白。建议先 base_shift1 看断链情况再逐步加到 3不要一步到位。提示如果两条真实事件的时间间隔很小而 base_shift 又偏大它们会被连成一条链。此时宁可接受部分断链也不要过连接。4.3 幅值阈值与链长筛选拦噪声、保事件幅值阈值作用在每层候选扫描之后row[i] th 才保留。th 不能全信号统一因为大尺度系数整体偏大正确做法是逐尺度计算th_j median(|coef_j|) 3·MAD_j其中 MAD 为该尺度系数的绝对中位差。用中位数和 MAD 而不是均值和标准差是为了抵抗脉冲尖峰对统计量的拉偏。链长筛选指标 min_chain_len 表示候选点连续出现在多少层尺度上。信噪比高时取 56信噪比低时取 3低于 3 基本等于没筛。再加一道跨尺度斜率校验对链上 log2|W| 对 log2 s 做线性回归斜率 α0 的链即使长度达标也优先剔除。参数推荐范围偏低后果偏高后果s_min28噪声多误检多漏掉窄脉冲s_maxN/8N/4高频信息丢失边界污染严重num_scales1632链断裂耗时、过连接base_shift15链断裂相邻事件错连thmed 3MAD假检测多弱事件漏检min_chain_len36噪声混入弱事件漏检4.4 边界长度与延拓方式第四个要调的量CWT 在信号两端因数据不足产生假极值。常见做法是先用 reflect 对称延拓信号跑完 CWT 后裁掉左右各 int(s_max) 的区域。延拓长度不是固定值它随 s_max 线性增长s_max128 时两端各有约 128 个点不可信。如果目标突变恰好落在前 128 点内无论如何调参都会看到一条贴边生长的假链那不是程序写错了。注意链长统计必须在裁边之后进行。边界假链链长往往相当可观先统计后裁边会把假链混进结果。5. 用模极大值程序定位轴承故障冲击并校准阈值参数5.1 从原始采集数据到冲击时刻表实际数据处理流程我一般按六步走去均值按先验频段做带通滤波。注意滤波阶数选择线性相位 FIR避免突变点被相位失真平移。根据目标冲击宽度定 s_min根据信号长度定 s_max调用 build_scales 生成尺度序列。调用 modmax_scan 得到 mod 和 cand。逐尺度按 med 3MAD 过滤候选点。调用 link_across_scales再把相邻层配对用字典合并成链统计每链出现的尺度次数。保留链长 ≥ min_chain_len 的链对链上位置按尺度加权平均权重取 1/s小尺度权重更高得到最终时间戳。第 5 步是链管理合并规则是若当前层候选 p 与上一层候选 q 配对成功且 q 已属于某链则 p 也归属该链若 q 尚未归属任何链则新建一条链。这是典型的并查集问题不需要递归。5.2 一组已知突变位置的验证样例拿到程序后先别直接上实测数据构造一个已知答案的信号fs1000时长 1 秒0.300s 处加阶跃0.600s 处加窄脉冲再叠加 SNR 约 10dB 的高斯白噪声。运行完整流程后把检出时刻与真实时刻比较误差不超过 5ms即 5 个采样点才算命中。记录形式参考下面模板试验真实时刻检出时刻偏差阶跃0.300s0.301s1ms窄脉冲0.600s0.598s-2ms同时统计误检数检出的链里没有任何真实事件对应的链占总链数比例。误检率超过 10% 时优先提高 min_chain_len漏检率超过 10% 时优先放宽幅值阈值而不是链长。调参顺序固定为“先尺度范围再链长最后幅值阈值”否则两个指标来回跳。5.3 进阶技巧先画多尺度乘积图再看模极大值线多尺度乘积不是模极大值程序的替代而是它的体检工具。思路是真实突变在每个尺度上都超过噪声水平把连续 k 个尺度的“过阈值程度”相乘噪声被随机压低真实事件处乘积接近 1。def multiscale_product(mod, scales, k3, th_fac3): flags [] for row in mod: mad 1.4826 * np.median(np.abs(row - np.median(row))) th th_fac * mad flags.append((np.abs(row) - th).clip(0) / th) flags np.array(flags) prod np.ones(flags.shape[1]) for j in range(len(scales) - k 1): window flags[j:j k] prod np.maximum(prod, window.prod(axis0)) return prod逻辑说明每层先把超过 3 倍 MAD 的部分映射成 01 的值再做滑窗乘积最后对所有窗口取最大值。真实事件的乘积曲线峰值尖锐噪声区域基本贴零。实现时注意 prod 用 np.maximum 保留历史窗口最大值避免多个小于 1 的数连续相乘把值压到接近 0。拿到乘积曲线后再看模极大值图曲线峰值位置对应的白色折线应该整齐跨越大半个尺度范围而噪声链集中在乘积接近 0 的时间带。我调参时先看这幅曲线再决定 min_chain_len 取 3 还是 5——峰值清晰就敢取 6峰值被噪声埋掉就先提高 s_min 过滤毛刺重新扫描。把这两个工具叠在一起小波模极大值程序才算真正从“能出图”变成“能出结论”。本文还有配套的精品资源点击获取