Add time based 2.5D WFS plane wave driving function.

Major improvements to time domain.

examples/time_domain.py

@@ -1,33 +1,66 @@
 """
 Create some examples in the time domain.
 
-Simulate and plot impulse behavior for Wafe Field Synthesis.
+Simulate and plot impulse behavior for Wave Field Synthesis.
 
 """
 
+import numpy as np
 import matplotlib.pyplot as plt
 import sfs
 
+# simulation parameters
 grid = sfs.util.xyz_grid([-3, 3], [-3, 3], 0, spacing=0.01)
 my_cmap = 'YlOrRd'
 N = 56  # number of secondary sources
 R = 1.5  # radius of spherical/circular array
-x0, nx0, a0 = sfs.array.circular(N, R)  # get secondary source positions
-xs = [2, 2, 0]
+x0, n0, a0 = sfs.array.circular(N, R)  # get secondary source positions
 
-# compute driving function
-d_delay, d_weight, d_line = sfs.time.drivingfunction.wfs_25d_ps(xs, x0, nx0)
+# create dirac
+signal = [1]
 
+# POINT SOURCE
+xs = [2, 2, 0]  # position of virtual source
+t = 0.008
+# compute driving signals
+d_delay, d_weight = sfs.time.drivingfunction.wfs_25d_point(x0, n0, xs)
+d, t_offset = sfs.time.drivingfunction.driving_signals(d_delay, d_weight,
+                                                       signal)
+t -= t_offset
 # test soundfield
-a = sfs.mono.drivingfunction.source_selection_point(nx0, x0, xs)
+a = sfs.mono.drivingfunction.source_selection_point(n0, x0, xs)
 twin = sfs.tapering.tukey(a, .3)
-p = sfs.time.soundfield.synthesize_p(d_line, twin * a0, x0, grid, 200)
+p = sfs.time.soundfield.p_array(x0, d, twin * a0, t, grid)
 p = p * 100  # scale absolute amplitude
 
 plt.figure(figsize=(10, 10))
 sfs.plot.level(p, grid, cmap=my_cmap)
-sfs.plot.loudspeaker_2d(x0, nx0, twin)
+sfs.plot.loudspeaker_2d(x0, n0, twin)
 plt.grid()
 sfs.plot.virtualsource_2d(xs)
-plt.title('impuls_ps_wfs_25d')
-plt.savefig('impuls_ps_wfs_25d' + '.png')
+plt.title('impulse_ps_wfs_25d')
+plt.savefig('impulse_ps_wfs_25d.png')
+
+# PLANE WAVE
+pw_angle = 30  # traveling direction of plane wave
+npw = sfs.util.direction_vector(np.radians(pw_angle))
+t = -0.001
+
+# compute driving signals
+d_delay, d_weight = sfs.time.drivingfunction.wfs_25d_plane(x0, n0, npw)
+d, t_offset = sfs.time.drivingfunction.driving_signals(d_delay,
+                                                       d_weight, signal)
+t -= t_offset
+# test soundfield
+a = sfs.mono.drivingfunction.source_selection_plane(n0, npw)
+twin = sfs.tapering.tukey(a, .3)
+p = sfs.time.soundfield.p_array(x0, d, twin * a0, t, grid)
+
+plt.figure(figsize=(10, 10))
+sfs.plot.level(p, grid, cmap=my_cmap)
+sfs.plot.loudspeaker_2d(x0, n0, twin)
+plt.grid()
+sfs.plot.virtualsource_2d([0, 0], npw, type='plane')
+plt.title('impulse_pw_wfs_25d')
+plt.savefig('impulse_pw_wfs_25d.png')
+# plt.show()

sfs/time/drivingfunction.py

@@ -1,64 +1,199 @@
-"""Compute time based driving functions for various systems."""
+"""Compute time based driving functions for various systems.
 
+.. include:: math-definitions.rst
+
+"""
+from __future__ import division
 import numpy as np
 from numpy.core.umath_tests import inner1d  # element-wise inner product
 from .. import defs
 from .. import util
 
 
-def wfs_25d_ps(xs, x0, nx0, xref=[0, 0, 0], fs=None, c=None):
-    """Point source by 2.5-dimensional WFS.
+def wfs_25d_plane(x0, n0, n=[0, 1, 0], xref=[0, 0, 0], c=None):
+    r"""Plane wave model by 2.5-dimensional WFS.
+
+    Parameters
+    ----------
+    x0 : (N, 3) array_like
+        Sequence of secondary source positions.
+    n0 : (N, 3) array_like
+        Sequence of secondary source orientations.
+    n : (3,) array_like, optional
+        Normal vector (propagation direction) of synthesized plane wave.
+    xref : (3,) array_like, optional
+        Reference position
+    c : float, optional
+        Speed of sound
+
+    Returns
+    -------
+    delays : (N,) numpy.ndarray
+        Delays of secondary sources in seconds.
+    weights: (N,) numpy.ndarray
+        Weights of secondary sources.
+
+    Notes
+    -----
+    2.5D correction factor
 
-    ::
+    .. math::
 
-         2.5D correction factor
-                ______________
-         g0 = \| 2pi |xref-x0|
+        g_0 = \sqrt{2 \pi |x_\mathrm{ref} - x_0|}
 
+    d using a plane wave as source model
 
-         d_2.5D using a point source as source model
+    .. math::
 
-                               g0  (x0-xs) nx0
-         d_2.5D(x0,t) = h(t) * --- ------------- delta(t-|x0-xs|/c)
-                               2pi |x0-xs|^(3/2)
+        d_{2.5D}(x_0,t) = h(t)
+        2 g_0 \scalarprod{n}{n_0}
+        \dirac{t - \frac{1}{c} \scalarprod{n}{x_0}}
 
-         see Wierstorf et al. (2015), eq.(#d:wfs:ps:2.5D)
+    with wfs(2.5D) prefilter h(t), which is not implemented yet.
+
+    References
+    ----------
+    See http://sfstoolbox.org/en/latest/#equation-d.wfs.pw.2.5D
 
     """
     if c is None:
         c = defs.c
-    if fs is None:
-        fs = defs.fs
     x0 = util.asarray_of_rows(x0)
-    nx0 = util.asarray_of_rows(nx0)
-    xs = util.asarray_1d(xs)
-    g0 = np.sqrt(2*np.pi*np.linalg.norm(xref-x0, axis=1))
-    ds = x0 - xs
-    r = np.linalg.norm(x0-xs, axis=1)
-    delay = 1/c * r
-    weight = g0/(2*np.pi) * inner1d(ds, nx0) / r**(3/2)
-    sig = np.ones(len(x0))  # dirac
-    line = delayline(sig, delay*fs, weight)
-    return delay, weight, line
+    n0 = util.asarray_of_rows(n0)
+    n = util.asarray_1d(n)
+    xref = util.asarray_1d(xref)
+    g0 = np.sqrt(2 * np.pi * np.linalg.norm(xref - x0, axis=1))
+    delays = inner1d(n, x0) / c
+    weights = 2 * g0 * inner1d(n, n0)
+    return delays, weights
+
+
+def wfs_25d_point(x0, n0, xs, xref=[0, 0, 0], c=None):
+    r"""Point source by 2.5-dimensional WFS.
+
+    Parameters
+    ----------
+    x0 : (N, 3) array_like
+        Sequence of secondary source positions.
+    n0 : (N, 3) array_like
+        Sequence of secondary source orientations.
+    xs : (3,) array_like
+        Virtual source position.
+    xref : (3,) array_like, optional
+        Reference position
+    c : float, optional
+        Speed of sound
+
+    Returns
+    -------
+    delays : (N,) numpy.ndarray
+        Delays of secondary sources in seconds.
+    weights: (N,) numpy.ndarray
+        Weights of secondary sources.
+
+    Notes
+    -----
+    2.5D correction factor
 
+    .. math::
 
-def delayline(sig, delay, weight):
-    """Delayline with weights, no fractional delay yet.
+         g_0 = \sqrt{2 \pi |x_\mathrm{ref} - x_0|}
 
-    Delay in samples.
-    sig (m x n) where columns n are channels and rows m the corresponding
-    signals.
-    Weight is array with weighting for each channel.
+
+    d using a point source as source model
+
+    .. math::
+
+         d_{2.5D}(x_0,t) = h(t)
+         \frac{g_0  \scalarprod{(x_0 - x_s)}{n_0}}
+         {2\pi |x_0 - x_s|^{3/2}}
+         \dirac{t - \frac{|x_0 - x_s|}{c}}
+
+    with wfs(2.5D) prefilter h(t), which is not implemented yet.
+
+    References
+    ----------
+    See http://sfstoolbox.org/en/latest/#equation-d.wfs.ps.2.5D
 
     """
-    delay = np.rint(delay).astype(int)  # no fractional delay, samples
-    extention = np.max(delay)  # avoid cutting signal
-    channels = len(weight)  # extract no. of channels
-    sig_length = sig.ndim
-    out = np.zeros([sig_length+extention, channels])  # create shape
+    if c is None:
+        c = defs.c
+    x0 = util.asarray_of_rows(x0)
+    n0 = util.asarray_of_rows(n0)
+    xs = util.asarray_1d(xs)
+    xref = util.asarray_1d(xref)
+    g0 = np.sqrt(2 * np.pi * np.linalg.norm(xref - x0, axis=1))
+    ds = x0 - xs
+    r = np.linalg.norm(ds, axis=1)
+    delays = r/c
+    weights = g0 * inner1d(ds, n0) / (2 * np.pi * r**(3/2))
+    return delays, weights
+
+
+def driving_signals(delays, weights, signal, fs=None):
+    """Get driving signals per secondary source.
+
+    Returned signals are the delayed and weighted mono input signal
+    (with N samples) per channel (C).
+
+    Parameters
+    ----------
+    delays : (C,) array_like
+        Delay in seconds for each channel, negative values allowed.
+    weights : (C,) array_like
+        Amplitude weighting factor for each channel.
+    signal : (N,) array_like
+        Excitation signal (mono) which gets weighted and delayed.
+    fs: int, optional
+        Sampling frequency in Hertz.
+
+    Returns
+    -------
+    driving_signals : (N, C) numpy.ndarray
+        Driving signal per channel (column represents channel).
+    t_offset : float
+        Simulation point in time offset (seconds).
 
-    for channel in range(channels):
-        cdelay = delay[channel]  # get channel delay
-        out[cdelay:cdelay+sig_length, channel] = sig[channel] * weight[channel]
+    """
+    delays = util.asarray_1d(delays)
+    weights = util.asarray_1d(weights)
+    d, t_offset = apply_delays(signal, delays, fs)
+    return d * weights, t_offset
+
+
+def apply_delays(signal, delays, fs=None):
+    """Apply delays for every channel.
+
+    A mono input signal gets delayed for each channel individually. The
+    simultation point in time is shifted by the smallest delay provided,
+    which allows negative delays as well.
+
+    Parameters
+    ----------
+    signal : (N,) array_like
+        Mono excitation signal (with N samples) which gets delayed.
+    delays : (C,) array_like
+        Delay in seconds for each channel (C), negative values allowed.
+    fs: int, optional
+        Sampling frequency in Hertz.
+
+    Returns
+    -------
+    out : (N, C) numpy.ndarray
+        Output signals (column represents channel).
+    t_offset : float
+        Simulation point in time offset (seconds).
 
-    return out
+    """
+    if fs is None:
+        fs = defs.fs
+    signal = util.asarray_1d(signal)
+    delays = util.asarray_1d(delays)
+
+    delays_samples = np.rint(fs * delays).astype(int)
+    offset_samples = delays_samples.min()
+    delays_samples -= offset_samples
+    out = np.zeros((delays_samples.max() + len(signal), len(delays_samples)))
+    for channel, cdelay in enumerate(delays_samples):
+        out[cdelay:cdelay + len(signal), channel] = signal
+    return out, offset_samples / fs