diff --git a/README.rst b/README.rst index 48967571..520d4029 100644 --- a/README.rst +++ b/README.rst @@ -43,7 +43,7 @@ General purpose shared utilities for the diffpy libraries. The diffpy.utils package provides functions for extracting array data from variously formatted text files, an interpolation function based on the Whittaker-Shannon formula that can be used to resample a PDF or other profile -function over a new grid, `Diffraction_objects` for conveniently manipulating +function over a new grid, `DiffractionObject` for conveniently manipulating diffraction data, and some wx GUI utilities used by the PDFgui program. diff --git a/src/diffpy/utils/__init__.py b/src/diffpy/utils/__init__.py index 41a4fa54..c9909a8a 100644 --- a/src/diffpy/utils/__init__.py +++ b/src/diffpy/utils/__init__.py @@ -20,5 +20,3 @@ # silence the pyflakes syntax checker assert __version__ or True - -# End of file diff --git a/src/diffpy/utils/parsers/__init__.py b/src/diffpy/utils/parsers/__init__.py index cd95dd63..bab9943c 100644 --- a/src/diffpy/utils/parsers/__init__.py +++ b/src/diffpy/utils/parsers/__init__.py @@ -15,5 +15,3 @@ """Various utilities related to data parsing and manipulation. """ - -# End of file diff --git a/src/diffpy/utils/parsers/loaddata.py b/src/diffpy/utils/parsers/loaddata.py index 870e6016..39c4b163 100644 --- a/src/diffpy/utils/parsers/loaddata.py +++ b/src/diffpy/utils/parsers/loaddata.py @@ -344,6 +344,3 @@ def isfloat(s): except ValueError: pass return False - - -# End of file diff --git a/src/diffpy/utils/resampler.py b/src/diffpy/utils/resampler.py index e88c4f0b..5f6062d2 100644 --- a/src/diffpy/utils/resampler.py +++ b/src/diffpy/utils/resampler.py @@ -15,25 +15,26 @@ """Various utilities related to data parsing and manipulation.""" -import numpy +import numpy as np -# NOTE - this should be faster than resample below and conforms more closely to -# numpy.interp. I'm keeping resample for legacy reasons. def wsinterp(x, xp, fp, left=None, right=None): """One-dimensional Whittaker-Shannon interpolation. - This uses the Whittaker-Shannon interpolation formula to interpolate the value of fp (array), - which is defined over xp (array), at x (array or float). + Reconstruct a continuous signal from discrete data points by utilizing sinc functions + as interpolation kernels. This function interpolates the values of fp (array), + which are defined over xp (array), at new points x (array or float). + The implementation is based on E. T. Whittaker's 1915 paper + (https://doi.org/10.1017/S0370164600017806). Parameters ---------- x: ndarray - Desired range for interpolation. + The x values at which interpolation is computed. xp: ndarray - Defined range for fp. + The array of known x values. fp: ndarray - Function to be interpolated. + The array of y values associated with xp. left: float If given, set fp for x < xp[0] to left. Otherwise, if left is None (default) or not given, set fp for x < xp[0] to fp evaluated at xp[-1]. @@ -43,24 +44,23 @@ def wsinterp(x, xp, fp, left=None, right=None): Returns ------- - float: - If input x is a scalar (not an array), return the interpolated value at x. - ndarray: - If input x is an array, return the interpolated array with dimensions of x. + ndarray or float + The interpolated values at points x. Returns a single float if x is a scalar, + otherwise returns a numpy.ndarray. """ - scalar = numpy.isscalar(x) + scalar = np.isscalar(x) if scalar: - x = numpy.array(x) + x = np.array(x) x.resize(1) # shape = (nxp, nx), nxp copies of x data span axis 1 - u = numpy.resize(x, (len(xp), len(x))) + u = np.resize(x, (len(xp), len(x))) # Must take transpose of u for proper broadcasting with xp. # shape = (nx, nxp), v(xp) data spans axis 1 v = (xp - u.T) / (xp[1] - xp[0]) # shape = (nx, nxp), m(v) data spans axis 1 - m = fp * numpy.sinc(v) + m = fp * np.sinc(v) # Sum over m(v) (axis 1) - fp_at_x = numpy.sum(m, axis=1) + fp_at_x = np.sum(m, axis=1) # Enforce left and right if left is None: @@ -100,36 +100,33 @@ def resample(r, s, dr): dr0 = r[1] - r[0] # Constant timestep if dr0 < dr: - rnew = numpy.arange(r[0], r[-1] + 0.5 * dr, dr) - snew = numpy.interp(rnew, r, s) + rnew = np.arange(r[0], r[-1] + 0.5 * dr, dr) + snew = np.interp(rnew, r, s) return rnew, snew elif dr0 > dr: # Tried to pad the end of s to dampen, but nothing works. # m = (s[-1] - s[-2]) / dr0 # b = (s[-2] * r[-1] - s[-1] * r[-2]) / dr0 - # rpad = r[-1] + numpy.arange(1, len(s))*dr0 + # rpad = r[-1] + np.arange(1, len(s))*dr0 # spad = rpad * m + b - # spad = numpy.concatenate([s,spad]) - # rnew = numpy.arange(0, rpad[-1], dr) - # snew = numpy.zeros_like(rnew) + # spad = np.concatenate([s,spad]) + # rnew = np.arange(0, rpad[-1], dr) + # snew = np.zeros_like(rnew) # Accommodate for the fact that r[0] might not be 0 # u = (rnew-r[0]) / dr0 # for n in range(len(spad)): - # snew += spad[n] * numpy.sinc(u - n) + # snew += spad[n] * np.sinc(u - n) - # sel = numpy.logical_and(rnew >= r[0], rnew <= r[-1]) + # sel = np.logical_and(rnew >= r[0], rnew <= r[-1]) - rnew = numpy.arange(0, r[-1], dr) - snew = numpy.zeros_like(rnew) + rnew = np.arange(0, r[-1], dr) + snew = np.zeros_like(rnew) u = (rnew - r[0]) / dr0 for n in range(len(s)): - snew += s[n] * numpy.sinc(u - n) - sel = numpy.logical_and(rnew >= r[0], rnew <= r[-1]) + snew += s[n] * np.sinc(u - n) + sel = np.logical_and(rnew >= r[0], rnew <= r[-1]) return rnew[sel], snew[sel] # If we got here, then no resampling is required return r.copy(), s.copy() - - -# End of file diff --git a/src/diffpy/utils/version.py b/src/diffpy/utils/version.py index e42fed78..0bc397e4 100644 --- a/src/diffpy/utils/version.py +++ b/src/diffpy/utils/version.py @@ -22,5 +22,3 @@ from importlib.metadata import version __version__ = version("diffpy.utils") - -# End of file diff --git a/src/diffpy/utils/wx/__init__.py b/src/diffpy/utils/wx/__init__.py index c2d0617b..3f7417ef 100644 --- a/src/diffpy/utils/wx/__init__.py +++ b/src/diffpy/utils/wx/__init__.py @@ -15,5 +15,3 @@ """Utilities related wx Python GUIs. """ - -# End of file diff --git a/src/diffpy/utils/wx/gridutils.py b/src/diffpy/utils/wx/gridutils.py index 1a25bfb3..c459ac30 100644 --- a/src/diffpy/utils/wx/gridutils.py +++ b/src/diffpy/utils/wx/gridutils.py @@ -163,6 +163,3 @@ def _indicesToBlocks(indices): i0 = i rv = [tuple(ij) for ij in rngs] return rv - - -# End of file diff --git a/tests/test_resample.py b/tests/test_resample.py index deddcca2..784932bb 100644 --- a/tests/test_resample.py +++ b/tests/test_resample.py @@ -1,3 +1,5 @@ +import random + import numpy as np import pytest @@ -5,26 +7,26 @@ def test_wsinterp(): - import random - - # Check known points are unchanged by interpolation # FIXME: if another SW interp function exists, run comparisons for interpolated points + # Sampling rate ssr = 44100**-1 # Standard sampling rate for human-hearable frequencies - t = ssr + + # Creating a symmetric set of sample points around zero. n = 5 - xp = np.array([i * t for i in range(-n, n + 1, 1)]) - x = np.array([i * t for i in range(-n - 1, n + 2, 1)]) + xp = np.array([i * ssr for i in range(-n, n + 1, 1)]) + x = np.array([i * ssr for i in range(-n - 1, n + 2, 1)]) + assert len(xp) == 11 and len(x) == 13 - # Interpolate a few random datasets + # Generate a new set of fp values across 10 trial runs trials = 10 - for trial in range(trials): - fp = np.array([random.random() * ssr for i in range(-n, n + 1, 1)]) + + for _ in range(trials): + # Create random function values (fp) at the points defined in xp above + fp = np.array([random.random() * ssr for _ in range(-n, n + 1, 1)]) fp_at_x = wsinterp(x, xp, fp) + + # Check that the known points are unchanged by interpolation assert np.allclose(fp_at_x[1:-1], fp) for i in range(len(x)): assert fp_at_x[i] == pytest.approx(wsinterp(x[i], xp, fp)) - - -if __name__ == "__main__": - test_wsinterp()