📚 课程目录第4章基于二次多项式插值的峰值定位

4. 基于二次多项式插值的峰值定位

任务说明

布里渊光时域反射（BOTDR）系统中，通过洛伦兹或高斯拟合可估计布里渊频移，但逐点拟合计算量大。本案例实现基于二次多项式插值的峰值定位方法：先对离散布里渊谱进行粗搜索找到最大幅值点，再利用该点及其左右邻点构造二次多项式，通过解析求导获得亚像素级峰值位置，从而在少量采样点下提高频移分辨率。

核心算法

设离散布里渊谱为 \((x_i, y_i)\)，其中 \(x_i\) 为频率，\(y_i\) 为幅值。粗搜索得到最大值点索引 \(k\)，取三点 \((x_{k-1}, y_{k-1})\)、\((x_k, y_k)\)、\((x_{k+1}, y_{k+1})\)。构造二次多项式 \(p(x) = a x^2 + b x + c\)，其峰值位置为：

\(x_{peak} = x_k + \dfrac{\Delta x}{2} \cdot \dfrac{y_{k-1} - y_{k+1}}{y_{k-1} - 2y_k + y_{k+1}}\)

其中 \(\Delta x = x_{k+1} - x_k\) 为频率采样间隔。该式由三点插值二次函数求导为零直接导出，无需矩阵求逆。

代码实现


import numpy as np

def quadratic_peak_interpolation(freq, amp):
    """
    基于二次多项式插值的峰值定位
    参数:
        freq: 1D array, 频率轴 (GHz)
        amp:  1D array, 对应幅值 (任意线性单位)
    返回:
        peak_freq: float, 亚像素峰值频率 (GHz)
        peak_idx:  int,   粗搜索最大值索引
    """
    if len(freq) != len(amp) or len(freq) < 3:
        raise ValueError("freq 和 amp 长度需相等且至少为3")
    
    # 1. 粗搜索：最大幅值点
    k = int(np.argmax(amp))
    
    # 边界处理：若最大值在端点，无法插值，直接返回该点
    if k == 0 or k == len(freq) - 1:
        return freq[k], k
    
    # 2. 取三点
    x0, x1, x2 = freq[k-1], freq[k], freq[k+1]
    y0, y1, y2 = amp[k-1], amp[k], amp[k+1]
    
    # 3. 二次插值峰值偏移量（以 x1 为参考）
    dx = x2 - x1  # 等间隔采样时 = x1 - x0
    denom = y0 - 2*y1 + y2
    
    if abs(denom) < 1e-12:
        # 退化情况：三点共线或数值噪声，返回粗峰值
        return x1, k
    
    delta = 0.5 * dx * (y0 - y2) / denom
    
    # 4. 亚像素峰值位置
    peak_freq = x1 + delta
    return peak_freq, k

# ========== 示例：构造一个离散布里渊谱并定位 ==========
if __name__ == "__main__":
    # 真实布里渊频移设为 10.850 GHz
    true_peak = 10.850
    # 频率采样：间隔 5 MHz，覆盖 10.82 ~ 10.88 GHz
    freq = np.arange(10.820, 10.881, 0.005)
    
    # 洛伦兹型布里渊谱（线宽 30 MHz），加入少量噪声
    gamma = 0.030  # GHz
    amp = 1.0 / (1.0 + ((freq - true_peak) / gamma)**2)
    rng = np.random.default_rng(42)
    amp += rng.normal(0, 0.005, size=len(freq))
    
    # 粗搜索峰值（直接取最大点）
    coarse_idx = int(np.argmax(amp))
    coarse_freq = freq[coarse_idx]
    
    # 二次插值精定位
    peak_freq, peak_idx = quadratic_peak_interpolation(freq, amp)
    
    print(f"真实峰值:      {true_peak:.6f} GHz")
    print(f"粗搜索峰值:    {coarse_freq:.6f} GHz (索引 {coarse_idx})")
    print(f"插值精定位:    {peak_freq:.6f} GHz (索引 {peak_idx})")
    print(f"粗搜索误差:    {abs(coarse_freq - true_peak)*1e6:.1f} MHz")
    print(f"插值后误差:    {abs(peak_freq - true_peak)*1e6:.1f} MHz")

落地与验证

运行方式：将代码保存为 peak_interp.py，在装有 NumPy 的 Python 环境中执行 python peak_interp.py。

验证方法：观察输出中“插值后误差”是否显著小于“粗搜索误差”。由于频率采样间隔为 5 MHz，粗搜索最大误差理论上可达 ±2.5 MHz；二次插值后误差通常可降至亚 MHz 量级。若代码正确，插值结果应更接近真实值 10.850 GHz。

判定标准：程序无异常退出，且插值后误差小于粗搜索误差即视为通过。注意：由于示例中加入了随机噪声，每次运行的具体误差数值会略有波动，但插值精度提升的趋势应稳定成立。