Skip to content
Merged
Show file tree
Hide file tree
Changes from 5 commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion README.rst
Original file line number Diff line number Diff line change
Expand Up @@ -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.

Expand Down
2 changes: 0 additions & 2 deletions src/diffpy/utils/__init__.py
Original file line number Diff line number Diff line change
Expand Up @@ -20,5 +20,3 @@

# silence the pyflakes syntax checker
assert __version__ or True

# End of file
2 changes: 0 additions & 2 deletions src/diffpy/utils/parsers/__init__.py
Original file line number Diff line number Diff line change
Expand Up @@ -15,5 +15,3 @@

"""Various utilities related to data parsing and manipulation.
"""

# End of file
3 changes: 0 additions & 3 deletions src/diffpy/utils/parsers/loaddata.py
Original file line number Diff line number Diff line change
Expand Up @@ -344,6 +344,3 @@ def isfloat(s):
except ValueError:
pass
return False


# End of file
61 changes: 29 additions & 32 deletions src/diffpy/utils/resampler.py
Original file line number Diff line number Diff line change
Expand Up @@ -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].
Expand All @@ -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:
Expand Down Expand Up @@ -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
2 changes: 0 additions & 2 deletions src/diffpy/utils/version.py
Original file line number Diff line number Diff line change
Expand Up @@ -22,5 +22,3 @@
from importlib.metadata import version

__version__ = version("diffpy.utils")

# End of file
2 changes: 0 additions & 2 deletions src/diffpy/utils/wx/__init__.py
Original file line number Diff line number Diff line change
Expand Up @@ -15,5 +15,3 @@

"""Utilities related wx Python GUIs.
"""

# End of file
3 changes: 0 additions & 3 deletions src/diffpy/utils/wx/gridutils.py
Original file line number Diff line number Diff line change
Expand Up @@ -163,6 +163,3 @@ def _indicesToBlocks(indices):
i0 = i
rv = [tuple(ij) for ij in rngs]
return rv


# End of file
28 changes: 15 additions & 13 deletions tests/test_resample.py
Original file line number Diff line number Diff line change
@@ -1,30 +1,32 @@
import random

import numpy as np
import pytest

from diffpy.utils.resampler import wsinterp


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()
Loading