|
| 1 | +import numpy as np |
| 2 | +from warnings import warn |
| 3 | +from scipy import sparse |
| 4 | + |
| 5 | +from pynumdiff.utils import utility |
| 6 | + |
| 7 | + |
| 8 | +def spectraldiff(x, dt, params=None, options=None, high_freq_cutoff=None, even_extension=True, pad_to_zero_dxdt=True): |
| 9 | + """Take a derivative in the fourier domain, with high frequency attentuation. |
| 10 | +
|
| 11 | + :param np.array[float] x: data to differentiate |
| 12 | + :param float dt: step size |
| 13 | + :param list[float] or float params: (**deprecated**, prefer :code:`high_freq_cutoff`) |
| 14 | + :param dict options: (**deprecated**, prefer :code:`even_extension` |
| 15 | + and :code:`pad_to_zero_dxdt`) a dictionary consisting of {'even_extension': (bool), 'pad_to_zero_dxdt': (bool)} |
| 16 | + :param float high_freq_cutoff: The high frequency cutoff as a multiple of the Nyquist frequency: Should be between 0 |
| 17 | + and 1. Frequencies below this threshold will be kept, and above will be zeroed. |
| 18 | + :param bool even_extension: if True, extend the data with an even extension so signal starts and ends at the same value. |
| 19 | + :param bool pad_to_zero_dxdt: if True, extend the data with extra regions that smoothly force the derivative to |
| 20 | + zero before taking FFT. |
| 21 | +
|
| 22 | + :return: tuple[np.array, np.array] of\n |
| 23 | + - **x_hat** -- estimated (smoothed) x |
| 24 | + - **dxdt_hat** -- estimated derivative of x |
| 25 | + """ |
| 26 | + if params != None: # Warning to support old interface for a while. Remove these lines along with params in a future release. |
| 27 | + warn("`params` and `options` parameters will be removed in a future version. Use `high_freq_cutoff`, " + |
| 28 | + "`even_extension`, and `pad_to_zero_dxdt` instead.", DeprecationWarning) |
| 29 | + high_freq_cutoff = params[0] if isinstance(params, list) else params |
| 30 | + if options != None: |
| 31 | + if 'even_extension' in options: even_extension = options['even_extension'] |
| 32 | + if 'pad_to_zero_dxdt' in options: pad_to_zero_dxdt = options['pad_to_zero_dxdt'] |
| 33 | + elif high_freq_cutoff == None: |
| 34 | + raise ValueError("`high_freq_cutoff` must be given.") |
| 35 | + |
| 36 | + L = len(x) |
| 37 | + |
| 38 | + # make derivative go to zero at ends (optional) |
| 39 | + if pad_to_zero_dxdt: |
| 40 | + padding = 100 |
| 41 | + pre = x[0]*np.ones(padding) # extend the edges |
| 42 | + post = x[-1]*np.ones(padding) |
| 43 | + x = np.hstack((pre, x, post)) |
| 44 | + kernel = utility.mean_kernel(padding//2) |
| 45 | + x_hat = utility.convolutional_smoother(x, kernel) # smooth the edges in |
| 46 | + x_hat[padding:-padding] = x[padding:-padding] # replace middle with original signal |
| 47 | + x = x_hat |
| 48 | + else: |
| 49 | + padding = 0 |
| 50 | + |
| 51 | + # Do even extension (optional) |
| 52 | + if even_extension is True: |
| 53 | + x = np.hstack((x, x[::-1])) |
| 54 | + |
| 55 | + # If odd, make N even, and pad x |
| 56 | + N = len(x) |
| 57 | + |
| 58 | + # Define the frequency range. |
| 59 | + k = np.concatenate((np.arange(N//2 + 1), np.arange(-N//2 + 1, 0))) |
| 60 | + if N % 2 == 0: k[N//2] = 0 # odd derivatives get the Nyquist element zeroed out |
| 61 | + omega = k*2*np.pi/(dt*N) # turn wavenumbers into frequencies in radians/s |
| 62 | + |
| 63 | + # Frequency based smoothing: remove signals with a frequency higher than high_freq_cutoff |
| 64 | + discrete_cutoff = int(high_freq_cutoff*N/2) # Nyquist is at N/2 location, and we're cutting off as a fraction of that |
| 65 | + omega[discrete_cutoff:N-discrete_cutoff] = 0 |
| 66 | + |
| 67 | + # Derivative = 90 deg phase shift |
| 68 | + dxdt_hat = np.real(np.fft.ifft(1.0j * omega * np.fft.fft(x))) |
| 69 | + dxdt_hat = dxdt_hat[padding:L+padding] |
| 70 | + |
| 71 | + # Integrate to get x_hat |
| 72 | + x_hat = utility.integrate_dxdt_hat(dxdt_hat, dt) |
| 73 | + x0 = utility.estimate_integration_constant(x[padding:L+padding], x_hat) |
| 74 | + x_hat = x_hat + x0 |
| 75 | + |
| 76 | + return x_hat, dxdt_hat |
| 77 | + |
| 78 | + |
| 79 | +def rbfdiff(x, _t, sigma=1, lmbd=0.01): |
| 80 | + """Find smoothed function and derivative estimates by fitting noisy data with radial-basis-functions. Naively, |
| 81 | + fill a matrix with basis function samples and solve a linear inverse problem against the data, but truncate tiny |
| 82 | + values to make columns sparse. Each basis function "hill" is topped with a "tower" of height :code:`lmbd` to reach |
| 83 | + noisy data samples, and the final smoothed reconstruction is found by razing these and only keeping the hills. |
| 84 | +
|
| 85 | + :param np.array[float] x: data to differentiate |
| 86 | + :param float or array[float] _t: This function supports variable step size. This parameter is either the constant |
| 87 | + :math:`\\Delta t` if given as a single float, or data locations if given as an array of same length as :code:`x`. |
| 88 | + :param float sigma: controls width of radial basis functions |
| 89 | + :param float lmbd: controls smoothness |
| 90 | +
|
| 91 | + :return: tuple[np.array, np.array] of\n |
| 92 | + - **x_hat** -- estimated (smoothed) x |
| 93 | + - **dxdt_hat** -- estimated derivative of x |
| 94 | + """ |
| 95 | + if np.isscalar(_t): |
| 96 | + t = np.arange(len(x))*_t |
| 97 | + else: # support variable step size for this function |
| 98 | + if len(x) != len(_t): raise ValueError("If `_t` is given as array-like, must have same length as `x`.") |
| 99 | + t = _t |
| 100 | + |
| 101 | + # The below does the approximate equivalent of this code, but sparsely in O(N sigma^2), since the rbf falls off rapidly |
| 102 | + # t_i, t_j = np.meshgrid(t,t) |
| 103 | + # r = t_j - t_i # radius |
| 104 | + # rbf = np.exp(-(r**2) / (2 * sigma**2)) # radial basis function kernel, O(N^2) entries |
| 105 | + # drbfdt = -(r / sigma**2) * rbf # derivative of kernel |
| 106 | + # rbf_regularized = rbf + lmbd*np.eye(len(t)) |
| 107 | + # alpha = np.linalg.solve(rbf_regularized, x) # O(N^3) |
| 108 | + |
| 109 | + cutoff = np.sqrt(-2 * sigma**2 * np.log(1e-4)) |
| 110 | + rows, cols, vals, dvals = [], [], [], [] |
| 111 | + for n in range(len(t)): |
| 112 | + # Only consider points within a cutoff. Gaussian drops below eps at distance ~ sqrt(-2*sigma^2 log eps) |
| 113 | + l = np.searchsorted(t, t[n] - cutoff) # O(log N) to find indices of points within cutoff |
| 114 | + r = np.searchsorted(t, t[n] + cutoff) # finds index where new value should be inserted |
| 115 | + for j in range(l, r): # width of this is dependent on sigma. [l, r) is correct inclusion/exclusion |
| 116 | + radius = t[n] - t[j] |
| 117 | + v = np.exp(-radius**2 / (2 * sigma**2)) |
| 118 | + dv = -radius / sigma**2 * v # take derivative of radial basis function, because d/dt coef*f(t) = coef*df/dt |
| 119 | + rows.append(n); cols.append(j); vals.append(v); dvals.append(dv) |
| 120 | + |
| 121 | + rbf = sparse.csr_matrix((vals, (rows, cols)), shape=(len(t), len(t))) # Build sparse kernels, O(N sigma) entries |
| 122 | + drbfdt = sparse.csr_matrix((dvals, (rows, cols)), shape=(len(t), len(t))) |
| 123 | + rbf_regularized = rbf + lmbd*sparse.eye(len(t), format="csr") # identity matrix gives a little extra height at the centers |
| 124 | + alpha = sparse.linalg.spsolve(rbf_regularized, x) # solve sparse system targeting the noisy data, O(N sigma^2) |
| 125 | + |
| 126 | + return rbf @ alpha, drbfdt @ alpha # find samples of reconstructions using the smooth bases |
0 commit comments