Savitzky–Golay滤波器与IIR或FIR线性滤波器

#1 楼

由于现有答案和评论中的讨论主要集中在Savitzky-Golay滤波器实际上是什么（非常有用）上，因此我将通过提供一些有关如何实际选择平滑滤波器的信息来尝试添加到现有答案中，据我了解，问题实际上是关于什么的。

首先，我想重复讨论中从其他答案中得出的清楚的内容：一方面，将问题引入线性和时不变（LTI）FIR / IIR滤波器，另一方面将其转化为Savitzky-Golay滤波器是一种误导。 Savitkzy-Golay滤波器只是根据特定条件（局部多项式逼近）设计的标准FIR滤波器。所以问题中提到的所有滤波器都是LTI滤波器。

剩下的问题是如何选择平滑滤波器。如果计算复杂性和/或存储是一个问题，IIR滤波器可能比FIR滤波器更可取，因为IIR滤波器通常会以比FIR滤波器低得多的滤波器阶数实现相当的噪声抑制（即阻带衰减）。但是请注意，如果需要实时处理，IIR滤波器的一个可能的缺点是它们不能具有精确的线性相位响应。因此，所需信号将遭受一些相位失真。对于离线处理，即使使用IIR滤波器，也可以通过应用零相位滤波来避免相位失真。

除了上一段中讨论的考虑因素外，主要还是要考虑设计准则，而滤波器是FIR还是IIR则不是那么重要，因为FIR滤波器可以任意精度近似任何（稳定的）IIR滤波器，并且FIR滤波器可以由IIR滤波器近似，即使后者可能要困难得多。适当的设计标准显然取决于数据和噪声的属性。在进行平滑处理时，我们通常会假设数据已充分过采样（即平滑）。如果噪声主要具有高频成分，即，如果数据和噪声之间几乎没有频谱重叠，则我们希望最大程度地降低阻带衰减或最小化阻带能量，同时尽可能保留所需信号。在这种情况下，我们可以选择使用Parks-McClellan算法根据极小极大准则设计的线性相位FIR滤波器。我们还可以通过选择最小二乘法来最小化阻带能量（即，最小化阻带中的噪声功率）。通过选择受约束的最小二乘设计，可以在两个标准（最小二乘和最小二乘）之间进行混合，这可以最小化阻带能量，同时限制通带中的最大近似误差。

如果噪声频谱与信号频谱明显重叠，需要采取更谨慎的方法，并且蛮力衰减将无法很好地工作，因为您会留下过多的噪声（通过选择过高的截止频率），或者会使所需的信号失真太多。在这种情况下，Savitzky-Golay（S-G）滤波器可能是一个不错的选择。要付出的代价是中等的阻带衰减，但是优点之一是很好地保留了某些信号特性。这与以下事实有关：SG滤波器具有平坦的通带响应，即，

$$ \ frac {d ^ kH（e ^ {j \ omega}）} {d \ omega ^ k} {\ huge |} _ {\ omega = 0} = 0 \ quad k = 1,2，\ ldots，r \ tag {1} $$

其中$ r $是近似多项式的阶数，而$ H（e ^ {j \ omega}）$是滤波器的频率响应。属性$（1）$保证输入信号的前$ r $矩保留在输出中，这意味着所需信号中的峰的宽度和高度都得到了很好的保留。

当然，上面讨论的两种类型的平滑滤波器（高阻带衰减和SG）之间也存在折衷。我们可以设计一个具有一定平坦度的FIR滤波器，其\ω= 0，并使用剩余的自由度来最大化阻带衰减或最小化阻带能量。对于FIR滤波器，所产生的设计问题非常简单（且凸），并且可以使用几个软件包中可用的常规优化例程来获得给定应用程序的最佳滤波器。

对SG滤波器理论感兴趣的人，以下是我可以推荐的最相关的参考文献：




Savitzky和Golay的原始论文

罗纳德·W·舍弗（Ronald W. Schafer）的教程文件
Orfanidis的DSP书籍的8.3.5节

#2 楼

像其他任何东西一样，有时某些工具比其他工具更好。

移动平均（MA）滤波器可用于平滑数据，并且是FIR。它们几乎是您可以想到的最简单的过滤器，并且只要您不打算对任何突然的跳跃或多项式趋势建模，它们就可以很好地完成许多任务。请记住，尽管它们本质上只是低通滤波器，所以当您在信号中关注的数据是低频且相当窄的带宽时，它们会发挥最佳作用。

Savitzky-Golay（SG）滤波器是一组特殊的FIR滤波器，当卷积沿信号滑动时，它们基本上适合您的时间序列的多项式。 SG滤波器可用于信号，在信号中您所关心的不一定是低频和相当窄的频带。

我认为您会发现，如果您阅读Wikipedia页面，您已经彻底链接了，它非常详细地说明了SG和MA滤波器之间的区别。不过请记住：最后，它们都是做类似事情的工具，这取决于您确定何时使用正确的工具

编辑：

由于似乎有些困惑，认为SG过滤器在某种程度上是“自适应”的，所以我提供了一个简单的MATLAB示例。正如Dan所指出的那样，可以使它们具有适应性，但是它们的基本实现常常不是。通过检查代码，您将看到这只是一个矩阵，并进行了一些特殊处理。在传统意义上，关于此滤波器的任何内容都不是“自适应的”，您只是在选择多项式拟合，以及通过卷积将多项式拟合到信号上的长度。 SG是必不可少的FIR。我下面的脚本产生了这个图：


看这个图，您可以看到MA和SG本质上完成了相同的事情，但是有一些重要的区别：


MA在抑制噪声方面做得很好，但是在捕获信号中的瞬态跳变方面做得很差。我们可以通过增加滤波器的长度来抑制更多的噪声，但是在瞬变中它的表现会更差。在任何瞬态附近，这种效果都会被视为“拖影”，您应该可以在所示的图中看到。
SG在捕获信号瞬态方面做得很好，但是做得并不好抑制噪声（至少与相同大小的MA相比）。我们可以通过增加帧长来改善非瞬态附近的噪声抑制，但这会在任何瞬态附近引入类似于Gibb现象的振铃。

为了让您更好地了解这些滤波器的工作原理，我鼓励您在此处使用代码并进行操作，并查看SG过滤器的所有各个部分如何工作。

MATLAB示例的代码：

% Generate a signal "s" that has square waves, and scale it with a
% polynomial of order 5
up = 1*ones(1,100);
down = zeros(1,150);
s = [down down up up down up down up down up up up down down down down down];
n = numel(s);
nn = linspace(0,4,numel(s));
s = s .* (nn .^5);
sn = (s + 4*randn(size(s))).';

% smooth it with SMA of length 15
sz = 15;
h = 1/sz * ones(1,sz);
sn_sma = conv(sn,h,'same');

% smooth it with sgolay of frame length 15
m = (sz-1)/2;
% look up the SG matrix for this order and size
B = sgolay(5, sz);

% compute the steady state response for the signal, i.e. everywhere that
% isnt the first or last "frame" for the SG filter
steady = conv(sn, B(m+1,:), 'same');
% handle the transiet portion at the start of the signal
y_st   = B(1:m,:)*sn(1:sz);
% handle the transient portion at the end of the signal
y_en   = B(sz -m+1 : sz, :) * sn(n - sz+1:n);

% combine our results
sn_sg  = steady;
sn_sg(1:m) = y_st;
sn_sg(n-m+1:n) = y_en;

% plots
figure(1);
hold off;
plot(sn,'Color',[0.75 0.75 0.75]);
hold on;
plot(sn_sma,'b');
plot(sn_sg,'r');

legend('Noisy Signal','MA Smoothing','SG Smoothing, order 5','Location','NorthWest');

#3 楼

注意

我先前的回答（在此编辑之前）将Savitzky-Golay（SG）滤波器表示为非线性的，随时间变化的输入数据相关项是错误的，原因是对Savitzky如何过早的误解-Golay（SG）过滤器根据提供的Wiki链接计算其输出。因此，现在我针对那些人的利益对其进行了纠正，他们还将看到FIR-LTI过滤如何实现S-G过滤器。感谢@MattL。为了纠正他的问题，他提供了重要的联系，并在调查该问题期间表现出了他的耐心（我从未表现出过）。虽然老实说，我更喜欢冗长的反对，但这显然没有必要。另外
请注意，正确的答案是另一个，这个答案只是为了进一步澄清SG滤波器的LTI属性。

现在，当有人（从未使用过）之前的那些滤波器面对SG滤波器的定义是对给定数据的低阶LSE多项式拟合，他/她将立即得出结论，即那些滤波器是与数据相关的，非线性的和随时间（移位）变化的自适应滤波器。

然而，多项式拟合过程可以由SG自己巧妙地解释，因此它可以实现完全独立于数据的时不变线性滤波，因此使SG成为固定的LTI-FIR滤波器。 br />
下面是MattL提供的链接中的最短摘要。有关似乎缺少的任何详细信息，请查阅原始文档，或要求进行澄清。但是我不想在这里重新生成整个文档。

现在考虑$ 2M + 1 $输入数据值$ x [-M]，x [-M + 1]，... ，x [0]，x [1]，...，x [M] $，它们以$ n = 0 $为中心，并且我们希望将阶次为$ N $的多项式$ p [n] $拟合到其中，其中$ n = -M，-M + 1，...，-1,0,1，... M $是整数时间索引：

$$ p [n] = \ sum_ {k = 0} ^ {N} a_k n ^ k = a_0 + a_1 n + a_2 n ^ 2 + ... + a_N n ^ N $$

经典的LSE多项式拟合过程计算这些系数$ a_k $，以找出使误差平方和最小的最优$ N ^ {th} $阶多项式$ p [n] $。 mathcal {E} = \ sum _ {-M} ^ {M}（p [n] -x [n]）^ 2 $$

在给定的数据向量$ x = \ left [x [-M]，x [-M + 1]，...，x [0]，x [1]，...，x [M] \ right] ^ {T} $。

通过将$ \ mathcal {E} $的导数设置为零，可以获得这些最优多项式系数$ a_k $：

$$ \ frac {\ partial \ mathcal {E}} {\ partial a_i} = 0 ~~~，~~~ \ text {for} ~~~ i = 0,1，..，N \ tag {1} $$

现在供那些熟悉的人使用使用LSE polyfit过程，我将简单地写出定义最佳系数集的矩阵方程式（来自链接）：

$$ \ boxed {a =（A ^ {T} A）^ {-1} A ^ {T} x = H x} \ tag {2} $$，其中$ x $是$（2M + 1）\乘以1 $的输入数据向量，$ H $是LSE polyfit矩阵， $ 2M + 1 $ x $ N $矩阵$ A $是时间瞬时矩阵（整数瞬时$ n $的幂）；即，请注意，$ A $和因此的$ H $都是输入数据值的独立变量，因为$ A $由下式给出：

$$ A =
\ begin {bmatrix}
\ alpha_ {n，i}
\ end {bmatrix}
=
\ begin {bmatrix}
（-M）^ 0＆（-M）^ 1＆ ...＆（-M）^ N \\
（-M + 1）^ 0＆（-M + 1）^ 1＆...＆（-M + 1）^ N \\
＆... \\
（0）^ 0＆（0）^ 1＆...＆（0）^ N \\
＆... \\
（ M）^ 0＆（M）^ 1＆...＆（M）^ N \\
\ end {bmatrix}

现在让我们往后靠一会儿，在这里讨论要点。

正如式（2）清楚地表明的那样，即使$ A $和$ H $与输入数据无关，并且仅取决于时间索引$ n $，最佳LSE多项式系数$ a_k $还是取决于输入数据。此外，它们还取决于窗口$ M $的大小和多项式$ N $的阶数。此外，随着窗口沿输入数据$ x [n] $滑动，系数$ a_k $应该重新计算（更新），因此也将与时间有关。这就是在下面摘录的$ 2 ^ {nd} $页面上的链接中定义SG过滤器的方式：

...此（LSE polyfit）可以在输入的每个样本中重复进行，每次产生一个新的多项式和一个新的输出序列y [n] ...

那么我们如何克服这个令人困惑的惊喜呢？通过将以下内容定义和定义SG过滤器输出：

顺序为$ N $的SG过滤器，对于每次采样$ n $，接受输入集$ x [n] $并产生单个输出样本$ y [n] $定义为以$ n = 0 $求值的多项式$ p [n] $；即，

$$ \ boxed {y [n] = y [0] = \ sum_ {m = 0} ^ {N} a_m n ^ m = a_0} $$

也就是说，对于$ x [n] $（以$ n = d $为中心）的$ 2M + 1 $个样本的每个输入集，SG滤波器生成输出，表示为$ y [n] $，等于与样本$ x [n] $的特定窗口相关联的最佳LSE多项式$ p [n] $的单个系数$ a_0 $。顺便说一下，当窗口沿输入数据长度滑动时，每次$ n = d $时，都会根据样本$ x [dM]，x [d-M + 1]的窗口来计算新的输出样本$ y [d] $ ]，...，x [d-1]，x [d]，x [d + 1]，... x [d + M] $。这是一个鼻孔过滤器。

现在是时候显示系数$ a_0 $是作为窗口中输入信号值$ x [n] $的线性组合而获得的，因此滤波器输出$ y [n] $是输入值的线性组合$ x [n] $。这就是通过FIR滤波器对（LTI）卷积的定义。时间$ n $的输出是其输入$ x [n] $和滤波器系数$ h [n] $的线性组合。但是，此S-G滤波器的滤波器系数是多少？让我们来看一下。

再次考虑$ a_k $的计算：

$$ a = H x $$

$$
\ begin {bmatrix}
a_0 \\
a_1 \\
\ vdots \\
a_N \\
\ end {bmatrix}
=
/> \ begin {bmatrix}
h（0,0）＆h（0,1）＆...＆h（0,2M）\\
h（1,0）＆h（ 1,1）＆...＆h（1,2M）\\
＆... \\
h（N，0）＆h（0,1）＆...＆h （0,2M）\\
\ end {bmatrix}
\ cdot
\ begin {bmatrix}
x [-M] \\
x [-M +1] \\
... \\
x [M] \\
\ end {bmatrix}
$$

可以看到单个系数$ a_0 $由$ H $第一行的点积与输入数据矢量$ x $给出；即，

$$ a_0 = H（0，n）\ cdot x = \ sum H（0，k）x [k] = H（0，-n）\ star x [n] $$

其中，在最后一个等式中，我们通过考虑SG滤波器的脉冲响应为$$ \ boxed {h [n] = H（ 0，-n）} $$，

更具体地说，它是阶为$ N $的SG滤波器的脉冲响应，窗口长度为$ 2M + 1 $。

对于一个长度为$ L $的输入$ x [n] $，通过一个LTI卷积可以得到窗口大小为$ 2M + 1 $的N阶SG滤波器的完整输出$ y [n] $。脉冲响应$ h_N [n] $

$$ \ boxed {y [n] = x [n] \ star h_N [n]} $$

COMMENT

多项式系数$ a_k $取决于输入数据的事实并不能防止滤波器成为LTI FIR。因为可以定义脉冲响应$ h [n] $来表示要从输入样本的线性组合中计算出的输出$ y [n] $。输入样本$ x $的线性组合固有地由定义$ p [n] $的最佳系数$ a_k $的矩阵乘积$ a = H x $隐含，因此$ a_k $的任何线性组合也将导致： FIR LTI过滤器$ h [n] $表示LSE多项式拟合过程。

MATLAB / OCTVE代码

以下简单的MATLAB / OCTAVE可用于计算那些SG滤波器脉冲响应$ h [n] $（请注意，它的内置SG设计器可能会产生由link-pdf概述的$ h [n] $的不同集合）

% Savitzky-Golay Filter
% 
clc; clear all; close all;

N = 3;                      % a0,a1,a2,a3 : 3rd order polynomial
M = 4;                      % x[-M],..x[M] . 2M + 1 data

A = zeros(2*M+1,N+1);
for n = -M:M
    A(n+M+1,:) = n.^[0:N];
end

H = (A'*A)^(-1)* A';        % LSE fit matrix

h = H(1,:);                 % S-G filter impulse response (nancausal symmetric FIR)

figure,subplot(2,1,1)
stem([-M:M],h);
title(['Impulse response h[n] of Savitzky-Golay filter of order N = ' num2str(N), ' and window size 2M+1 =  ' , num2str(2*M+1)]);

subplot(2,1,2)
plot(linspace(-1,1,1024), abs(fftshift(fft(h,1024))));
title('Frequency response magnitude of h[n]');

输出为：



希望这可以澄清问题。