|
15 | 15 |
|
16 | 16 | """Various utilities related to data parsing and manipulation."""
|
17 | 17 |
|
18 |
| -import numpy |
| 18 | +import numpy as np |
19 | 19 |
|
20 | 20 |
|
21 |
| -# NOTE - this should be faster than resample below and conforms more closely to |
22 |
| -# numpy.interp. I'm keeping resample for legacy reasons. |
23 | 21 | def wsinterp(x, xp, fp, left=None, right=None):
|
24 | 22 | """One-dimensional Whittaker-Shannon interpolation.
|
25 | 23 |
|
26 |
| - This uses the Whittaker-Shannon interpolation formula to interpolate the value of fp (array), |
27 |
| - which is defined over xp (array), at x (array or float). |
| 24 | + Reconstruct a continuous signal from discrete data points by utilizing sinc functions |
| 25 | + as interpolation kernels. This function interpolates the values of fp (array), |
| 26 | + which are defined over xp (array), at new points x (array or float). |
28 | 27 |
|
29 | 28 | Parameters
|
30 | 29 | ----------
|
31 | 30 | x: ndarray
|
32 |
| - Desired range for interpolation. |
| 31 | + The x values at which interpolation is computed. |
33 | 32 | xp: ndarray
|
34 |
| - Defined range for fp. |
| 33 | + The array of known x values. |
35 | 34 | fp: ndarray
|
36 |
| - Function to be interpolated. |
| 35 | + The array of y values associated xp. |
37 | 36 | left: float
|
38 | 37 | If given, set fp for x < xp[0] to left. Otherwise, if left is None (default) or not given,
|
39 | 38 | set fp for x < xp[0] to fp evaluated at xp[-1].
|
40 | 39 | right: float
|
41 | 40 | If given, set fp for x > xp[-1] to right. Otherwise, if right is None (default) or not given, set fp for
|
42 |
| - x > xp[-1] to fp evaluated at xp[-1]. |
| 41 | + x > xp[-1] to fp evaluated at xp[-1] |
43 | 42 |
|
44 | 43 | Returns
|
45 | 44 | -------
|
46 |
| - float: |
47 |
| - If input x is a scalar (not an array), return the interpolated value at x. |
48 |
| - ndarray: |
49 |
| - If input x is an array, return the interpolated array with dimensions of x. |
| 45 | + ndarray or float |
| 46 | + Interpolated values at points x. Returns a single float if x is a scalar, |
| 47 | + otherwise returns a numpy.ndarray. |
50 | 48 | """
|
51 |
| - scalar = numpy.isscalar(x) |
| 49 | + scalar = np.isscalar(x) |
52 | 50 | if scalar:
|
53 |
| - x = numpy.array(x) |
| 51 | + x = np.array(x) |
54 | 52 | x.resize(1)
|
55 | 53 | # shape = (nxp, nx), nxp copies of x data span axis 1
|
56 |
| - u = numpy.resize(x, (len(xp), len(x))) |
| 54 | + u = np.resize(x, (len(xp), len(x))) |
57 | 55 | # Must take transpose of u for proper broadcasting with xp.
|
58 | 56 | # shape = (nx, nxp), v(xp) data spans axis 1
|
59 | 57 | v = (xp - u.T) / (xp[1] - xp[0])
|
60 | 58 | # shape = (nx, nxp), m(v) data spans axis 1
|
61 |
| - m = fp * numpy.sinc(v) |
| 59 | + m = fp * np.sinc(v) |
62 | 60 | # Sum over m(v) (axis 1)
|
63 |
| - fp_at_x = numpy.sum(m, axis=1) |
| 61 | + fp_at_x = np.sum(m, axis=1) |
64 | 62 |
|
65 | 63 | # Enforce left and right
|
66 | 64 | if left is None:
|
@@ -100,36 +98,33 @@ def resample(r, s, dr):
|
100 | 98 | dr0 = r[1] - r[0] # Constant timestep
|
101 | 99 |
|
102 | 100 | if dr0 < dr:
|
103 |
| - rnew = numpy.arange(r[0], r[-1] + 0.5 * dr, dr) |
104 |
| - snew = numpy.interp(rnew, r, s) |
| 101 | + rnew = np.arange(r[0], r[-1] + 0.5 * dr, dr) |
| 102 | + snew = np.interp(rnew, r, s) |
105 | 103 | return rnew, snew
|
106 | 104 |
|
107 | 105 | elif dr0 > dr:
|
108 | 106 | # Tried to pad the end of s to dampen, but nothing works.
|
109 | 107 | # m = (s[-1] - s[-2]) / dr0
|
110 | 108 | # b = (s[-2] * r[-1] - s[-1] * r[-2]) / dr0
|
111 |
| - # rpad = r[-1] + numpy.arange(1, len(s))*dr0 |
| 109 | + # rpad = r[-1] + np.arange(1, len(s))*dr0 |
112 | 110 | # spad = rpad * m + b
|
113 |
| - # spad = numpy.concatenate([s,spad]) |
114 |
| - # rnew = numpy.arange(0, rpad[-1], dr) |
115 |
| - # snew = numpy.zeros_like(rnew) |
| 111 | + # spad = np.concatenate([s,spad]) |
| 112 | + # rnew = np.arange(0, rpad[-1], dr) |
| 113 | + # snew = np.zeros_like(rnew) |
116 | 114 | # Accommodate for the fact that r[0] might not be 0
|
117 | 115 | # u = (rnew-r[0]) / dr0
|
118 | 116 | # for n in range(len(spad)):
|
119 |
| - # snew += spad[n] * numpy.sinc(u - n) |
| 117 | + # snew += spad[n] * np.sinc(u - n) |
120 | 118 |
|
121 |
| - # sel = numpy.logical_and(rnew >= r[0], rnew <= r[-1]) |
| 119 | + # sel = np.logical_and(rnew >= r[0], rnew <= r[-1]) |
122 | 120 |
|
123 |
| - rnew = numpy.arange(0, r[-1], dr) |
124 |
| - snew = numpy.zeros_like(rnew) |
| 121 | + rnew = np.arange(0, r[-1], dr) |
| 122 | + snew = np.zeros_like(rnew) |
125 | 123 | u = (rnew - r[0]) / dr0
|
126 | 124 | for n in range(len(s)):
|
127 |
| - snew += s[n] * numpy.sinc(u - n) |
128 |
| - sel = numpy.logical_and(rnew >= r[0], rnew <= r[-1]) |
| 125 | + snew += s[n] * np.sinc(u - n) |
| 126 | + sel = np.logical_and(rnew >= r[0], rnew <= r[-1]) |
129 | 127 | return rnew[sel], snew[sel]
|
130 | 128 |
|
131 | 129 | # If we got here, then no resampling is required
|
132 | 130 | return r.copy(), s.copy()
|
133 |
| - |
134 |
| - |
135 |
| -# End of file |
0 commit comments