常见量化算法系列 - Savitzky-Golay 滤波
-
-
阿凡 2025，大牛市快来了 V：flj99103681
0 人点赞了该文章 · 269 浏览
常见量化算法系列 Savitzky-Golay 滤波：保峰的"精修磨皮"文 / 阿凡 |
如果说高斯平滑是"美图秀秀一键磨皮"，那 Savitzky-Golay（简称 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₂]，使得多项式曲线与窗口内实际数据点的误差平方和最小。
这就是最小二乘法的标准形式：
最小二乘正规方程 (VT · V) · a = VT · 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 为例：
|
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)] |
| ♦ |
| 运算流程回顾 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) # 二阶导（加速度） |
类比开车： |
| 为什么求导有用？ 一阶导数过零点 → 价格达到局部最高/最低点 二阶导数过零点 → 趋势拐点（由涨转跌、由跌转涨） |
杀 手 级 功 能 同一个范德蒙矩阵，取不同行就能得到平滑、一阶导、二阶导——一个算法，多种用途。 |
| ♦ |
三层拆解回顾：
第一层：5 点 2 阶卷积核 归一化系数：[-3, 12, 17, 12, -3] ÷ 35 中心权重最高（17/35），最外侧两个点是负权重（-3/35）——这是保峰的关键。 |
第二层：窗口内 5 个采样值 y₋₂ y₋₁ yᵢ yᵢ₊₁ yᵢ₊₂ |
第三层：加权求和 |
最终公式 yismooth = (-3·yi-2 + 12·yi-1 + 17·yi + 12·yi+1 - 3·yi+2) ÷ 35 |
| 系数提前固定：只由窗口长度和多项式阶数决定，和数据无关。 ⚡ 高效：预计算核 + 滑动卷积，适合实时处理。 🛡️ 负权重保峰：最外侧的 -3 抵消峰衰减，维持极值高度。 |
| ♦ |
输入带噪数据 → 设定参数 → 计算卷积核 → 滑动窗口 → 边界处理 → 输出平滑序列
| 1 | 输入一维带噪时序数据 y[1:N] |
| 2 | 设定：窗口长度 L=2m+1、多项式阶数 n（常用 2/3 阶） |
| 3 | 构建范德蒙矩阵 → 最小二乘求卷积核 |
| 4 | 滑动窗口遍历：① 截取局部数据 ② 卷积核加权求和 |
| 5 | 边界处理（最近邻填充） |
| 6 | 输出完整平滑序列 |
| ♦ |
|
|
|
| ♦ |
第三部分
量化策略中的用途
| 场景 | 用法 |
| 趋势跟踪 | SG 平滑后的价格作为趋势基准，比均线更灵敏 |
| 拐点检测 | 二阶导数过零点 = 潜在买卖信号 |
| 指标预处理 | 先 SG 平滑再算 RSI/MACD，减少假信号 |
| 多周期分析 | 不同 window 对应不同时间尺度 |
SG 滤波 vs 加权平均
| 特性 | 加权平均 | SG 滤波 |
| 原理 | 加权平均 | 局部多项式拟合 |
| 波峰保留 | 会削低 | ✅ 基本保留 |
| 波谷保留 | 会填高 | ✅ 基本保留 |
| 能否求导 | 需额外处理 | ✅ 内置求导 |
| 计算速度 | 快 | 稍慢 |
| 适用场景 | 简单去噪 | 精细分析、拐点检测 |
| ♦ |
第四部分
注意事项
window 必须为奇数：偶数没有明确中心点 |
poly_order < window：否则方程欠定（未知数比方程多） |
window 过大 → 过度平滑，丢失局部信号 window 过小 → 去噪效果不足 |
不是预测工具：后验平滑，让历史数据更清晰 |
推荐起始参数： |
最 后 的 话 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}%") |
5 天前