savitzky.py
Go to the documentation of this file.
```00001 #!/usr/bin/env python
00002 #
00003 # Copyright 2017 Fraunhofer Institute for Manufacturing Engineering and Automation (IPA)
00004 #
00006 # you may not use this file except in compliance with the License.
00007 # You may obtain a copy of the License at
00008 #
00010 #
00011 # Unless required by applicable law or agreed to in writing, software
00013 # WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
00014 # See the License for the specific language governing permissions and
00015 # limitations under the License.
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
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()
```

cob_voltage_control
Author(s): Alexander Bubeck
autogenerated on Sat Jun 8 2019 21:02:33