|
| 1 | +import numpy as np |
| 2 | +import scipy |
| 3 | +from warnings import warn |
| 4 | + |
| 5 | +from pynumdiff.utils import utility |
| 6 | + |
| 7 | + |
| 8 | +def splinediff(x, dt, params=None, options={}, order=3, s=None, num_iterations=1): |
| 9 | + """Find smoothed data and derivative estimates by fitting a smoothing spline to the data with |
| 10 | + scipy.interpolate.UnivariateSpline. |
| 11 | +
|
| 12 | + :param np.array[float] x: data to differentiate |
| 13 | + :param float dt: step size |
| 14 | + :param list params: (**deprecated**, prefer :code:`order`, :code:`cutoff_freq`, and :code:`num_iterations`) |
| 15 | + :param dict options: (**deprecated**, prefer :code:`num_iterations`) an empty dictionary or {'iterate': (bool)} |
| 16 | + :param int order: polynomial order of the spline. A kth order spline can be differentiated k times. |
| 17 | + :param float s: positive smoothing factor used to choose the number of knots. Number of knots will be increased |
| 18 | + until the smoothing condition is satisfied: :math:`\\sum_t (x[t] - \\text{spline}[t])^2 \\leq s` |
| 19 | + :param int num_iterations: how many times to apply smoothing |
| 20 | +
|
| 21 | + :return: tuple[np.array, np.array] of\n |
| 22 | + - **x_hat** -- estimated (smoothed) x |
| 23 | + - **dxdt_hat** -- estimated derivative of x |
| 24 | + """ |
| 25 | + if params != None: # Warning to support old interface for a while. Remove these lines along with params in a future release. |
| 26 | + warn("`params` and `options` parameters will be removed in a future version. Use `order`, `s`, and " + |
| 27 | + "`num_iterations` instead.", DeprecationWarning) |
| 28 | + order, s = params[0:2] |
| 29 | + if 'iterate' in options and options['iterate']: |
| 30 | + num_iterations = params[2] |
| 31 | + |
| 32 | + t = np.arange(len(x))*dt |
| 33 | + |
| 34 | + x_hat = x |
| 35 | + for _ in range(num_iterations): |
| 36 | + spline = scipy.interpolate.UnivariateSpline(t, x_hat, k=order, s=s) |
| 37 | + x_hat = spline(t) |
| 38 | + |
| 39 | + dspline = spline.derivative() |
| 40 | + dxdt_hat = dspline(t) |
| 41 | + |
| 42 | + return x_hat, dxdt_hat |
| 43 | + |
| 44 | + |
| 45 | +def polydiff(x, dt, params=None, options=None, poly_order=None, window_size=None, step_size=1, |
| 46 | + kernel='friedrichs'): |
| 47 | + """Fit polynomials to the data, and differentiate the polynomials. |
| 48 | +
|
| 49 | + :param np.array[float] x: data to differentiate |
| 50 | + :param float dt: step size |
| 51 | + :param list[int] params: (**deprecated**, prefer :code:`poly_order` and :code:`window_size`) |
| 52 | + :param dict options: (**deprecated**, prefer :code:`step_size` and :code:`kernel`) |
| 53 | + a dictionary consisting of {'sliding': (bool), 'step_size': (int), 'kernel_name': (str)} |
| 54 | + :param int poly_order: order of the polynomial |
| 55 | + :param int window_size: size of the sliding window, if not given no sliding |
| 56 | + :param int step_size: step size for sliding |
| 57 | + :param str kernel: name of kernel to use for weighting and smoothing windows ('gaussian' or 'friedrichs') |
| 58 | +
|
| 59 | + :return: tuple[np.array, np.array] of\n |
| 60 | + - **x_hat** -- estimated (smoothed) x |
| 61 | + - **dxdt_hat** -- estimated derivative of x |
| 62 | + """ |
| 63 | + if params != None: |
| 64 | + warn("`params` and `options` parameters will be removed in a future version. Use `poly_order` " + |
| 65 | + "and `window_size` instead.", DeprecationWarning) |
| 66 | + poly_order = params[0] |
| 67 | + if len(params) > 1: window_size = params[1] |
| 68 | + if options != None: |
| 69 | + if 'sliding' in options and not options['sliding']: window_size = None |
| 70 | + if 'step_size' in options: step_size = options['step_size'] |
| 71 | + if 'kernel_name' in options: kernel = options['kernel_name'] |
| 72 | + elif poly_order == None or window_size == None: |
| 73 | + raise ValueError("`poly_order` and `window_size` must be given.") |
| 74 | + |
| 75 | + if window_size < poly_order*3: |
| 76 | + window_size = poly_order*3+1 |
| 77 | + if window_size % 2 == 0: |
| 78 | + window_size += 1 |
| 79 | + warn("Kernel window size should be odd. Added 1 to length.") |
| 80 | + |
| 81 | + def _polydiff(x, dt, poly_order, weights=None): |
| 82 | + t = np.arange(len(x))*dt |
| 83 | + |
| 84 | + r = np.polyfit(t, x, poly_order, w=weights) # polyfit returns highest order first |
| 85 | + dr = np.polyder(r) # power rule already implemented for us |
| 86 | + |
| 87 | + dxdt_hat = np.polyval(dr, t) # evaluate the derivative and original polynomials at points t |
| 88 | + x_hat = np.polyval(r, t) # smoothed x |
| 89 | + |
| 90 | + return x_hat, dxdt_hat |
| 91 | + |
| 92 | + if not window_size: |
| 93 | + return _polydiff(x, dt, poly_order) |
| 94 | + |
| 95 | + kernel = {'gaussian':utility.gaussian_kernel, 'friedrichs':utility.friedrichs_kernel}[kernel](window_size) |
| 96 | + return utility.slide_function(_polydiff, x, dt, kernel, poly_order, stride=step_size, pass_weights=True) |
| 97 | + |
| 98 | + |
| 99 | +def savgoldiff(x, dt, params=None, options=None, poly_order=None, window_size=None, smoothing_win=None): |
| 100 | + """Use the Savitzky-Golay to smooth the data and calculate the first derivative. It uses |
| 101 | + scipy.signal.savgol_filter. The Savitzky-Golay is very similar to the sliding polynomial fit, |
| 102 | + but slightly noisier, and much faster. |
| 103 | +
|
| 104 | + :param np.array[float] x: data to differentiate |
| 105 | + :param float dt: step size |
| 106 | + :param list params: (**deprecated**, prefer :code:`poly_order`, :code:`window_size`, and :code:`smoothing_win`) |
| 107 | + :param dict options: (**deprecated**) |
| 108 | + :param int poly_order: order of the polynomial |
| 109 | + :param int window_size: size of the sliding window, must be odd (if not, 1 is added) |
| 110 | + :param int smoothing_win: size of the window used for gaussian smoothing, a good default is |
| 111 | + window_size, but smaller for high frequnecy data |
| 112 | +
|
| 113 | + :return: tuple[np.array, np.array] of\n |
| 114 | + - **x_hat** -- estimated (smoothed) x |
| 115 | + - **dxdt_hat** -- estimated derivative of x |
| 116 | + """ |
| 117 | + if params != None: # Warning to support old interface for a while. Remove these lines along with params in a future release. |
| 118 | + warn("`params` and `options` parameters will be removed in a future version. Use `poly_order`, " + |
| 119 | + "`window_size`, and `smoothing_win` instead.", DeprecationWarning) |
| 120 | + poly_order, window_size, smoothing_win = params |
| 121 | + elif poly_order == None or window_size == None or smoothing_win == None: |
| 122 | + raise ValueError("`poly_order`, `window_size`, and `smoothing_win` must be given.") |
| 123 | + |
| 124 | + window_size = np.clip(window_size, poly_order + 1, len(x)-1) |
| 125 | + if window_size % 2 == 0: |
| 126 | + window_size += 1 # window_size needs to be odd |
| 127 | + warn("Kernel window size should be odd. Added 1 to length.") |
| 128 | + smoothing_win = min(smoothing_win, len(x)-1) |
| 129 | + |
| 130 | + dxdt_hat = scipy.signal.savgol_filter(x, window_size, poly_order, deriv=1)/dt |
| 131 | + |
| 132 | + kernel = utility.gaussian_kernel(smoothing_win) |
| 133 | + dxdt_hat = utility.convolutional_smoother(dxdt_hat, kernel) |
| 134 | + |
| 135 | + x_hat = utility.integrate_dxdt_hat(dxdt_hat, dt) |
| 136 | + x0 = utility.estimate_integration_constant(x, x_hat) |
| 137 | + x_hat = x_hat + x0 |
| 138 | + |
| 139 | + return x_hat, dxdt_hat |
0 commit comments