小波模极大值定位信号突变点：原理与MATLAB实现

简介：这是一套基于MATLAB的小波模极大值信号处理示例，面向需要学习小波分析、信号突变点检测与特征提取的科研人员或高年级学生。资源围绕cgau连续高斯小波展开，通过qiyidian.m脚本演示了从信号小波分解（wavedec）、模极大值定位、阈值筛选到利用waverec重建的完整流程，脚本注释清晰、结构紧凑，能让读者直观理解小波系数绝对值最大处与信号突变点的对应关系，掌握模极大值方法的参数设置思路，适合作为小波分析入门实践及后续算法改造的蓝本。压缩包共1个文件，仅1个m脚本，大小563B，代码量小、便于阅读修改。目前已有227人学习下载；读者拿到后可直接运行或逐行调试，并可结合地震、医学图像、故障诊断等场景进一步迁移应用，是快速上手小波模极大值的实用范例。

1. 小波模极大值能定位到哪些信号突变点

工程师拿到振动波形、心电数据或电力暂态信号时，第一件事往往不是滤波，而是问：跳变发生在哪个时刻？有多陡？傅里叶变换只能给出频率成分，小波变换却能同时给出时间和尺度，小波模极大值（WTMM）则更进一步，只盯住绝对值最大、尺度间能连成线的极值点，这些点对应信号奇异位置。qiyidian.rar里的qiyidian.m脚本，演示了从信号分解到奇异点定位的完整链路。这套方法适合故障诊断、生物医学信号处理、结构健康监测场景的工程师，特别是拿着现场波形却不知道突变点落在哪个采样点的那类需求。后面的章节会从选基原理、MATLAB实现、噪声判据到亚采样级时间戳依次展开。

2. 小波基函数与模极大值的数学定义

模极大值方法的第一步是选小波基。很多新手直接用默认的db3跑DWT，虽然也能得到系数，但突变点在低频段的定位会明显偏移。造成偏移的原因不是DWT本身出错，而是离散变换的尺度间隔太粗，无法连续追踪极值点随尺度的移动轨迹。这一章从原理出发，把选基和解算路径说清楚。

2.1 CWT与DWT在极值追踪上的关键差异

离散小波变换通过wavedec把信号分解为不同频带的系数。它的计算效率极高，但尺度按2的幂次跳变，相邻尺度之间的系数位置相差两倍，不方便追踪极值点沿尺度变化的轨迹。连续小波变换（CWT）则把尺度变量离散成连续数组（比如1、1.5、2、2.5...），每个尺度都保留完整的时间轴信息，模极大值点可以在尺度方向上相连成线，也就是所谓的模极大值脊线。

在MATLAB中，R2016b之前的cwt函数直接支持小波名称和尺度数组，之后版本改成了基于频率的参数输入。下面两种写法在工程中都很常见：

scales = 1:0.5:32; coefs_old = cwt(signal, scales, 'cgau4'); [wt, freqs] = cwt(signal, 'amor', fs);

第一段里的'cgau4'是连续高斯小波的第4阶导数。导数阶数对应小波的消失矩：cgau1有一个符号变化，适合识别阶跃；cgau4有四个符号变化，适合识别更细微的奇点。第二段返回的wt是复Morlet小波系数，模值和相位同时存在，适合时频分析，但在模极大值检测任务里相位信息通常用不上。如果信号段较长，优先用第二段写法，频带划分更细；但一个尺度一个尺度追踪极大值时，第一段写法更直观。

2.2 模极大值的判定条件与实现

模极大值的数学条件：固定某个尺度scales = s，对平移参数b求偏导并令其为零——等效于在小波系数绝对值曲线上寻找局部最大值。用代码实现时，问题被转换成计算绝对值的离散局部极大值：

% 假设coefs是尺度个数 x 信号长度的二维矩阵 num_scales = size(coefs, 1); maxpos_cell = cell(num_scales, 1); for i = 1:num_scales abs_curve = abs(coefs(i, :)); tf = islocalmax(abs_curve); % 用峰高阈值剔除弱峰 thr = 0.05 * max(abs_curve); tf = tf & (abs_curve >= thr); maxpos_cell{i} = find(tf); end

这段代码在每个尺度上独立找到所有局部极大值点。islocalmax默认把两端边界视为非极大值，这使得信号起始和结束位置的突变点会丢失。如果需要保留边界处的极值，可以在处理前给信号两端各延长一段数据，计算完再截掉对应位置的极值点。延长的长度建议不小于小波有效支撑长度的一半，这样边界效应不会污染目标区间。

求完极大值后，还需要把不同尺度上位置相近的点连成脊线。下面给出一个最简的贪婪匹配逻辑，实际脚本可以在此基础上增加幅值约束，避免把幅值衰减过快的点强行连入脊线：

ridge_len = num_scales; ridge = zeros(1, ridge_len); ridge(end) = maxpos_cell{end}(1); % 从最大尺度开始 for k = num_scales-1:-1:1 cand = maxpos_cell{k}; [~, idx] = min(abs(cand - ridge(k+1))); if abs(cand(idx) - ridge(k+1)) <= 5 ridge(k) = cand(idx); else ridge(k) = NaN; % 断点 end end

位置容差设置为5个采样点，对1000 Hz采样率的信号来说大约是5毫秒，足以滤掉大部分偶然接近的噪声峰。容差太大会把两条平行脊线误并成一条，取值建议以信号最小时间尺度的1/3为上限。如果多个尺度上的极值点同时满足容差条件，优先选择幅值更接近前一个尺度极大值的候选，这样能保证脊线的幅值连续。

2.3 小波基参数对结果的影响

实践中的参数选择可以按信号类型做一个对照：

选择依据是突变形状越复杂，需要越高阶的小波，计算量也随之增大。尺度步长从1改成0.5时，极值点的定位精度大约提升30%，但计算时间接近翻倍。做故障报警场景，粗尺度就够了；做特征分类任务则需要在细尺度上把相邻两个突变点分开，这时步长和容差都要收紧。还有一条经验：先跑一次不加阈值的结果，看看极大值数量级，再根据目标数量反推最小峰高，比拍脑袋设阈值容易收敛。

3. qiyidian.m的完整流程：wavedec分解到模极大值定位

qiyidian.m这个命名取自“奇异点”谐音，脚本核心任务就是定位奇异点。分析脚本功能结构，可以分成三个环节：先选择小波基并配置分解层数，再完成多层小波分解，最后逐层提取细节系数的模极大值并按层合并成奇异点集合。

3.1 从wavedec到waverec的信号分解与重构

小波分解的函数是wavedec：

load leleccum; % 内置含噪信号 signal = leleccum(1:1024); N = 3; wname = 'db4'; [C, L] = wavedec(signal, N, wname); a3 = wrcoef('a', C, L, wname, 3); d1 = wrcoef('d', C, L, wname, 1); d2 = wrcoef('d', C, L, wname, 2); d3 = wrcoef('d', C, L, wname, 3);

wavedec返回值C是所有系数的拼接向量，L是一个长度N+2的数组，记录各层系数个数。wrcoef通过重构算法把细节系数恢复成与原始信号等长的序列，这是为了在相同时间轴上比较不同层的极值位置。若直接用detcoef会得到长度逐层折半的系数段，时间坐标需要换算成乘2的层数次方，定位误差也随层数累积。

如果要用模极大值重构信号，通常不会直接用waverec，而是先对系数做屏蔽再重构。屏蔽操作的代码如下：

C_shrunk = C; nonMaxPos = setdiff(1:numel(C), targetPosInC); C_shrunk(nonMaxPos) = 0; sig_recon = waverec(C_shrunk, L, wname);

这里的targetPosInC是模极大值位置映射到C向量中的索引。这样重构出的信号仅保留突变部分，相当于完成了去噪和压缩，重构信号与原始信号的差值就代表被滤除的平稳分量。做这一步的前提是必须清楚C向量中各层系数的索引范围，wavedec的系数排布顺序是近似系数在最前、细节系数按层递增排列，需要仔细对齐。下表给出了3层分解时C向量结构的参考信息：

3.2 极值点搜索与时间段合并

逐层搜索极值的循环代码在上一章给出，这里讨论搜索之后如何处理。不同尺度搜索出的极值点数量不同，尺度越小数量越多。为了得到稳定的奇异点，我把三层极值点合并成时间段，条件是它们在时间轴上相邻不超过3个采样点：

merged_points = []; for k = 1:min([numel(d1_pos), numel(d2_pos), numel(d3_pos)]) base = d1_pos(k); if any(abs(d2_pos - base) <= 3) && any(abs(d3_pos - base) <= 3) merged_points(end+1) = base; end end

这个条件比直接intersect宽松，因为各层系数经过重构后虽然时间轴对齐，但边缘振铃会使极值点偏移一到两个采样点。合并三段位置后得到的点集是候选奇异点，数量通常只有极值点总数的1/5不到，这也是模极大值方法用于信号压缩时的原理之一。处理低频漂移明显的信号时，合并窗口要放宽到5~7个点，否则基线漂移会产生额外极值点，让合并结果偏多。

3.3 脚本复现和相位对齐验证

拿到qiyidian.m文件后，先不要急着运行。打开脚本找到小波变换语句，确认使用的是cwt还是wavedec，再检查时域波形图绘制是否使用subplot。添加输出代码将极大值位置写到变量mmp_pos，然后比较其与原始信号波形的对应关系：

highlight_signal = zeros(size(signal)); highlight_signal(merged_points) = signal(merged_points); plot(signal); hold on; stem(merged_points, signal(merged_points), 'r');

如果红点都落在波形跳变处，说明极值定位有效。若大量红点出现在平滑区域，则说明阈值过低或小波阶数选择不当，应从层数和阈值两方面调参。如果红点整体延迟了数个采样点，问题多半出在cwt的边界填充方式上，将信号翻转后再跑一次对比即可确认。整个脚本跑完后，把每一层wrcoef输出保存成mat文件，便于对比不同层数对定位精度的影响。

4. Lipschitz指数与阈值策略：从噪声中筛选真实模极大值

模极大值集合中，真实奇异点的模极大值随尺度增大保持稳定，而噪声产生的模极大值随尺度增大快速衰减。这一点在时域表现为噪声峰被平滑掉，在小波域表现为系数幅值急剧下降。Lipschitz指数把这个视觉观察量化成可计算的数值，是区分真实突变和噪声伪峰的比较严格的判据。

4.1 基于小波系数的Lipschitz指数估计

信号在x0处的Lipschitz指数alpha满足：

log2|Wf(2^j, x0)| <= log2K + j*alpha

实际操作时，选取某一条脊线上j从1到J的系数绝对值，计算log2之后做线性回归，回归斜率就是alpha的估计值。MATLAB中的实现步骤：

jvals = 1:J; logcoef = log2(abs(ridge_coefs(1:J))); pfit = polyfit(jvals, logcoef, 1); alpha = pfit(1);

polyfit返回的一次项系数pfit(1)就是斜率。若alpha在0附近，说明该点是阶跃或陡坡；若alpha为负，说明它是一个白噪声主导的孤立值，应从特征集合里剔除。alpha大于0.5时，该点接近光滑区域的缓变段，并非严格意义上的奇异点，在特征工程里需要单独标记。回归时要注意：如果脊线上某几个尺度出现NaN（脊线断裂），对应的点不能参与拟合，需要先用rmmissing清洗。

4.2 无标度区间的选择：拟合范围决定指数精度

回归效果取决于参与拟合的尺度范围。尺度太小的系数受信号原始波形影响大，尺度太大的系数受两个不同奇异点之间的相干叠加影响大。惯例做法是取最大尺度的前1/2或1/3区间参与拟合。

fit_scales = 1:floor(max_scales/2); P = polyfit(fit_scales, logcoef(fit_scales), 1); resid = sum(abs(logcoef(fit_scales) - polyval(P, fit_scales)).^2);

残差resid衡量拟合的线性程度。如果一个奇异点真正是单点突变，残差很小；如果是多个奇点距离过近，相邻脊线的系数叠加会让残差变大。这时应将该点标记为“复合奇点”，保留位置但弃用其alpha值进行特征建模。还有一点需要提醒：对于只有3层分解的信号，拟合只有3个点，回归结果非常不稳定，建议至少做5层分解再估计alpha。

4.3 自适应阈值调参表

我在真实噪声环境中的调参过程通常参考下面这张表，脚本调试时可以直接套用：

工具栏的maxlocw函数可以自动搜索极大值位置，但各版本输入输出格式略有不同，使用前先查帮助文档，确认是否需要传入阈值参数。maxlocw适合生产环境批量处理，qiyidian.m中的极值搜索代码则适合出图和分析。两种方式输出结果一致时，再确认参数没有问题。

5. 模极大值脊线串联与亚采样级奇异点时间戳

最后一个具体技巧是：将散落在各尺度的模极大值点连成脊线，并利用插值把奇异点定位精度提高到亚采样级别。这是qiyidian.m这类脚本落地到数据采集系统时最常被追问的功能。

5.1 从粗尺度到细尺度的脊线回溯匹配

脊线匹配从最大尺度开始向小尺度回溯。在最大尺度选择幅值最强的若干个极大值作为种子点，然后逐层向下，在每个较小尺度上查找与种子点最近的极大值位置。位置差距不超过设定容差则连入脊线，超过则不连接：

seed_pos = maxpos{num_scales}(1:3); ridge_bank = cell(numel(seed_pos), 1); for r = 1:numel(seed_pos) cur = seed_pos(r); chain = cur; for s = num_scales-1:-1:1 [mind, idx] = min(abs(maxpos{s} - cur)); if mind <= tol cur = maxpos{s}(idx); chain = [chain, cur]; else break; end end ridge_bank{r} = chain; end

匹配结束后，脊线长度不一致是正常现象。短脊线多是噪声，长度超过总尺度数2/3的脊线锚定真实奇异点。对该类脊线，取最小尺度端的位置作为奇异点候选，并在原始信号上做相位校正。这里有个容易忽略的地方：种子点选在前几个极值点时，如果信号里同时存在两个幅值接近的突变点，种子顺序会影响后续匹配，建议按幅值降序排列后再取前几个种子。

5.2 三点插值把定位精度推进到亚采样级

定位精度的提升来自于对极值点附近的系数曲线做插值。选择极小值附近三个相邻点做三次样条插值，取插值曲线的峰值坐标：

idx0 = coarse_pos; % 极值点的整数坐标 seg = abs(coefs(min_scale, idx0-1:idx0+1)); xq = -1:0.1:1; vq = interp1(-1:1, seg, xq, 'spline'); [~, jmax] = max(vq); fine_shift = xq(jmax); fine_pos = idx0 + fine_shift;

fine_shift可正可负，指示极值点真正的顶点位于整采样点的哪一侧。对长信号来说，亚采样级定位能显著减少连续事件之间的时间间隔测量误差，尤其在脉宽本身小于两个采样周期的场合，这种细化直接决定测到的脉冲起始时刻是否可重复。插值前最好检查seg的单调特性，如果三个点里有两个值相等，峰值位置落在中间点，插值结果稳定。调整xq步长时注意，步长缩小一半，插值点数就翻倍，超声回波信号常用0.05步长以获得ns级分辨率。

5.3 对位检查与特征表输出

工程中每完成一轮提取，都建议做一次对位检查。用细化的时间戳序列与原始信号差分结果的过零点对比，检查是否一致。如果差异超过用户定义的时差容忍区间，多半是阈值过低引入了伪极大值。此步完成后，把时间戳、尺度、alpha值写入文件：

T = table(event_id, fine_pos, alpha, 'VariableNames', {'EventID', 'TimeStamp', 'Alpha'}); writetable(T, 'singularity_features.csv');

生成的特征表可以直接进入下游分类器或统计模型。EventID要保持稳定，同一位置在不同分解参数下应映射到相同的ID，这需要在前端做一次去重合并。将时间戳按采样率换算成秒并保留4位小数，可以把亚采样级定位的精度完整保留在CSV里，避免后续读取时被截断。

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