diff --git a/sfs/time/source.py b/sfs/time/source.py
index 6bce271f..f71e1bbe 100644
--- a/sfs/time/source.py
+++ b/sfs/time/source.py
@@ -8,11 +8,36 @@
 
 from __future__ import division
 import numpy as np
+from scipy.interpolate import interp1d
 from .. import util
 from .. import defs
 
 
-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
@@ -32,6 +57,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
     -------
@@ -58,6 +86,14 @@ def point(xs, signal, observation_time, grid, c=None):
     weights = 1 / (4 * np.pi * r)
     delays = r / c
     base_time = observation_time - signal_offset
-    return weights * np.interp(base_time - delays,
-                               np.arange(len(data)) / samplerate,
-                               data, left=0, right=0)
+    p_interpolant = interpolator(np.arange(len(data)), data)
+    return p_interpolant((base_time - delays) * samplerate) * weights
+
+
+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