常见量化算法系列 - Savitzky-Golay 滤波

: 阿凡 2025，大牛市快来了 V：flj99103681

0 People liked this article · 409 views

常见量化算法系列

Savitzky-Golay 滤波：保峰的"精修磨皮"

文 / 阿凡

如果说高斯平滑是"美图秀秀一键磨皮"，那 Savitzky-Golay（简称 SG）滤波就是"专业修图师精修"。

高斯平滑：波峰被削低，波谷被填高（像磨皮过度，五官扁平了）
SG 滤波：保留波峰波谷的高度和宽度（只磨瑕疵，五官依然立体）

为什么？因为它用的不是"加权平均"，而是局部多项式拟合。

💡 提示：本文所有代码均可直接复制运行，完整源码在文章末尾。

第一部分

从零手搓：SG 滤波完整实现

不调用任何现成库，我们用 NumPy 从零开始，一步一步实现 SG 滤波。整个过程分为 5 个步骤。

♦

Step 1：构建范德蒙矩阵

SG 滤波的核心是"在局部窗口内拟合一条多项式曲线"。要拟合曲线，首先要构造范德蒙矩阵（Vandermonde matrix）。

假设窗口大小为 5，窗口中心的索引为 0，则窗口内各点的相对位置为：

k = [-2, -1, 0, 1, 2]

用 2 阶多项式（y = a₀ + a₁·k + a₂·k²）去拟合，范德蒙矩阵的每一列就是 k 的 0 次、1 次、2 次幂：

PYTHON · Step 1：范德蒙矩阵

import numpy as np

# 窗口大小 5，中心点索引为 0
window = 5
half = window // 2          # half = 2
k = np.arange(-half, half + 1)   # [-2, -1, 0, 1, 2]

# 用 2 阶多项式拟合
poly_order = 2

# 构建范德蒙矩阵：每一列 = k 的 j 次幂
V = np.vander(k, poly_order + 1, increasing=True)
# V 的形状: (5, 3)
# V = [[ 1, -2, 4],     ← k=-2 时的 [k⁰, k¹, k²]
#       [ 1, -1, 1],     ← k=-1 时的 [k⁰, k¹, k²]
#       [ 1, 0, 0],     ← k= 0 时的 [k⁰, k¹, k²]
#       [ 1, 1, 1],     ← k= 1 时的 [k⁰, k¹, k²]
#       [ 1, 2, 4]]     ← k= 2 时的 [k⁰, k¹, k²]

范德蒙矩阵的本质：每一行代表一个采样点的"多项式基"，每一列代表一个幂次。它把"用多项式描述一组点"这件事，转化成了矩阵运算。

Step 2：最小二乘求卷积核

范德蒙矩阵 V 建好后，我们要解的问题是：找到一组多项式系数 a = [a₀, a₁, a₂]，使得多项式曲线与窗口内实际数据点的误差平方和最小。

这就是最小二乘法的标准形式：

最小二乘正规方程

(V^T · V) · a = V^T · y

但我们不需要真的去解某一个窗口的 a！SG 滤波最巧妙的地方在于：

核心洞察

拟合曲线在窗口中心（k=0）的值，恰好等于原始数据的线性组合。这个线性组合的系数只跟"窗口多大"和"用几阶多项式"有关，跟具体数据无关！

所以我们可以提前算好这个系数数组（卷积核），之后每次只需要做加权求和即可。

PYTHON · Step 2：求卷积核

# 求解 (V^T @ V)^(-1) @ V^T
# 这个矩阵的每一行，就是对应导数阶数的卷积核
A = np.linalg.pinv(V.T @ V) @ V.T

# 取第 0 行（对应平滑操作，deriv=0）作为卷积核
coeffs = A[0]
# coeffs ≈ [-0.0857, 0.3429, 0.4857, 0.3429, -0.0857]
# 乘以 35 后就是整数：[-3, 12, 17, 12, -3]
# 归一化（让系数和为 1）
coeffs = coeffs / np.sum(coeffs)

为什么只取第 0 行？

因为我们要求的是平滑值（deriv=0），即拟合曲线在窗口中心的函数值。

如果想求一阶导数（deriv=1），就取 A[1]；求二阶导数（deriv=2），就取 A[2]。同一个矩阵，搞定所有事情。

Step 3：边界填充

卷积核宽度为 5，在信号两端会"悬空"。我们需要在两端补一些数据，让窗口能完整滑动。

PYTHON · Step 3：边界填充

# 窗口半宽
pad = window // 2 # pad = 2

# 用最近邻值填充两端（最稳妥的方式）
x_padded = np.pad(x, pad, mode='nearest')

# 举例：x = [10, 20, 30, 40, 50]
# x_padded = [10, 10, 10, 20, 30, 40, 50, 50, 50]
# ↑补2个原始数据 ↑补2个

Step 4：卷积操作

有了卷积核和填充后的数据，剩下的就是一个滑动加权求和——这正是卷积的定义。

PYTHON · Step 4：卷积

# 对填充后的信号做卷积，valid 模式保证输出长度正确
y = np.correlate(x_padded, coeffs, mode='valid')

# 验证输出长度：
# len(valid) = len(padded) - len(coeffs) + 1
# = (n + 2*pad) - (2*pad + 1) + 1 = n ✅

手动画一个卷积过程

以窗口 = 5，核 = [-3, 12, 17, 12, -3] / 35 为例：

窗口数据: [100, 105, 110, 108, 112]
卷积核: [-3, 12, 17, 12, -3] 35

y = (-3×100 + 12×105 + 17×110 + 12×108 - 3×112) ÷ 35
= (-300 + 1260 + 1870 + 1296 - 336) ÷ 35
= 3790 ÷ 35 ≈ 108.29

注意最外侧的 -3 是负权重——它抵消了中间权重过大导致的"削峰"效应，这就是保峰的数学本质。

Step 5：打包成函数

把上面 4 步合并成一个完整的函数：

PYTHON · 手写 SG 滤波器

def manual_savgol(x, window=5, poly_order=2):
    """手写 Savitzky-Golay 平滑滤波器"""
    x = np.asarray(x, dtype=float)
    half = window // 2

    # Step 1: 构建范德蒙矩阵
    k = np.arange(-half, half + 1, dtype=float)
    V = np.vander(k, poly_order + 1, increasing=True)

    # Step 2: 求卷积核（最小二乘）
    A = np.linalg.pinv(V.T @ V) @ V.T
    coeffs = A[0] # deriv=0 取平滑核
    coeffs = coeffs / np.sum(coeffs) # 归一化

    # Step 3: 边界填充
    x_padded = np.pad(x, half, mode='nearest')

    # Step 4: 卷积
    y = np.correlate(x_padded, coeffs, mode='valid')
    return y[:len(x)]

♦

2 S-G 滤波滑动窗口局部拟合原理

运算流程回顾

1. 滑动窗口截取一段局部数据（图中 5 个黑点）

2. 用范德蒙矩阵 + 最小二乘拟合低阶多项式（绿色虚线）

3. 取窗口中心拟合值作为平滑结果（红色圆点）

4. 窗口向右滑动 1 个点，重复计算

核心：不做简单平均，而是局部拟合曲线——多项式贴合原始波形趋势，不会削峰。

♦

第二部分

进阶：用 SG 滤波求导数

还记得我们之前说的吗？矩阵 A 的每一行对应不同阶数的导数。刚才只用了 A[0]（平滑），现在把 A[1] 和 A[2] 也利用起来。

PYTHON · 手写 SG 求导

def manual_savgol_deriv(x, window=5, poly_order=2, deriv=1):
    """手写 SG 滤波求导"""
    x = np.asarray(x, dtype=float)
    half = window // 2

    # Step 1: 范德蒙矩阵（同上）
    k = np.arange(-half, half + 1, dtype=float)
    V = np.vander(k, poly_order + 1, increasing=True)

    # Step 2: 取 deriv 对应的行
    A = np.linalg.pinv(V.T @ V) @ V.T
    coeffs = A[deriv]

    # Step 3: 求导核需要乘以阶乘
    import math
    coeffs = coeffs * math.factorial(deriv)

    # Step 4 & 5: 填充 + 卷积（同上）
    x_padded = np.pad(x, half, mode='nearest')
    y = np.correlate(x_padded, coeffs, mode='valid')
    return y[:len(x)]

# 使用示例
velocity = manual_savgol_deriv(prices, deriv=1)    # 一阶导（速度）
accel     = manual_savgol_deriv(prices, deriv=2)    # 二阶导（加速度）

类比开车：
价格 = 高速上的位置
一阶导数 = 车速表上的速度
二阶导数 = 踩油门还是刹车的感觉

为什么求导有用？

一阶导数过零点 → 价格达到局部最高/最低点

二阶导数过零点 → 趋势拐点（由涨转跌、由跌转涨）

杀手级功能

同一个范德蒙矩阵，取不同行就能得到平滑、一阶导、二阶导——一个算法，多种用途。

♦

3 S-G 卷积核线性运算拆解

三层拆解回顾：

第一层：5 点 2 阶卷积核

归一化系数：[-3, 12, 17, 12, -3] ÷ 35

中心权重最高（17/35），最外侧两个点是负权重（-3/35）——这是保峰的关键。

第二层：窗口内 5 个采样值

y₋₂ y₋₁ yᵢ yᵢ₊₁ yᵢ₊₂

第三层：加权求和

最终公式

y_i^smooth = (-3·y_i-2 + 12·y_i-1 + 17·y_i + 12·y_i+1 - 3·y_i+2) ÷ 35

系数提前固定：只由窗口长度和多项式阶数决定，和数据无关。
⚡ 高效：预计算核 + 滑动卷积，适合实时处理。
🛡️ 负权重保峰：最外侧的 -3 抵消峰衰减，维持极值高度。

♦

4 SG 滤波完整流程

输入带噪数据 → 设定参数 → 计算卷积核 → 滑动窗口 → 边界处理 → 输出平滑序列

1	输入一维带噪时序数据 y[1:N]
2	设定：窗口长度 L=2m+1、多项式阶数 n（常用 2/3 阶）
3	构建范德蒙矩阵 → 最小二乘求卷积核
4	滑动窗口遍历：① 截取局部数据 ② 卷积核加权求和
5	边界处理（最近邻填充）
6	输出完整平滑序列

♦

5 SG 滤波三大核心特性

1. 降噪

加权平均抑制随机噪声，提升信噪比。

2. 保形

局部多项式拟合，保留峰值高度、拐点位置——SG 滤波的最大优势。

3. 高效

预计算核系数，非线性拟合 → 线性卷积，计算开销极低。

♦

第三部分

量化策略中的用途

场景	用法
趋势跟踪	SG 平滑后的价格作为趋势基准，比均线更灵敏
拐点检测	二阶导数过零点 = 潜在买卖信号
指标预处理	先 SG 平滑再算 RSI/MACD，减少假信号
多周期分析	不同 window 对应不同时间尺度

SG 滤波 vs 加权平均

特性	加权平均	SG 滤波
原理	加权平均	局部多项式拟合
波峰保留	会削低	✅ 基本保留
波谷保留	会填高	✅ 基本保留
能否求导	需额外处理	✅ 内置求导
计算速度	快	稍慢
适用场景	简单去噪	精细分析、拐点检测

♦

第四部分

注意事项

window 必须为奇数：偶数没有明确中心点

poly_order < window：否则方程欠定（未知数比方程多）

window 过大 → 过度平滑，丢失局部信号

window 过小 → 去噪效果不足

不是预测工具：后验平滑，让历史数据更清晰

推荐起始参数：window=11, poly_order=3

最后的话

SG 滤波就像一位精修师：

高斯平滑把照片整体磨了一层，五官和皮肤一起变模糊。
SG 滤波先看清五官轮廓，再精准磨掉皮肤瑕疵。
结果：噪音少了，趋势信号还在，波峰波谷不变形。

一句话：想要保留信号形状的平滑，选 SG 滤波。

附录

完整手搓代码

PYTHON · 完全手搓 · 不依赖 scipy

"""
Savitzky-Golay 滤波器 —— 纯 NumPy 手写实现
============================================
不调用 scipy.signal.savgol_filter，从零推导每一步。

核心思想：
  在局部窗口内，不是简单"取平均"，
  而是先"画一条最合适的多项式曲线"，
  再取曲线中心点的值作为平滑结果。
"""

import numpy as np
import math

# ===================================================================
# 第一部分：核心工具函数
# ===================================================================

def _build_coeffs(window, poly_order, deriv=0):
    """
    构建 SG 滤波的卷积核。

    推导过程：
      1. 窗口内有 N=window 个点，相对位置 k = [-m, ..., 0, ..., m]
      2. 用 n 阶多项式 y(k) = a0 + a1*k + ... + an*k^n 去拟合
      3. 最小二乘：min Σ[y_true(k) - y_fit(k)]²
      4. 求导令为 0 → (V^T·V)·a = V^T·y（正规方程）
      5. 解出 a 后，中心点 k=0 的拟合值 = a0
         而 a0 恰好是原始数据的线性组合！
      6. 所以系数 c[k] = (V^T·V)^{-1}·V^T 的第 deriv 行
    """
    if window % 2 == 0:
        raise ValueError("窗口长度必须是奇数（要有明确的中心点）")
    if poly_order >= window:
        raise ValueError("多项式阶数必须 < 窗口长度")

    half = (window - 1) // 2
    k = np.arange(-half, half + 1, dtype=float)

    # 范德蒙矩阵：V[i,j] = k[i]^j
    V = np.vander(k, poly_order + 1, increasing=True)

    # 最小二乘解：直接对 V 做伪逆，避免 (V^T @ V) 条件数恶化
    A = np.linalg.pinv(V)

    # 取对应导数阶数的行
    coeffs = A[deriv]

    # 平滑操作：归一化使系数和为 1
    if deriv == 0:
        coeffs = coeffs / np.sum(coeffs)
    else:
        # 导数操作：乘以阶乘
        coeffs = coeffs * math.factorial(deriv)

    return coeffs

# ===================================================================
# 第二部分：SG 平滑滤波器
# ===================================================================

def savgol_smooth(x, window=11, poly_order=3):
    """
    Savitzky-Golay 平滑滤波器（纯手写）。

    参数：
        x           : 输入的一维信号数组
        window      : 滑动窗口长度（奇数）
        poly_order : 多项式阶数（2 或 3 最常用）

    返回：
        y : 平滑后的信号，长度与输入相同
    """
    x = np.asarray(x, dtype=float)
    if x.ndim != 1:
        raise ValueError("只支持一维信号")

    # Step 1: 构建卷积核
    coeffs = _build_coeffs(window, poly_order, deriv=0)

    # Step 2: 边界填充（最近邻）
    pad = window // 2
    x_padded = np.pad(x, pad, mode='nearest')

    # Step 3: 相关性计算（⚠️ 用 correlate 不用 convolve）
    # convolve 会翻转核，导致求导时符号相反
    y = np.correlate(x_padded, coeffs, mode='valid')

    return y

# ===================================================================
# 第三部分：SG 求导
# ===================================================================

def savgol_deriv(x, window=11, poly_order=3, deriv=1, delta=1.0):
    """
    Savitzky-Golay 求导。

    deriv=1 : 一阶导数（速度/斜率）
    deriv=2 : 二阶导数（加速度/曲率）
    delta   : 横坐标步长（默认 1.0）
    """
    x = np.asarray(x, dtype=float)

    # Step 1: 构建导数卷积核
    coeffs = _build_coeffs(window, poly_order, deriv=deriv)

    # Step 2: 边界填充
    pad = window // 2
    x_padded = np.pad(x, pad, mode='nearest')

    # Step 3: 相关性计算（⚠️ 用 correlate，不用 convolve）
    y = np.correlate(x_padded, coeffs, mode='valid')

    # 若横坐标步长不为 1，除以 delta^deriv
    y = y / (delta ** deriv)

    return y

# ===================================================================
# 第四部分：可视化验证
# ===================================================================

if __name__ == "__main__":
    import matplotlib.pyplot as plt

    # 构造测试信号
    np.random.seed(42)
    t = np.linspace(0, 4 * np.pi, 200)

    # 真实信号：正弦波 + 两个尖峰
    true_signal = (
        2.0 * np.sin(t)
        + 3.0 * np.exp(-((t - 5)**2) / 0.5)
        + 2.5 * np.exp(-((t - 8)**2) / 0.3)
    )

    # 加噪音
    noise = np.random.normal(0, 0.4, size=t.shape)
    noisy_signal = true_signal + noise

    # 对比：普通滑动平均 vs SG 滤波
    w = 15
    pad = w // 2
    ma = np.correlate(
        np.pad(noisy_signal, pad, mode='nearest'),
        np.ones(w) / w,
        mode='valid'
    )[:len(t)]

    sg = savgol_smooth(noisy_signal, window=w, poly_order=3)

    # 画图
    fig, axes = plt.subplots(2, 1, figsize=(12, 8), sharex=True)

    ax1 = axes[0]
    ax1.plot(t, noisy_signal, 'lightgray', label='带噪信号', alpha=0.8)
    ax1.plot(t, true_signal, 'black', line, linewidth=2, label='真实信号')
    ax1.plot(t, ma, 'red', linewidth=2, label='滑动平均（峰被削平）')
    ax1.set_title('滑动平均的问题：峰被削平，拐点变圆')
    ax1.legend(); ax1.grid(True, alpha=0.3)

    ax2 = axes[1]
    ax2.plot(t, noisy_signal, 'lightgray', label='带噪信号', alpha=0.8)
    ax2.plot(t, true_signal, 'black', line, linewidth=2, label='真实信号')
    ax2.plot(t, sg, 'blue', linewidth=2, label='SG 滤波（保峰保谷）')
    ax2.set_title('SG 滤波的优势：保峰高、保拐点、保波形')
    ax2.legend(); ax2.grid(True, alpha=0.3)
    ax2.set_xlabel('采样点')

    plt.tight_layout()
    plt.savefig('savgol_handmade.png', dpi=150)
    print("✅ 对比图已保存为 savgol_handmade.png")

    # 定量对比
    mae_ma = np.mean(np.abs(ma - true_signal))
    mae_sg = np.mean(np.abs(sg - true_signal))
    print(f"\n平均绝对误差 (MAE)：")
    print(f" 滑动平均 : {mae_ma:.4f}")
    print(f" SG 滤波 : {mae_sg:.4f}")
    print(f" 提升幅度 : {(1 - mae_sg/mae_ma)*100:.1f}%")