Change interpolation API, reimplement sinc interpolation, add time signal zeropadding

Author: Till Rettberg

sfs/td/source.py

@@ -30,7 +30,31 @@
 from .. import util as _util
 
 
-def point(xs, signal, observation_time, grid, c=None, interpolator_kind='linear'):
+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
@@ -50,6 +74,9 @@ def point(xs, signal, observation_time, grid, c=None, interpolator_kind='linear'
         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
     -------
@@ -86,46 +113,19 @@ def point(xs, signal, observation_time, grid, c=None, interpolator_kind='linear'
         weights = 1 / (4 * _np.pi * r)
     delays = r / c
     base_time = observation_time - signal_offset
-    if interpolator_kind == 'sinc':
-        p = _sinc_interp(data, _np.arange(len(data)),
-                         _np.array((base_time - delays) * samplerate))
-    else:
-        interpolator = _interp1d(_np.arange(len(data)), data,
-                                kind=interpolator_kind, bounds_error=False,
-                                fill_value=0)
-        p = interpolator((base_time - delays) * samplerate)
+    p = interpolator(_np.arange(len(data)), data)((base_time - delays) * samplerate)
     # weights can be +-infinity
     with _np.errstate(invalid='ignore'):
         return weights * p
 
 
-def _sinc_interp(x, s, u):
-    """
-    Ideal sinc interpolation of a signal
-    adapted from https://gist.github.com/endolith/1297227
-
-    Parameters
-    ----------
-    x : (N,) array_like
-        Signal to be interpolated.
-    s : (N,) array_like
-        Sampling instants of signal.
-    u : (N,) array_like
-        Sampling instants after interpolation.
-
-    Returns
-    -------
-    numpy.ndarray
-        Interpolated signal
-    """
-
-    # sampling period
-    T = s[1] - s[0]
-    # perform sinc interpolation
-    sincM = _np.tile(u, (len(s), 1)) - _np.tile(s[:, _np.newaxis], (1, len(u)))
-    y = _np.dot(x, _np.sinc(sincM/T))
+def sinc_interpolator(x, y):
+    x = _util.asarray_1d(x)
+    y = _util.asarray_1d(y)
 
-    return 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,