自定义峰值查找算法漏检振荡峰值致阻尼率计算异常排查

问题根因 #

• 自定义峰值检测函数存在逻辑缺陷
  现有手写find_peaks函数处理序列首个点（i==0）的分支执行完成后未设置跳转，会继续执行后续局部峰值判定逻辑，此时访问y[i-1]等价于访问序列末尾元素y[-1]，直接干扰靠前位置的峰值判定，造成前序峰值漏检。同时手写逻辑未做噪声容错、边界补全，对非单调振荡序列适配性差。
• 拟合模型与信号特征不匹配
  现有指数模型形式为A * np.exp(-lbda * t)，仅能描述从t=0开始单调衰减的过程，无法匹配振幅先上升、后阻尼衰减至零的信号特征，即便全量峰值识别正确，直接拟合全部峰值也会产生结果偏差。

修复方案 #

• 替换手写峰值检测逻辑为scipy.signal.find_peaks，该实现自带完善的边界判定、噪声过滤能力，可通过高度、突出度等参数筛选有效峰值，避免前序峰值漏检。
• 识别峰值序列的最大值点，仅截取最大值点之后的衰减段峰值开展指数拟合，同时适配两类场景：常规小振幅振荡场景下最大值点为首个峰值，截取逻辑不会丢失有效数据；含初始增幅段的场景下可避开上升段数据，保证拟合数据符合指数衰减假设。

参考实现代码 #

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
from scipy.signal import find_peaks

def exp_decay(t, A, lbda):
    return A * np.exp(-lbda * t)

Tarr=[Tplot0,Tplot1,Tplot2,Tplot3,Tplot4]
Zarr=[Zrel0,Zrel1,Zrel2,Zrel3,Zrel4]
Alpharr=[alpha0,alpha1,alpha2,alpha3,alpha4]

for x0,y0,alf in zip(Tarr,Zarr,Alpharr):
    # 检测有效峰值，可按需传入height、prominence参数过滤噪声伪峰
    peak_idx, _ = find_peaks(y0)
    peak_x0 = x0[peak_idx]
    peak_y0 = y0[peak_idx]
    
    # 截取最大峰值后的衰减段用于拟合
    max_peak_idx = np.argmax(peak_y0)
    fit_x = peak_x0[max_peak_idx:]
    fit_y = peak_y0[max_peak_idx:]
    
    popt, pcov = curve_fit(exp_decay, fit_x, fit_y)
    print(*[f"{val:.2f} +/- {err:.2f} __ " for val, err in zip(popt, np.sqrt(np.diag(pcov)))])
    print("decay rate = ", popt[1])
    
    plt.plot(x0, y0)
    plt.plot(fit_x, exp_decay(fit_x, *popt), '-r', label=f"{popt[0]:.2f}*exp(-{popt[1]:.2f}t)")
    plt.xlabel("time")
    plt.ylabel("Zrel")
    plt.title(f"n = {n} ; alpha = {alf}")
    plt.legend()
    plt.show()

适配说明 #

• 若信号存在明显噪声，可给find_peaks传入prominence参数设置峰值最小突出度，过滤噪声导致的伪峰值，提升检测准确率。
• 该逻辑对两类测试场景均兼容：常规小振幅振荡场景下拟合逻辑与原有效果一致，含初始增幅段的场景下可准确识别全部峰值、避开非衰减段数据，输出正确阻尼率。

内容的提问来源于stack exchange，提问作者IAmOneWithTheScientist