diff --git a/doc/source/examples/resampleexample.rst b/doc/source/examples/resampleexample.rst
index f7ba5e18..3166e42b 100644
--- a/doc/source/examples/resampleexample.rst
+++ b/doc/source/examples/resampleexample.rst
@@ -58,7 +58,7 @@ given enough datapoints.
      nickel_resample = wsinterp(grid, nickel_grid, nickel_func)
      target_resample = wsinterp(grid, target_grid, target_func)
 
-   We can now plot the difference to see that these two functions are quite similar.:
+   We can now plot the difference to see that these two functions are quite similar.::
 
      plt.plot(grid, target_resample)
      plt.plot(grid, nickel_resample)
@@ -78,3 +78,10 @@ given enough datapoints.
    In the case of our dataset, our band-limit is ``qmax=25.0`` and our function spans :math:`r \in (0.0, 60.0)`.
    Thus, our original grid requires :math:`25.0 * 60.0 / \pi < 478`. Since our grid has :math:`601` datapoints, our
    reconstruction was perfect as shown from the comparison between ``Nickel.gr`` and ``NiTarget.gr``.
+
+   This computation is implemented in the function ``nsinterp``.::
+
+     from diffpy.utils.resampler import nsinterp
+     qmin = 0
+     qmax = 25
+     nickel_resample = (nickel_grid, nickel_func, qmin, qmax)
diff --git a/news/nsinterp.rst b/news/nsinterp.rst
new file mode 100644
index 00000000..9b716b87
--- /dev/null
+++ b/news/nsinterp.rst
@@ -0,0 +1,23 @@
+**Added:**
+
+* Function nsinterp for automatic interpolation onto the Nyquist-Shannon grid.
+
+**Changed:**
+
+* <news item>
+
+**Deprecated:**
+
+* <news item>
+
+**Removed:**
+
+* <news item>
+
+**Fixed:**
+
+* <news item>
+
+**Security:**
+
+* <news item>
diff --git a/src/diffpy/utils/resampler.py b/src/diffpy/utils/resampler.py
index 5f6062d2..890f29fb 100644
--- a/src/diffpy/utils/resampler.py
+++ b/src/diffpy/utils/resampler.py
@@ -77,6 +77,46 @@ def wsinterp(x, xp, fp, left=None, right=None):
     return fp_at_x
 
 
+def nsinterp(xp, fp, qmin=0, qmax=25, left=None, right=None):
+    """One-dimensional Whittaker-Shannon interpolation onto the Nyquist-Shannon grid.
+
+    Takes a band-limited function fp and original grid xp and resamples fp on the NS grid.
+    Uses the minimum number of points N required by the Nyquist sampling theorem.
+    N = (qmax-qmin)(rmax-rmin)/pi, where rmin and rmax are the ends of the real-space ranges.
+    fp must be finite, and the user inputs qmin and qmax of the frequency-domain.
+
+    Parameters
+    ----------
+    xp: ndarray
+        The array of known x values.
+    fp: ndarray
+        The array of y values associated with xp.
+    qmin: float
+        The lower band limit in the frequency domain.
+    qmax: float
+        The upper band limit in the frequency domain.
+
+    Returns
+    -------
+    x: ndarray
+        The Nyquist-Shannon grid computed for the given qmin and qmax.
+    fp_at_x: ndarray
+        The interpolated values at points x. Returns a single float if x is a scalar,
+        otherwise returns a numpy.ndarray.
+    """
+    # Ensure numpy array
+    xp = np.array(xp)
+    rmin = np.min(xp)
+    rmax = np.max(xp)
+
+    nspoints = int(np.round((qmax - qmin) * (rmax - rmin) / np.pi))
+
+    x = np.linspace(rmin, rmax, nspoints)
+    fp_at_x = wsinterp(x, xp, fp)
+
+    return x, fp_at_x
+
+
 def resample(r, s, dr):
     """Resample a PDF on a new grid.
 
diff --git a/tests/test_resample.py b/tests/test_resample.py
index 784932bb..6e6294ed 100644
--- a/tests/test_resample.py
+++ b/tests/test_resample.py
@@ -3,7 +3,7 @@
 import numpy as np
 import pytest
 
-from diffpy.utils.resampler import wsinterp
+from diffpy.utils.resampler import nsinterp, wsinterp
 
 
 def test_wsinterp():
@@ -30,3 +30,23 @@ def test_wsinterp():
         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))
+
+
+def test_nsinterp():
+    # Create a cosine function cos(2x) for x \in [0, 3pi]
+    xp = np.linspace(0, 3 * np.pi, 100)
+    fp = np.cos(2 * xp)
+
+    # Want to resample onto the grid {0, pi, 2pi, 3pi}
+    x = np.array([0, np.pi, 2 * np.pi, 3 * np.pi])
+
+    # Get wsinterp result
+    ws_f = wsinterp(x, xp, fp)
+
+    # Use nsinterp with qmin-qmax=4/3
+    qmin = np.random.uniform()
+    qmax = qmin + 4 / 3
+    ns_x, ns_f = nsinterp(xp, fp, qmin, qmax)
+
+    assert np.allclose(x, ns_x)
+    assert np.allclose(ws_f, ns_f)