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

原创已于 2025-08-11 14:30:22 修改 · 883 阅读

6 ·

CC 4.0 BY-SA版权

文章标签：

#python #人工智能 #算法

于 2025-08-11 14:19:25 首次发布

python 同时被 2 个专栏收录

58 篇文章

订阅专栏

论文算法复现

1 篇文章

订阅专栏

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

在信号处理领域，从噪声信号中准确地识别出峰值是一项重要的挑战。最近，我研究了一篇题为《一种在噪声周期和准周期信号中自动检测峰值的高效算法》（原文名称

An Efficient Algorithm for Automatic Peak Detection in Noisy

Periodic and Quasi-Periodic Signals

）的论文，该论文提出了一种名为**自动多尺度峰值检测（AMPD）**的创新方法。

AMPD算法的独特之处在于它无需用户手动设置任何参数，就能在嘈杂的周期性和准周期性信号中高效地工作。这篇博客文章将带你深入了解AMPD算法的四个核心步骤，并展示如何用Python代码实现它。

算法步骤详解

AMPD算法的核心思想是利用信号在不同尺度下的局部极大值分布特性来识别真正的峰值。以下是其详细步骤：

线性去趋势: 在分析信号之前，算法首先对信号进行线性去趋势处理。这一步是为了消除信号的长期趋势，确保分析的重点集中在周期性振荡上。简单来说，就是拟合一条直线并从信号中减去它。
计算局部极大值刻度图（LMS）: 这是算法最关键的一步。它使用一个移动窗口，其长度从2到N/2逐步增加，N是信号的长度。对于每个窗口，算法会检查信号中的每个点是否是该窗口内的局部极大值。如果不是，则在对应的LMS矩阵位置上赋予一个随机数加上一个常数$\alpha$（通常设为1）；如果是，则赋予0。这个过程会生成一个LMS矩阵M，其中的每一行都代表了一个特定的尺度，而每一列则代表了信号中的一个时间点。矩阵中0的位置正是该尺度下局部极大值出现的地方。
对LMS进行逐行求和与重塑: 接下来，算法对LMS矩阵M进行逐行求和，得到一个向量$\gamma$。这个向量的最小值所对应的尺度，即为信号中周期性特征最显著的尺度$\lambda$。然后，算法会重塑LMS矩阵，只保留尺度1到$\lambda之间的行，生成一个新的矩阵M_r$。这一步有效地滤除了与信号主要周期性不相关的噪声和高频振荡。
峰值检测: 最后一步是计算重塑后的矩阵M_r的逐列标准差，得到一个标准差向量$\sigma$。由于矩阵M_r中的峰值点在所有尺度下都是局部极大值，因此其对应的列在重塑后仍然全部为0。这意味着在这些列中，标准差将为0。所以，算法通过寻找标准差向量中值为0的索引，就能精准地定位出信号中的所有峰值。

Python代码实现与可视化

以下是我使用Python实现的AMPD算法。代码清晰地展示了上述四个步骤，并利用matplotlib库对算法的中间过程进行了可视化，帮助我们更直观地理解它的工作原理。

import numpy as np
import matplotlib.pyplot as plt
from matplotlib import rcParams

# 设置matplotlib中文支持
rcParams['font.sans-serif'] = ['SimHei']  # 设置中文字体支持
rcParams['axes.unicode_minus'] = False  # 解决负号显示问题


def ampd_algorithm(x):
    """
    AMPD（自动多尺度峰值检测）算法的实现。
    本实现遵循论文《一种在噪声周期和准周期信号中自动检测峰值的高效算法》中概述的步骤。

    Args:
        x (np.ndarray): 一维均匀采样信号。

    Returns:
        np.ndarray: 对应于检测到的峰值的索引数组。
    """
    N = len(x)

    # 论文中L定义为floor(N/2)-1。我们使用整数除法 N // 2 - 1。
    L = N // 2 - 1

    # 1. 线性去趋势
    # 论文中首先指定对信号进行线性去趋势。
    t = np.arange(N)
    p = np.polyfit(t, x, 1)
    x_detrended = x - np.polyval(p, t)

    # 2. 计算局部极大值刻度图 (LMS)
    # 这是一个L x N的矩阵。
    M = np.zeros((L, N))
    alpha = 1.0

    # 遍历每个尺度k
    for k in range(1, L + 1):
        # 使用移动窗口方法，窗口长度为2k。
        # 论文中的索引是1-based，Python中的是0-based。
        # 循环i从k+1到N-k-1（0-based）以检查局部极大值。
        for i in range(k + 1, N - k):
            # 检查该点是否为窗口中的局部极大值。
            # 论文中更可靠的索引解释是将中心点与两侧k步的点进行比较。
            if x_detrended[i] > x_detrended[i - k] and x_detrended[i] > x_detrended[i + k]:
                M[k - 1, i] = 0
            else:
                # 如果不是局部极大值，则分配一个随机数+alpha。
                M[k - 1, i] = np.random.uniform(0, 1) + alpha

        # 处理论文中指定的边界条件：
        # “对于 i=1,...,k+1，以及 i=N-k+2,...,N，将值 r+alpha 分配给 m_{k,i}。”
        # 这意味着前k+1和后 N-(N-k+1)+1 = k 个元素。
        r_alpha = np.random.uniform(0, 1) + alpha
        M[k - 1, :k] = r_alpha
        M[k - 1, N - k:] = r_alpha

    # 3. 对LMS进行逐行求和
    # gamma向量包含有关局部极大值随尺度的分布信息。
    gamma = np.sum(M, axis=1)

    # gamma的全局最小值代表包含最多局部极大值的尺度。
    # lambda_val是尺度索引（1-based），因此我们在np.argmin结果上加1。
    lambda_val = np.argmin(gamma) + 1

    # 通过移除k > lambda的行来重塑LMS矩阵。
    Mr = M[:lambda_val, :]

    # 4. 峰值检测
    # 计算重塑后的LMS矩阵的逐列标准差。
    sigma = np.std(Mr, axis=0)

    # 峰值是标准差为零的索引。
    peaks = np.where(sigma == 0)[0]

    return peaks


# --- 基于论文（图2）的模拟数据生成 ---
fs = 800  # 采样频率 (sampling frequency)
N = 3000  # 采样点数 (number of samples)
t = np.arange(N) / fs
f1, f2, f3 = 10, 70, 5  # 频率 (frequencies)
a, b, c, d = 1, 1, 0.5, 0.1  # 振幅和噪声因子 (amplitudes and noise factor)

# 信号由三个正弦波和高斯白噪声组成
# The signal is composed of three sine waves and Gaussian white noise
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 算法
# Apply the AMPD algorithm
detected_peaks = ampd_algorithm(x)

# --- 绘制结果 (Plot the results) ---
# 为了更好地可视化，我们绘制去趋势后的信号。
# To visualize better, we plot the detrended signal.
t_detrended = np.arange(N)
p_detrended = np.polyfit(t_detrended, x, 1)
x_detrended = x - np.polyval(p_detrended, t_detrended)

plt.figure(figsize=(14, 8))
plt.plot(t, x_detrended, color='black', label='模拟信号')
plt.plot(t[detected_peaks], x_detrended[detected_peaks], 'ro', markersize=8, label='检测到的峰值')
plt.title('AMPD 算法应用于模拟信号', fontsize=16)
plt.xlabel('时间 [s]', fontsize=12)
plt.ylabel('振幅 [a.u.]', fontsize=12)
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.show()

print("检测到的峰值索引:", detected_peaks)
print(f"检测到的峰值总数: {len(detected_peaks)}")

还有一个版本代码如下：

import numpy as np
import matplotlib.pyplot as plt
import matplotlib
import scipy.linalg as la
from matplotlib import rcParams

# 设置matplotlib中文支持
rcParams['font.sans-serif'] = ['SimHei']  # 设置中文字体支持
rcParams['axes.unicode_minus'] = False  # 解决负号显示问题


def least_squares(data):
    """
    最小二乘线性去趋势
    - params:
        - data: 输入数据
    - return:
        - a: 截距 (y = a + b*x)
        - b: 斜率 (y = a + b*x)
    """
    x = np.array(range(len(data)))
    A = np.vstack([x ** 0, x ** 1])
    sol, r, rank, s = la.lstsq(A.T, np.array(data))
    return sol[0], sol[1]  # 修正：返回两个标量值


def automatic_peak_detection(data):
    # 最小二乘，线性去趋势
    a, b = least_squares(data)
    x = np.array(range(len(data)))
    trend = a + b * x
    detrended_data = data - trend  # 使用去趋势后的数据

    fig, ((ax0, ax1), (ax2, ax3)) = plt.subplots(2, 2, figsize=(12, 10))

    # 滑动窗口
    alpha = 1.
    # L行N列
    rows = len(data) // 2 - 1
    m = np.random.rand(rows, len(data)) * alpha + 1

    for k in range(1, len(data) // 2 + 1):
        left = k
        right = len(data) - k
        for i in range(left, right):
            # 检测局部最大值
            if detrended_data[i] > detrended_data[i - k] and detrended_data[i] > detrended_data[i + k]:
                m[k - 1, i] = 0

    ax0.imshow(m, cmap="gray", aspect='auto')
    ax0.set_title("局部最大值标度图(LMS)")
    ax0.set_xlabel("数据点索引")
    ax0.set_ylabel("尺度(k)")

    # 对m矩阵逐行求和
    gamma = np.sum(m, axis=1)
    ax1.plot(range(1, len(gamma) + 1), gamma, 'b-')
    ax1.set_title("LMS行和分布")
    ax1.set_xlabel("尺度(k)")
    ax1.set_ylabel("行和值")

    # 找出最小行和对应的尺度
    lambda_ = np.argmin(gamma)
    local_maxima = m[:lambda_ + 1, :]

    ax2.imshow(local_maxima, cmap="gray", aspect='auto')
    ax2.set_title(f"裁剪后的LMS(λ={lambda_ + 1})")
    ax2.set_xlabel("数据点索引")
    ax2.set_ylabel("尺度(1-λ)")

    # 计算列标准差
    sigma = np.std(local_maxima, axis=0)
    ax3.plot(sigma, 'g-')
    ax3.set_title("列标准差分布")
    ax3.set_xlabel("数据点索引")
    ax3.set_ylabel("标准差")
    ax3.axhline(y=0, color='r', linestyle='--', alpha=0.5)

    plt.tight_layout()

    # 峰值检测
    peaks = []
    for i in range(len(data)):
        if sigma[i] == 0:
            peaks.append(i)
    return peaks, detrended_data


# 构造含噪声的数据
np.random.seed(42)  # 固定随机种子以便复现结果
x = np.arange(1000)
base_signal = np.sin(0.2 * x) + 2 * np.sin(0.4 * x)
noise = 0.8 * np.random.randn(len(x))  # 生成一维噪声
y = base_signal + noise

# 执行峰值检测
peaks, detrended_y = automatic_peak_detection(y)

# 绘制最终结果
plt.figure(figsize=(14, 6))
plt.plot(x, y, 'b-', label='原始信号', alpha=0.7)
plt.plot(x, detrended_y, 'g-', label='去趋势信号', alpha=0.7)
plt.scatter(peaks, y[peaks], color='red', s=40, zorder=5, label='检测峰值')
plt.xlim(0, len(y))
plt.title("多尺度自动峰值检测结果")
plt.xlabel("数据点索引")
plt.ylabel("信号值")
plt.legend()
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

# 输出统计信息
print(f"检测到峰值数量: {len(peaks)}")
print(f"峰值位置索引: {peaks}")