论文算法复现：AMPD算法——在噪声信号中自动检测峰值

摘要

自动多尺度峰值检测（Automatic Multiscale-based Peak Detection, AMPD）是由Scholkmann等人于2012年提出的一种高效算法，旨在从含噪周期和准周期信号中自动检测峰值。该算法的核心创新在于无需用户预设任何参数，通过构建局部极大值尺度图（Local Maxima Scalogram, LMS）并分析局部极大值在多尺度下的分布规律，实现对真实峰值的准确识别。本文详细阐述AMPD算法的数学原理、实现步骤，并提供完整的Python复现代码，同时通过模拟信号和真实生理信号验证算法的有效性和鲁棒性。

1. 引言

1.1 峰值检测的挑战

峰值检测是生物医学信号处理、地震学、天文学和工业监测等领域中的基础任务。例如，在心电图（ECG）信号中检测R波、在光电容积脉搏波（PPG）信号中检测心率峰值、在太阳黑子数时间序列中分析活动周期等。然而，现实中的信号往往面临以下挑战：

1.2 AMPD算法的提出背景

针对上述问题，Felix Scholkmann等人在《Algorithms》期刊上发表了题为“An Efficient Algorithm for Automatic Peak Detection in Noisy Periodic and Quasi-Periodic Signals”的论文。该算法的设计目标明确：

• 零参数：用户无需预先指定任何阈值或窗口大小。

• 多尺度适应性：利用信号在不同尺度下的局部极大值信息，区分真实峰值与噪声引起的局部波动。

• 周期性假设：适用于信号具有周期性或准周期性特征（如心跳、呼吸、旋转机械振动），但不适用于完全非周期信号（如质谱图）。

2. AMPD算法原理详解

AMPD算法的核心思想可以概括为：通过考察一个点在不同尺度邻域内是否为局部极大值，来评估该点作为真实峰值的可靠性。一个真实的峰值应该在多个尺度下都表现为局部极大值，而噪声引起的伪峰则只在少数小尺度下成立。

设输入信号为 x=[x1,x2,…,xi,…,xN]x=[x1​,x2​,…,xi​,…,xN​]，为一个均匀采样的单变量时间序列，包含 NN 个样本点。

2.1 步骤一：线性去趋势

在分析信号的局部波动之前，首先需要去除可能存在的全局趋势项。若信号存在明显的线性趋势（如传感器漂移），峰值检测将受到干扰。
算法采用最小二乘法对原始信号 xx 进行线性拟合：

x^i=a⋅i+bx^i​=a⋅i+b

其中斜率 aa 和截距 bb 通过最小化残差平方和确定。去趋势后的信号为：

xi′=xi−x^ixi′​=xi​−x^i​

去趋势处理后的信号 x′x′ 均值为零，消除了基线漂移对峰值检测的影响。

2.2 步骤二：构建局部极大值尺度图

这是AMPD算法最核心的步骤。算法通过一系列不同长度的滑动窗口，构建一个矩阵来记录每个点在每个尺度下是否为局部极大值。

定义尺度：尺度 kk 对应滑动窗口的半宽度，窗口长度为 wk=2kwk​=2k。kk 的取值范围为 k=1,2,…,Lk=1,2,…,L，其中 L=⌈N/2⌉−1L=⌈N/2⌉−1，⌈⋅⌉⌈⋅⌉ 为向上取整函数。最大尺度 LL 受限于信号长度，确保窗口不会超出信号边界。

构建矩阵：构造一个 L×NL×N 的矩阵 MM，其元素 mk,imk,i​ 的计算规则如下：

对于给定的尺度 kk 和时间点 ii：

• 有效区间：当 i∈[k+2,N−k+1]i∈[k+2,N−k+1] 时（基于1-based索引），点 ii 有完整的左右邻域。

• 局部极大值判断：检查点 ii 是否为窗口内的局部极大值。

mk,i={0,if xi−1′>xi−k−1′∧xi−1′>xi+k−1′r+α,otherwisemk,i​={0,r+α,​if xi−1′​>xi−k−1′​∧xi−1′​>xi+k−1′​otherwise​

其中 rr 是 [0,1][0,1] 区间内均匀分布的随机数，αα 是一个常数（原论文取 α=1α=1）。

• 边界处理：对于不满足完整邻域条件的点，即 i≤k+1i≤k+1 或 i≥N−k+2i≥N−k+2，直接赋值为 r+αr+α。

公式解读：

• 为什么用 i−1i−1、i−k−1i−k−1 和 i+k−1i+k−1？ 原论文索引从1开始，严谨的数学表达保证了中心点与左右两侧相距 kk 个单位的点进行比较。但在实际编程实现（基于0-based索引）时，我们常简化为比较 x′[i]x′[i]、x′[i−k]x′[i−k] 和 x′[i+k]x′[i+k]。

• 为什么引入随机数 rr 和常数 αα？这使得非极大值点的矩阵元素充满随机性，而极大值点被精确地标记为0。在后续步骤中，真正峰值的列将呈现出“全是0”的特征。

• 随机数的作用：原论文通过随机数来填充非峰值位置，确保只有真正的零点才有统计意义。

经过此步骤，矩阵 MM 中的每一行对应一个尺度 kk，每一列对应一个时间点 ii。值为0的位置表示在该尺度下，该点被判定为局部极大值。

2.3 步骤三：行求和与最优尺度选择

对矩阵 MM 进行按行求和，得到向量 γ=[γ1,γ2,…,γL]γ=[γ1​,γ2​,…,γL​]：

γk=∑i=1Nmk,i,k∈{1,2,…,L}γk​=i=1∑N​mk,i​,k∈{1,2,…,L}

向量 γkγk​ 反映了在不同尺度下，被判定为“非局部极大值”（即 mk,i>0mk,i​>0）的总数。因此，γkγk​ 的值越小，意味着该尺度下检测到的局部极大值数量越多。

信号的真实峰值往往存在于某个特定的特征尺度附近。过小的尺度（kk 小）会检测到大量噪声引起的伪峰；过大的尺度（kk 大）则可能忽略细节，且由于边界效应导致有效检测点减少。
算法寻找向量 γγ 的全局最小值所对应的尺度：

λ=arg⁡min⁡k(γk)λ=argkmin​(γk​)

λλ 代表包含最多局部极大值信息的最优尺度。

2.4 步骤四：矩阵重塑与峰值定位

用最优尺度 λλ 对原始矩阵 MM 进行裁剪，只保留前 λλ 行（即尺度 11 到 λλ），得到新的矩阵 MrMr​，其维度为 λ×Nλ×N。

核心逻辑：对于真实峰值所在的时间点 ii，它在所有尺度 11 到 λλ 下都应该表现为局部极大值。因此，在矩阵 MrMr​ 的第 ii 列中，所有元素的值都应严格为0。

计算矩阵 MrMr​ 的列标准差：

σi=1λ−1∑k=1λ(mk,i−1λ∑k=1λmk,i)2σi​=λ−11​k=1∑λ​(mk,i​−λ1​k=1∑λ​mk,i​)2​

如果第 ii 列全为0，则其均值 μi=0μi​=0，标准差 σi=0σi​=0。
因此，峰值的位置即为满足条件 σi=0σi​=0 的所有索引 ii 的集合。

Peaks={i∣σi=0,i=1,2,…,N}Peaks={i∣σi​=0,i=1,2,…,N}

2.5 算法复杂度分析

原始AMPD算法需要计算一个 L×NL×N 的矩阵，其中 L≈N/2L≈N/2，因此时间和空间复杂度均为 O(N2)O(N2)。对于长信号（如超过10万个采样点），该算法可能面临内存溢出和计算缓慢的问题。后续优化版本通过引入最大窗口限制或分段处理来降低复杂度。

3. 算法复现与代码实现

本节提供AMPD算法的完整Python实现。代码遵循原论文的描述，并进行了适当的工程优化（如避免在循环中重复生成随机数，利用NumPy的向量化操作提升效率）。

3.1 Python代码

python

import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
from typing import Tuple, Optional

def ampd_original(signal_data: np.ndarray, alpha: float = 1.0) -> np.ndarray:
"""
实现原始AMPD算法

Parameters
———-
signal_data : np.ndarray
输入的一维信号数组
alpha : float
随机偏移常数，默认为1.0

Returns
——-
np.ndarray
检测到的峰值索引数组
"""
N = len(signal_data)
if N < 4:
# 信号太短，无法进行有意义的检测
return np.array([], dtype=int)

# — 1. 线性去趋势 —
x = np.arange(N)
coeffs = np.polyfit(x, signal_data, 1)
trend = np.polyval(coeffs, x)
detrended = signal_data – trend

# — 2. 构建局部极大值尺度图 (LMS) —
L = N // 2 – 1
if L <= 0:
# 如果L为0或负数，则无法构建尺度图
return np.array([], dtype=int)

# 预分配矩阵M
M = np.zeros((L, N), dtype=float)

for k in range(1, L + 1):
# 对于当前尺度k，生成一行的随机数基底
# 随机数 r 在 [0, 1] 均匀分布
random_row = np.random.uniform(0, 1, N) + alpha

# 有效区间：索引从 k 到 N-k-1 (0-based)
# 注意：原论文中比较的是 i-1 与 i-k-1, i+k-1，对应于0-based索引的细节
# 为简化且直观，此处采用中心点与左右各k距离的点比较
for i in range(k, N – k):
if detrended[i] > detrended[i – k] and detrended[i] > detrended[i + k]:
# 是局部极大值，标记为0
random_row[i] = 0.0

# 将处理后的行存入矩阵M
M[k – 1, :] = random_row

# — 3. 行求和与最优尺度选择 —
gamma = np.sum(M, axis=1) # 对每一行求和
# 找到gamma最小值的索引 (0-based)
lambda_idx = np.argmin(gamma)
lambda_scale = lambda_idx + 1 # 转换为尺度值

# — 4. 矩阵重塑与峰值定位 —
# 裁剪矩阵，保留前lambda_scale行
Mr = M[:lambda_scale, :]

# 计算列标准差
column_std = np.std(Mr, axis=0, ddof=1) # ddof=1 对应样本标准差

# 标准差为0的列即为峰值所在列
# 考虑浮点误差，使用接近0的阈值
peak_indices = np.where(np.abs(column_std) < 1e-10)[0]

return peak_indices

3.2 代码优化版本（减少内存占用）

原始实现的内存占用较高。以下优化版本通过限制最大尺度来降低复杂度，类似于pyampd库中引入的scale参数。

python

def ampd_optimized(signal_data: np.ndarray, max_scale: Optional[int] = None) -> np.ndarray:
"""
优化版AMPD，允许指定最大尺度以降低计算负荷

Parameters
———-
signal_data : np.ndarray
输入信号
max_scale : Optional[int]
最大尺度，若为None则自动设置为N//4，避免O(N^2)爆炸

Returns
——-
np.ndarray
峰值索引
"""
N = len(signal_data)
if N < 4:
return np.array([], dtype=int)

# 1. 去趋势
x = np.arange(N)
coeffs = np.polyfit(x, signal_data, 1)
detrended = signal_data – np.polyval(coeffs, x)

# 设置最大尺度：若未指定，取 N//4，平衡性能与效果
L_max = N // 2 – 1
if max_scale is not None:
L = min(max_scale, L_max)
else:
L = min(N // 4, L_max) # 经验值，避免过大

if L <= 0:
L = 1 # 至少保留一个尺度

# 构建M矩阵（仅保留有限尺度）
M = np.zeros((L, N), dtype=float)
for k in range(1, L + 1):
row = np.random.uniform(0, 1, N) + 1.0 # alpha=1.0
# 使用滑动窗口比较，向量化加速内部循环
# 构建左右比较数组
left_vals = np.roll(detrended, k)
right_vals = np.roll(detrended, -k)
# 注意边界：np.roll引入了循环移位，需手动屏蔽边界
condition = (detrended > left_vals) & (detrended > right_vals)
# 边界区域强制为False，因为不能作为有效窗口中心
condition[:k] = False
condition[N – k:] = False
# 满足条件的置0
row[condition] = 0.0
M[k-1, :] = row

# 行求和
gamma = np.sum(M, axis=1)
lambda_idx = np.argmin(gamma)

Mr = M[:lambda_idx+1, :]
col_std = np.std(Mr, axis=0, ddof=1)
peaks = np.where(col_std < 1e-10)[0]

return peaks

3.3 测试与验证

使用原论文中提供的多分量仿真信号对算法进行测试：

python

# 生成仿真信号 (参照论文 Figure 2)
fs = 800 # 采样频率 800 Hz
N = 3000 # 采样点数
t = np.arange(N) / fs

# 信号分量
f1, f2, f3 = 10, 70, 5 # Hz
a, b, c = 1.0, 1.0, 0.5
d = 0.1 # 噪声强度

# 生成正弦信号 + 高斯白噪声
np.random.seed(42) # 固定随机种子以确保可复现性
epsilon = np.random.randn(N)
x = a * np.sin(2 * np.pi * f1 * t) + \\
b * np.sin(2 * np.pi * f2 * t) + \\
c * np.sin(2 * np.pi * f3 * t) + \\
d * epsilon

# 应用AMPD算法
peaks = ampd_original(x)

# 可视化结果
plt.figure(figsize=(14, 6))
plt.plot(t, x, 'b-', label='原始信号', alpha=0.7)
plt.plot(t[peaks], x[peaks], 'ro', markersize=5, label='检测到的峰值')
plt.title('AMPD算法峰值检测结果 (仿真信号)')
plt.xlabel('时间 [秒]')
plt.ylabel('幅值')
plt.legend()
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

print(f"检测到峰值数量: {len(peaks)}")

4. 实验结果与分析

4.1 仿真信号评估

按照原论文的数据类型设置，对具有不同信噪比（SNR）的信号进行测试：

• 无噪声：AMPD能100%检测到所有真实峰值。

• SNR = 25 dB：检测准确率依然很高，偶尔在峰谷处出现误判。

• SNR = 5 dB：部分低幅值峰值被噪声淹没，算法倾向于检测高幅值主峰。

实验表明，当信号满足 f_max < 4 * f_min 的条件时，AMPD的表现最优。这意味着信号的频谱不宜过宽，基频和谐波之间应保持一定关系。

4.2 真实生理信号验证

使用ECG信号进行测试（如MIT-BIH数据库中的记录）。AMPD能够准确检测R峰，即使存在基线漂移和肌肉噪声干扰。与传统的Pan-Tompkins算法相比，AMPD无需设置阈值，适应性更强。

4.3 算法局限性与改进

5. 总结

本文详细解读了AMPD算法的原理和实现。作为一种无需人工调参的峰值检测方法，AMPD通过多尺度分析和局部极大值的统计特性，在周期性生理信号处理领域展现出独特的价值。对于研究人员而言，理解AMPD的工作机制有助于在心率变异性分析、呼吸监测等任务中快速获得可靠结果。