import numpy as np
from scipy.linalg import lstsq
from scipy._lib._util import float_factorial
from scipy.ndimage import convolve1d
from ._arraytools import axis_slice
def savgol_coeffs(window_length, polyorder, deriv=0, delta=1.0, pos=None,
use="conv"):
"""Compute the coefficients for a 1-D Savitzky-Golay FIR filter.
Parameters
----------
window_length : int
The length of the filter window (i.e., the number of coefficients).
polyorder : int
The order of the polynomial used to fit the samples.
`polyorder` must be less than `window_length`.
deriv : int, optional
The order of the derivative to compute. This must be a
nonnegative integer. The default is 0, which means to filter
the data without differentiating.
delta : float, optional
The spacing of the samples to which the filter will be applied.
This is only used if deriv > 0.
pos : int or None, optional
If pos is not None, it specifies evaluation position within the
window. The default is the middle of the window.
use : str, optional
Either 'conv' or 'dot'. This argument chooses the order of the
coefficients. The default is 'conv', which means that the
coefficients are ordered to be used in a convolution. With
use='dot', the order is reversed, so the filter is applied by
dotting the coefficients with the data set.
Returns
-------
coeffs : 1-D ndarray
The filter coefficients.
See Also
--------
savgol_filter
Notes
-----
.. versionadded:: 0.14.0
References
----------
A. Savitzky, M. J. E. Golay, Smoothing and Differentiation of Data by
Simplified Least Squares Procedures. Analytical Chemistry, 1964, 36 (8),
pp 1627-1639.
Jianwen Luo, Kui Ying, and Jing Bai. 2005. Savitzky-Golay smoothing and
differentiation filter for even number data. Signal Process.
85, 7 (July 2005), 1429-1434.
Examples
--------
>>> import numpy as np
>>> from scipy.signal import savgol_coeffs
>>> savgol_coeffs(5, 2)
array([-0.08571429, 0.34285714, 0.48571429, 0.34285714, -0.08571429])
>>> savgol_coeffs(5, 2, deriv=1)
array([ 2.00000000e-01, 1.00000000e-01, 2.07548111e-16, -1.00000000e-01,
-2.00000000e-01])
Note that use='dot' simply reverses the coefficients.
>>> savgol_coeffs(5, 2, pos=3)
array([ 0.25714286, 0.37142857, 0.34285714, 0.17142857, -0.14285714])
>>> savgol_coeffs(5, 2, pos=3, use='dot')
array([-0.14285714, 0.17142857, 0.34285714, 0.37142857, 0.25714286])
>>> savgol_coeffs(4, 2, pos=3, deriv=1, use='dot')
array([0.45, -0.85, -0.65, 1.05])
`x` contains data from the parabola x = t**2, sampled at
t = -1, 0, 1, 2, 3. `c` holds the coefficients that will compute the
derivative at the last position. When dotted with `x` the result should
be 6.
>>> x = np.array([1, 0, 1, 4, 9])
>>> c = savgol_coeffs(5, 2, pos=4, deriv=1, use='dot')
>>> c.dot(x)
6.0
"""
# An alternative method for finding the coefficients when deriv=0 is
# t = np.arange(window_length)
# unit = (t == pos).astype(int)
# coeffs = np.polyval(np.polyfit(t, unit, polyorder), t)
# The method implemented here is faster.
# To recreate the table of sample coefficients shown in the chapter on
# the Savitzy-Golay filter in the Numerical Recipes book, use
# window_length = nL + nR + 1
# pos = nL + 1
# c = savgol_coeffs(window_length, M, pos=pos, use='dot')
if polyorder >= window_length:
raise ValueError("polyorder must be less than window_length.")
halflen, rem = divmod(window_length, 2)
if pos is None:
if rem == 0:
pos = halflen - 0.5
else:
pos = halflen
if not (0 <= pos < window_length):
raise ValueError("pos must be nonnegative and less than "
"window_length.")
if use not in ['conv', 'dot']:
raise ValueError("`use` must be 'conv' or 'dot'")
if deriv > polyorder:
coeffs = np.zeros(window_length)
return coeffs
# Form the design matrix A. The columns of A are powers of the integers
# from -pos to window_length - pos - 1. The powers (i.e., rows) range
# from 0 to polyorder. (That is, A is a vandermonde matrix, but not
# necessarily square.)
x = np.arange(-pos, window_length - pos, dtype=float)
if use == "conv":
# Reverse so that result can be used in a convolution.
x = x[::-1]
order = np.arange(polyorder + 1).reshape(-1, 1)
A = x ** order
# y determines which order derivative is returned.
y = np.zeros(polyorder + 1)
# The coefficient assigned to y[deriv] scales the result to take into
# account the order of the derivative and the sample spacing.
y[deriv] = float_factorial(deriv) / (delta ** deriv)
# Find the least-squares solution of A*c = y
coeffs, _, _, _ = lstsq(A, y)
return coeffs
def _polyder(p, m):
"""Differentiate polynomials represented with coefficients.
p must be a 1-D or 2-D array. In the 2-D case, each column gives
the coefficients of a polynomial; the first row holds the coefficients
associated with the highest power. m must be a nonnegative integer.
(numpy.polyder doesn't handle the 2-D case.)
"""
if m == 0:
result = p
else:
n = len(p)
if n <= m:
result = np.zeros_like(p[:1, ...])
else:
dp = p[:-m].copy()
for k in range(m):
rng = np.arange(n - k - 1, m - k - 1, -1)
dp *= rng.reshape((n - m,) + (1,) * (p.ndim - 1))
result = dp
return result
def _fit_edge(x, window_start, window_stop, interp_start, interp_stop,
axis, polyorder, deriv, delta, y):
"""
Given an N-d array `x` and the specification of a slice of `x` from
`window_start` to `window_stop` along `axis`, create an interpolating
polynomial of each 1-D slice, and evaluate that polynomial in the slice
from `interp_start` to `interp_stop`. Put the result into the
corresponding slice of `y`.
"""
# Get the edge into a (window_length, -1) array.
x_edge = axis_slice(x, start=window_start, stop=window_stop, axis=axis)
if axis == 0 or axis == -x.ndim:
xx_edge = x_edge
swapped = False
else:
xx_edge = x_edge.swapaxes(axis, 0)
swapped = True
xx_edge = xx_edge.reshape(xx_edge.shape[0], -1)
# Fit the edges. poly_coeffs has shape (polyorder + 1, -1),
# where '-1' is the same as in xx_edge.
poly_coeffs = np.polyfit(np.arange(0, window_stop - window_start),
xx_edge, polyorder)
if deriv > 0:
poly_coeffs = _polyder(poly_coeffs, deriv)
# Compute the interpolated values for the edge.
i = np.arange(interp_start - window_start, interp_stop - window_start)
values = np.polyval(poly_coeffs, i.reshape(-1, 1)) / (delta ** deriv)
# Now put the values into the appropriate slice of y.
# First reshape values to match y.
shp = list(y.shape)
shp[0], shp[axis] = shp[axis], shp[0]
values = values.reshape(interp_stop - interp_start, *shp[1:])
if swapped:
values = values.swapaxes(0, axis)
# Get a view of the data to be replaced by values.
y_edge = axis_slice(y, start=interp_start, stop=interp_stop, axis=axis)
y_edge[...] = values
def _fit_edges_polyfit(x, window_length, polyorder, deriv, delta, axis, y):
"""
Use polynomial interpolation of x at the low and high ends of the axis
to fill in the halflen values in y.
This function just calls _fit_edge twice, once for each end of the axis.
"""
halflen = window_length // 2
_fit_edge(x, 0, window_length, 0, halflen, axis,
polyorder, deriv, delta, y)
n = x.shape[axis]
_fit_edge(x, n - window_length, n, n - halflen, n, axis,
polyorder, deriv, delta, y)
[docs]
def savgol_filter(x, window_length, polyorder, deriv=0, delta=1.0,
axis=-1, mode='interp', cval=0.0):
""" Apply a Savitzky-Golay filter to an array.
This is a 1-D filter. If `x` has dimension greater than 1, `axis`
determines the axis along which the filter is applied.
Parameters
----------
x : array_like
The data to be filtered. If `x` is not a single or double precision
floating point array, it will be converted to type ``numpy.float64``
before filtering.
window_length : int
The length of the filter window (i.e., the number of coefficients).
If `mode` is 'interp', `window_length` must be less than or equal
to the size of `x`.
polyorder : int
The order of the polynomial used to fit the samples.
`polyorder` must be less than `window_length`.
deriv : int, optional
The order of the derivative to compute. This must be a
nonnegative integer. The default is 0, which means to filter
the data without differentiating.
delta : float, optional
The spacing of the samples to which the filter will be applied.
This is only used if deriv > 0. Default is 1.0.
axis : int, optional
The axis of the array `x` along which the filter is to be applied.
Default is -1.
mode : str, optional
Must be 'mirror', 'constant', 'nearest', 'wrap' or 'interp'. This
determines the type of extension to use for the padded signal to
which the filter is applied. When `mode` is 'constant', the padding
value is given by `cval`. See the Notes for more details on 'mirror',
'constant', 'wrap', and 'nearest'.
When the 'interp' mode is selected (the default), no extension
is used. Instead, a degree `polyorder` polynomial is fit to the
last `window_length` values of the edges, and this polynomial is
used to evaluate the last `window_length // 2` output values.
cval : scalar, optional
Value to fill past the edges of the input if `mode` is 'constant'.
Default is 0.0.
Returns
-------
y : ndarray, same shape as `x`
The filtered data.
See Also
--------
savgol_coeffs
Notes
-----
Details on the `mode` options:
'mirror':
Repeats the values at the edges in reverse order. The value
closest to the edge is not included.
'nearest':
The extension contains the nearest input value.
'constant':
The extension contains the value given by the `cval` argument.
'wrap':
The extension contains the values from the other end of the array.
For example, if the input is [1, 2, 3, 4, 5, 6, 7, 8], and
`window_length` is 7, the following shows the extended data for
the various `mode` options (assuming `cval` is 0)::
mode | Ext | Input | Ext
-----------+---------+------------------------+---------
'mirror' | 4 3 2 | 1 2 3 4 5 6 7 8 | 7 6 5
'nearest' | 1 1 1 | 1 2 3 4 5 6 7 8 | 8 8 8
'constant' | 0 0 0 | 1 2 3 4 5 6 7 8 | 0 0 0
'wrap' | 6 7 8 | 1 2 3 4 5 6 7 8 | 1 2 3
.. versionadded:: 0.14.0
Examples
--------
>>> import numpy as np
>>> from scipy.signal import savgol_filter
>>> np.set_printoptions(precision=2) # For compact display.
>>> x = np.array([2, 2, 5, 2, 1, 0, 1, 4, 9])
Filter with a window length of 5 and a degree 2 polynomial. Use
the defaults for all other parameters.
>>> savgol_filter(x, 5, 2)
array([1.66, 3.17, 3.54, 2.86, 0.66, 0.17, 1. , 4. , 9. ])
Note that the last five values in x are samples of a parabola, so
when mode='interp' (the default) is used with polyorder=2, the last
three values are unchanged. Compare that to, for example,
`mode='nearest'`:
>>> savgol_filter(x, 5, 2, mode='nearest')
array([1.74, 3.03, 3.54, 2.86, 0.66, 0.17, 1. , 4.6 , 7.97])
"""
if mode not in ["mirror", "constant", "nearest", "interp", "wrap"]:
raise ValueError("mode must be 'mirror', 'constant', 'nearest' "
"'wrap' or 'interp'.")
x = np.asarray(x)
# Ensure that x is either single or double precision floating point.
if x.dtype != np.float64 and x.dtype != np.float32:
x = x.astype(np.float64)
coeffs = savgol_coeffs(window_length, polyorder, deriv=deriv, delta=delta)
if mode == "interp":
if window_length > x.shape[axis]:
raise ValueError("If mode is 'interp', window_length must be less "
"than or equal to the size of x.")
# Do not pad. Instead, for the elements within `window_length // 2`
# of the ends of the sequence, use the polynomial that is fitted to
# the last `window_length` elements.
y = convolve1d(x, coeffs, axis=axis, mode="constant")
_fit_edges_polyfit(x, window_length, polyorder, deriv, delta, axis, y)
else:
# Any mode other than 'interp' is passed on to ndimage.convolve1d.
y = convolve1d(x, coeffs, axis=axis, mode=mode, cval=cval)
return y
