diff --git a/news/ws.rst b/news/ws.rst
new file mode 100644
index 00000000..59526e7b
--- /dev/null
+++ b/news/ws.rst
@@ -0,0 +1,23 @@
+**Added:**
+
+* Added Whittaker-Shannon interpolation option for grid_interpolation.
+
+**Changed:**
+
+* Use WS interpolation for Nyquist grid.
+
+**Deprecated:**
+
+* <news item>
+
+**Removed:**
+
+* <news item>
+
+**Fixed:**
+
+* <news item>
+
+**Security:**
+
+* <news item>
diff --git a/src/diffpy/pdfgui/control/fitdataset.py b/src/diffpy/pdfgui/control/fitdataset.py
index 607b14c4..559e9dbd 100644
--- a/src/diffpy/pdfgui/control/fitdataset.py
+++ b/src/diffpy/pdfgui/control/fitdataset.py
@@ -22,6 +22,7 @@
 from diffpy.pdfgui.control.controlerrors import ControlStatusError
 from diffpy.pdfgui.control.parameter import Parameter
 from diffpy.pdfgui.control.pdfdataset import PDFDataSet
+from diffpy.utils.resampler import wsinterp
 
 
 class FitDataSet(PDFDataSet):
@@ -610,13 +611,14 @@ def _updateRcalcSampling(self):
         if frmax - (rcalcfirst + nrcalc * frstep) > frstep * 1e-8:
             nrcalc += 1
         newrcalc = rcalcfirst + frstep * numpy.arange(nrcalc + 1)
+        tp = self.getFitSamplingType()
         # Gcalc:
         if len(self._Gcalc) > 0:
-            newGcalc = grid_interpolation(self._rcalc, self._Gcalc, newrcalc)
+            newGcalc = grid_interpolation(self._rcalc, self._Gcalc, newrcalc, tp=tp)
             self._Gcalc = list(newGcalc)
         # dGcalc
         if len(self._dGcalc) > 0:
-            newdGcalc = grid_interpolation(self._rcalc, self._dGcalc, newrcalc)
+            newdGcalc = grid_interpolation(self._rcalc, self._dGcalc, newrcalc, tp=tp)
             self._dGcalc = list(newdGcalc)
         # invalidate Gtrunc and dGtrunc
         self._Gtrunc = []
@@ -709,7 +711,8 @@ def _set_dGcalc(self, value):
     def _get_Gtrunc(self):
         self._updateRcalcSampling()
         if not self._Gtrunc:
-            newGtrunc = grid_interpolation(self.robs, self.Gobs, self.rcalc)
+            tp = self.getFitSamplingType()
+            newGtrunc = grid_interpolation(self.robs, self.Gobs, self.rcalc, tp=tp)
             self._Gtrunc = list(newGtrunc)
         return self._Gtrunc
 
@@ -724,13 +727,15 @@ def _set_Gtrunc(self, value):
     def _get_dGtrunc(self):
         self._updateRcalcSampling()
         if not self._dGtrunc:
+            tp = self.getFitSamplingType()
             # use sum to avoid index error for empty arrays
             newdGtrunc = grid_interpolation(
                 self.robs,
                 self.dGobs,
                 self.rcalc,
-                youtleft=sum(self.dGobs[:1]),
-                youtright=sum(self.dGobs[-1:]),
+                left=sum(self.dGobs[:1]),
+                right=sum(self.dGobs[-1:]),
+                tp=tp,
             )
             self._dGtrunc = list(newdGtrunc)
         return self._dGtrunc
@@ -775,19 +780,8 @@ def _set_crw(self, value):
 ##############################################################################
 # helper functions
 ##############################################################################
+def _linear_interpolation(x0, y0, x1, youtleft, youtright):
 
-
-def grid_interpolation(x0, y0, x1, youtleft=0.0, youtright=0.0):
-    """Linear interpolation of x0, y0 values to a new grid x1.
-
-    x0       -- original x-grid, must be equally spaced
-    y0       -- original y values
-    x1       -- new x-grid, it can have any spacing
-    youtleft -- value for interpolated y1 for x1 below the x0 range
-    youtright -- value for interpolated y1 for x1 above the x0 range
-
-    Return numpy.array of interpolated y1 values.
-    """
     x0 = numpy.asarray(x0, copy=None, dtype=float)
     y0 = numpy.asarray(y0, copy=None, dtype=float)
     n0 = len(x0)
@@ -820,6 +814,47 @@ def grid_interpolation(x0, y0, x1, youtleft=0.0, youtright=0.0):
     return y1
 
 
+def grid_interpolation(x0, y0, x1, left=None, right=None, tp=None):
+    """Interpolate values from one grid onto another using either linear or
+    Whittaker–Shannon interpolation.
+
+    Parameters
+    ----------
+    x0 : array_like
+        Original x-grid, must be equally spaced.
+    y0 : array_like
+        Original values defined on x0.
+    x1 : array_like
+        New x-grid upon which to interpolate.
+    tp : {'data', 'Nyquist', 'custom'}, optional
+        Corresponding fit sampling type. Use Whittaker–Shannon interpolation
+        for Nyquist resampling and linear interpolation otherwise.
+        If not provided, linear interpolation is used.
+    left : float, optional
+        Value for interpolated y1 for x1 below the x0 range.
+        Default: if tp='Nyquist' then y1[0] is used. Otherwise 0.0 is used.
+    right : float, optional
+        Value for interpolated y1 for x1 above the x0 range.
+        Default: if tp='Nyquist' then y1[-1] is used. Otherwise 0.0 is used.
+
+    Returns
+    -------
+    numpy.ndarray
+        Array of interpolated values on the new grid x1.
+
+    Notes
+    -----
+    When tp='Nyquist', the function calls :func:`wsinterp` to perform Whittaker–Shannon interpolation.
+    Otherwise it uses the internal :func:`_linear_interpolation` routine.
+    """
+    if tp == "Nyquist":
+        return wsinterp(x1, x0, y0, left, right)
+    else:
+        left = 0.0 if left is None else left
+        right = 0.0 if right is None else right
+        return _linear_interpolation(x0, y0, x1, left, right)
+
+
 # simple test code
 if __name__ == "__main__":
     FitDataSet("name")
diff --git a/tests/test_fitdataset.py b/tests/test_fitdataset.py
index 7ec0d1ff..2e8c0655 100644
--- a/tests/test_fitdataset.py
+++ b/tests/test_fitdataset.py
@@ -39,7 +39,7 @@ def test_grid_interpolation(self):
         y0 = numpy.sin(x0)
         x1 = [-6, x0[0], -0.2, x0[-1], 37]
         y1a = fds.grid_interpolation(x0, y0, x1)
-        y1b = fds.grid_interpolation(x0, y0, x1, youtleft=637, youtright=638)
+        y1b = fds.grid_interpolation(x0, y0, x1, left=637, right=638)
         # outside values
         self.assertEqual(0.0, y1a[0])
         self.assertEqual(637, y1b[0])