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