00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018 import numpy as np
00019 from math import factorial
00020 import pylab
00021
00022 class savitzky_golay():
00023 r"""Smooth (and optionally differentiate) data with a Savitzky-Golay filter.
00024 The Savitzky-Golay filter removes high frequency noise from data.
00025 It has the advantage of preserving the original shape and
00026 features of the signal better than other types of filtering
00027 approaches, such as moving averages techniques.
00028 Parameters
00029 ----------
00030 y : array_like, shape (N,)
00031 the values of the time history of the signal.
00032 window_size : int
00033 the length of the window. Must be an odd integer number.
00034 order : int
00035 the order of the polynomial used in the filtering.
00036 Must be less then `window_size` - 1.
00037 deriv: int
00038 the order of the derivative to compute (default = 0 means only smoothing)
00039 Returns
00040 -------
00041 ys : ndarray, shape (N)
00042 the smoothed signal (or it's n-th derivative).
00043 Notes
00044 -----
00045 The Savitzky-Golay is a type of low-pass filter, particularly
00046 suited for smoothing noisy data. The main idea behind this
00047 approach is to make for each point a least-square fit with a
00048 polynomial of high order over a odd-sized window centered at
00049 the point.
00050 Examples
00051 --------
00052 t = np.linspace(-4, 4, 500)
00053 y = np.exp( -t**2 ) + np.random.normal(0, 0.05, t.shape)
00054 ysg = savitzky_golay(y, window_size=31, order=4)
00055 import matplotlib.pyplot as plt
00056 plt.plot(t, y, label='Noisy signal')
00057 plt.plot(t, np.exp(-t**2), 'k', lw=1.5, label='Original signal')
00058 plt.plot(t, ysg, 'r', label='Filtered signal')
00059 plt.legend()
00060 plt.show()
00061 References
00062 ----------
00063 .. [1] A. Savitzky, M. J. E. Golay, Smoothing and Differentiation of
00064 Data by Simplified Least Squares Procedures. Analytical
00065 Chemistry, 1964, 36 (8), pp 1627-1639.
00066 .. [2] Numerical Recipes 3rd Edition: The Art of Scientific Computing
00067 W.H. Press, S.A. Teukolsky, W.T. Vetterling, B.P. Flannery
00068 Cambridge University Press ISBN-13: 9780521880688
00069 """
00070
00071 def __init__(self, window_size, order, deriv=0, rate=1):
00072
00073 self.window_size = window_size
00074 self.order = order
00075 self.deriv = deriv
00076 self.rate = rate
00077
00078 try:
00079 self.window_size = np.abs(np.int(self.window_size))
00080 self.order = np.abs(np.int(self.order))
00081 except ValueError, msg:
00082 raise ValueError("window_size and order have to be of type int")
00083 if self.window_size % 2 != 1 or self.window_size < 1:
00084 raise TypeError("window_size size must be a positive odd number")
00085 if self.window_size < self.order + 2:
00086 raise TypeError("window_size is too small for the polynomials order")
00087
00088 def filter(self, y):
00089
00090 self.order_range = range(self.order+1)
00091 half_window = (self.window_size -1) // 2
00092 # precompute coefficients
00093 b = np.mat([[k**i for i in self.order_range] for k in range(-half_window, half_window+1)])
00094 m = np.linalg.pinv(b).A[self.deriv] * self.rate**self.deriv * factorial(self.deriv)
00095 # pad the signal at the extremes with
00096 # values taken from the signal itself
00097 firstvals = y[0] - np.abs( y[1:half_window+1][::-1] - y[0] )
00098 lastvals = y[-1] + np.abs(y[-half_window-1:-1][::-1] - y[-1])
00099 y = np.concatenate((firstvals, y, lastvals))
00100 return np.convolve( m[::-1], y, mode='valid')
00101
00102 if __name__ == "__main__":
00103
00104 # create some sample twoD data
00105 x = np.linspace(-3,3,100)
00106 Z = np.exp( -(x))
00107
00108 # add noise
00109 Zn = Z + np.random.normal( 0, 0.2, Z.shape )
00110
00111 sg = savitzky_golay(window_size=29, order=4)
00112
00113 l = sg.filter(Zn)
00114
00115 pylab.plot(l)
00116
00117 pylab.plot(Zn)
00118
00119 pylab.show()