savitzky.py
Go to the documentation of this file.
1 #!/usr/bin/env python
2 #
3 # Copyright 2017 Fraunhofer Institute for Manufacturing Engineering and Automation (IPA)
4 #
6 # you may not use this file except in compliance with the License.
7 # You may obtain a copy of the License at
8 #
10 #
11 # Unless required by applicable law or agreed to in writing, software
13 # WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 # See the License for the specific language governing permissions and
15 # limitations under the License.
16
17
18 import numpy as np
19 from math import factorial
20 import pylab
21
23  r"""Smooth (and optionally differentiate) data with a Savitzky-Golay filter.
24  The Savitzky-Golay filter removes high frequency noise from data.
25  It has the advantage of preserving the original shape and
26  features of the signal better than other types of filtering
27  approaches, such as moving averages techniques.
28  Parameters
29  ----------
30  y : array_like, shape (N,)
31  the values of the time history of the signal.
32  window_size : int
33  the length of the window. Must be an odd integer number.
34  order : int
35  the order of the polynomial used in the filtering.
36  Must be less then `window_size` - 1.
37  deriv: int
38  the order of the derivative to compute (default = 0 means only smoothing)
39  Returns
40  -------
41  ys : ndarray, shape (N)
42  the smoothed signal (or it's n-th derivative).
43  Notes
44  -----
45  The Savitzky-Golay is a type of low-pass filter, particularly
46  suited for smoothing noisy data. The main idea behind this
47  approach is to make for each point a least-square fit with a
48  polynomial of high order over a odd-sized window centered at
49  the point.
50  Examples
51  --------
52  t = np.linspace(-4, 4, 500)
53  y = np.exp( -t**2 ) + np.random.normal(0, 0.05, t.shape)
54  ysg = savitzky_golay(y, window_size=31, order=4)
55  import matplotlib.pyplot as plt
56  plt.plot(t, y, label='Noisy signal')
57  plt.plot(t, np.exp(-t**2), 'k', lw=1.5, label='Original signal')
58  plt.plot(t, ysg, 'r', label='Filtered signal')
59  plt.legend()
60  plt.show()
61  References
62  ----------
63  ..  A. Savitzky, M. J. E. Golay, Smoothing and Differentiation of
64  Data by Simplified Least Squares Procedures. Analytical
65  Chemistry, 1964, 36 (8), pp 1627-1639.
66  ..  Numerical Recipes 3rd Edition: The Art of Scientific Computing
67  W.H. Press, S.A. Teukolsky, W.T. Vetterling, B.P. Flannery
68  Cambridge University Press ISBN-13: 9780521880688
69  """
70
71  def __init__(self, window_size, order, deriv=0, rate=1):
72
73  self.window_size = window_size
74  self.order = order
75  self.deriv = deriv
76  self.rate = rate
77
78  try:
79  self.window_size = np.abs(np.int(self.window_size))
80  self.order = np.abs(np.int(self.order))
81  except ValueError as msg:
82  raise ValueError("ValueError: window_size and order have to be of type int - {}".format(msg))
83  if self.window_size % 2 != 1 or self.window_size < 1:
84  raise TypeError("window_size size must be a positive odd number")
85  if self.window_size < self.order + 2:
86  raise TypeError("window_size is too small for the polynomials order")
87
88  def filter(self, y):
89
90  self.order_range = list(range(self.order+1))
91  half_window = (self.window_size -1) // 2
92  # precompute coefficients
93  b = np.mat([[k**i for i in self.order_range] for k in range(-half_window, half_window+1)])
94  m = np.linalg.pinv(b).A[self.deriv] * self.rate**self.deriv * factorial(self.deriv)
95  # pad the signal at the extremes with
96  # values taken from the signal itself
97  firstvals = y - np.abs( y[1:half_window+1][::-1] - y )
98  lastvals = y[-1] + np.abs(y[-half_window-1:-1][::-1] - y[-1])
99  y = np.concatenate((firstvals, y, lastvals))
100  return np.convolve( m[::-1], y, mode='valid')
101
102 if __name__ == "__main__":
103
104  # create some sample twoD data
105  x = np.linspace(-3,3,100)
106  Z = np.exp( np.negative(x))
107
109  Zn = Z + np.random.normal( 0, 0.2, Z.shape )
110
111  sg = savitzky_golay(window_size=29, order=4)
112
113  l = sg.filter(Zn)
114
115  pylab.plot(l)
116
117  pylab.plot(Zn)
118
119  pylab.show()

cob_voltage_control
Author(s): Alexander Bubeck
autogenerated on Wed Apr 7 2021 02:11:57