Skip to content

Instantly share code, notes, and snippets.

@denis-bz
Last active December 19, 2015 12:49
Show Gist options
  • Save denis-bz/5957279 to your computer and use it in GitHub Desktop.
Save denis-bz/5957279 to your computer and use it in GitHub Desktop.
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.
""" 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
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment