Calculating Savitsky-Golay filter coefficients - normalization

Question

Using the matrix equations found on Wikipedia I wrote some Python code to compute the Savitsky-Golay coefficients for $m$ data-points, using a polynomial of degree $k$. Here it is:

from numpy import *
from numpy.linalg import *

k = 4       # Degree of polynomial to use for regression
m = 25      # Window size, that is, the number of points. Must be an odd number. 

z = range(-(m-1)/2, (m-1)/2+1)  # z is the distance from the central point. 

jacobian, b = [], []
for i in z:
  for j in range(k+1):
    b.append(i**j)
  jacobian.append(b)
  b=[]

J = matrix(jacobian)

C = (J.T*J).I*J.T

I do realize that this code is in no way optimized and it could've been done in a much cleaner way, however I'm not a (professional) programmer nor a mathematician and would be really happy if I could get it to work at all.

My problem is that the coefficients computed become rather small, using the parameters above ($k=4$, $m=25$) on the order of $\pm 10^{-6}$. If I "go one step back" and multiply with the determinant I get really huge numbers ($\pm 10^{22}$).

I tried finding the greatest common divisor of these but got the result 2048 which leaves me with huge numbers (and it can't be right, result of numerical errors due to small/huge numbers or something?). I would like to have the coefficients in the form they're usually found in tables, with a normalization constant for each derivative, and an integer coefficient for each data point. Since I'm interested in the point of maximum curvature in a data series I want to calculate its 3rd derivative, I've been unable to find a table with these coefficients.

Answer 1

The SciPy library for scientific computing in Python contains functions for Savitzky-Golay filtering in its scipy.signal module since version 0.14.0, specifically you'll most likely be interested in the scipy.signal.savgol_coeffs function if you're only interested in the coefficients, or the scipy.signal.savgol_filter function which provides a nice interface for filtering alone.

When implementing some DSP method in Python, it's usually best to look into SciPy or scikits-* libraries rather than reinvent the wheel.

Answer 2

I do not think that the method you are using is the one you want. The method you chose is good if you have very few measurements, but in your case you have too many and there are much easier methods to calculate curvature if that is the case.

Curvature is defined as $k(x) = f''(x)/|1+f'(x)^2|^{3/2}$. To calculate f'(x) and f''(x) there some easy formulas:

$f'(x) = ( f(x+h)-f(x-h) )/(2h)$ <------- Use this equation first

$f''(x) = ( f'(x+h)-f'(x-h) )/(2h)$, <------ Then use this equation with your calculated $f'$ values

where $f(x)$ is your measurement on position $x$, $h$ is the distance to another measurement(in your case a small integer value). Use these equations with different h and take a mean value the get a good approximation of $f'(x)$ and $f''(x)$. Read up on taylor series expansion if you want an explanation of these equations.

This should be sufficient to calculate the curvature on a given point.

By the way: This row "C = (J.TJ).IJ.T" is really not a good way to do what you were trying to do. This method is used to explain the theory, but not to be used in code. The C value is nice to have but it is really slow to calculate/ not always possible to calculate. This problem is very common(the method is called "least squares"). Because of this there are a lot of built in solutions to get what wiki calls C * y:

I don't know exactly how do this in python but in matlab C * y is calculated as J\y