第十二篇：智能峰值/谷值检测算法详解

在前几篇中，我们专注于语音信号的预处理和预测延拓，但实际应用中我们常需要从波形中自动定位极值点——例如基音周期中寻找波峰、频谱中识别谐波峰值、或心跳信号中检测 R 波峰值。MATLAB 自带的 findpeaks 函数（需 Signal Processing Toolbox）功能强大，但有时我们需要更轻量、更可控的实现。findpeakm.m 提供了完全自包含的峰值/谷值检测算法，支持二次插值亚样本精度、宽度容差去重，并且无需任何工具箱依赖。本篇将深入剖析其核心算法，从一阶差分定位到插值优化，再到局部峰值去重，带你全面掌握这个实用工具。

---

1. 峰值检测的核心思想

峰值检测的本质是在离散序列中找到局部极大值点。一个点 $$x(k)$$ 是峰值，当且仅当其左侧相邻点小于或等于它，且右侧相邻点也小于或等于它（对于严格的峰值，两侧严格小于）。但实际信号常伴有噪声，直接比较容易出现虚假极值。

基本流程：

1. 计算相邻点的差分（一阶导数近似）；
2. 找出差分由正变负的位置 —— 这些位置就是峰值候选点；
3. 若信号包含谷值，则可将信号取反后再执行上述流程；
4. 可选：对候选点进行二次插值，获得亚样本级精度的位置和幅值；
5. 可选：根据距离容差剔除相距太近的峰值（保留幅值较大者）。

findpeakm.m 完整实现了上述逻辑，并支持 'q' 插值模式和 'v' 谷值模式。

---

2. 代码结构总览

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. 差分法定位峰值——核心算法详解

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 的行为）。

---

4. 二次插值 —— 实现亚样本精度

当信号峰值不是恰好落在采样点上时，我们可以用抛物线拟合三个相邻点（左、候选、右）来精确估计峰值的真实位置和幅值。

设候选点索引为 $$k$$，其幅值为 $$x_k$$，左右点为 $$x_{k-1}$$ 和 $$x_{k+1}$$。构造二次多项式：

$$f(t)=a{t}^2+bt+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 才进行插值。

---

5. 邻近峰值去重（宽度容差 w）

当两个峰值之间的距离小于等于 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。

---

6. 谷值检测模式

若 m 中包含 'v'，则函数先将信号取反（x = -x），执行完上述峰值检测后，再将幅值取反回来（v = -v）。这样谷值就变成了峰值，复用了同一套逻辑。

---

7. 实战示例：检测正弦波中的峰值

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。二次插值使得位置误差远小于采样间隔。

---

8. 与其他函数的关系

• 在语音分析中，峰值检测可用于基音周期提取（检测相邻波峰距离）。
• 在频谱分析中，检测谐波峰值可进行共振峰估计。
• 与 findSegment 结合，可在每个有效语音段内单独检测峰值，避免静音段的伪峰。

---

9. 本讲小结

• 我们深入剖析了 findpeakm.m 的核心算法：差分定位、二次插值、邻近去重。
• 理解了谷值模式、平坦区处理、亚样本精度等高级特性。
• 通过示例展示了其简单易用的接口和强大的绘图功能。

峰值检测是信号分析中极为常见的基础操作，掌握这个工具将使你在处理各种波形特征提取时游刃有余。下一篇，我们将从峰值检测回到语音端点分割，学习如何将帧级 VAD 标记转换为连续的语音段结构体。

---

📥 所有代码均已打包，点击下方链接免费获取：
下载链接

---

下篇预告：语音端点检测后处理——连续有话段分割。我们将进入 findSegment.m，了解如何将离散的 0/1 标签聚合成有意义的语音片段，并计算每个片段的时长，为后续的语音识别或特征提取准备结构化数据。敬请期待！