# Savitzky-Golay filters & Julia

I always have found that presentations of the Savitzky-Golay filters were over tricky. I present here a simple derivation of these formula and a possible implementation in Julia.

## Derivation of the formula

Given a polynomial of degree $d$ with values

$p(x_i)=\sum\limits_{j=0}^d c_j x_i^j$

defined for a window of size $2n+1$, $\{x_i,\ i\in\llbracket -n,n \rrbracket \}$

We want to find the value of its k-order derivative $p^{(k)}(x_0)$ in the middle of the window assuming that the $c_j$ are founded solving a least-squares problem:

$\min\limits_{\mathbf{c}} \frac{1}{2} \| \mathbf{V} \mathbf{c} - \mathbf{y} \|_2^2$

where $\{y_i, i \in\llbracket -n,n \rrbracket \}$ is our signal values and $\mathbf{V}$ is the Vandermonde matrix:

$\mathbf{V}= \left( \begin{array}{c|c|c} \vdots & \vdots & \vdots \\ 1 & x_i^{(j-1)} & x_i^d \\ \vdots & \vdots & \vdots \end{array} \right)$

using the normal equation

$\mathbf{c}=(\mathbf{V}^t.\mathbf{V})^{-1}.\mathbf{V}^t.\mathbf{y}$

and a QR decomposition, $\mathbf{V}=\mathbf{Q}.\mathbf{R}$ we get

$\mathbf{c}=\mathbf{R}^{-1}.\mathbf{Q}^t.\mathbf{y}$

now we can express all the polynomial values $p(x_i)$ in a vector $\mathbf{p}=\mathbf{V}.\mathbf{c}$. Lets rewrite this in matrix form:

$\underbrace{\left( \begin{array}{c} p(x_{-n}) \\ \vdots \\ p(x_{0}) \\ \vdots \\ p(x_{+n}) \end{array} \right)}\limits_{\mathbf{p}}=\underbrace{ \left( \begin{array}{c|c|c} \vdots & \vdots & \vdots \\ 1 & x_i^{(j-1)} & x_i^d \\ \vdots & \vdots & \vdots \end{array} \right)}\limits_{\mathbf{V}}.\underbrace{\left( \begin{array}{c} c_0 \\ \vdots \\ c_n \end{array} \right)}\limits_{\mathbf{c}}$

Now the “big trick” is to write the Taylor series and to remember that this formula is exact for polynomial functions:

$\forall i,\ P(x_i) = \sum\limits_{j=0}^d \frac{x_i^j}{j!} P^{(j)}(x_0)$

Lets rewrite this in matrix form:

$\underbrace{ \left( \begin{array}{c} p(x_{-n}) \\ \vdots \\ p(x_{0}) \\ \vdots \\ p(x_{n}) \\ \end{array} \right) }_{\mathbf{p}} = \underbrace{ \left( \begin{array}{c|c|c} \vdots & \vdots & \vdots \\ 1 & \frac{x_i^{(j-1)}}{(j-1)!} & \frac{x_i^{d}}{d!} \\ \vdots & \vdots & \vdots \end{array} \right) }_{\mathbf{T}} \underbrace{ \left( \begin{array}{c} P^{(0)}(x_0) \\ \vdots \\ P^{(k)}(x_0) \\ \vdots \\ P^{(d)}(x_0) \\ \end{array} \right) }_{\mathbf{p^\delta}}$

With a good eye we see that $\mathbf{V}=\mathbf{T}.\mathbf{D}$ where $\mathbf{D}$ is a diagonal matrix:

$\underbrace{ \left( \begin{array}{c|c|c} \vdots & \vdots & \vdots \\ 1 & x_i^{(j-1)} & x_i^d \\ \vdots & \vdots & \vdots \end{array} \right)}\limits_{\mathbf{V}} = \underbrace{ \left( \begin{array}{c|c|c} \vdots & \vdots & \vdots \\ 1 & \frac{x_i^{(j-1)}}{(j-1)!} & \frac{x_i^{d}}{d!} \\ \vdots & \vdots & \vdots \end{array} \right) }_{\mathbf{T}}.\underbrace{\left( \begin{array}{ccc} 1 & & \\ & (j-1)! & \\ & & d! \end{array} \right)}\limits_{\mathbf{D}}$

That’s all, we only have to group pieces:

$\mathbf{V}.\mathbf{c}=\mathbf{P}=\mathbf{T}.\mathbf{p^\delta}=\mathbf{V}.\mathbf{D}^{-1}.\mathbf{p^\delta}$

With the QR decomposition $\mathbf{V}=\mathbf{Q}.\mathbf{R}$ and $\mathbf{c}=\mathbf{R}^{-1}.\mathbf{Q}^t.\mathbf{y}$ we get:

$\mathbf{Q}.\mathbf{Q}^t.\mathbf{y}=\mathbf{Q}.\mathbf{R}.\mathbf{D}^{-1}.\mathbf{p^\delta}$

using the fact that $\mathbf{Q}^t.\mathbf{Q}=\mathbf{I}$ we get:

$\mathbf{Q}^t.\mathbf{y}=\mathbf{R}.\mathbf{D}^{-1}.\mathbf{p^\delta}$

hence we have:

$\boxed{ \mathbf{p^\delta} = \mathbf{D}.\mathbf{R}^{-1}.\mathbf{Q}^t.\mathbf{y} }$

which is our final formula.

## Symbolic computation to check that it works

We can use mathematica to do a symbolic computation using $\mathbf{p^\delta} = \mathbf{D}.\mathbf{R}^{-1}.\mathbf{Q}^t.\mathbf{y}$.

For a window width of $2n+1=7$ points and a polynomial of degree $d=2$ we get:

n = 3; d = 2;
V = Table[If[j != 0, i^j, 1], {i, -n, n}, {j, 0, d}];
{Qt, R} = QRDecomposition[V];
DD = DiagonalMatrix[Table[Factorial[i], {i, 0, d}]];
DD.Inverse[R].Qt // TeXForm


$\left( \begin{array}{ccccccc} -\frac{2}{21} & \frac{1}{7} & \frac{2}{7} & \frac{1}{3} & \frac{2}{7} & \frac{1}{7} & -\frac{2}{21} \\ -\frac{3}{28} & -\frac{1}{14} & -\frac{1}{28} & 0 & \frac{1}{28} & \frac{1}{14} & \frac{3}{28} \\ \frac{5}{42} & 0 & -\frac{1}{14} & -\frac{2}{21} & -\frac{1}{14} & 0 & \frac{5}{42} \end{array} \right)$

The first row is the smoothing $P(0)$ filter, the second one the smoothed first order derivative $P'(0)$ filter and the last one the smoothed second order derivative $P''(0)$ filter.

## A Julia implementation

Here I present a direct implementation in julia.

We first initialize a Vandermonde matrix:

function vandermonde(halfWindow::Int, polyDeg::Int,T::Type=Float64)

@assert halfWindow>=0
@assert polyDeg>=0

x=T[i for i in -halfWindow:halfWindow]

n = polyDeg+1
m = length(x)

V = Array{T}(m, n)

for i = 1:m
V[i,1] = T(1)
end
for j = 2:n
for i = 1:m
V[i,j] = x[i] * V[i,j-1]
end
end

return V
end


We then compute the filter coefficients with the presented formula.

Attention we return the transposed matrix because in Julia it is more convenient to use SG[:,1] which is a julia vector. SG[1,:] would be a (1,n) julia matrix.

function SG(halfWindow::Int, polyDeg::Int,T::Type=Float64)

@assert 2*halfWindow>polyDeg

V=vandermonde(halfWindow,polyDeg,T)
Q,R=qr(V)
SG=R\Q'

for i in 1:size(SG,1)
SG[i,:]*=factorial(i-1)
end

# CAVEAT: returns the transposed matrix

return SG'
end


The final step to use the filter is to provide function to do the cross-correlation.

I do not want to talk too much about this subroutine because in a future post I will show how to efficiently compute such kind of convolution. Here we use a FFT, but with a short filter it is much more efficient (around 10 times faster) to use a direct computation. I will show how to implement

$\gamma[k]=\sum\limits_i\alpha[i]\beta[k+\lambda i],\ with\ \lambda\in\mathbb{Z}^*$

which can be used to compute discrete and stationary wavelet transform for instance.

One last thing, here we use constant padding to manage the boundaries.

function apply_filter{T}(filter::StridedVector{T},signal::StridedVector{T})

@assert isodd(length(filter))

halfWindow = round(Int,(length(filter)-1)/2)

[signal[1]*ones(halfWindow);
signal;
signal[end]*ones(halfWindow)]

return filter_cross_signal[2*halfWindow+1:end-2*halfWindow]
end


Finally I have included a small usage example. To see well the effect of Savitzky-Golay filters, I have over smoothed with a large window width $2.n+1$, $n=20$

using Winston

sg=SG(20,3) # halt-window, polynomal degree

#________________

smoothed=apply_filter(sg[:,1],s)

plot(s,"r")
oplot(smoothed)
title("Smoothed")
savefig("smoothed.png")

#________________

smoothed_d1=apply_filter(sg[:,2],s)

plot(smoothed_d1)
title("Smoothed derivative")
savefig("smoothed_d1.png")

#________________

smoothed_d2=apply_filter(sg[:,3],s)

plot(smoothed_d2)
title("Smoothed 2-derivative")
savefig("smoothed_d2.png")


Here is the resulting plots:

## Final word

To reproduce these figures you can find the complete code on github.