Use DelayedSignal in wfs_25d_fir_prefilter

Author: Till Rettberg

doc/readthedocs-environment.yml

@@ -1,6 +1,6 @@
 name: SFS-readthedocs
 dependencies:
-  - python>=3
+  - python==3.5
   - sphinx>=1.3
   - sphinx_rtd_theme
   - numpy

examples/time_domain.py

@@ -15,18 +15,18 @@
 N = 56  # number of secondary sources
 R = 1.5  # radius of spherical/circular array
 x0, n0, a0 = sfs.array.circular(N, R)  # get secondary source positions
+fs = 44100  # sampling rate
 
-# create dirac
-signal = [1]
+# unit impulse
+signal = [1], fs
 
 # 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
+d = sfs.time.drivingfunction.driving_signals(d_delay, d_weight, signal)
+
 # test soundfield
 a = sfs.mono.drivingfunction.source_selection_point(n0, x0, xs)
 twin = sfs.tapering.tukey(a, .3)
@@ -48,9 +48,8 @@
 
 # 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
+d = sfs.time.drivingfunction.driving_signals(d_delay, d_weight, signal)
+
 # test soundfield
 a = sfs.mono.drivingfunction.source_selection_plane(n0, npw)
 twin = sfs.tapering.tukey(a, .3)
@@ -63,4 +62,27 @@
 sfs.plot.virtualsource_2d([0, 0], npw, type='plane')
 plt.title('impulse_pw_wfs_25d')
 plt.savefig('impulse_pw_wfs_25d.png')
+
+# FOCUSED SOURCE
+xs = np.r_[0.5, 0.5, 0]  # position of virtual source
+xref = np.r_[0, 0, 0]
+nfs = sfs.util.normalize_vector(xref - xs)  # main n of fsource
+t = 0.003  # compute driving signals
+d_delay, d_weight = sfs.time.drivingfunction.wfs_25d_focused(x0, n0, xs)
+d = sfs.time.drivingfunction.driving_signals(d_delay, d_weight, signal)
+
+# test soundfield
+a = sfs.mono.drivingfunction.source_selection_focused(nfs, x0, xs)
+twin = sfs.tapering.tukey(a, .3)
+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, n0, twin)
+plt.grid()
+sfs.plot.virtualsource_2d(xs)
+plt.title('impulse_fs_wfs_25d')
+plt.savefig('impulse_fs_wfs_25d.png')
+
 # plt.show()

sfs/defs.py

@@ -8,6 +8,3 @@
 
 # tolerance used for secondary source selection
 selection_tolerance = 1e-6
-
-# sampling frequency
-fs = 44100

sfs/mono/drivingfunction.py

@@ -86,7 +86,7 @@ def _wfs_plane(omega, x0, n0, n=[0, 1, 0], c=None):
     """
     x0 = util.asarray_of_rows(x0)
     n0 = util.asarray_of_rows(n0)
-    n = util.asarray_1d(n)
+    n = util.normalize_vector(n)
     k = util.wavenumber(omega, c)
     return 2j * k * np.inner(n, n0) * np.exp(-1j * k * np.inner(n, x0))
 
@@ -106,7 +106,7 @@ def wfs_25d_plane(omega, x0, n0, n=[0, 1, 0], xref=[0, 0, 0], c=None,
     """
     x0 = util.asarray_of_rows(x0)
     n0 = util.asarray_of_rows(n0)
-    n = util.asarray_1d(n)
+    n = util.normalize_vector(n)
     xref = util.asarray_1d(xref)
     k = util.wavenumber(omega, c)
     return wfs_25d_preeq(omega, omalias, c) * \
@@ -179,7 +179,7 @@ def wfs_25d_preeq(omega, omalias, c):
 def delay_3d_plane(omega, x0, n0, n=[0, 1, 0], c=None):
     """Plane wave by simple delay of secondary sources."""
     x0 = util.asarray_of_rows(x0)
-    n = util.asarray_1d(n)
+    n = util.normalize_vector(n)
     k = util.wavenumber(omega, c)
     return np.exp(-1j * k * np.inner(n, x0))
 
@@ -191,7 +191,7 @@ def source_selection_plane(n0, n):
 
     """
     n0 = util.asarray_of_rows(n0)
-    n = util.asarray_1d(n)
+    n = util.normalize_vector(n)
     return np.inner(n, n0) >= defs.selection_tolerance
 
 
@@ -225,7 +225,7 @@ def source_selection_focused(ns, x0, xs):
     """
     x0 = util.asarray_of_rows(x0)
     xs = util.asarray_1d(xs)
-    ns = util.asarray_1d(ns)
+    ns = util.normalize_vector(ns)
     ds = xs - x0
     return inner1d(ns, ds) >= defs.selection_tolerance
 
@@ -251,7 +251,7 @@ def nfchoa_2d_plane(omega, x0, r0, n=[0, 1, 0], max_order=None, c=None):
     """
     x0 = util.asarray_of_rows(x0)
     k = util.wavenumber(omega, c)
-    n = util.asarray_1d(n)
+    n = util.normalize_vector(n)
     phi, _, r = util.cart2sph(*n)
     phi0 = util.cart2sph(*x0.T)[0]
     M = _max_order_circular_harmonics(len(x0), max_order)
@@ -305,7 +305,7 @@ def nfchoa_25d_plane(omega, x0, r0, n=[0, 1, 0], max_order=None, c=None):
     """
     x0 = util.asarray_of_rows(x0)
     k = util.wavenumber(omega, c)
-    n = util.asarray_1d(n)
+    n = util.normalize_vector(n)
     phi, _, r = util.cart2sph(*n)
     phi0 = util.cart2sph(*x0.T)[0]
     M = _max_order_circular_harmonics(len(x0), max_order)
@@ -345,7 +345,7 @@ def sdm_2d_plane(omega, x0, n0, n=[0, 1, 0], c=None):
     """
     x0 = util.asarray_of_rows(x0)
     n0 = util.asarray_of_rows(n0)
-    n = util.asarray_1d(n)
+    n = util.normalize_vector(n)
     k = util.wavenumber(omega, c)
     return k * n[1] * np.exp(-1j * k * n[0] * x0[:, 0])
 
@@ -361,7 +361,7 @@ def sdm_25d_plane(omega, x0, n0, n=[0, 1, 0], xref=[0, 0, 0], c=None):
     """
     x0 = util.asarray_of_rows(x0)
     n0 = util.asarray_of_rows(n0)
-    n = util.asarray_1d(n)
+    n = util.normalize_vector(n)
     xref = util.asarray_1d(xref)
     k = util.wavenumber(omega, c)
     return 4j * np.exp(-1j*k*n[1]*xref[1]) / hankel2(0, k*n[1]*xref[1]) * \
@@ -420,6 +420,7 @@ def esa_edge_2d_plane(omega, x0, n=[0, 1, 0], alpha=3/2*np.pi, Nc=None, c=None):
 
     """
     x0 = np.asarray(x0)
+    n = util.normalize_vector(n)
     k = util.wavenumber(omega, c)
     phi_s = np.arctan2(n[1], n[0]) + np.pi
     L = x0.shape[0]
@@ -477,6 +478,7 @@ def esa_edge_dipole_2d_plane(omega, x0, n=[0, 1, 0], alpha=3/2*np.pi, Nc=None, c
 
     """
     x0 = np.asarray(x0)
+    n = util.normalize_vector(n)
     k = util.wavenumber(omega, c)
     phi_s = np.arctan2(n[1], n[0]) + np.pi
     L = x0.shape[0]

sfs/mono/source.py

@@ -489,7 +489,7 @@ def plane(omega, x0, n0, grid, c=None):
     """
     k = util.wavenumber(omega, c)
     x0 = util.asarray_1d(x0)
-    n0 = util.asarray_1d(n0)
+    n0 = util.normalize_vector(n0)
     grid = util.as_xyz_components(grid)
     return np.exp(-1j * k * np.inner(grid - x0, n0))