1import numpy as np 2from scipy.linalg import lstsq 3from scipy._lib._util import float_factorial 4from scipy.ndimage import convolve1d 5from ._arraytools import axis_slice 6 7 8def savgol_coeffs(window_length, polyorder, deriv=0, delta=1.0, pos=None, 9 use="conv"): 10 """Compute the coefficients for a 1-D Savitzky-Golay FIR filter. 11 12 Parameters 13 ---------- 14 window_length : int 15 The length of the filter window (i.e., the number of coefficients). 16 `window_length` must be an odd positive integer. 17 polyorder : int 18 The order of the polynomial used to fit the samples. 19 `polyorder` must be less than `window_length`. 20 deriv : int, optional 21 The order of the derivative to compute. This must be a 22 nonnegative integer. The default is 0, which means to filter 23 the data without differentiating. 24 delta : float, optional 25 The spacing of the samples to which the filter will be applied. 26 This is only used if deriv > 0. 27 pos : int or None, optional 28 If pos is not None, it specifies evaluation position within the 29 window. The default is the middle of the window. 30 use : str, optional 31 Either 'conv' or 'dot'. This argument chooses the order of the 32 coefficients. The default is 'conv', which means that the 33 coefficients are ordered to be used in a convolution. With 34 use='dot', the order is reversed, so the filter is applied by 35 dotting the coefficients with the data set. 36 37 Returns 38 ------- 39 coeffs : 1-D ndarray 40 The filter coefficients. 41 42 References 43 ---------- 44 A. Savitzky, M. J. E. Golay, Smoothing and Differentiation of Data by 45 Simplified Least Squares Procedures. Analytical Chemistry, 1964, 36 (8), 46 pp 1627-1639. 47 48 See Also 49 -------- 50 savgol_filter 51 52 Notes 53 ----- 54 55 .. versionadded:: 0.14.0 56 57 Examples 58 -------- 59 >>> from scipy.signal import savgol_coeffs 60 >>> savgol_coeffs(5, 2) 61 array([-0.08571429, 0.34285714, 0.48571429, 0.34285714, -0.08571429]) 62 >>> savgol_coeffs(5, 2, deriv=1) 63 array([ 2.00000000e-01, 1.00000000e-01, 2.07548111e-16, -1.00000000e-01, 64 -2.00000000e-01]) 65 66 Note that use='dot' simply reverses the coefficients. 67 68 >>> savgol_coeffs(5, 2, pos=3) 69 array([ 0.25714286, 0.37142857, 0.34285714, 0.17142857, -0.14285714]) 70 >>> savgol_coeffs(5, 2, pos=3, use='dot') 71 array([-0.14285714, 0.17142857, 0.34285714, 0.37142857, 0.25714286]) 72 73 `x` contains data from the parabola x = t**2, sampled at 74 t = -1, 0, 1, 2, 3. `c` holds the coefficients that will compute the 75 derivative at the last position. When dotted with `x` the result should 76 be 6. 77 78 >>> x = np.array([1, 0, 1, 4, 9]) 79 >>> c = savgol_coeffs(5, 2, pos=4, deriv=1, use='dot') 80 >>> c.dot(x) 81 6.0 82 """ 83 84 # An alternative method for finding the coefficients when deriv=0 is 85 # t = np.arange(window_length) 86 # unit = (t == pos).astype(int) 87 # coeffs = np.polyval(np.polyfit(t, unit, polyorder), t) 88 # The method implemented here is faster. 89 90 # To recreate the table of sample coefficients shown in the chapter on 91 # the Savitzy-Golay filter in the Numerical Recipes book, use 92 # window_length = nL + nR + 1 93 # pos = nL + 1 94 # c = savgol_coeffs(window_length, M, pos=pos, use='dot') 95 96 if polyorder >= window_length: 97 raise ValueError("polyorder must be less than window_length.") 98 99 halflen, rem = divmod(window_length, 2) 100 101 if rem == 0: 102 raise ValueError("window_length must be odd.") 103 104 if pos is None: 105 pos = halflen 106 107 if not (0 <= pos < window_length): 108 raise ValueError("pos must be nonnegative and less than " 109 "window_length.") 110 111 if use not in ['conv', 'dot']: 112 raise ValueError("`use` must be 'conv' or 'dot'") 113 114 if deriv > polyorder: 115 coeffs = np.zeros(window_length) 116 return coeffs 117 118 # Form the design matrix A. The columns of A are powers of the integers 119 # from -pos to window_length - pos - 1. The powers (i.e., rows) range 120 # from 0 to polyorder. (That is, A is a vandermonde matrix, but not 121 # necessarily square.) 122 x = np.arange(-pos, window_length - pos, dtype=float) 123 if use == "conv": 124 # Reverse so that result can be used in a convolution. 125 x = x[::-1] 126 127 order = np.arange(polyorder + 1).reshape(-1, 1) 128 A = x ** order 129 130 # y determines which order derivative is returned. 131 y = np.zeros(polyorder + 1) 132 # The coefficient assigned to y[deriv] scales the result to take into 133 # account the order of the derivative and the sample spacing. 134 y[deriv] = float_factorial(deriv) / (delta ** deriv) 135 136 # Find the least-squares solution of A*c = y 137 coeffs, _, _, _ = lstsq(A, y) 138 139 return coeffs 140 141 142def _polyder(p, m): 143 """Differentiate polynomials represented with coefficients. 144 145 p must be a 1-D or 2-D array. In the 2-D case, each column gives 146 the coefficients of a polynomial; the first row holds the coefficients 147 associated with the highest power. m must be a nonnegative integer. 148 (numpy.polyder doesn't handle the 2-D case.) 149 """ 150 151 if m == 0: 152 result = p 153 else: 154 n = len(p) 155 if n <= m: 156 result = np.zeros_like(p[:1, ...]) 157 else: 158 dp = p[:-m].copy() 159 for k in range(m): 160 rng = np.arange(n - k - 1, m - k - 1, -1) 161 dp *= rng.reshape((n - m,) + (1,) * (p.ndim - 1)) 162 result = dp 163 return result 164 165 166def _fit_edge(x, window_start, window_stop, interp_start, interp_stop, 167 axis, polyorder, deriv, delta, y): 168 """ 169 Given an N-d array `x` and the specification of a slice of `x` from 170 `window_start` to `window_stop` along `axis`, create an interpolating 171 polynomial of each 1-D slice, and evaluate that polynomial in the slice 172 from `interp_start` to `interp_stop`. Put the result into the 173 corresponding slice of `y`. 174 """ 175 176 # Get the edge into a (window_length, -1) array. 177 x_edge = axis_slice(x, start=window_start, stop=window_stop, axis=axis) 178 if axis == 0 or axis == -x.ndim: 179 xx_edge = x_edge 180 swapped = False 181 else: 182 xx_edge = x_edge.swapaxes(axis, 0) 183 swapped = True 184 xx_edge = xx_edge.reshape(xx_edge.shape[0], -1) 185 186 # Fit the edges. poly_coeffs has shape (polyorder + 1, -1), 187 # where '-1' is the same as in xx_edge. 188 poly_coeffs = np.polyfit(np.arange(0, window_stop - window_start), 189 xx_edge, polyorder) 190 191 if deriv > 0: 192 poly_coeffs = _polyder(poly_coeffs, deriv) 193 194 # Compute the interpolated values for the edge. 195 i = np.arange(interp_start - window_start, interp_stop - window_start) 196 values = np.polyval(poly_coeffs, i.reshape(-1, 1)) / (delta ** deriv) 197 198 # Now put the values into the appropriate slice of y. 199 # First reshape values to match y. 200 shp = list(y.shape) 201 shp[0], shp[axis] = shp[axis], shp[0] 202 values = values.reshape(interp_stop - interp_start, *shp[1:]) 203 if swapped: 204 values = values.swapaxes(0, axis) 205 # Get a view of the data to be replaced by values. 206 y_edge = axis_slice(y, start=interp_start, stop=interp_stop, axis=axis) 207 y_edge[...] = values 208 209 210def _fit_edges_polyfit(x, window_length, polyorder, deriv, delta, axis, y): 211 """ 212 Use polynomial interpolation of x at the low and high ends of the axis 213 to fill in the halflen values in y. 214 215 This function just calls _fit_edge twice, once for each end of the axis. 216 """ 217 halflen = window_length // 2 218 _fit_edge(x, 0, window_length, 0, halflen, axis, 219 polyorder, deriv, delta, y) 220 n = x.shape[axis] 221 _fit_edge(x, n - window_length, n, n - halflen, n, axis, 222 polyorder, deriv, delta, y) 223 224 225def savgol_filter(x, window_length, polyorder, deriv=0, delta=1.0, 226 axis=-1, mode='interp', cval=0.0): 227 """ Apply a Savitzky-Golay filter to an array. 228 229 This is a 1-D filter. If `x` has dimension greater than 1, `axis` 230 determines the axis along which the filter is applied. 231 232 Parameters 233 ---------- 234 x : array_like 235 The data to be filtered. If `x` is not a single or double precision 236 floating point array, it will be converted to type ``numpy.float64`` 237 before filtering. 238 window_length : int 239 The length of the filter window (i.e., the number of coefficients). 240 `window_length` must be a positive odd integer. If `mode` is 'interp', 241 `window_length` must be less than or equal to the size of `x`. 242 polyorder : int 243 The order of the polynomial used to fit the samples. 244 `polyorder` must be less than `window_length`. 245 deriv : int, optional 246 The order of the derivative to compute. This must be a 247 nonnegative integer. The default is 0, which means to filter 248 the data without differentiating. 249 delta : float, optional 250 The spacing of the samples to which the filter will be applied. 251 This is only used if deriv > 0. Default is 1.0. 252 axis : int, optional 253 The axis of the array `x` along which the filter is to be applied. 254 Default is -1. 255 mode : str, optional 256 Must be 'mirror', 'constant', 'nearest', 'wrap' or 'interp'. This 257 determines the type of extension to use for the padded signal to 258 which the filter is applied. When `mode` is 'constant', the padding 259 value is given by `cval`. See the Notes for more details on 'mirror', 260 'constant', 'wrap', and 'nearest'. 261 When the 'interp' mode is selected (the default), no extension 262 is used. Instead, a degree `polyorder` polynomial is fit to the 263 last `window_length` values of the edges, and this polynomial is 264 used to evaluate the last `window_length // 2` output values. 265 cval : scalar, optional 266 Value to fill past the edges of the input if `mode` is 'constant'. 267 Default is 0.0. 268 269 Returns 270 ------- 271 y : ndarray, same shape as `x` 272 The filtered data. 273 274 See Also 275 -------- 276 savgol_coeffs 277 278 Notes 279 ----- 280 Details on the `mode` options: 281 282 'mirror': 283 Repeats the values at the edges in reverse order. The value 284 closest to the edge is not included. 285 'nearest': 286 The extension contains the nearest input value. 287 'constant': 288 The extension contains the value given by the `cval` argument. 289 'wrap': 290 The extension contains the values from the other end of the array. 291 292 For example, if the input is [1, 2, 3, 4, 5, 6, 7, 8], and 293 `window_length` is 7, the following shows the extended data for 294 the various `mode` options (assuming `cval` is 0):: 295 296 mode | Ext | Input | Ext 297 -----------+---------+------------------------+--------- 298 'mirror' | 4 3 2 | 1 2 3 4 5 6 7 8 | 7 6 5 299 'nearest' | 1 1 1 | 1 2 3 4 5 6 7 8 | 8 8 8 300 'constant' | 0 0 0 | 1 2 3 4 5 6 7 8 | 0 0 0 301 'wrap' | 6 7 8 | 1 2 3 4 5 6 7 8 | 1 2 3 302 303 .. versionadded:: 0.14.0 304 305 Examples 306 -------- 307 >>> from scipy.signal import savgol_filter 308 >>> np.set_printoptions(precision=2) # For compact display. 309 >>> x = np.array([2, 2, 5, 2, 1, 0, 1, 4, 9]) 310 311 Filter with a window length of 5 and a degree 2 polynomial. Use 312 the defaults for all other parameters. 313 314 >>> savgol_filter(x, 5, 2) 315 array([1.66, 3.17, 3.54, 2.86, 0.66, 0.17, 1. , 4. , 9. ]) 316 317 Note that the last five values in x are samples of a parabola, so 318 when mode='interp' (the default) is used with polyorder=2, the last 319 three values are unchanged. Compare that to, for example, 320 `mode='nearest'`: 321 322 >>> savgol_filter(x, 5, 2, mode='nearest') 323 array([1.74, 3.03, 3.54, 2.86, 0.66, 0.17, 1. , 4.6 , 7.97]) 324 325 """ 326 if mode not in ["mirror", "constant", "nearest", "interp", "wrap"]: 327 raise ValueError("mode must be 'mirror', 'constant', 'nearest' " 328 "'wrap' or 'interp'.") 329 330 x = np.asarray(x) 331 # Ensure that x is either single or double precision floating point. 332 if x.dtype != np.float64 and x.dtype != np.float32: 333 x = x.astype(np.float64) 334 335 coeffs = savgol_coeffs(window_length, polyorder, deriv=deriv, delta=delta) 336 337 if mode == "interp": 338 if window_length > x.size: 339 raise ValueError("If mode is 'interp', window_length must be less " 340 "than or equal to the size of x.") 341 342 # Do not pad. Instead, for the elements within `window_length // 2` 343 # of the ends of the sequence, use the polynomial that is fitted to 344 # the last `window_length` elements. 345 y = convolve1d(x, coeffs, axis=axis, mode="constant") 346 _fit_edges_polyfit(x, window_length, polyorder, deriv, delta, axis, y) 347 else: 348 # Any mode other than 'interp' is passed on to ndimage.convolve1d. 349 y = convolve1d(x, coeffs, axis=axis, mode=mode, cval=cval) 350 351 return y 352