Compare Savitzky-Golay filter f( x + xnoise ) ~ Savgol( f( x )) + ynoise for f(x) = sin( 2 pi freq x ), to try to elucidate Press et al., Numerical Recipes p. 772.

savgol-noise.py

""" Compare Savitzky-Golay filter f( x + xnoise )  ~  Savgol( f( x )) + ynoise
for f(x) = sin( 2 pi freq x )
to try to elucidate Press et al., Numerical Recipes p. 772:
    "for irregularly sampled data ... one can pretend that the data points *are* equally spaced
    ... if change in f across the full width ..." 

(Simpler: noise amplification of linear filters:
    |convolve( filter, x )| <= |filter| * |x|, || = rms
    |flat filter| = 1 / sqrt(N)
    |Savgol| ~ 2 * that

"""

from __future__ import division
import sys
import numpy as np

__version__ = "2013-07-09 jul denis"

def savgol_coef( n, d=4 ):
    """ NR pp. 766-772 """
    # n 9 d 4: d 4: [ 3 -13 7 31  42  31 7 -13 3]
    # noise amplification rms .65  > 1/3
    assert n % 2 == 1, n
    x = np.arange( - (n//2), n//2 + 1. )
    A = np.array([ x**k for k in range( d+1 )]) .T  # Ajk = basisk( xj )
    # print A
    return np.linalg.pinv( A )[0]

#...............................................................................
Ns = [9, 17, 33]
d = 4
nx = 100
nrand = 100
plot = 0
seed = 0

exec( "\n".join( sys.argv[1:] ))  # run this.py n= ...  from sh or ipython
np.set_printoptions( 2, threshold=100, edgeitems=10, suppress=True )
np.random.seed(seed)

title = """
Savitzky-Golay filter f( x + xnoise )  ~  Savgol( f( x )) + ynoise, NR p. 772:
    ratio ynoise.std / (xnoise.std / sqrt N)
    f(x) = sin( 2 pi freq x )\
"""
print title

x = np.arange( nx + 0.  )
xnoise = np.random.uniform( -.5, .5, size=(nrand, nx) )
xnoise_std = xnoise.std()  # 1 / sqrt 12  ~ .3

def f( x ):
    y = np.sin( 2*np.pi * freq * x )
    return y

#...............................................................................
    # err = ynoise - ynonoise for sin( 2 pi freq x ) --
for freq in np.arange( 0, .51, .1 ):
    print "\nfreq %.2g --" % freq
    for N in Ns:
        coefs = savgol_coef( N, d )
        ynonoise = np.convolve( f(x), coefs, mode="valid" )
        ynoise = np.array([ np.convolve( f( x + noise ), coefs, mode="valid" )
                        for noise in xnoise ])
        err = (ynoise - ynonoise) / (xnoise_std / np.sqrt(N))
        errstd, errmax = err.std(), err.__abs__().max()
        print "N %2d: err / (sigma / sqrt N)  av %.2g  max %.2g " % ( N, errstd, errmax )

        # freq .1 .. .5 -> ~ 1 1.6 2.4 3 4.8