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:
Hope this clarifies the issue.
