如何不预拟合衰减正弦即可精准绘制阻尼函数衰减包络

需求说明 #

本次计算目标为阻尼函数衰减率：通过对阻尼振荡的幅值拟合指数曲线，总结其衰减规律，效果参考标准阻尼正弦波的衰减包络。待处理信号存在无规则快速局部振荡、振荡特性非均匀的特点，初始因账号权限限制未上传波形，后续已补充实际信号图。

现有方案的精度问题 #

目前常用的两类实现思路均存在明显短板：

• 阻尼正弦函数整体拟合法：标准阻尼正弦模型仅能匹配均匀振荡场景，无法适配信号的非均匀振荡特性，对带高频局部扰动的信号拟合效果极差，结果可信度低
• 分区间取最大值拟合法：计算精度完全依赖人工设定的区间划分长度，目前无通用的区间长度选取标准：区间过大易丢失幅值变化细节，区间过小会将局部高频振荡的峰值误判为全局包络幅值，结果稳定性差

精度更优的实现路径（Python全流程可实现） #

核心思路是绕开预设振荡模型、跳过人工区间划分环节，先自适应提取信号真实上包络，再基于包络计算衰减率，具体步骤如下：

1. 高频噪声预处理：采用零相位低通滤波滤除信号中无规则的快速局部振荡，优先选择filtfilt类零相位滤波算法，避免普通滤波带来的相位偏移导致峰值位置错位
2. 自适应包络提取：
   • 调用峰值检测算法自动识别滤波后信号的所有局部极大值点，通过峰值最小间距参数自动过滤不符合全局衰减趋势的伪峰值，无需手动划分区间
   • 对提取到的有效峰值点做三次样条插值，得到连续平滑的上包络曲线
3. 衰减率拟合：对得到的上包络曲线取自然对数，将指数拟合问题转化为线性最小二乘拟合问题，通过拟合得到的线性斜率即可直接换算得到指数衰减率，线性拟合的抗噪性远高于直接做非线性指数拟合。

核心实现代码 #

import numpy as np
from scipy.signal import find_peaks, filtfilt, butter
from scipy.interpolate import CubicSpline

def lowpass_filter(raw_signal: np.ndarray, sample_rate: float, cutoff_freq: float) -> np.ndarray:
    """零相位低通滤波，滤除高频局部振荡"""
    b, a = butter(N=3, Wn=cutoff_freq, fs=sample_rate, btype='low')
    return filtfilt(b, a, raw_signal)

def extract_upper_envelope(time_seq: np.ndarray, filtered_signal: np.ndarray, min_peak_dist: int = 5) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
    """自适应提取信号上包络，min_peak_dist根据采样率调整，过滤相邻过近的伪峰值"""
    peak_idx, _ = find_peaks(filtered_signal, distance=min_peak_dist)
    peak_time = time_seq[peak_idx]
    peak_amp = filtered_signal[peak_idx]
    spline_interp = CubicSpline(peak_time, peak_amp)
    upper_envelope = spline_interp(time_seq)
    return upper_envelope, peak_time, peak_amp

def calc_decay_rate(time_seq: np.ndarray, envelope: np.ndarray) -> float:
    """基于包络计算指数衰减率，模型为 A*exp(-alpha*t)"""
    log_envelope = np.log(envelope)
    # 线性拟合 ln(A) - alpha*t = ln(envelope)，一次项系数为 -alpha
    fit_coeff = np.polyfit(time_seq, log_envelope, deg=1)
    return -fit_coeff[0]

# 调用示例
# filtered_sig = lowpass_filter(raw_sig, sample_rate=fs, cutoff_freq=high_freq_cutoff)
# upper_env, _, _ = extract_upper_envelope(t_seq, filtered_sig, min_peak_dist=peak_spacing)
# decay_alpha = calc_decay_rate(t_seq, upper_env)

方案优势 #

• 无需预设信号振荡模型，对非均匀振荡的阻尼信号适配性强，不受局部高频扰动干扰
• 无需人工划分幅值提取区间，峰值检测为自适应过程，从根源上避免了区间长度选择带来的人为误差
• 全流程依赖Python生态中numpy、scipy的成熟内置函数实现，数值稳定性有保障，计算效率高

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