Savitzky-Golay滤波

一、直观理解：在滑动窗口里做局部多项式拟合

Savitzky-Golay 滤波（简称 SG 滤波）可以理解为：不是简单把窗口里的点平均掉，而是在每个局部窗口内拟合一条低阶多项式，再取该多项式在窗口中心的值作为平滑结果。

对比很直观：

• 移动平均假设“窗口内信号近似是常数”，所以用均值；
• SG 滤波假设“窗口内信号近似是低阶曲线”，所以用多项式拟合。

因此，当信号本身有峰、谷、上升沿时，SG 往往比同宽度的移动平均更能保留峰高和局部斜率，同时仍能抑制随机噪声。

一个常用类比是：用直尺去贴合一段轻微弯曲的绳子。若只用平均高度，峰会被削平；若允许贴合一条抛物线，峰顶附近的形状就能保留得更好。

SG 滤波特别常见于：

• 光谱、色谱等带尖峰的实验曲线；
• 需要后续计算一阶/二阶导数的信号；
• 希望去噪但又不想把峰位置和峰面积破坏太多的场景。

二、适用场景与优缺点

类型	说明
适用场景	等间隔采样的一维信号；噪声近似随机且频率高于信号变化；需要保留峰形、拐点或导数信息；移动平均削峰过重时的替代方案。
核心优点	在同类线性平滑中对多项式状信号的保真较好；可同时得到平滑值与数值导数；参数（窗口、阶数）含义清楚。
局限与缺点	窗口与阶数选不好会欠平滑或过拟合噪声；两端边界效应明显；对突变台阶、脉冲噪声不如中位数类方法稳健；非等间隔数据需先整理。

使用前建议确认：

1. 数据基本等间隔，缺失已插补或已剔除；
2. 噪声主要是高频抖动，而不是系统性偏差；
3. 真正关心的结构可以用低阶多项式在局部近似；
4. 峰宽与采样间隔的关系已知，便于设定窗口；
5. 若用于监督学习特征，滤波参数只在训练阶段选定。

若序列目标是“只要大趋势、峰形不重要”，移动平均往往更简单；若目标是“去噪且尽量保峰/保导”，再优先考虑 SG。

三、核心原理：局部最小二乘多项式

设窗口长度为奇数 $w=2m+1$，多项式阶数为 $d$，且必须满足 $d<w$。

对位于中心时刻 $t$ 的窗口观测

$$\mathbf y_t=\big(y_{t-m},\ldots,y_t,\ldots,y_{t+m}\big)^\mathsf{T}$$

在局部坐标 $u=-m,\ldots,m$ 上拟合

$$p(u)=a_0+a_1 u+\cdots+a_d u^d$$

使残差平方和最小：

$$\min_{a_0,\ldots,a_d}\sum_{u=-m}^{m}\big(y_{t+u}-p(u)\big)^2$$

平滑值取中心点：

$$\hat y_t=p(0)=a_0$$

若需要一阶导，则

$$\widehat{y'_t}=p'(0)=a_1$$

（实际应用中还要除以采样间隔 $\Delta t$，把“每采样点”换算成物理单位。）

矩阵形式与卷积实现

令设计矩阵 $A$ 的元素为 $A_{u,j}=u^j$（$u=-m,\ldots,m$，$j=0,\ldots,d$），则

$$\hat{\mathbf a}=(A^\mathsf{T}A)^{-1}A^\mathsf{T}\mathbf y_t$$

由于每个窗口都在做同类最小二乘，SG 等价于用一组固定卷积系数对原序列做滤波。这也是它计算快、结果稳定的原因。

与移动平均的关系

当 $d=0$ 时，局部“多项式”就是常数，SG 退化为普通移动平均。$d$ 提高后，模型能跟踪更弯曲的局部形状，但也更容易把噪声当成弯曲拟合进去。

四、关键决策：窗口与阶数

1. 窗口长度 $w$

• $w$ 必须为奇数，保证有唯一中心点；
• $w$ 应覆盖若干个噪声相关长度，但尽量小于主要峰的半宽；
• 采样率变化后，同一 $w$ 对应的物理平滑尺度会变，不能照搬旧参数。

经验上可从较小窗口试起（如 5、7、9、11），观察峰高保留与残差平稳性。

2. 多项式阶数 $d$

• 常用 $d=2$ 或 $d=3$；
• $d$ 太低：峰顶被削、弯曲跟不上；
• $d$ 太高：接近“插值噪声”，平滑名存实亡；
• 必须保证 $d\le w-1$，实务中常取 $d\ll w$。

3. 导数阶数

若目标是求导数，不要先随便平滑再做差分；应直接使用 SG 的导数输出，并对 $\Delta t$ 做正确换算。导数会放大高频噪声，窗口往往需要比只做平滑时更大一些。

4. 边界策略

序列两端无法构成完整对称窗口。MATLAB sgolayfilt 会对端部采用特殊处理，但端部结果仍比中段更不可靠。若端部恰好是研究重点，应报告其不确定性，或改用更合适的边界模型。

五、MATLAB 示例

% 构造带峰信号：两个高斯峰 + 噪声
rng(2);
t = linspace(0, 10, 401)';
yTrue = exp(-(t-3.5).^2/0.15) + 0.7*exp(-(t-6.5).^2/0.25);
y = yTrue + 0.08*randn(size(t));

% SG 参数
polyOrder = 3;     % 多项式阶数 d
frameLen  = 15;    % 窗口长度 w（奇数，且 > d）

ySG = sgolayfilt(y, polyOrder, frameLen);

% 对照：同窗口移动平均（更容易削峰）
yMA = movmean(y, frameLen);

figure;
plot(t, y, ".", "DisplayName", "含噪观测");
hold on;
plot(t, yTrue, "k-", "LineWidth", 1.2, "DisplayName", "真值");
plot(t, ySG, "LineWidth", 1.8, "DisplayName", "SG 滤波");
plot(t, yMA, "LineWidth", 1.5, "DisplayName", "移动平均");
xlabel("t");
ylabel("信号");
title("Savitzky-Golay 与移动平均的保峰对比");
legend("Location", "best");
grid on;

% 峰高保留程度（以第一个峰附近最大值近似）
idx = t > 2.5 & t < 4.5;
fprintf("真值峰高     = %.4f\n", max(yTrue(idx)));
fprintf("SG 峰高      = %.4f\n", max(ySG(idx)));
fprintf("移动平均峰高 = %.4f\n", max(yMA(idx)));

% 如需一阶导数：用 sgolay 生成系数后卷积
[b, g] = sgolay(polyOrder, frameLen);
dt = t(2) - t(1);
% g 的第 2 列对应一阶导卷积核（还要按阶乘与 dt 换算）
yDeriv = conv(y, factorial(1) * g(:, 2) / dt^1, "same");

figure;
plot(t, yDeriv, "LineWidth", 1.5);
xlabel("t");
ylabel("dy/dt（SG 估计）");
title("SG 数值一阶导");
grid on;

说明：sgolayfilt 需要 Signal Processing Toolbox。若环境受限，可按第三节最小二乘公式对每个窗口求解 $a_0$，核心思想相同。

六、如何解释结果

SG 的输出首先是平滑后的序列，有时还包含导数序列。解释时应围绕“噪声被抑制了多少、形状被改变了多少”。

建议报告：

1. 采样间隔与序列长度；
2. 窗口长度、多项式阶数，以及它们对应的物理时间尺度；
3. 平滑前后噪声水平（如残差标准差）与峰高/峰位变化；
4. 是否计算了导数，导数单位如何换算；
5. 边界段是否裁剪或降权；
6. 与移动平均或其他滤波的对比结论；
7. 滤波结果用于可视化、特征还是机理分析。

若峰面积、峰位置是评价指标，必须定量给出滤波前后的偏差，不能只附一张“更光滑”的图。

七、常见误区

1. 把窗口写成与采样率无关的“魔法数”

同样 frameLen=15，在 1 Hz 与 100 Hz 采样下物理含义完全不同。应写成“约等于 xxx 秒/分钟的局部拟合窗”。

2. 阶数盲目取高

阶数接近窗口长度时，SG 几乎穿过所有点，失去滤波意义，只是在做昂贵的恒等变换。

3. 把滤波曲线当作原始测量

任何线性滤波都会改变噪声谱，也可能轻微平移峰位、改变峰高。讨论测量不确定度时要用原始数据或明确传播滤波影响。

4. 忽略边界

端部点的有效信息少于中部，赛题若关注起步段/收尾段，不要把边界平滑值与中段等量齐观。

5. 对脉冲噪声硬套 SG

SG 本质是最小二乘，对离群脉冲敏感。先识别并处理脉冲，或改用稳健平滑，通常更合适。

八、小结

Savitzky-Golay 滤波的核心思想可以概括为：

在每个滑动窗口内用低阶多项式做局部最小二乘拟合，以中心拟合值作为平滑结果，从而在去噪的同时更好保留峰形与局部斜率。

它是移动平均的自然升级：当信号局部更像曲线而不是水平线时，SG 通常更合适。使用关键在于让窗口匹配峰宽与噪声尺度，并让多项式阶数保持克制。只要把参数、边界和保峰效果写清楚，SG 就是数学建模中非常稳妥的一维信号预处理工具。