Savitzky-Golay smoothing filter for not equally spaced data
Because of the way the Savitzky-Golay filter is derived (i.e as local least-squares polynomial fits), there's a natural generalization to nonuniform sampling--it's just much more computationally expensive.

Savitzky-Golay Filters in General

For the standard filter, the idea is to fit a polynomial to a local set of samples [using least squares], then replace the center sample with the value of the polynomial at the center index (i.e. at 0). That means the standard SG filter coefficients can be generated by inverting a Vandermonde matrix of sample indicies. For example, to generate a local parabolic fit across five samples $y_0\dots y_4$ (with local indicies -2,-1,0,1,2), the system of design equations $Ac = y$ would be as follows:

$$
\begin{bmatrix}-2^0 &amp; -2^1 &amp; -2^2 \\ -1^0 &amp; -1^1 &amp; -1^2 \\ 0^0 &amp; 0^1 &amp; 0^2 \\ 1^0 &amp; 1^1 &amp; 1^2 \\ 2^0 &amp; 2^1 &amp; 2^2 \end{bmatrix} \begin{bmatrix} c_0 \\ c_1 \\ c_2 \end{bmatrix} = \begin{bmatrix} y_0 \\ y_1 \\ y_2 \\ y_3 \\ y_4 \end{bmatrix}.
$$

In the above, the $c_0 \dots  c_2$ are the unknown coefficients of the least squares polynomial $c_0 + c_1x + c_2x^2$. Since the value of the polynomial at $x = 0$ is just $c_0$, computing the pseudoinverse of the design matrix (i.e. $c = (A^TA)^{-1}A^T y\space $) will yield the SG filter coefficients in the top row. In this case, they would be

$$
\begin{bmatrix}c_0 \\ c_1 \\ c_2 \end{bmatrix} =
\begin{bmatrix}
-3 &amp; 12 &amp; 17 &amp; 12 &amp; -3 \\
-7 &amp; -4 &amp; 0 &amp; 4 &amp; 7 \\
 5 &amp; -3 &amp; -5 &amp; -3 &amp; 5 \\
\end{bmatrix}
\begin{bmatrix} y_0 \\ y_1 \\ y_2 \\ y_3 \\ y_4 \end{bmatrix}.
$$

Note that since the derivative of $c_0 + c_1x + c_2x^2$ is $c_1 + 2c_2x$, the second row of the matrix (which evaluates $c_1$) will be a smoothed derivative filter. The same argument applies for successive rows--they give smoothed higher-order derivatives. Note that I scaled the matrix by 35 so the first row would match the smoothing coefficients given on Wikipedia (above). The derivative filters each differ by other scaling factors.

Nonuniform Sampling

When the samples are evenly spaced, the filter coefficients are translation-invariant, so the result is just an FIR filter. For nonuniform samples, the coefficients will differ based on the local sample spacing, so the design matrix will need to be constructed and inverted at each sample. If the nonuniform sample times are $x_n$, and we construct local coordinates $t_n$ with each center sample time fixed at $0$, i.e.

$$
\begin{align}
t_{-2} &amp; = x_{-2} - x_0 \\
t_{-1} &amp; = x_{-1} - x_0 \\ 
t_0 &amp; = x_0 - x_0 \\ 
t_1 &amp; = x_1 - x_0 \\
t_2 &amp; = x_2 - x_0
\end{align}
$$

then each design matrix will be of the following form:

$$
A = 
\begin{bmatrix}
 t_{-2}^0 &amp;  t_{-2}^1 &amp; t_{-2}^2 \\
 t_{-1}^0 &amp;  t_{-1}^1 &amp; t_{-1}^2 \\
 t_0^0 &amp;  t_0^1 &amp; t_0^2 \\
 t_1^0 &amp;  t_1^1 &amp; t_1^2 \\
 t_2^0 &amp;  t_2^1 &amp; t_2^2
\end{bmatrix} =
\begin{bmatrix}
 1 &amp;  t_{-2} &amp; t_{-2}^2 \\
 1 &amp;  t_{-1} &amp; t_{-1}^2 \\
 1 &amp;  0 &amp; 0 \\
 1 &amp;  t_1 &amp; t_1^2 \\
 1 &amp;  t_2 &amp; t_2^2
\end{bmatrix}.
$$

The first row of the pseudoinverse of $A$ dotted with the local sample values will yield $c_0$, the smoothed value at that sample.