diff --git a/sfs/td/source.py b/sfs/td/source.py
index 338dfdb7..af25673b 100644
--- a/sfs/td/source.py
+++ b/sfs/td/source.py
@@ -24,12 +24,37 @@
 
 """
 import numpy as _np
+from scipy.interpolate import interp1d as _interp1d
 
 from .. import default as _default
 from .. import util as _util
 
 
-def point(xs, signal, observation_time, grid, c=None):
+def linear_interpolator(x, y):
+    """1d linear interpolator with zero-padding.
+
+    Parameters
+    ----------
+    x : (N,) array_like
+        Sampling points.
+    y : (N,) array_like
+        Values at sampling points.
+
+    Returns
+    -------
+    function
+        Piecewise linear interpolant.
+
+    """
+    x = _util.asarray_1d(x)
+    y = _util.asarray_1d(y)
+    x = _np.concatenate([_np.array([min(x)-1]), x, _np.array([max(x)+1])])
+    y = _np.concatenate([_np.array([0]), y, _np.array([0])])
+    return _interp1d(x, y, bounds_error=False, fill_value=0)
+
+
+def point(xs, signal, observation_time, grid, c=None,
+          interpolator=linear_interpolator):
     r"""Source model for a point source: 3D Green's function.
 
     Calculates the scalar sound pressure field for a given point in
@@ -49,6 +74,9 @@ def point(xs, signal, observation_time, grid, c=None):
         See `sfs.util.xyz_grid()`.
     c : float, optional
         Speed of sound.
+    interpolator : function, optional
+        A function which constructs and returns a 1d interpolator.
+        see: linear_interpolator, sinc_interpolator
 
     Returns
     -------
@@ -75,6 +103,7 @@ def point(xs, signal, observation_time, grid, c=None):
     xs = _util.asarray_1d(xs)
     data, samplerate, signal_offset = _util.as_delayed_signal(signal)
     data = _util.asarray_1d(data)
+    observation_time = _util.asarray_1d(observation_time)
     grid = _util.as_xyz_components(grid)
     if c is None:
         c = _default.c
@@ -84,12 +113,19 @@ def point(xs, signal, observation_time, grid, c=None):
         weights = 1 / (4 * _np.pi * r)
     delays = r / c
     base_time = observation_time - signal_offset
-    points_at_time = _np.interp(base_time - delays,
-                               _np.arange(len(data)) / samplerate,
-                               data, left=0, right=0)
+    p = interpolator(_np.arange(len(data)), data)((base_time - delays) * samplerate)
     # weights can be +-infinity
     with _np.errstate(invalid='ignore'):
-        return weights * points_at_time
+        return weights * p
+
+
+def sinc_interpolator(x, y):
+    x = _util.asarray_1d(x)
+    y = _util.asarray_1d(y)
+
+    def f(xnew):
+        return sum([y[i] * _np.sinc(xnew - x[i]) for i in range(len(x))])
+    return f
 
 
 def point_image_sources(x0, signal, observation_time, grid, L, max_order,