9
$\begingroup$
  • A traditional IIR / FIR filter (lowpass to remove the high freq oscillations), e.g. moving average,

  • or a Savitzky-Golay filter

can all be useful to smoothen a signal, such as an envelope signal:

enter image description here

For which application would a Savitzky-Golay filter be more interesting than a classical lowpass?

What makes it different than a standard filter, and what does it add in comparison to standard filters?

Does it adapt itself to the input data?

Is it better for transient preserving?


Have you ever been in an engineering situation one day when you decided "Let's use a SG filter instead of moving average or another FIR lowpass! It's better because this and this and this..."? Then this question is for you!

$\endgroup$
3
$\begingroup$

Since the discussion in the existing answers and comments has mainly focused on what Savitzky-Golay filters actually are (which was very useful), I will try to add to the existing answers by providing some information on how to actually choose a smoothing filter, which is, to my understanding, what the question is actually about.

First of all, I'd like to repeat what has become clear in the discussion spawning from the other answers: the categorization of smoothing filters in the question into linear and time-invariant (LTI) FIR/IIR filters on the one hand, and Savitzky-Golay filters on the other hand is misleading. A Savitkzy-Golay filter is just a standard FIR filter designed according to a specific criterion (local polynomial approximation). So all filters mentioned in the question are LTI filters.

The remaining question is how to choose a smoothing filter. If computational complexity and/or memory are an issue, IIR filters may be preferable over FIR filters, because they typically will achieve comparable noise suppression (i.e., stopband attenuation) with a much lower filter order than FIR filters. But note that if real-time processing is necessary, one possible disadvantage of IIR filters is that they cannot have an exactly linear phase response. So the desired signal will suffer some phase distortions. For offline-processing, the phase distortions can be avoided, even with IIR filters, by applying zero-phase filtering.

Apart from the considerations discussed in the previous paragraph, it is mainly the design criterion that matters, not so much if the filter is FIR or IIR, because any (stable) IIR filter can be approximated with arbitrary accuracy by an FIR filter, and any FIR filter can be approximated by an IIR filter, even though the latter can be much more difficult. The appropriate design criterion obviously depends on the properties of the data and the noise. When it comes to smoothing we usually assume sufficiently oversampled (i.e., smooth) data. If the noise has mainly high frequency components, i.e., if there is little spectral overlap between the data and the noise, we want to maximize stop band attenuation, or minimize stop band energy, while preserving the desired signal as well as possible. In this case we could choose a linear-phase FIR filter designed according to a minimax criterion using the Parks-McClellan algorithm. We could also minimize the stop band energy (i.e., minimize noise power in the stop band) by choosing a least squares method. A mix between the two criteria (minimax and least squares) is possible by choosing a constrained least squares design, which minimizes stop band energy while constraining the maximum approximation error in the pass band.

If the noise spectrum significantly overlaps with the signal spectrum, a more careful approach is required, and brute-force attenuation will not work well because either you leave too much noise (by choosing the cut-off frequency too high) or you distort the desired signal too much. In this case Savitzky-Golay (S-G) filters may be a good choice. The price to be paid is mediocre stopband attenuation, but one advantage is that some signal properties are preserved very well. This has to do with the fact that S-G filters have a flat pass band response, i.e.,

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

where $r$ is the order of the approximating polynomial and $H(e^{j\omega})$ is the filter's frequency response. Property $(1)$ guarantees that the first $r$ moments of the input signal are preserved in the output, which means that the width and height of peaks in the desired signal are well preserved.

Of course there's also a compromise between the two types of smoothing filters discussed above (high stopband attenuation and S-G). We could design an FIR filter with a certain degree of flatness at $\omega=0$ and use the remaining degrees of freedom to maximize stop band attenuation, or minimize stopband energy. In the case of FIR filters, the resulting design problem is sufficiently simple (and convex), and general optimization routines available in several software packages can be used to obtain the optimal filter for the given application.

For the ones interested in the theory of S-G filters, the most relevant references I can recommend are the following:

$\endgroup$
2
$\begingroup$

As with anything, sometimes certain tools are better than others.

Moving average (MA) filters can be used to smooth data, and are FIR. They're pretty much the simplest filter you can come up with, and they work well for a lot of tasks as long as you aren't trying to model any sudden jumps or polynomial trends. Keep in mind though that these are essentially just low-pass filter, so they work best when the data you care about in the signal is low frequency and fairly narrow-band.

Savitzky-Golay (SG) filters are a special group of FIR filters which essentially fit a polynomial to your time series as the convolution slides along the signal. SG filters are useful for signals where the things you care about aren't necessarily low frequency and fairly narrow band.

I think you'll find that if you read the Wikipedia page you've linked pretty thoroughly, it explains the difference between SG and MA filters in quite rigorous detail. Keep in mind though: in the end, they're both just tools to do similar things, it's up to you to find out when to use the right tool

EDIT:

Since there seems to be some confusion that SG filters are "adaptive" in some way at a basic level, I've included a simple MATLAB example. As Dan as pointed out, these can be made adaptive, but their basic implementation often is not. By inspecting the code, you'll see that this is simply just a matrix look up with some special handling. Nothing about this filter is "adaptive" in the traditional sense, you simply are picking a polynomial fit and the length over which that polynomial will be fit onto the signal via convolution; SG is essential FIR. The script I have below produces this plot: SG filter comparison with MA

Looking at this figure, you can see that the MA and the SG essentially accomplish the same thing, but with some important distinctions:

  1. The MA does a great job of suppressing the noise, but it does a poor job capturing the transient jumps in the signal. We can suppress even more noise by increasing the length of the filter, but then it will perform even worse on the transients; this effect will be seen as "smearing" near any transients, which you should be able to see in the figure shown.
  2. The SG does a great job of capturing the transients of the signal, but does a not so great job of suppressing the noise (at least in comparison to a MA of the same size). We can improve noise suppression near non-transients by increasing the frame length, but this will introduce ringing analogous to Gibb's phenomenon near any transients.

For you to get a better understanding of how these filters work, I'd encourage you to take the code here and manipulate it, and see how all the individual parts of the SG filter work.

Code for the MATLAB example:

% 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');
$\endgroup$
  • 1
    $\begingroup$ The point seems to be (by seeing the other answer) that SG filter is a "fully data dependent nonlinear time-varying filter whose coefficients are re-computed for each short segment of its input". $\endgroup$ – g6kxjv1ozn Sep 27 '18 at 6:59
  • 1
    $\begingroup$ The SG filter, to my understanding from implementing it several times, is not an adaptive filter at all, especially in comparison to your average adaptive filter like an LMS or RLS. I would completely disagree with the statement that the filter weights are time varying. SG filters are essentially a table lookup, you filter with values from the table to compute a transient response, and then you go back and handle the edge cases at the start/end of the signal. I'll edit my post with a MATLAB example to show this to you. $\endgroup$ – matthewjpollard Sep 27 '18 at 11:49
  • 2
    $\begingroup$ @matthewjpollard To note, I don't personally have significant experience using this filter, but to me the SG filter as best implemented does appear very much to be an adaptive filter implementation, with time varying coefficients. The way you applied the filter in your code is not (as you treated the entire sequence as the data "sub-set"), but the way specifically as described at Wikipedia en.wikipedia.org/wiki/Savitzky%E2%80%93Golay_filter and in the paper itself by Savitzky and Golay is indeed adaptive: pdfs.semanticscholar.org/4830/… $\endgroup$ – Dan Boschen Sep 27 '18 at 12:50
  • 2
    $\begingroup$ @matthewjpollard In your real time systems is your data ever continuously streaming, so that you are recomputing coefficients over shorter intervals or do you always work in small blocks of data? $\endgroup$ – Dan Boschen Sep 27 '18 at 13:33
  • 2
    $\begingroup$ Thanks Matt. So perhaps we could associate what your doing as adaptive/time varying in the sense that coefficients are computed for each data collect (same coefficients within a data collect however with proper treatment of start and finish but varying from one collect to the next- if I understand correctly). Thanks for sharing your code and example application. $\endgroup$ – Dan Boschen Sep 27 '18 at 17:19
2
$\begingroup$

NOTE

my previous answer (before this edit) denoting the Savitzky-Golay (S-G) filter as a nonlinear, time-varying input data dependent was wrong, due to a premature mis-interpretation of how a Savitzky-Golay (S-G) filter computes its output according to the wiki link provided. So now I'm correcting it for the benefit of those who would also see how S-G filters are implementable by FIR-LTI filtering. Thanks to @MattL. for his correction, the great link he has provided and the patience he had (that I could have never shown) during my investigation of the issue. Though I'ld honestly prefer a more verbose objection which is clearly not necessary nevertheless. Also please note that the correct answer is the other one, this one is just for additional clarification of LTI property of S-G filters.

Now it's not surprising that when someone (who has never used those filters before) faces the definition of the S-G filter as a low order LSE polynomial fit to given data he/she would immediately jump into the conclusion that those are data dependent, nonlinear and time (shift) varying, adaptive filters.

Yet, the polynomial fitting procedure is cleverly interpreted by S-G themselves, such that it enables a completely data independent, time-invariant, linear filtering possible, hence making S-G as a fixed LTI-FIR filter.

The below is a shortest summary from the link provided by MattL. For any details that seems to be missing, please consult the original document, or ask to clarify. But I would not like to re-produce the whole document here.

Now consider $2M+1$ input data values $x[-M],x[-M+1],...,x[0],x[1],...,x[M]$ which are centered around $n=0$ and into which we want to fit a polynomial $p[n]$ of order $N$, with $n = -M,-M+1,...,-1,0,1,...M$ being the integer time indices:

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

The classical LSE polynomial fit procedure computes those coefficients $a_k$ to find out the optimal $N^{th}$ order polynomial $p[n]$ that minimises the sum of error squares

$$ \mathcal{E} = \sum_{-M}^{M} (p[n]-x[n])^2 $$

over the given data vector $x = \left[ x[-M],x[-M+1],...,x[0],x[1],...,x[M] \right]^{T}$.

Those optimal polynomial coefficients $a_k$ are obtained by setting the derivative of $\mathcal{E}$ to zero:

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

Now for those who are familiar with LSE polyfit procedure, I'll simply write the resulting matrix equation (from the link) that defines the optimal coefficient set:

$$ \boxed{ a = (A^{T}A)^{-1} A^{T} x = H x } \tag{2} $$ where $x$ is the $(2M+1) \times 1$ input data vector, $H$ is the LSE polyfit matrix and the $2M+1$ by $N$ matrix $A$ is the time instant matrix (powers of integer time instants $n$) ; i.e., Note that both $A$ and consequently $H$ are indpendent of input data values as $A$ is given by:

$$ 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} $$

Now lets lean back for a moment and discuss a point here.

As eq.(2) clearly indicates, eventhough $A$ and $H$ are independent of input data and depend only on the time indices $n$, the optimal LSE polynomial coefficients $a_k$ are input data dependent. Furthermore, they also depend on the size of the window $M$ and order of the polynomial $N$. Moreover, as the window slides along the input data $x[n]$, the coefficient $a_k$ should be re computed (updated) and hence will be time dependent as well. That's exactly how the S-G filter is defined in the link at the $2^{nd}$ page excerpted below:

...This (the LSE polyfit) can be repeated at each sample of the input, each time producing a new polynomial and a new value of the output sequence y[n]...

So how do we overcome this puzzling surprise? By interpreting and defining the S-G filter output to be the following:

The S-G filter of order $N$, for each time sample $n$, accepts an input set $x[n]$ and produces a single output sample $y[n]$ defined to be the polynomial $p[n]$ evaluated at $n=0$; i.e,

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

That is, for each input set of $2M+1$ samples of $x[n]$ (centered about $n=d$) the S-G filter produces the output, denoted by $y[n]$, which is equivalent to the single coefficient $a_0$ of the optimal LSE polynomial $p[n]$ associated with that particular window of samples $x[n]$. Incidentally as the window slides along the input data length, at every time $n=d$ a new output sample $y[d]$ is computed according to the window of samples $x[d-M],x[d-M+1],...,x[d-1],x[d],x[d+1],...x[d+M]$. Here this is a nancusal filter.

Now is the time to show that the coefficient $a_0$ is obtained as a linear combination of input signal values $x[n]$ in the window, and the filter output $y[n]$ is therefore a linear combination of input values $x[n]$. And this is the very definition of the (LTI) convolution via an FIR filter; the output at time $n$ is the linear combination of its input $x[n]$ and filter coefficients $h[n]$. But then, what are the filter coefficents for this S-G filter ? Lets see.

Again consider the computation of $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} $$

from which we can see that the single coefficent $a_0$ is given by the dot product of the first row of $H$ with the input data vector $x$; i.e.,

$$ a_0 = H(0,n) \cdot x = \sum H(0,k) x[k] = H(0,-n) \star x[n] $$

where in the last equality, we have interpreted the dot product sum as the convolution sum by considering the impulse response of the S-G filter to be $$\boxed{ h[n] = H(0,-n) }$$,

more specifically it's the impulse response of S-G filter of order $N$ with a window length of $2M+1$.

And the complete output $y[n]$ of the N-th order S-G filter witha window size of $2M+1$, for an input $x[n]$ of length $L$ is obtained by a single LTI convolution with the impulse response $h_N[n]$

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

COMMENT

The fact that polynomial coefficients $a_k$ depend on the input data, does not prevent the filter from being an LTI FIR. Because an impulse response $h[n]$ can be defined to represent the output $y[n]$ to be computed from a linear combinations of input samples. The linear combinations of input samples $x$ are inherently implied by the matrix product $a = H x$ that defines the optimal coefficients $a_k$ of $p[n]$, hence any linear combination of $a_k$ would also result in an FIR LTI filter $h[n]$ to represent the LSE polynomial fit procedure.

MATLAB / OCTVE CODE

The following simple MATLAB/OCTAVE can be used to compute those S-G filter impulse responses $h[n]$ (Note that its built in S-G designer may produce a different set of $h[n]$ as outlined by linked-pdf)

% 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]');

The output is:

enter image description here

Hope this clarifies the issue.

$\endgroup$
  • 2
    $\begingroup$ @Fat32 I think it is because it was a long back and forth list of comments, so to keep the board clean they typically move it "to chat". It's all still there, just not cluttering up the main page. That is why the system suggests moving it to chat when the back and forth gets to long. Not to fret, everyone still loves you. $\endgroup$ – Dan Boschen Sep 28 '18 at 19:11
  • 1
    $\begingroup$ @g6kxjv1ozn I'm comign to that point... please wait... $\endgroup$ – Fat32 Sep 28 '18 at 22:51
  • 2
    $\begingroup$ @Fat32 Great work! I read through it but will need to read through it and it is written very clearly, I will just need to follow through with pencil and paper step by step to totally see it as you now do. Thanks for putting this all here. $\endgroup$ – Dan Boschen Sep 29 '18 at 1:36
  • 4
    $\begingroup$ @DanBoschen: The polynomial is indeed not needed, it's just one way to describe S-G filters. They could also be described as filters minimizing the noise reduction ratio (i.e., the energy of the frequency response) subject to the frequency response being $1$ at DC, and additional flatness constraints at $\omega=0$. Cf. the section on S-G filters in Orfanidis' book linked to in my answer. $\endgroup$ – Matt L. Sep 29 '18 at 13:04
  • 2
    $\begingroup$ @DanBoschen Yes Dan the code shows just the filter impulse response h[n]... Let me repeat my self again; eventhough the polynomial coefficeints $a_k$ do certainly (linearly) depend on (local) input {x[n-M],...,x[n+M]}, the filter impulse response $h[n]$ doesnt. Our confusion stemmed from the fact that we thought the filter coefficients were a function of polynomial coefficients which would then make the filter as time varying and nonlinear, but as it turns out, that's not the case. $\endgroup$ – Fat32 Sep 30 '18 at 10:02

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.