Savitzky-Golay filters & Julia

By: Picaud Vincent

Re-posted from: https://pixorblog.wordpress.com/2016/07/13/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)
    
    padded_signal = 
        [signal[1]*ones(halfWindow);
         signal;
         signal[end]*ones(halfWindow)]

    filter_cross_signal = conv(filter[end:-1:1], padded_signal)

    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

s=readdlm("signal.txt")[:,1]

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:

smoothed.png smoothed_d1.png smoothed_d2.png

Final word

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