第十二篇：智能峰值/谷值检测算法详解
在前几篇中，我们专注于语音信号的预处理和预测延拓，但实际应用中我们常需要从波形中自动定位极值点——例如基音周期中寻找波峰、频谱中识别谐波峰值、或心跳信号中检测 R 波峰值。MATLAB 自带的 findpeaks 函数（需 Signal Processing Toolbox）功能强大，但有时我们需要更轻量、更可控的实现。findpeakm.m 提供了完全自包含的峰值/谷值检测算法，支持二次插值亚样本精度、宽度容差去重，并且无需任何工具箱依赖。本篇将深入剖析其核心算法，从一阶差分定位到插值优化，再到局部峰值去重，带你全面掌握这个实用工具。
峰值检测的本质是在离散序列中找到局部极大值点。一个点 x(k)
是峰值，当且仅当其左侧相邻点小于或等于它，且右侧相邻点也小于或等于它（对于严格的峰值，两侧严格小于）。但实际信号常伴有噪声，直接比较容易出现虚假极值。
基本流程：
findpeakm.m 完整实现了上述逻辑，并支持 ‘q’ 插值模式和 ‘v’ 谷值模式。
function [k, v] = findpeakm(x, m, w)
% 输入：x – 信号向量
% m – 模式字符串：'q' 表示二次插值，'v' 表示寻找谷值
% w – 宽度容差（样本数），若两个峰值距离 ≤ w，则剔除较低的
% 输出：k – 峰值位置（若 'q' 模式则为浮点数）
% v – 峰值幅值
主要步骤：
· 将信号转为列向量，若为谷值模式则取反；
· 计算差分 dx = x(2:end) – x(1:end-1)；
· 找到上升段和下降段的索引（dx>0 为上升，dx<0 为下降）；
· 利用上升/下降段的起始和结束位置，确定峰值所在位置（处理平坦区）；
· 若为 ‘q’ 模式，用二次插值修正位置和幅值；
· 若指定了 w，则进行邻近峰值去重；
· 若为谷值模式，将幅值取反回来；
· 若无输出参数，则自动绘图。
3.1 计算差分并标记上升/下降
dx = x(2:end) – x(1:end–1);
r = find(dx > 0); % 上升段的起始索引（指向上方点的索引）
f = find(dx < 0); % 下降段的起始索引
dx(i) = x(i+1) – x(i)。若 dx(i) > 0，说明从 i 到 i+1 是上升趋势；若 dx(i) < 0，则为下降。
3.2 计算相对于上升和下降的“时间距离”
这段代码是 findpeakm 的精髓，它通过累计“自上次上升/下降以来的样本数”来定位峰值。
dr = r;
dr(2:end) = r(2:end) – r(1:end–1);
rc = repmat(1, nx, 1);
rc(r+1) = 1 – dr;
rc(1) = 0;
rs = cumsum(rc); % rs 向量：每个样本点距离最近一次上升点的样本数
类似地，fs 计算距离最近一次下降点的样本数。
3.3 确定峰值候选位置
峰值应满足：
· 它离最近一次上升点很近（rs < fs）：说明这个点正处于上升之后；
· 它离最近一次下降点也很近（fq < rq）：说明它即将转为下降；
· 并且 floor((fq – rs)/2) == 0：这个条件确保了峰值位于平坦区的中心（若存在平坦段）。
k = find((rs < fs) & (fq < rq) & (floor((fq – rs)/2) == 0));
v = x(k);
在没有平坦区时，fq – rs 通常为 1，此时 floor(0.5)=0，满足条件。如果出现一个平台（plateau），例如 [1, 2, 2, 1]，该条件会将峰值定位在平台的中心位置（第 2 或第 3 个点，取决于 floor 的行为）。
当信号峰值不是恰好落在采样点上时，我们可以用抛物线拟合三个相邻点（左、候选、右）来精确估计峰值的真实位置和幅值。
设候选点索引为 k
，其幅值为 x_k
，左右点为 x_{k-1}
和 x_{k+1}
。构造二次多项式：
f(t) = a t^2 + b t + c
令 t=0
对应候选点 k
，则：
· f(0) = c = x_k
· f(-1) = a – b + c = x_{k-1}
· f(1) = a + b + c = x_{k+1}
解得：
a = \\frac{x_{k-1} + x_{k+1}}{2} – x_k
b = \\frac{x_{k+1} – x_{k-1}}{2}
极值点位置（相对于候选点）为 t_{\\max} = -\\frac{b}{2a}
，对应的幅值为：
f(t_{\\max}) = x_k – \\frac{b^2}{4a}
当 a > 0
时，抛物线开口向上，那是极小值，但峰值处应为 a < 0
（开口向下）。若 a \\approx 0
，说明为平坦区，则取中心。
代码实现：
if any(m=='q')
b = 0.5 * (x(k+1) – x(k–1));
a = x(k) – b – x(k–1);
j = (a > 0); % 通常 a<0，此处 j 用于区分平坦区
v(j) = x(k(j)) + 0.25 * b(j).^2 ./ a(j);
k(j) = k(j) + 0.5 * b(j) ./ a(j);
k(~j) = k(~j) + (fq(k(~j)) – rs(k(~j))) / 2; % 平坦区取中心
end
注意：这里 a > 0 是异常情况（极小值），但实际峰值处 a 应为负。代码用 a 判断是否为平坦区（a 接近 0），如果 a>0 则强制按极小值修正，但通常不会发生。更稳健的实现应检查 a < -eps 才进行插值。
当两个峰值之间的距离小于等于 w 时，我们只保留幅值较大的那个。这是一个非极大值抑制（NMS）过程。
if nargin > 2
j = find(k(2:end) – k(1:end–1) <= w);
while any(j)
j = j + (v(j) >= v(j+1)); % 若前一个更高，则删除后一个；否则删除前一个
k(j) = [];
v(j) = [];
j = find(k(2:end) – k(1:end–1) <= w);
end
end
这段代码非常巧妙：j 是相邻峰值距离不足 w 的位置索引，然后根据幅值比较，决定删除前一个还是后一个（j 指向较低者）。循环直至所有距离均大于 w。
若 m 中包含 ‘v’，则函数先将信号取反（x = -x），执行完上述峰值检测后，再将幅值取反回来（v = -v）。这样谷值就变成了峰值，复用了同一套逻辑。
t = 0:0.01:1;
x = sin(2*pi*5*t) + 0.1*randn(size(t));
[k, v] = findpeakm(x, 'q', 0.1); % 二次插值，去重宽度 0.1 样本
% 绘制结果
findpeakm(x, 'q', 0.1); % 无输出时自动绘图
你会看到峰值位置被精确标记，且幅值接近 1。二次插值使得位置误差远小于采样间隔。
· 在语音分析中，峰值检测可用于基音周期提取（检测相邻波峰距离）。
· 在频谱分析中，检测谐波峰值可进行共振峰估计。
· 与 findSegment 结合，可在每个有效语音段内单独检测峰值，避免静音段的伪峰。
· 我们深入剖析了 findpeakm.m 的核心算法：差分定位、二次插值、邻近去重。
· 理解了谷值模式、平坦区处理、亚样本精度等高级特性。
· 通过示例展示了其简单易用的接口和强大的绘图功能。
峰值检测是信号分析中极为常见的基础操作，掌握这个工具将使你在处理各种波形特征提取时游刃有余。下一篇，我们将从峰值检测回到语音端点分割，学习如何将帧级 VAD 标记转换为连续的语音段结构体。
📥 所有代码均已打包，点击下方链接免费获取：
下载链接
下篇预告：语音端点检测后处理——连续有话段分割。我们将进入 findSegment.m，了解如何将离散的 0/1 标签聚合成有意义的语音片段，并计算每个片段的时长，为后续的语音识别或特征提取准备结构化数据。敬请期待！
思考题：如果信号中有一个很宽的平坦峰值（如方波的顶部），findpeakm 会定位在平台中心。若你想要检测平台两侧的边缘，应如何修改算法？欢迎评论区讨论。
