From e16d1813aa59a695b25cea71ea1bbe6678c65b01 Mon Sep 17 00:00:00 2001 From: Frank Schultz Date: Tue, 12 Mar 2019 10:06:54 +0100 Subject: [PATCH 01/15] new mono drivingfunction handling --- doc/examples/sound_field_synthesis.py | 22 +- sfs/mono/__init__.py | 13 +- sfs/mono/esaedge.py | 400 ++++++++++++++++++ sfs/mono/nfchoa.py | 259 ++++++++++++ sfs/mono/sdm.py | 284 +++++++++++++ sfs/mono/soundfigure.py | 2 +- sfs/mono/wfs.py | 582 ++++++++++++++++++++++++++ 7 files changed, 1548 insertions(+), 14 deletions(-) create mode 100644 sfs/mono/esaedge.py create mode 100644 sfs/mono/nfchoa.py create mode 100644 sfs/mono/sdm.py create mode 100644 sfs/mono/wfs.py diff --git a/doc/examples/sound_field_synthesis.py b/doc/examples/sound_field_synthesis.py index 733ca2b..81cd8be 100644 --- a/doc/examples/sound_field_synthesis.py +++ b/doc/examples/sound_field_synthesis.py @@ -44,22 +44,22 @@ # === compute driving function and determine active secondary sources === -#d, selection, secondary_source = sfs.mono.drivingfunction.delay_3d_plane(omega, array.x, array.n, npw) +#d, selection, secondary_source = sfs.mono.wfs.plane_3d_delay(omega, array.x, array.n, npw) -#d, selection, secondary_source = sfs.mono.drivingfunction.wfs_2d_line(omega, array.x, array.n, xs) +#d, selection, secondary_source = sfs.mono.wfs.line_2d(omega, array.x, array.n, xs) -#d, selection, secondary_source = sfs.mono.drivingfunction.wfs_2d_plane(omega, array.x, array.n, npw) -d, selection, secondary_source = sfs.mono.drivingfunction.wfs_25d_plane(omega, array.x, array.n, npw, xref) -#d, selection, secondary_source = sfs.mono.drivingfunction.wfs_3d_plane(omega, array.x, array.n, npw) +#d, selection, secondary_source = sfs.mono.wfs.plane_2d(omega, array.x, array.n, npw) +d, selection, secondary_source = sfs.mono.wfs.plane_25d(omega, array.x, array.n, npw, xref) +#d, selection, secondary_source = sfs.mono.wfs.plane_3d(omega, array.x, array.n, npw) -#d, selection, secondary_source = sfs.mono.drivingfunction.wfs_2d_point(omega, array.x, array.n, xs) -#d, selection, secondary_source = sfs.mono.drivingfunction.wfs_25d_point(omega, array.x, array.n, xs) -#d, selection, secondary_source = sfs.mono.drivingfunction.wfs_3d_point(omega, array.x, array.n, xs) +#d, selection, secondary_source = sfs.mono.wfs.point_2d(omega, array.x, array.n, xs) +#d, selection, secondary_source = sfs.mono.wfs.point_25d(omega, array.x, array.n, xs) +#d, selection, secondary_source = sfs.mono.wfs.point_3d(omega, array.x, array.n, xs) -#d, selection, secondary_source = sfs.mono.drivingfunction.nfchoa_2d_plane(omega, array.x, R, npw) +#d, selection, secondary_source = sfs.mono.nfchoa.plane_2d(omega, array.x, R, npw) -#d, selection, secondary_source = sfs.mono.drivingfunction.nfchoa_25d_point(omega, array.x, R, xs) -#d, selection, secondary_source = sfs.mono.drivingfunction.nfchoa_25d_plane(omega, array.x, R, npw) +#d, selection, secondary_source = sfs.mono.nfchoa.point_25d(omega, array.x, R, xs) +#d, selection, secondary_source = sfs.mono.nfchoa.plane_25d(omega, array.x, R, npw) # === compute tapering window === diff --git a/sfs/mono/__init__.py b/sfs/mono/__init__.py index e33502e..0fa4f6f 100644 --- a/sfs/mono/__init__.py +++ b/sfs/mono/__init__.py @@ -3,14 +3,23 @@ .. autosummary:: :toctree: - drivingfunction + esaedge + nfchoa + sdm + wfs + source soundfigure """ import numpy as _np -from . import drivingfunction + +from . import esaedge +from . import nfchoa +from . import sdm +from . import wfs + from . import source from . import soundfigure diff --git a/sfs/mono/esaedge.py b/sfs/mono/esaedge.py new file mode 100644 index 0000000..7c84618 --- /dev/null +++ b/sfs/mono/esaedge.py @@ -0,0 +1,400 @@ +"""Compute driving functions for various systems. + +.. include:: math-definitions.rst + +.. plot:: + :context: reset + + import matplotlib.pyplot as plt + import numpy as np + import sfs + + plt.rcParams['figure.figsize'] = 6, 6 + + xs = -1.5, 1.5, 0 + xs_focused = -0.5, 0.5, 0 + # normal vector for plane wave: + npw = sfs.util.direction_vector(np.radians(-45)) + # normal vector for focused source: + ns_focused = sfs.util.direction_vector(np.radians(-45)) + f = 300 # Hz + omega = 2 * np.pi * f + R = 1.5 # Radius of circular loudspeaker array + + grid = sfs.util.xyz_grid([-2, 2], [-2, 2], 0, spacing=0.02) + + array = sfs.array.circular(N=32, R=R) + + def plot(d, selection, secondary_source): + p = sfs.mono.synthesize(d, selection, array, secondary_source, grid=grid) + sfs.plot.soundfield(p, grid) + sfs.plot.loudspeaker_2d(array.x, array.n, selection * array.a, size=0.15) + +""" + +import numpy as np +from numpy.core.umath_tests import inner1d # element-wise inner product +from scipy.special import jn, hankel2 +from .. import util +from .. import default +from . import source as _source + + +def plane_2d(omega, x0, n=[0, 1, 0], alpha=3/2*np.pi, Nc=None, + c=None): + r"""Driving function for 2-dimensional plane wave with edge ESA. + + Driving function for a virtual plane wave using the 2-dimensional ESA + for an edge-shaped secondary source distribution consisting of + monopole line sources. + + Parameters + ---------- + omega : float + Angular frequency. + x0 : int(N, 3) array_like + Sequence of secondary source positions. + n : (3,) array_like, optional + Normal vector of synthesized plane wave. + alpha : float, optional + Outer angle of edge. + Nc : int, optional + Number of elements for series expansion of driving function. Estimated + if not given. + c : float, optional + Speed of sound + + Returns + ------- + d : (N,) numpy.ndarray + Complex weights of secondary sources. + selection : (N,) numpy.ndarray + Boolean array containing ``True`` or ``False`` depending on + whether the corresponding secondary source is "active" or not. + secondary_source_function : callable + A function that can be used to create the sound field of a + single secondary source. See `sfs.mono.synthesize()`. + + Notes + ----- + One leg of the secondary sources has to be located on the x-axis (y0=0), + the edge at the origin. + + Derived from :cite:`Spors2016` + + """ + 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] + + r = np.linalg.norm(x0, axis=1) + phi = np.arctan2(x0[:, 1], x0[:, 0]) + phi = np.where(phi < 0, phi+2*np.pi, phi) + + if Nc is None: + Nc = np.ceil(2 * k * np.max(r) * alpha/np.pi) + + epsilon = np.ones(Nc) # weights for series expansion + epsilon[0] = 2 + + d = np.zeros(L, dtype=complex) + for m in np.arange(Nc): + nu = m*np.pi/alpha + d = d + 1/epsilon[m] * np.exp(1j*nu*np.pi/2) * np.sin(nu*phi_s) \ + * np.cos(nu*phi) * nu/r * jn(nu, k*r) + + d[phi > 0] = -d[phi > 0] + + selection = util.source_selection_all(len(x0)) + return 4*np.pi/alpha * d, selection, secondary_source_line(omega, c) + + +def plane_2d_with_dipole_ssd(omega, x0, n=[0, 1, 0], alpha=3/2*np.pi, Nc=None, + c=None): + r"""Driving function for 2-dimensional plane wave with edge dipole ESA. + + Driving function for a virtual plane wave using the 2-dimensional ESA + for an edge-shaped secondary source distribution consisting of + dipole line sources. + + Parameters + ---------- + omega : float + Angular frequency. + x0 : int(N, 3) array_like + Sequence of secondary source positions. + n : (3,) array_like, optional + Normal vector of synthesized plane wave. + alpha : float, optional + Outer angle of edge. + Nc : int, optional + Number of elements for series expansion of driving function. Estimated + if not given. + c : float, optional + Speed of sound + + Returns + ------- + d : (N,) numpy.ndarray + Complex weights of secondary sources. + selection : (N,) numpy.ndarray + Boolean array containing ``True`` or ``False`` depending on + whether the corresponding secondary source is "active" or not. + secondary_source_function : callable + A function that can be used to create the sound field of a + single secondary source. See `sfs.mono.synthesize()`. + + Notes + ----- + One leg of the secondary sources has to be located on the x-axis (y0=0), + the edge at the origin. + + Derived from :cite:`Spors2016` + + """ + 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] + + r = np.linalg.norm(x0, axis=1) + phi = np.arctan2(x0[:, 1], x0[:, 0]) + phi = np.where(phi < 0, phi+2*np.pi, phi) + + if Nc is None: + Nc = np.ceil(2 * k * np.max(r) * alpha/np.pi) + + epsilon = np.ones(Nc) # weights for series expansion + epsilon[0] = 2 + + d = np.zeros(L, dtype=complex) + for m in np.arange(Nc): + nu = m*np.pi/alpha + d = d + 1/epsilon[m] * np.exp(1j*nu*np.pi/2) * np.cos(nu*phi_s) \ + * np.cos(nu*phi) * jn(nu, k*r) + + return 4*np.pi/alpha * d + + +def line_2d(omega, x0, xs, alpha=3/2*np.pi, Nc=None, c=None): + r"""Driving function for 2-dimensional line source with edge ESA. + + Driving function for a virtual line source using the 2-dimensional ESA + for an edge-shaped secondary source distribution consisting of line + sources. + + Parameters + ---------- + omega : float + Angular frequency. + x0 : int(N, 3) array_like + Sequence of secondary source positions. + xs : (3,) array_like + Position of synthesized line source. + alpha : float, optional + Outer angle of edge. + Nc : int, optional + Number of elements for series expansion of driving function. Estimated + if not given. + c : float, optional + Speed of sound + + Returns + ------- + d : (N,) numpy.ndarray + Complex weights of secondary sources. + selection : (N,) numpy.ndarray + Boolean array containing ``True`` or ``False`` depending on + whether the corresponding secondary source is "active" or not. + secondary_source_function : callable + A function that can be used to create the sound field of a + single secondary source. See `sfs.mono.synthesize()`. + + Notes + ----- + One leg of the secondary sources has to be located on the x-axis (y0=0), + the edge at the origin. + + Derived from :cite:`Spors2016` + + """ + x0 = np.asarray(x0) + k = util.wavenumber(omega, c) + phi_s = np.arctan2(xs[1], xs[0]) + if phi_s < 0: + phi_s = phi_s + 2*np.pi + r_s = np.linalg.norm(xs) + L = x0.shape[0] + + r = np.linalg.norm(x0, axis=1) + phi = np.arctan2(x0[:, 1], x0[:, 0]) + phi = np.where(phi < 0, phi+2*np.pi, phi) + + if Nc is None: + Nc = np.ceil(2 * k * np.max(r) * alpha/np.pi) + + epsilon = np.ones(Nc) # weights for series expansion + epsilon[0] = 2 + + d = np.zeros(L, dtype=complex) + idx = (r <= r_s) + for m in np.arange(Nc): + nu = m*np.pi/alpha + f = 1/epsilon[m] * np.sin(nu*phi_s) * np.cos(nu*phi) * nu/r + d[idx] = d[idx] + f[idx] * jn(nu, k*r[idx]) * hankel2(nu, k*r_s) + d[~idx] = d[~idx] + f[~idx] * jn(nu, k*r_s) * hankel2(nu, k*r[~idx]) + + d[phi > 0] = -d[phi > 0] + + selection = util.source_selection_all(len(x0)) + return -1j*np.pi/alpha * d, selection, secondary_source_line(omega, c) + + +def line_2d_with_dipole_ssd(omega, x0, xs, alpha=3/2*np.pi, Nc=None, c=None): + r"""Driving function for 2-dimensional line source with edge dipole ESA. + + Driving function for a virtual line source using the 2-dimensional ESA + for an edge-shaped secondary source distribution consisting of dipole line + sources. + + Parameters + ---------- + omega : float + Angular frequency. + x0 : (N, 3) array_like + Sequence of secondary source positions. + xs : (3,) array_like + Position of synthesized line source. + alpha : float, optional + Outer angle of edge. + Nc : int, optional + Number of elements for series expansion of driving function. Estimated + if not given. + c : float, optional + Speed of sound + + Returns + ------- + d : (N,) numpy.ndarray + Complex weights of secondary sources. + selection : (N,) numpy.ndarray + Boolean array containing ``True`` or ``False`` depending on + whether the corresponding secondary source is "active" or not. + secondary_source_function : callable + A function that can be used to create the sound field of a + single secondary source. See `sfs.mono.synthesize()`. + + Notes + ----- + One leg of the secondary sources has to be located on the x-axis (y0=0), + the edge at the origin. + + Derived from :cite:`Spors2016` + + """ + x0 = np.asarray(x0) + k = util.wavenumber(omega, c) + phi_s = np.arctan2(xs[1], xs[0]) + if phi_s < 0: + phi_s = phi_s + 2*np.pi + r_s = np.linalg.norm(xs) + L = x0.shape[0] + + r = np.linalg.norm(x0, axis=1) + phi = np.arctan2(x0[:, 1], x0[:, 0]) + phi = np.where(phi < 0, phi+2*np.pi, phi) + + if Nc is None: + Nc = np.ceil(2 * k * np.max(r) * alpha/np.pi) + + epsilon = np.ones(Nc) # weights for series expansion + epsilon[0] = 2 + + d = np.zeros(L, dtype=complex) + idx = (r <= r_s) + for m in np.arange(Nc): + nu = m*np.pi/alpha + f = 1/epsilon[m] * np.cos(nu*phi_s) * np.cos(nu*phi) + d[idx] = d[idx] + f[idx] * jn(nu, k*r[idx]) * hankel2(nu, k*r_s) + d[~idx] = d[~idx] + f[~idx] * jn(nu, k*r_s) * hankel2(nu, k*r[~idx]) + + return -1j*np.pi/alpha * d + + +def point_25d(omega, x0, xs, xref=[2, -2, 0], alpha=3/2*np.pi, + Nc=None, c=None): + r"""Driving function for 2.5-dimensional point source with edge ESA. + + Driving function for a virtual point source using the 2.5-dimensional + ESA for an edge-shaped secondary source distribution consisting of point + sources. + + Parameters + ---------- + omega : float + Angular frequency. + x0 : int(N, 3) array_like + Sequence of secondary source positions. + xs : (3,) array_like + Position of synthesized line source. + xref: (3,) array_like or float + Reference position or reference distance + alpha : float, optional + Outer angle of edge. + Nc : int, optional + Number of elements for series expansion of driving function. Estimated + if not given. + c : float, optional + Speed of sound + + Returns + ------- + d : (N,) numpy.ndarray + Complex weights of secondary sources. + selection : (N,) numpy.ndarray + Boolean array containing ``True`` or ``False`` depending on + whether the corresponding secondary source is "active" or not. + secondary_source_function : callable + A function that can be used to create the sound field of a + single secondary source. See `sfs.mono.synthesize()`. + + Notes + ----- + One leg of the secondary sources has to be located on the x-axis (y0=0), + the edge at the origin. + + Derived from :cite:`Spors2016` + + """ + x0 = np.asarray(x0) + xs = np.asarray(xs) + xref = np.asarray(xref) + + if np.isscalar(xref): + a = np.linalg.norm(xref)/np.linalg.norm(xref-xs) + else: + a = np.linalg.norm(xref-x0, axis=1)/np.linalg.norm(xref-xs) + + d, selection, _ = line_2d(omega, x0, xs, alpha=alpha, Nc=Nc, c=c) + return 1j*np.sqrt(a) * d, selection, secondary_source_point(omega, c) + + +def secondary_source_point(omega, c): + """Create a point source for use in `sfs.mono.synthesize()`.""" + + def secondary_source(position, _, grid): + return _source.point(omega, position, grid, c) + + return secondary_source + + +def secondary_source_line(omega, c): + """Create a line source for use in `sfs.mono.synthesize()`.""" + + def secondary_source(position, _, grid): + return _source.line(omega, position, grid, c) + + return secondary_source diff --git a/sfs/mono/nfchoa.py b/sfs/mono/nfchoa.py new file mode 100644 index 0000000..950eae9 --- /dev/null +++ b/sfs/mono/nfchoa.py @@ -0,0 +1,259 @@ +"""Compute driving functions for various systems. + +.. include:: math-definitions.rst + +.. plot:: + :context: reset + + import matplotlib.pyplot as plt + import numpy as np + import sfs + + plt.rcParams['figure.figsize'] = 6, 6 + + xs = -1.5, 1.5, 0 + xs_focused = -0.5, 0.5, 0 + # normal vector for plane wave: + npw = sfs.util.direction_vector(np.radians(-45)) + # normal vector for focused source: + ns_focused = sfs.util.direction_vector(np.radians(-45)) + f = 300 # Hz + omega = 2 * np.pi * f + R = 1.5 # Radius of circular loudspeaker array + + grid = sfs.util.xyz_grid([-2, 2], [-2, 2], 0, spacing=0.02) + + array = sfs.array.circular(N=32, R=R) + + def plot(d, selection, secondary_source): + p = sfs.mono.synthesize(d, selection, array, secondary_source, grid=grid) + sfs.plot.soundfield(p, grid) + sfs.plot.loudspeaker_2d(array.x, array.n, selection * array.a, size=0.15) + +""" + +import numpy as np +from numpy.core.umath_tests import inner1d # element-wise inner product +from scipy.special import jn, hankel2 +from .. import util +from .. import default +from . import source as _source + + +def plane_2d(omega, x0, r0, n=[0, 1, 0], max_order=None, c=None): + r"""Driving function for 2-dimensional NFC-HOA for a virtual plane wave. + + Parameters + ---------- + omega : float + Angular frequency of plane wave. + x0 : (N, 3) array_like + Sequence of secondary source positions. + r0 : float + Radius of circular secondary source distribution. + n : (3,) array_like, optional + Normal vector (traveling direction) of plane wave. + max_order : float, optional + Maximum order of circular harmonics used for the calculation. + c : float, optional + Speed of sound. + + Returns + ------- + d : (N,) numpy.ndarray + Complex weights of secondary sources. + selection : (N,) numpy.ndarray + Boolean array containing only ``True`` indicating that + all secondary source are "active" for NFC-HOA. + secondary_source_function : callable + A function that can be used to create the sound field of a + single secondary source. See `sfs.mono.synthesize()`. + + Notes + ----- + .. math:: + + D(\phi_0, \omega) = + -\frac{2\i}{\pi r_0} + \sum_{m=-M}^M + \frac{\i^{-m}}{\Hankel{2}{m}{\wc r_0}} + \e{\i m (\phi_0 - \phi_\text{pw})} + + See http://sfstoolbox.org/#equation-D.nfchoa.pw.2D. + + Examples + -------- + .. plot:: + :context: close-figs + + d, selection, secondary_source = sfs.mono.nfchoa.plane_2d( + omega, array.x, R, npw) + plot(d, selection, secondary_source) + + """ + if max_order is None: + max_order = util.max_order_circular_harmonics(len(x0)) + + x0 = util.asarray_of_rows(x0) + k = util.wavenumber(omega, c) + n = util.normalize_vector(n) + phi, _, r = util.cart2sph(*n) + phi0 = util.cart2sph(*x0.T)[0] + d = 0 + for m in range(-max_order, max_order + 1): + d += 1j**-m / hankel2(m, k * r0) * np.exp(1j * m * (phi0 - phi)) + selection = util.source_selection_all(len(x0)) + return -2j / (np.pi*r0) * d, selection, secondary_source_point(omega, c) + + +def point_25d(omega, x0, r0, xs, max_order=None, c=None): + r"""Driving function for 2.5-dimensional NFC-HOA for a virtual point source. + + Parameters + ---------- + omega : float + Angular frequency of point source. + x0 : (N, 3) array_like + Sequence of secondary source positions. + r0 : float + Radius of circular secondary source distribution. + xs : (3,) array_like + Position of point source. + max_order : float, optional + Maximum order of circular harmonics used for the calculation. + c : float, optional + Speed of sound. + + Returns + ------- + d : (N,) numpy.ndarray + Complex weights of secondary sources. + selection : (N,) numpy.ndarray + Boolean array containing only ``True`` indicating that + all secondary source are "active" for NFC-HOA. + secondary_source_function : callable + A function that can be used to create the sound field of a + single secondary source. See `sfs.mono.synthesize()`. + + Notes + ----- + .. math:: + + D(\phi_0, \omega) = + \frac{1}{2 \pi r_0} + \sum_{m=-M}^M + \frac{\hankel{2}{|m|}{\wc r}}{\hankel{2}{|m|}{\wc r_0}} + \e{\i m (\phi_0 - \phi)} + + See http://sfstoolbox.org/#equation-D.nfchoa.ps.2.5D. + + Examples + -------- + .. plot:: + :context: close-figs + + d, selection, secondary_source = sfs.mono.nfchoa.point_25d( + omega, array.x, R, xs) + plot(d, selection, secondary_source) + + """ + if max_order is None: + max_order = util.max_order_circular_harmonics(len(x0)) + + x0 = util.asarray_of_rows(x0) + k = util.wavenumber(omega, c) + xs = util.asarray_1d(xs) + phi, _, r = util.cart2sph(*xs) + phi0 = util.cart2sph(*x0.T)[0] + hr = util.spherical_hn2(range(0, max_order + 1), k * r) + hr0 = util.spherical_hn2(range(0, max_order + 1), k * r0) + d = 0 + for m in range(-max_order, max_order + 1): + d += hr[abs(m)] / hr0[abs(m)] * np.exp(1j * m * (phi0 - phi)) + selection = util.source_selection_all(len(x0)) + return d / (2 * np.pi * r0), selection, secondary_source_point(omega, c) + + +def plane_25d(omega, x0, r0, n=[0, 1, 0], max_order=None, c=None): + r"""Driving function for 2.5-dimensional NFC-HOA for a virtual plane wave. + + Parameters + ---------- + omega : float + Angular frequency of point source. + x0 : (N, 3) array_like + Sequence of secondary source positions. + r0 : float + Radius of circular secondary source distribution. + n : (3,) array_like, optional + Normal vector (traveling direction) of plane wave. + max_order : float, optional + Maximum order of circular harmonics used for the calculation. + c : float, optional + Speed of sound. + + Returns + ------- + d : (N,) numpy.ndarray + Complex weights of secondary sources. + selection : (N,) numpy.ndarray + Boolean array containing only ``True`` indicating that + all secondary source are "active" for NFC-HOA. + secondary_source_function : callable + A function that can be used to create the sound field of a + single secondary source. See `sfs.mono.synthesize()`. + + Notes + ----- + .. math:: + + D(\phi_0, \omega) = + \frac{2\i}{r_0} + \sum_{m=-M}^M + \frac{\i^{-|m|}}{\wc \hankel{2}{|m|}{\wc r_0}} + \e{\i m (\phi_0 - \phi_\text{pw})} + + See http://sfstoolbox.org/#equation-D.nfchoa.pw.2.5D. + + Examples + -------- + .. plot:: + :context: close-figs + + d, selection, secondary_source = sfs.mono.nfchoa.plane_25d( + omega, array.x, R, npw) + plot(d, selection, secondary_source) + + """ + if max_order is None: + max_order = util.max_order_circular_harmonics(len(x0)) + + x0 = util.asarray_of_rows(x0) + k = util.wavenumber(omega, c) + n = util.normalize_vector(n) + phi, _, r = util.cart2sph(*n) + phi0 = util.cart2sph(*x0.T)[0] + d = 0 + hn2 = util.spherical_hn2(range(0, max_order + 1), k * r0) + for m in range(-max_order, max_order + 1): + d += (-1j)**abs(m) / (k * hn2[abs(m)]) * np.exp(1j * m * (phi0 - phi)) + selection = util.source_selection_all(len(x0)) + return 2*1j / r0 * d, selection, secondary_source_point(omega, c) + + +def secondary_source_point(omega, c): + """Create a point source for use in `sfs.mono.synthesize()`.""" + + def secondary_source(position, _, grid): + return _source.point(omega, position, grid, c) + + return secondary_source + + +def secondary_source_line(omega, c): + """Create a line source for use in `sfs.mono.synthesize()`.""" + + def secondary_source(position, _, grid): + return _source.line(omega, position, grid, c) + + return secondary_source diff --git a/sfs/mono/sdm.py b/sfs/mono/sdm.py new file mode 100644 index 0000000..50e3879 --- /dev/null +++ b/sfs/mono/sdm.py @@ -0,0 +1,284 @@ +"""Compute driving functions for various systems. + +.. include:: math-definitions.rst + +.. plot:: + :context: reset + + import matplotlib.pyplot as plt + import numpy as np + import sfs + + plt.rcParams['figure.figsize'] = 6, 6 + + xs = -1.5, 1.5, 0 + xs_focused = -0.5, 0.5, 0 + # normal vector for plane wave: + npw = sfs.util.direction_vector(np.radians(-45)) + # normal vector for focused source: + ns_focused = sfs.util.direction_vector(np.radians(-45)) + f = 300 # Hz + omega = 2 * np.pi * f + R = 1.5 # Radius of circular loudspeaker array + + grid = sfs.util.xyz_grid([-2, 2], [-2, 2], 0, spacing=0.02) + + array = sfs.array.circular(N=32, R=R) + + def plot(d, selection, secondary_source): + p = sfs.mono.synthesize(d, selection, array, secondary_source, grid=grid) + sfs.plot.soundfield(p, grid) + sfs.plot.loudspeaker_2d(array.x, array.n, selection * array.a, size=0.15) + +""" + +import numpy as np +from numpy.core.umath_tests import inner1d # element-wise inner product +from scipy.special import jn, hankel2 +from .. import util +from .. import default +from . import source as _source + + +def line_2d(omega, x0, n0, xs, c=None): + r"""Driving function for 2-dimensional SDM for a virtual line source. + + Parameters + ---------- + omega : float + Angular frequency of line source. + x0 : (N, 3) array_like + Sequence of secondary source positions. + n0 : (N, 3) array_like + Sequence of normal vectors of secondary sources. + xs : (3,) array_like + Position of line source. + c : float, optional + Speed of sound. + + Returns + ------- + d : (N,) numpy.ndarray + Complex weights of secondary sources. + selection : (N,) numpy.ndarray + Boolean array containing ``True`` or ``False`` depending on + whether the corresponding secondary source is "active" or not. + secondary_source_function : callable + A function that can be used to create the sound field of a + single secondary source. See `sfs.mono.synthesize()`. + + Notes + ----- + The secondary sources have to be located on the x-axis (y0=0). + Derived from :cite:`Spors2009`, Eq.(9), Eq.(4). + + Examples + -------- + .. plot:: + :context: close-figs + + array = sfs.array.linear(32, 0.2, orientation=[0, -1, 0]) + d, selection, secondary_source = sfs.mono.sdm.line_2d( + omega, array.x, array.n, xs) + plot(d, selection, secondary_source) + + """ + x0 = util.asarray_of_rows(x0) + n0 = util.asarray_of_rows(n0) + xs = util.asarray_1d(xs) + k = util.wavenumber(omega, c) + ds = x0 - xs + r = np.linalg.norm(ds, axis=1) + d = - 1j/2 * k * xs[1] / r * hankel2(1, k * r) + selection = util.source_selection_all(len(x0)) + return d, selection, secondary_source_line(omega, c) + + +def plane_2d(omega, x0, n0, n=[0, 1, 0], c=None): + r"""Driving function for 2-dimensional SDM for a virtual plane wave. + + Parameters + ---------- + omega : float + Angular frequency of plane wave. + x0 : (N, 3) array_like + Sequence of secondary source positions. + n0 : (N, 3) array_like + Sequence of normal vectors of secondary sources. + n: (3,) array_like, optional + Normal vector (traveling direction) of plane wave. + c : float, optional + Speed of sound. + + Returns + ------- + d : (N,) numpy.ndarray + Complex weights of secondary sources. + selection : (N,) numpy.ndarray + Boolean array containing ``True`` or ``False`` depending on + whether the corresponding secondary source is "active" or not. + secondary_source_function : callable + A function that can be used to create the sound field of a + single secondary source. See `sfs.mono.synthesize()`. + + Notes + ----- + The secondary sources have to be located on the x-axis (y0=0). + Derived from :cite:`Ahrens2012`, Eq.(3.73), Eq.(C.5), Eq.(C.11): + + .. math:: + + D(\x_0,k) = k_\text{pw,y} \e{-\i k_\text{pw,x} x} + + Examples + -------- + .. plot:: + :context: close-figs + + array = sfs.array.linear(32, 0.2, orientation=[0, -1, 0]) + d, selection, secondary_source = sfs.mono.sdm.plane_2d( + omega, array.x, array.n, npw) + plot(d, selection, secondary_source) + + """ + x0 = util.asarray_of_rows(x0) + n0 = util.asarray_of_rows(n0) + n = util.normalize_vector(n) + k = util.wavenumber(omega, c) + d = k * n[1] * np.exp(-1j * k * n[0] * x0[:, 0]) + selection = util.source_selection_all(len(x0)) + return d, selection, secondary_source_line(omega, c) + + +def plane_25d(omega, x0, n0, n=[0, 1, 0], xref=[0, 0, 0], c=None): + r"""Driving function for 2.5-dimensional SDM for a virtual plane wave. + + Parameters + ---------- + omega : float + Angular frequency of plane wave. + x0 : (N, 3) array_like + Sequence of secondary source positions. + n0 : (N, 3) array_like + Sequence of normal vectors of secondary sources. + n: (3,) array_like, optional + Normal vector (traveling direction) of plane wave. + xref : (3,) array_like, optional + Reference point for synthesized sound field. + c : float, optional + Speed of sound. + + Returns + ------- + d : (N,) numpy.ndarray + Complex weights of secondary sources. + selection : (N,) numpy.ndarray + Boolean array containing ``True`` or ``False`` depending on + whether the corresponding secondary source is "active" or not. + secondary_source_function : callable + A function that can be used to create the sound field of a + single secondary source. See `sfs.mono.synthesize()`. + + Notes + ----- + The secondary sources have to be located on the x-axis (y0=0). + Eq.(3.79) from :cite:`Ahrens2012`. + + Examples + -------- + .. plot:: + :context: close-figs + + array = sfs.array.linear(32, 0.2, orientation=[0, -1, 0]) + d, selection, secondary_source = ( + sfs.mono.sdm._25d( + omega, array.x, array.n, npw, [0, -1, 0])) + plot(d, selection, secondary_source) + + """ + x0 = util.asarray_of_rows(x0) + n0 = util.asarray_of_rows(n0) + n = util.normalize_vector(n) + xref = util.asarray_1d(xref) + k = util.wavenumber(omega, c) + d = 4j * np.exp(-1j*k*n[1]*xref[1]) / hankel2(0, k*n[1]*xref[1]) * \ + np.exp(-1j*k*n[0]*x0[:, 0]) + selection = util.source_selection_all(len(x0)) + return d, selection, secondary_source_point(omega, c) + + +def point_25d(omega, x0, n0, xs, xref=[0, 0, 0], c=None): + r"""Driving function for 2.5-dimensional SDM for a virtual point source. + + Parameters + ---------- + omega : float + Angular frequency of point source. + x0 : (N, 3) array_like + Sequence of secondary source positions. + n0 : (N, 3) array_like + Sequence of normal vectors of secondary sources. + xs: (3,) array_like + Position of virtual point source. + xref : (3,) array_like, optional + Reference point for synthesized sound field. + c : float, optional + Speed of sound. + + Returns + ------- + d : (N,) numpy.ndarray + Complex weights of secondary sources. + selection : (N,) numpy.ndarray + Boolean array containing ``True`` or ``False`` depending on + whether the corresponding secondary source is "active" or not. + secondary_source_function : callable + A function that can be used to create the sound field of a + single secondary source. See `sfs.mono.synthesize()`. + + Notes + ----- + The secondary sources have to be located on the x-axis (y0=0). + Driving function from :cite:`Spors2010`, Eq.(24). + + Examples + -------- + .. plot:: + :context: close-figs + + array = sfs.array.linear(32, 0.2, orientation=[0, -1, 0]) + d, selection, secondary_source = ( + sfs.mono.sdm.point_25d( + omega, array.x, array.n, xs, [0, -1, 0])) + plot(d, selection, secondary_source) + + """ + x0 = util.asarray_of_rows(x0) + n0 = util.asarray_of_rows(n0) + xs = util.asarray_1d(xs) + xref = util.asarray_1d(xref) + k = util.wavenumber(omega, c) + ds = x0 - xs + r = np.linalg.norm(ds, axis=1) + d = 1/2 * 1j * k * np.sqrt(xref[1] / (xref[1] - xs[1])) * \ + xs[1] / r * hankel2(1, k * r) + selection = util.source_selection_all(len(x0)) + return d, selection, secondary_source_point(omega, c) + + +def secondary_source_point(omega, c): + """Create a point source for use in `sfs.mono.synthesize()`.""" + + def secondary_source(position, _, grid): + return _source.point(omega, position, grid, c) + + return secondary_source + + +def secondary_source_line(omega, c): + """Create a line source for use in `sfs.mono.synthesize()`.""" + + def secondary_source(position, _, grid): + return _source.line(omega, position, grid, c) + + return secondary_source diff --git a/sfs/mono/soundfigure.py b/sfs/mono/soundfigure.py index fcaf4f5..0c0cfb9 100644 --- a/sfs/mono/soundfigure.py +++ b/sfs/mono/soundfigure.py @@ -2,7 +2,7 @@ import numpy as np from .. import util -from .drivingfunction import wfs_3d_plane +from .wfs import plane_3d def wfs_3d_pw(omega, x0, n0, figure, npw=[0, 0, 1], c=None): diff --git a/sfs/mono/wfs.py b/sfs/mono/wfs.py new file mode 100644 index 0000000..d54f157 --- /dev/null +++ b/sfs/mono/wfs.py @@ -0,0 +1,582 @@ +"""Compute driving functions for various systems. + +.. include:: math-definitions.rst + +.. plot:: + :context: reset + + import matplotlib.pyplot as plt + import numpy as np + import sfs + + plt.rcParams['figure.figsize'] = 6, 6 + + xs = -1.5, 1.5, 0 + xs_focused = -0.5, 0.5, 0 + # normal vector for plane wave: + npw = sfs.util.direction_vector(np.radians(-45)) + # normal vector for focused source: + ns_focused = sfs.util.direction_vector(np.radians(-45)) + f = 300 # Hz + omega = 2 * np.pi * f + R = 1.5 # Radius of circular loudspeaker array + + grid = sfs.util.xyz_grid([-2, 2], [-2, 2], 0, spacing=0.02) + + array = sfs.array.circular(N=32, R=R) + + def plot(d, selection, secondary_source): + p = sfs.mono.synthesize(d, selection, array, secondary_source, grid=grid) + sfs.plot.soundfield(p, grid) + sfs.plot.loudspeaker_2d(array.x, array.n, selection * array.a, size=0.15) + +""" + +import numpy as np +from numpy.core.umath_tests import inner1d # element-wise inner product +from scipy.special import jn, hankel2 +from .. import util +from .. import default +from . import source as _source + + +def line_2d(omega, x0, n0, xs, c=None): + r"""Driving function for 2-dimensional WFS for a virtual line source. + + Parameters + ---------- + omega : float + Angular frequency of line source. + x0 : (N, 3) array_like + Sequence of secondary source positions. + n0 : (N, 3) array_like + Sequence of normal vectors of secondary sources. + xs : (3,) array_like + Position of virtual line source. + c : float, optional + Speed of sound. + + Returns + ------- + d : (N,) numpy.ndarray + Complex weights of secondary sources. + selection : (N,) numpy.ndarray + Boolean array containing ``True`` or ``False`` depending on + whether the corresponding secondary source is "active" or not. + secondary_source_function : callable + A function that can be used to create the sound field of a + single secondary source. See `sfs.mono.synthesize()`. + + Notes + ----- + .. math:: + + D(\x_0,\w) = \frac{\i}{2} \wc + \frac{\scalarprod{\x-\x_0}{\n_0}}{|\x-\x_\text{s}|} + \Hankel{2}{1}{\wc|\x-\x_\text{s}|} + + Examples + -------- + .. plot:: + :context: close-figs + + d, selection, secondary_source = sfs.mono.wfs.line_2d( + omega, array.x, array.n, xs) + plot(d, selection, secondary_source) + + """ + x0 = util.asarray_of_rows(x0) + n0 = util.asarray_of_rows(n0) + xs = util.asarray_1d(xs) + k = util.wavenumber(omega, c) + ds = x0 - xs + r = np.linalg.norm(ds, axis=1) + d = -1j/2 * k * inner1d(ds, n0) / r * hankel2(1, k * r) + selection = util.source_selection_line(n0, x0, xs) + return d, selection, secondary_source_line(omega, c) + + +def _point(omega, x0, n0, xs, c=None): + r"""Driving function for 2/3-dimensional WFS for a virtual point source. + + Parameters + ---------- + omega : float + Angular frequency of point source. + x0 : (N, 3) array_like + Sequence of secondary source positions. + n0 : (N, 3) array_like + Sequence of normal vectors of secondary sources. + xs : (3,) array_like + Position of virtual point source. + c : float, optional + Speed of sound. + + Returns + ------- + d : (N,) numpy.ndarray + Complex weights of secondary sources. + selection : (N,) numpy.ndarray + Boolean array containing ``True`` or ``False`` depending on + whether the corresponding secondary source is "active" or not. + secondary_source_function : callable + A function that can be used to create the sound field of a + single secondary source. See `sfs.mono.synthesize()`. + + Notes + ----- + .. math:: + + D(\x_0, \w) = \i\wc \frac{\scalarprod{\x_0-\x_\text{s}}{\n_0}} + {|\x_0-\x_\text{s}|^{\frac{3}{2}}} + \e{-\i\wc |\x_0-\x_\text{s}|} + + Examples + -------- + .. plot:: + :context: close-figs + + d, selection, secondary_source = sfs.mono.wfs.point_3d( + omega, array.x, array.n, xs) + plot(d, selection, secondary_source) + + """ + x0 = util.asarray_of_rows(x0) + n0 = util.asarray_of_rows(n0) + xs = util.asarray_1d(xs) + k = util.wavenumber(omega, c) + ds = x0 - xs + r = np.linalg.norm(ds, axis=1) + d = 1j * k * inner1d(ds, n0) / r ** (3 / 2) * np.exp(-1j * k * r) + selection = util.source_selection_point(n0, x0, xs) + return d, selection, secondary_source_point(omega, c) + + +point_2d = _point + + +def point_25d(omega, x0, n0, xs, xref=[0, 0, 0], c=None, omalias=None): + r"""Driving function for 2.5-dimensional WFS for a virtual point source. + + Parameters + ---------- + omega : float + Angular frequency of point source. + x0 : (N, 3) array_like + Sequence of secondary source positions. + n0 : (N, 3) array_like + Sequence of normal vectors of secondary sources. + xs : (3,) array_like + Position of virtual point source. + xref : (3,) array_like, optional + Reference point for synthesized sound field. + c : float, optional + Speed of sound. + + Returns + ------- + d : (N,) numpy.ndarray + Complex weights of secondary sources. + selection : (N,) numpy.ndarray + Boolean array containing ``True`` or ``False`` depending on + whether the corresponding secondary source is "active" or not. + secondary_source_function : callable + A function that can be used to create the sound field of a + single secondary source. See `sfs.mono.synthesize()`. + + Notes + ----- + .. math:: + + D(\x_0,\w) = \sqrt{\i\wc |\x_\text{ref}-\x_0|} + \frac{\scalarprod{\x_0-\x_\text{s}}{\n_0}} + {|\x_0-\x_\text{s}|^\frac{3}{2}} + \e{-\i\wc |\x_0-\x_\text{s}|} + + Examples + -------- + .. plot:: + :context: close-figs + + d, selection, secondary_source = sfs.mono.wfs.point_25d( + omega, array.x, array.n, xs) + plot(d, selection, secondary_source) + + """ + x0 = util.asarray_of_rows(x0) + n0 = util.asarray_of_rows(n0) + xs = util.asarray_1d(xs) + xref = util.asarray_1d(xref) + k = util.wavenumber(omega, c) + ds = x0 - xs + r = np.linalg.norm(ds, axis=1) + d = ( + preeq_25d(omega, omalias, c) * + np.sqrt(np.linalg.norm(xref - x0)) * inner1d(ds, n0) / + r ** (3 / 2) * np.exp(-1j * k * r)) + selection = util.source_selection_point(n0, x0, xs) + return d, selection, secondary_source_point(omega, c) + + +point_3d = _point + + +def _plane(omega, x0, n0, n=[0, 1, 0], c=None): + r"""Driving function for 2/3-dimensional WFS for a virtual plane wave. + + Parameters + ---------- + omega : float + Angular frequency of plane wave. + x0 : (N, 3) array_like + Sequence of secondary source positions. + n0 : (N, 3) array_like + Sequence of normal vectors of secondary sources. + n : (3,) array_like, optional + Normal vector (traveling direction) of plane wave. + c : float, optional + Speed of sound. + + Returns + ------- + d : (N,) numpy.ndarray + Complex weights of secondary sources. + selection : (N,) numpy.ndarray + Boolean array containing ``True`` or ``False`` depending on + whether the corresponding secondary source is "active" or not. + secondary_source_function : callable + A function that can be used to create the sound field of a + single secondary source. See `sfs.mono.synthesize()`. + + Notes + ----- + Eq.(17) from :cite:`Spors2008`: + + .. math:: + + D(\x_0,\w) = \i\wc \scalarprod{\n}{\n_0} + \e{-\i\wc\scalarprod{\n}{\x_0}} + + Examples + -------- + .. plot:: + :context: close-figs + + d, selection, secondary_source = sfs.mono.wfs.plane_3d( + omega, array.x, array.n, npw) + plot(d, selection, secondary_source) + + """ + x0 = util.asarray_of_rows(x0) + n0 = util.asarray_of_rows(n0) + n = util.normalize_vector(n) + k = util.wavenumber(omega, c) + d = 2j * k * np.inner(n, n0) * np.exp(-1j * k * np.inner(n, x0)) + selection = util.source_selection_plane(n0, n) + return d, selection, secondary_source_point(omega, c) + + +plane_2d = _plane + + +def plane_25d(omega, x0, n0, n=[0, 1, 0], xref=[0, 0, 0], c=None, + omalias=None): + r"""Driving function for 2.5-dimensional WFS for a virtual plane wave. + + Parameters + ---------- + omega : float + Angular frequency of plane wave. + x0 : (N, 3) array_like + Sequence of secondary source positions. + n0 : (N, 3) array_like + Sequence of normal vectors of secondary sources. + n : (3,) array_like, optional + Normal vector (traveling direction) of plane wave. + xref : (3,) array_like, optional + Reference point for synthesized sound field. + c : float, optional + Speed of sound. + omalias: float, optional + Angular frequency where spatial aliasing becomes prominent. + + Returns + ------- + d : (N,) numpy.ndarray + Complex weights of secondary sources. + selection : (N,) numpy.ndarray + Boolean array containing ``True`` or ``False`` depending on + whether the corresponding secondary source is "active" or not. + secondary_source_function : callable + A function that can be used to create the sound field of a + single secondary source. See `sfs.mono.synthesize()`. + + Notes + ----- + .. math:: + + D_\text{2.5D}(\x_0,\w) = \sqrt{\i\wc |\x_\text{ref}-\x_0|} + \scalarprod{\n}{\n_0} + \e{-\i\wc \scalarprod{\n}{\x_0}} + + Examples + -------- + .. plot:: + :context: close-figs + + d, selection, secondary_source = sfs.mono.wfs.plane_25d( + omega, array.x, array.n, npw) + plot(d, selection, secondary_source) + + """ + x0 = util.asarray_of_rows(x0) + n0 = util.asarray_of_rows(n0) + n = util.normalize_vector(n) + xref = util.asarray_1d(xref) + k = util.wavenumber(omega, c) + d = ( + preeq_25d(omega, omalias, c) * + np.sqrt(8*np.pi * np.linalg.norm(xref - x0, axis=-1)) * + np.inner(n, n0) * np.exp(-1j * k * np.inner(n, x0))) + selection = util.source_selection_plane(n0, n) + return d, selection, secondary_source_point(omega, c) + + +plane_3d = _plane + + +def _focused(omega, x0, n0, xs, ns, c=None): + r"""Driving function for 2/3-dimensional WFS for a focused source. + + Parameters + ---------- + omega : float + Angular frequency of focused source. + x0 : (N, 3) array_like + Sequence of secondary source positions. + n0 : (N, 3) array_like + Sequence of normal vectors of secondary sources. + xs : (3,) array_like + Position of focused source. + ns : (3,) array_like + Direction of focused source. + c : float, optional + Speed of sound. + + Returns + ------- + d : (N,) numpy.ndarray + Complex weights of secondary sources. + selection : (N,) numpy.ndarray + Boolean array containing ``True`` or ``False`` depending on + whether the corresponding secondary source is "active" or not. + secondary_source_function : callable + A function that can be used to create the sound field of a + single secondary source. See `sfs.mono.synthesize()`. + + Notes + ----- + .. math:: + + D(\x_0,\w) = \i\wc \frac{\scalarprod{\x_0-\x_\text{s}}{\n_0}} + {|\x_0-\x_\text{s}|^\frac{3}{2}} + \e{\i\wc |\x_0-\x_\text{s}|} + + Examples + -------- + .. plot:: + :context: close-figs + + d, selection, secondary_source = sfs.mono.wfs.focused_3d( + omega, array.x, array.n, xs_focused, ns_focused) + plot(d, selection, secondary_source) + + """ + x0 = util.asarray_of_rows(x0) + n0 = util.asarray_of_rows(n0) + xs = util.asarray_1d(xs) + k = util.wavenumber(omega, c) + ds = x0 - xs + r = np.linalg.norm(ds, axis=1) + d = 1j * k * inner1d(ds, n0) / r ** (3 / 2) * np.exp(1j * k * r) + selection = util.source_selection_focused(ns, x0, xs) + return d, selection, secondary_source_point(omega, c) + + +focused_2d = _focused + + +def focused_25d(omega, x0, n0, xs, ns, xref=[0, 0, 0], c=None, + omalias=None): + r"""Driving function for 2.5-dimensional WFS for a focused source. + + Parameters + ---------- + omega : float + Angular frequency of focused source. + x0 : (N, 3) array_like + Sequence of secondary source positions. + n0 : (N, 3) array_like + Sequence of normal vectors of secondary sources. + xs : (3,) array_like + Position of focused source. + ns : (3,) array_like + Direction of focused source. + xref : (3,) array_like, optional + Reference point for synthesized sound field. + c : float, optional + Speed of sound. + omalias: float, optional + Angular frequency where spatial aliasing becomes prominent. + + Returns + ------- + d : (N,) numpy.ndarray + Complex weights of secondary sources. + selection : (N,) numpy.ndarray + Boolean array containing ``True`` or ``False`` depending on + whether the corresponding secondary source is "active" or not. + secondary_source_function : callable + A function that can be used to create the sound field of a + single secondary source. See `sfs.mono.synthesize()`. + + Notes + ----- + .. math:: + + D(\x_0,\w) = \sqrt{\i\wc |\x_\text{ref}-\x_0|} + \frac{\scalarprod{\x_0-\x_\text{s}}{\n_0}} + {|\x_0-\x_\text{s}|^\frac{3}{2}} + \e{\i\wc |\x_0-\x_\text{s}|} + + Examples + -------- + .. plot:: + :context: close-figs + + d, selection, secondary_source = sfs.mono.wfs.focused_25d( + omega, array.x, array.n, xs_focused, ns_focused) + plot(d, selection, secondary_source) + + """ + x0 = util.asarray_of_rows(x0) + n0 = util.asarray_of_rows(n0) + xs = util.asarray_1d(xs) + xref = util.asarray_1d(xref) + k = util.wavenumber(omega, c) + ds = x0 - xs + r = np.linalg.norm(ds, axis=1) + d = ( + preeq_25d(omega, omalias, c) * + np.sqrt(np.linalg.norm(xref - x0)) * inner1d(ds, n0) / + r ** (3 / 2) * np.exp(1j * k * r)) + selection = util.source_selection_focused(ns, x0, xs) + return d, selection, secondary_source_point(omega, c) + + +focused_3d = _focused + + +def preeq_25d(omega, omalias, c): + r"""Pre-equalization filter for 2.5-dimensional WFS. + + Parameters + ---------- + omega : float + Angular frequency. + omalias: float + Angular frequency where spatial aliasing becomes prominent. + c : float + Speed of sound. + + Returns + ------- + float + Complex weight for given angular frequency. + + Notes + ----- + .. math:: + + H(\w) = \begin{cases} + \sqrt{\i \wc} & \text{for } \w \leq \w_\text{alias} \\ + \sqrt{\i \frac{\w_\text{alias}}{c}} & \text{for } \w > \w_\text{alias} + \end{cases} + + """ + if omalias is None: + return np.sqrt(1j * util.wavenumber(omega, c)) + else: + if omega <= omalias: + return np.sqrt(1j * util.wavenumber(omega, c)) + else: + return np.sqrt(1j * util.wavenumber(omalias, c)) + + +def plane_3d_delay(omega, x0, n0, n=[0, 1, 0], c=None): + r"""Delay-only driving function for a virtual plane wave. + + Parameters + ---------- + omega : float + Angular frequency of plane wave. + x0 : (N, 3) array_like + Sequence of secondary source positions. + n0 : (N, 3) array_like + Sequence of normal vectors of secondary sources. + n : (3,) array_like, optional + Normal vector (traveling direction) of plane wave. + c : float, optional + Speed of sound. + + Returns + ------- + d : (N,) numpy.ndarray + Complex weights of secondary sources. + selection : (N,) numpy.ndarray + Boolean array containing ``True`` or ``False`` depending on + whether the corresponding secondary source is "active" or not. + secondary_source_function : callable + A function that can be used to create the sound field of a + single secondary source. See `sfs.mono.synthesize()`. + + Notes + ----- + .. math:: + + D(\x_0,\w) = \e{-\i\wc\scalarprod{\n}{\x_0}} + + Examples + -------- + .. plot:: + :context: close-figs + + d, selection, secondary_source = sfs.mono.wfs.plane_3d_delay( + omega, array.x, array.n, npw) + plot(d, selection, secondary_source) + + """ + x0 = util.asarray_of_rows(x0) + n = util.normalize_vector(n) + k = util.wavenumber(omega, c) + d = np.exp(-1j * k * np.inner(n, x0)) + selection = util.source_selection_plane(n0, n) + return d, selection, secondary_source_point(omega, c) + + +def secondary_source_point(omega, c): + """Create a point source for use in `sfs.mono.synthesize()`.""" + + def secondary_source(position, _, grid): + return _source.point(omega, position, grid, c) + + return secondary_source + + +def secondary_source_line(omega, c): + """Create a line source for use in `sfs.mono.synthesize()`.""" + + def secondary_source(position, _, grid): + return _source.line(omega, position, grid, c) + + return secondary_source From bcc8239bfe1b8c8d45d73f38b8e31cc366d9676a Mon Sep 17 00:00:00 2001 From: Frank Schultz Date: Tue, 12 Mar 2019 19:55:07 +0100 Subject: [PATCH 02/15] __init.py__, soundfigure into wfs handling --- sfs/mono/__init__.py | 26 +- sfs/mono/drivingfunction.py | 1351 ----------------------------------- sfs/mono/esaedge.py | 20 +- sfs/mono/nfchoa.py | 20 +- sfs/mono/sdm.py | 29 +- sfs/mono/soundfigure.py | 46 -- sfs/mono/wfs.py | 56 +- 7 files changed, 67 insertions(+), 1481 deletions(-) delete mode 100644 sfs/mono/drivingfunction.py delete mode 100644 sfs/mono/soundfigure.py diff --git a/sfs/mono/__init__.py b/sfs/mono/__init__.py index 0fa4f6f..dd328f4 100644 --- a/sfs/mono/__init__.py +++ b/sfs/mono/__init__.py @@ -9,20 +9,38 @@ wfs source - soundfigure """ +from . import source as _source + + +def secondary_source_point(omega, c): + """Create a point source for use in `sfs.mono.synthesize()`.""" + + def secondary_source(position, _, grid): + return _source.point(omega, position, grid, c) + + return secondary_source + + +def secondary_source_line(omega, c): + """Create a line source for use in `sfs.mono.synthesize()`.""" + + def secondary_source(position, _, grid): + return _source.line(omega, position, grid, c) + + return secondary_source + + import numpy as _np +from . import source from . import esaedge from . import nfchoa from . import sdm from . import wfs -from . import source -from . import soundfigure - from .. import array as _array diff --git a/sfs/mono/drivingfunction.py b/sfs/mono/drivingfunction.py deleted file mode 100644 index 5dfb574..0000000 --- a/sfs/mono/drivingfunction.py +++ /dev/null @@ -1,1351 +0,0 @@ -"""Compute driving functions for various systems. - -.. include:: math-definitions.rst - -.. plot:: - :context: reset - - import matplotlib.pyplot as plt - import numpy as np - import sfs - - plt.rcParams['figure.figsize'] = 6, 6 - - xs = -1.5, 1.5, 0 - xs_focused = -0.5, 0.5, 0 - # normal vector for plane wave: - npw = sfs.util.direction_vector(np.radians(-45)) - # normal vector for focused source: - ns_focused = sfs.util.direction_vector(np.radians(-45)) - f = 300 # Hz - omega = 2 * np.pi * f - R = 1.5 # Radius of circular loudspeaker array - - grid = sfs.util.xyz_grid([-2, 2], [-2, 2], 0, spacing=0.02) - - array = sfs.array.circular(N=32, R=R) - - def plot(d, selection, secondary_source): - p = sfs.mono.synthesize(d, selection, array, secondary_source, grid=grid) - sfs.plot.soundfield(p, grid) - sfs.plot.loudspeaker_2d(array.x, array.n, selection * array.a, size=0.15) - -""" - -import numpy as np -from numpy.core.umath_tests import inner1d # element-wise inner product -from scipy.special import jn, hankel2 -from .. import util -from .. import default -from . import source as _source - - -def wfs_2d_line(omega, x0, n0, xs, c=None): - r"""Driving function for 2-dimensional WFS for a virtual line source. - - Parameters - ---------- - omega : float - Angular frequency of line source. - x0 : (N, 3) array_like - Sequence of secondary source positions. - n0 : (N, 3) array_like - Sequence of normal vectors of secondary sources. - xs : (3,) array_like - Position of virtual line source. - c : float, optional - Speed of sound. - - Returns - ------- - d : (N,) numpy.ndarray - Complex weights of secondary sources. - selection : (N,) numpy.ndarray - Boolean array containing ``True`` or ``False`` depending on - whether the corresponding secondary source is "active" or not. - secondary_source_function : callable - A function that can be used to create the sound field of a - single secondary source. See `sfs.mono.synthesize()`. - - Notes - ----- - .. math:: - - D(\x_0,\w) = \frac{\i}{2} \wc - \frac{\scalarprod{\x-\x_0}{\n_0}}{|\x-\x_\text{s}|} - \Hankel{2}{1}{\wc|\x-\x_\text{s}|} - - Examples - -------- - .. plot:: - :context: close-figs - - d, selection, secondary_source = sfs.mono.drivingfunction.wfs_2d_line( - omega, array.x, array.n, xs) - plot(d, selection, secondary_source) - - """ - x0 = util.asarray_of_rows(x0) - n0 = util.asarray_of_rows(n0) - xs = util.asarray_1d(xs) - k = util.wavenumber(omega, c) - ds = x0 - xs - r = np.linalg.norm(ds, axis=1) - d = -1j/2 * k * inner1d(ds, n0) / r * hankel2(1, k * r) - selection = util.source_selection_line(n0, x0, xs) - return d, selection, secondary_source_line(omega, c) - - -def _wfs_point(omega, x0, n0, xs, c=None): - r"""Driving function for 2/3-dimensional WFS for a virtual point source. - - Parameters - ---------- - omega : float - Angular frequency of point source. - x0 : (N, 3) array_like - Sequence of secondary source positions. - n0 : (N, 3) array_like - Sequence of normal vectors of secondary sources. - xs : (3,) array_like - Position of virtual point source. - c : float, optional - Speed of sound. - - Returns - ------- - d : (N,) numpy.ndarray - Complex weights of secondary sources. - selection : (N,) numpy.ndarray - Boolean array containing ``True`` or ``False`` depending on - whether the corresponding secondary source is "active" or not. - secondary_source_function : callable - A function that can be used to create the sound field of a - single secondary source. See `sfs.mono.synthesize()`. - - Notes - ----- - .. math:: - - D(\x_0, \w) = \i\wc \frac{\scalarprod{\x_0-\x_\text{s}}{\n_0}} - {|\x_0-\x_\text{s}|^{\frac{3}{2}}} - \e{-\i\wc |\x_0-\x_\text{s}|} - - Examples - -------- - .. plot:: - :context: close-figs - - d, selection, secondary_source = sfs.mono.drivingfunction.wfs_3d_point( - omega, array.x, array.n, xs) - plot(d, selection, secondary_source) - - """ - x0 = util.asarray_of_rows(x0) - n0 = util.asarray_of_rows(n0) - xs = util.asarray_1d(xs) - k = util.wavenumber(omega, c) - ds = x0 - xs - r = np.linalg.norm(ds, axis=1) - d = 1j * k * inner1d(ds, n0) / r ** (3 / 2) * np.exp(-1j * k * r) - selection = util.source_selection_point(n0, x0, xs) - return d, selection, secondary_source_point(omega, c) - - -wfs_2d_point = _wfs_point - - -def wfs_25d_point(omega, x0, n0, xs, xref=[0, 0, 0], c=None, omalias=None): - r"""Driving function for 2.5-dimensional WFS for a virtual point source. - - Parameters - ---------- - omega : float - Angular frequency of point source. - x0 : (N, 3) array_like - Sequence of secondary source positions. - n0 : (N, 3) array_like - Sequence of normal vectors of secondary sources. - xs : (3,) array_like - Position of virtual point source. - xref : (3,) array_like, optional - Reference point for synthesized sound field. - c : float, optional - Speed of sound. - - Returns - ------- - d : (N,) numpy.ndarray - Complex weights of secondary sources. - selection : (N,) numpy.ndarray - Boolean array containing ``True`` or ``False`` depending on - whether the corresponding secondary source is "active" or not. - secondary_source_function : callable - A function that can be used to create the sound field of a - single secondary source. See `sfs.mono.synthesize()`. - - Notes - ----- - .. math:: - - D(\x_0,\w) = \sqrt{\i\wc |\x_\text{ref}-\x_0|} - \frac{\scalarprod{\x_0-\x_\text{s}}{\n_0}} - {|\x_0-\x_\text{s}|^\frac{3}{2}} - \e{-\i\wc |\x_0-\x_\text{s}|} - - Examples - -------- - .. plot:: - :context: close-figs - - d, selection, secondary_source = sfs.mono.drivingfunction.wfs_25d_point( - omega, array.x, array.n, xs) - plot(d, selection, secondary_source) - - """ - x0 = util.asarray_of_rows(x0) - n0 = util.asarray_of_rows(n0) - xs = util.asarray_1d(xs) - xref = util.asarray_1d(xref) - k = util.wavenumber(omega, c) - ds = x0 - xs - r = np.linalg.norm(ds, axis=1) - d = ( - wfs_25d_preeq(omega, omalias, c) * - np.sqrt(np.linalg.norm(xref - x0)) * inner1d(ds, n0) / - r ** (3 / 2) * np.exp(-1j * k * r)) - selection = util.source_selection_point(n0, x0, xs) - return d, selection, secondary_source_point(omega, c) - - -wfs_3d_point = _wfs_point - - -def _wfs_plane(omega, x0, n0, n=[0, 1, 0], c=None): - r"""Driving function for 2/3-dimensional WFS for a virtual plane wave. - - Parameters - ---------- - omega : float - Angular frequency of plane wave. - x0 : (N, 3) array_like - Sequence of secondary source positions. - n0 : (N, 3) array_like - Sequence of normal vectors of secondary sources. - n : (3,) array_like, optional - Normal vector (traveling direction) of plane wave. - c : float, optional - Speed of sound. - - Returns - ------- - d : (N,) numpy.ndarray - Complex weights of secondary sources. - selection : (N,) numpy.ndarray - Boolean array containing ``True`` or ``False`` depending on - whether the corresponding secondary source is "active" or not. - secondary_source_function : callable - A function that can be used to create the sound field of a - single secondary source. See `sfs.mono.synthesize()`. - - Notes - ----- - Eq.(17) from :cite:`Spors2008`: - - .. math:: - - D(\x_0,\w) = \i\wc \scalarprod{\n}{\n_0} - \e{-\i\wc\scalarprod{\n}{\x_0}} - - Examples - -------- - .. plot:: - :context: close-figs - - d, selection, secondary_source = sfs.mono.drivingfunction.wfs_3d_plane( - omega, array.x, array.n, npw) - plot(d, selection, secondary_source) - - """ - x0 = util.asarray_of_rows(x0) - n0 = util.asarray_of_rows(n0) - n = util.normalize_vector(n) - k = util.wavenumber(omega, c) - d = 2j * k * np.inner(n, n0) * np.exp(-1j * k * np.inner(n, x0)) - selection = util.source_selection_plane(n0, n) - return d, selection, secondary_source_point(omega, c) - - -wfs_2d_plane = _wfs_plane - - -def wfs_25d_plane(omega, x0, n0, n=[0, 1, 0], xref=[0, 0, 0], c=None, - omalias=None): - r"""Driving function for 2.5-dimensional WFS for a virtual plane wave. - - Parameters - ---------- - omega : float - Angular frequency of plane wave. - x0 : (N, 3) array_like - Sequence of secondary source positions. - n0 : (N, 3) array_like - Sequence of normal vectors of secondary sources. - n : (3,) array_like, optional - Normal vector (traveling direction) of plane wave. - xref : (3,) array_like, optional - Reference point for synthesized sound field. - c : float, optional - Speed of sound. - omalias: float, optional - Angular frequency where spatial aliasing becomes prominent. - - Returns - ------- - d : (N,) numpy.ndarray - Complex weights of secondary sources. - selection : (N,) numpy.ndarray - Boolean array containing ``True`` or ``False`` depending on - whether the corresponding secondary source is "active" or not. - secondary_source_function : callable - A function that can be used to create the sound field of a - single secondary source. See `sfs.mono.synthesize()`. - - Notes - ----- - .. math:: - - D_\text{2.5D}(\x_0,\w) = \sqrt{\i\wc |\x_\text{ref}-\x_0|} - \scalarprod{\n}{\n_0} - \e{-\i\wc \scalarprod{\n}{\x_0}} - - Examples - -------- - .. plot:: - :context: close-figs - - d, selection, secondary_source = sfs.mono.drivingfunction.wfs_25d_plane( - omega, array.x, array.n, npw) - plot(d, selection, secondary_source) - - """ - x0 = util.asarray_of_rows(x0) - n0 = util.asarray_of_rows(n0) - n = util.normalize_vector(n) - xref = util.asarray_1d(xref) - k = util.wavenumber(omega, c) - d = ( - wfs_25d_preeq(omega, omalias, c) * - np.sqrt(8*np.pi * np.linalg.norm(xref - x0, axis=-1)) * - np.inner(n, n0) * np.exp(-1j * k * np.inner(n, x0))) - selection = util.source_selection_plane(n0, n) - return d, selection, secondary_source_point(omega, c) - - -wfs_3d_plane = _wfs_plane - - -def _wfs_focused(omega, x0, n0, xs, ns, c=None): - r"""Driving function for 2/3-dimensional WFS for a focused source. - - Parameters - ---------- - omega : float - Angular frequency of focused source. - x0 : (N, 3) array_like - Sequence of secondary source positions. - n0 : (N, 3) array_like - Sequence of normal vectors of secondary sources. - xs : (3,) array_like - Position of focused source. - ns : (3,) array_like - Direction of focused source. - c : float, optional - Speed of sound. - - Returns - ------- - d : (N,) numpy.ndarray - Complex weights of secondary sources. - selection : (N,) numpy.ndarray - Boolean array containing ``True`` or ``False`` depending on - whether the corresponding secondary source is "active" or not. - secondary_source_function : callable - A function that can be used to create the sound field of a - single secondary source. See `sfs.mono.synthesize()`. - - Notes - ----- - .. math:: - - D(\x_0,\w) = \i\wc \frac{\scalarprod{\x_0-\x_\text{s}}{\n_0}} - {|\x_0-\x_\text{s}|^\frac{3}{2}} - \e{\i\wc |\x_0-\x_\text{s}|} - - Examples - -------- - .. plot:: - :context: close-figs - - d, selection, secondary_source = sfs.mono.drivingfunction.wfs_3d_focused( - omega, array.x, array.n, xs_focused, ns_focused) - plot(d, selection, secondary_source) - - """ - x0 = util.asarray_of_rows(x0) - n0 = util.asarray_of_rows(n0) - xs = util.asarray_1d(xs) - k = util.wavenumber(omega, c) - ds = x0 - xs - r = np.linalg.norm(ds, axis=1) - d = 1j * k * inner1d(ds, n0) / r ** (3 / 2) * np.exp(1j * k * r) - selection = util.source_selection_focused(ns, x0, xs) - return d, selection, secondary_source_point(omega, c) - - -wfs_2d_focused = _wfs_focused - - -def wfs_25d_focused(omega, x0, n0, xs, ns, xref=[0, 0, 0], c=None, - omalias=None): - r"""Driving function for 2.5-dimensional WFS for a focused source. - - Parameters - ---------- - omega : float - Angular frequency of focused source. - x0 : (N, 3) array_like - Sequence of secondary source positions. - n0 : (N, 3) array_like - Sequence of normal vectors of secondary sources. - xs : (3,) array_like - Position of focused source. - ns : (3,) array_like - Direction of focused source. - xref : (3,) array_like, optional - Reference point for synthesized sound field. - c : float, optional - Speed of sound. - omalias: float, optional - Angular frequency where spatial aliasing becomes prominent. - - Returns - ------- - d : (N,) numpy.ndarray - Complex weights of secondary sources. - selection : (N,) numpy.ndarray - Boolean array containing ``True`` or ``False`` depending on - whether the corresponding secondary source is "active" or not. - secondary_source_function : callable - A function that can be used to create the sound field of a - single secondary source. See `sfs.mono.synthesize()`. - - Notes - ----- - .. math:: - - D(\x_0,\w) = \sqrt{\i\wc |\x_\text{ref}-\x_0|} - \frac{\scalarprod{\x_0-\x_\text{s}}{\n_0}} - {|\x_0-\x_\text{s}|^\frac{3}{2}} - \e{\i\wc |\x_0-\x_\text{s}|} - - Examples - -------- - .. plot:: - :context: close-figs - - d, selection, secondary_source = sfs.mono.drivingfunction.wfs_25d_focused( - omega, array.x, array.n, xs_focused, ns_focused) - plot(d, selection, secondary_source) - - """ - x0 = util.asarray_of_rows(x0) - n0 = util.asarray_of_rows(n0) - xs = util.asarray_1d(xs) - xref = util.asarray_1d(xref) - k = util.wavenumber(omega, c) - ds = x0 - xs - r = np.linalg.norm(ds, axis=1) - d = ( - wfs_25d_preeq(omega, omalias, c) * - np.sqrt(np.linalg.norm(xref - x0)) * inner1d(ds, n0) / - r ** (3 / 2) * np.exp(1j * k * r)) - selection = util.source_selection_focused(ns, x0, xs) - return d, selection, secondary_source_point(omega, c) - - -wfs_3d_focused = _wfs_focused - - -def wfs_25d_preeq(omega, omalias, c): - r"""Pre-equalization filter for 2.5-dimensional WFS. - - Parameters - ---------- - omega : float - Angular frequency. - omalias: float - Angular frequency where spatial aliasing becomes prominent. - c : float - Speed of sound. - - Returns - ------- - float - Complex weight for given angular frequency. - - Notes - ----- - .. math:: - - H(\w) = \begin{cases} - \sqrt{\i \wc} & \text{for } \w \leq \w_\text{alias} \\ - \sqrt{\i \frac{\w_\text{alias}}{c}} & \text{for } \w > \w_\text{alias} - \end{cases} - - """ - if omalias is None: - return np.sqrt(1j * util.wavenumber(omega, c)) - else: - if omega <= omalias: - return np.sqrt(1j * util.wavenumber(omega, c)) - else: - return np.sqrt(1j * util.wavenumber(omalias, c)) - - -def delay_3d_plane(omega, x0, n0, n=[0, 1, 0], c=None): - r"""Delay-only driving function for a virtual plane wave. - - Parameters - ---------- - omega : float - Angular frequency of plane wave. - x0 : (N, 3) array_like - Sequence of secondary source positions. - n0 : (N, 3) array_like - Sequence of normal vectors of secondary sources. - n : (3,) array_like, optional - Normal vector (traveling direction) of plane wave. - c : float, optional - Speed of sound. - - Returns - ------- - d : (N,) numpy.ndarray - Complex weights of secondary sources. - selection : (N,) numpy.ndarray - Boolean array containing ``True`` or ``False`` depending on - whether the corresponding secondary source is "active" or not. - secondary_source_function : callable - A function that can be used to create the sound field of a - single secondary source. See `sfs.mono.synthesize()`. - - Notes - ----- - .. math:: - - D(\x_0,\w) = \e{-\i\wc\scalarprod{\n}{\x_0}} - - Examples - -------- - .. plot:: - :context: close-figs - - d, selection, secondary_source = sfs.mono.drivingfunction.delay_3d_plane( - omega, array.x, array.n, npw) - plot(d, selection, secondary_source) - - """ - x0 = util.asarray_of_rows(x0) - n = util.normalize_vector(n) - k = util.wavenumber(omega, c) - d = np.exp(-1j * k * np.inner(n, x0)) - selection = util.source_selection_plane(n0, n) - return d, selection, secondary_source_point(omega, c) - - -def nfchoa_2d_plane(omega, x0, r0, n=[0, 1, 0], max_order=None, c=None): - r"""Driving function for 2-dimensional NFC-HOA for a virtual plane wave. - - Parameters - ---------- - omega : float - Angular frequency of plane wave. - x0 : (N, 3) array_like - Sequence of secondary source positions. - r0 : float - Radius of circular secondary source distribution. - n : (3,) array_like, optional - Normal vector (traveling direction) of plane wave. - max_order : float, optional - Maximum order of circular harmonics used for the calculation. - c : float, optional - Speed of sound. - - Returns - ------- - d : (N,) numpy.ndarray - Complex weights of secondary sources. - selection : (N,) numpy.ndarray - Boolean array containing only ``True`` indicating that - all secondary source are "active" for NFC-HOA. - secondary_source_function : callable - A function that can be used to create the sound field of a - single secondary source. See `sfs.mono.synthesize()`. - - Notes - ----- - .. math:: - - D(\phi_0, \omega) = - -\frac{2\i}{\pi r_0} - \sum_{m=-M}^M - \frac{\i^{-m}}{\Hankel{2}{m}{\wc r_0}} - \e{\i m (\phi_0 - \phi_\text{pw})} - - See http://sfstoolbox.org/#equation-D.nfchoa.pw.2D. - - Examples - -------- - .. plot:: - :context: close-figs - - d, selection, secondary_source = sfs.mono.drivingfunction.nfchoa_2d_plane( - omega, array.x, R, npw) - plot(d, selection, secondary_source) - - """ - if max_order is None: - max_order = util.max_order_circular_harmonics(len(x0)) - - x0 = util.asarray_of_rows(x0) - k = util.wavenumber(omega, c) - n = util.normalize_vector(n) - phi, _, r = util.cart2sph(*n) - phi0 = util.cart2sph(*x0.T)[0] - d = 0 - for m in range(-max_order, max_order + 1): - d += 1j**-m / hankel2(m, k * r0) * np.exp(1j * m * (phi0 - phi)) - selection = util.source_selection_all(len(x0)) - return -2j / (np.pi*r0) * d, selection, secondary_source_point(omega, c) - - -def nfchoa_25d_point(omega, x0, r0, xs, max_order=None, c=None): - r"""Driving function for 2.5-dimensional NFC-HOA for a virtual point source. - - Parameters - ---------- - omega : float - Angular frequency of point source. - x0 : (N, 3) array_like - Sequence of secondary source positions. - r0 : float - Radius of circular secondary source distribution. - xs : (3,) array_like - Position of point source. - max_order : float, optional - Maximum order of circular harmonics used for the calculation. - c : float, optional - Speed of sound. - - Returns - ------- - d : (N,) numpy.ndarray - Complex weights of secondary sources. - selection : (N,) numpy.ndarray - Boolean array containing only ``True`` indicating that - all secondary source are "active" for NFC-HOA. - secondary_source_function : callable - A function that can be used to create the sound field of a - single secondary source. See `sfs.mono.synthesize()`. - - Notes - ----- - .. math:: - - D(\phi_0, \omega) = - \frac{1}{2 \pi r_0} - \sum_{m=-M}^M - \frac{\hankel{2}{|m|}{\wc r}}{\hankel{2}{|m|}{\wc r_0}} - \e{\i m (\phi_0 - \phi)} - - See http://sfstoolbox.org/#equation-D.nfchoa.ps.2.5D. - - Examples - -------- - .. plot:: - :context: close-figs - - d, selection, secondary_source = sfs.mono.drivingfunction.nfchoa_25d_point( - omega, array.x, R, xs) - plot(d, selection, secondary_source) - - """ - if max_order is None: - max_order = util.max_order_circular_harmonics(len(x0)) - - x0 = util.asarray_of_rows(x0) - k = util.wavenumber(omega, c) - xs = util.asarray_1d(xs) - phi, _, r = util.cart2sph(*xs) - phi0 = util.cart2sph(*x0.T)[0] - hr = util.spherical_hn2(range(0, max_order + 1), k * r) - hr0 = util.spherical_hn2(range(0, max_order + 1), k * r0) - d = 0 - for m in range(-max_order, max_order + 1): - d += hr[abs(m)] / hr0[abs(m)] * np.exp(1j * m * (phi0 - phi)) - selection = util.source_selection_all(len(x0)) - return d / (2 * np.pi * r0), selection, secondary_source_point(omega, c) - - -def nfchoa_25d_plane(omega, x0, r0, n=[0, 1, 0], max_order=None, c=None): - r"""Driving function for 2.5-dimensional NFC-HOA for a virtual plane wave. - - Parameters - ---------- - omega : float - Angular frequency of point source. - x0 : (N, 3) array_like - Sequence of secondary source positions. - r0 : float - Radius of circular secondary source distribution. - n : (3,) array_like, optional - Normal vector (traveling direction) of plane wave. - max_order : float, optional - Maximum order of circular harmonics used for the calculation. - c : float, optional - Speed of sound. - - Returns - ------- - d : (N,) numpy.ndarray - Complex weights of secondary sources. - selection : (N,) numpy.ndarray - Boolean array containing only ``True`` indicating that - all secondary source are "active" for NFC-HOA. - secondary_source_function : callable - A function that can be used to create the sound field of a - single secondary source. See `sfs.mono.synthesize()`. - - Notes - ----- - .. math:: - - D(\phi_0, \omega) = - \frac{2\i}{r_0} - \sum_{m=-M}^M - \frac{\i^{-|m|}}{\wc \hankel{2}{|m|}{\wc r_0}} - \e{\i m (\phi_0 - \phi_\text{pw})} - - See http://sfstoolbox.org/#equation-D.nfchoa.pw.2.5D. - - Examples - -------- - .. plot:: - :context: close-figs - - d, selection, secondary_source = sfs.mono.drivingfunction.nfchoa_25d_plane( - omega, array.x, R, npw) - plot(d, selection, secondary_source) - - """ - if max_order is None: - max_order = util.max_order_circular_harmonics(len(x0)) - - x0 = util.asarray_of_rows(x0) - k = util.wavenumber(omega, c) - n = util.normalize_vector(n) - phi, _, r = util.cart2sph(*n) - phi0 = util.cart2sph(*x0.T)[0] - d = 0 - hn2 = util.spherical_hn2(range(0, max_order + 1), k * r0) - for m in range(-max_order, max_order + 1): - d += (-1j)**abs(m) / (k * hn2[abs(m)]) * np.exp(1j * m * (phi0 - phi)) - selection = util.source_selection_all(len(x0)) - return 2*1j / r0 * d, selection, secondary_source_point(omega, c) - - -def sdm_2d_line(omega, x0, n0, xs, c=None): - r"""Driving function for 2-dimensional SDM for a virtual line source. - - Parameters - ---------- - omega : float - Angular frequency of line source. - x0 : (N, 3) array_like - Sequence of secondary source positions. - n0 : (N, 3) array_like - Sequence of normal vectors of secondary sources. - xs : (3,) array_like - Position of line source. - c : float, optional - Speed of sound. - - Returns - ------- - d : (N,) numpy.ndarray - Complex weights of secondary sources. - selection : (N,) numpy.ndarray - Boolean array containing ``True`` or ``False`` depending on - whether the corresponding secondary source is "active" or not. - secondary_source_function : callable - A function that can be used to create the sound field of a - single secondary source. See `sfs.mono.synthesize()`. - - Notes - ----- - The secondary sources have to be located on the x-axis (y0=0). - Derived from :cite:`Spors2009`, Eq.(9), Eq.(4). - - Examples - -------- - .. plot:: - :context: close-figs - - array = sfs.array.linear(32, 0.2, orientation=[0, -1, 0]) - d, selection, secondary_source = sfs.mono.drivingfunction.sdm_2d_line( - omega, array.x, array.n, xs) - plot(d, selection, secondary_source) - - """ - x0 = util.asarray_of_rows(x0) - n0 = util.asarray_of_rows(n0) - xs = util.asarray_1d(xs) - k = util.wavenumber(omega, c) - ds = x0 - xs - r = np.linalg.norm(ds, axis=1) - d = - 1j/2 * k * xs[1] / r * hankel2(1, k * r) - selection = util.source_selection_all(len(x0)) - return d, selection, secondary_source_line(omega, c) - - -def sdm_2d_plane(omega, x0, n0, n=[0, 1, 0], c=None): - r"""Driving function for 2-dimensional SDM for a virtual plane wave. - - Parameters - ---------- - omega : float - Angular frequency of plane wave. - x0 : (N, 3) array_like - Sequence of secondary source positions. - n0 : (N, 3) array_like - Sequence of normal vectors of secondary sources. - n: (3,) array_like, optional - Normal vector (traveling direction) of plane wave. - c : float, optional - Speed of sound. - - Returns - ------- - d : (N,) numpy.ndarray - Complex weights of secondary sources. - selection : (N,) numpy.ndarray - Boolean array containing ``True`` or ``False`` depending on - whether the corresponding secondary source is "active" or not. - secondary_source_function : callable - A function that can be used to create the sound field of a - single secondary source. See `sfs.mono.synthesize()`. - - Notes - ----- - The secondary sources have to be located on the x-axis (y0=0). - Derived from :cite:`Ahrens2012`, Eq.(3.73), Eq.(C.5), Eq.(C.11): - - .. math:: - - D(\x_0,k) = k_\text{pw,y} \e{-\i k_\text{pw,x} x} - - Examples - -------- - .. plot:: - :context: close-figs - - array = sfs.array.linear(32, 0.2, orientation=[0, -1, 0]) - d, selection, secondary_source = sfs.mono.drivingfunction.sdm_2d_plane( - omega, array.x, array.n, npw) - plot(d, selection, secondary_source) - - """ - x0 = util.asarray_of_rows(x0) - n0 = util.asarray_of_rows(n0) - n = util.normalize_vector(n) - k = util.wavenumber(omega, c) - d = k * n[1] * np.exp(-1j * k * n[0] * x0[:, 0]) - selection = util.source_selection_all(len(x0)) - return d, selection, secondary_source_line(omega, c) - - -def sdm_25d_plane(omega, x0, n0, n=[0, 1, 0], xref=[0, 0, 0], c=None): - r"""Driving function for 2.5-dimensional SDM for a virtual plane wave. - - Parameters - ---------- - omega : float - Angular frequency of plane wave. - x0 : (N, 3) array_like - Sequence of secondary source positions. - n0 : (N, 3) array_like - Sequence of normal vectors of secondary sources. - n: (3,) array_like, optional - Normal vector (traveling direction) of plane wave. - xref : (3,) array_like, optional - Reference point for synthesized sound field. - c : float, optional - Speed of sound. - - Returns - ------- - d : (N,) numpy.ndarray - Complex weights of secondary sources. - selection : (N,) numpy.ndarray - Boolean array containing ``True`` or ``False`` depending on - whether the corresponding secondary source is "active" or not. - secondary_source_function : callable - A function that can be used to create the sound field of a - single secondary source. See `sfs.mono.synthesize()`. - - Notes - ----- - The secondary sources have to be located on the x-axis (y0=0). - Eq.(3.79) from :cite:`Ahrens2012`. - - Examples - -------- - .. plot:: - :context: close-figs - - array = sfs.array.linear(32, 0.2, orientation=[0, -1, 0]) - d, selection, secondary_source = ( - sfs.mono.drivingfunction.sdm_25d_plane( - omega, array.x, array.n, npw, [0, -1, 0])) - plot(d, selection, secondary_source) - - """ - x0 = util.asarray_of_rows(x0) - n0 = util.asarray_of_rows(n0) - n = util.normalize_vector(n) - xref = util.asarray_1d(xref) - k = util.wavenumber(omega, c) - d = 4j * np.exp(-1j*k*n[1]*xref[1]) / hankel2(0, k*n[1]*xref[1]) * \ - np.exp(-1j*k*n[0]*x0[:, 0]) - selection = util.source_selection_all(len(x0)) - return d, selection, secondary_source_point(omega, c) - - -def sdm_25d_point(omega, x0, n0, xs, xref=[0, 0, 0], c=None): - r"""Driving function for 2.5-dimensional SDM for a virtual point source. - - Parameters - ---------- - omega : float - Angular frequency of point source. - x0 : (N, 3) array_like - Sequence of secondary source positions. - n0 : (N, 3) array_like - Sequence of normal vectors of secondary sources. - xs: (3,) array_like - Position of virtual point source. - xref : (3,) array_like, optional - Reference point for synthesized sound field. - c : float, optional - Speed of sound. - - Returns - ------- - d : (N,) numpy.ndarray - Complex weights of secondary sources. - selection : (N,) numpy.ndarray - Boolean array containing ``True`` or ``False`` depending on - whether the corresponding secondary source is "active" or not. - secondary_source_function : callable - A function that can be used to create the sound field of a - single secondary source. See `sfs.mono.synthesize()`. - - Notes - ----- - The secondary sources have to be located on the x-axis (y0=0). - Driving function from :cite:`Spors2010`, Eq.(24). - - Examples - -------- - .. plot:: - :context: close-figs - - array = sfs.array.linear(32, 0.2, orientation=[0, -1, 0]) - d, selection, secondary_source = ( - sfs.mono.drivingfunction.sdm_25d_point( - omega, array.x, array.n, xs, [0, -1, 0])) - plot(d, selection, secondary_source) - - """ - x0 = util.asarray_of_rows(x0) - n0 = util.asarray_of_rows(n0) - xs = util.asarray_1d(xs) - xref = util.asarray_1d(xref) - k = util.wavenumber(omega, c) - ds = x0 - xs - r = np.linalg.norm(ds, axis=1) - d = 1/2 * 1j * k * np.sqrt(xref[1] / (xref[1] - xs[1])) * \ - xs[1] / r * hankel2(1, k * r) - selection = util.source_selection_all(len(x0)) - return d, selection, secondary_source_point(omega, c) - - -def esa_edge_2d_plane(omega, x0, n=[0, 1, 0], alpha=3/2*np.pi, Nc=None, - c=None): - r"""Driving function for 2-dimensional plane wave with edge ESA. - - Driving function for a virtual plane wave using the 2-dimensional ESA - for an edge-shaped secondary source distribution consisting of - monopole line sources. - - Parameters - ---------- - omega : float - Angular frequency. - x0 : int(N, 3) array_like - Sequence of secondary source positions. - n : (3,) array_like, optional - Normal vector of synthesized plane wave. - alpha : float, optional - Outer angle of edge. - Nc : int, optional - Number of elements for series expansion of driving function. Estimated - if not given. - c : float, optional - Speed of sound - - Returns - ------- - d : (N,) numpy.ndarray - Complex weights of secondary sources. - selection : (N,) numpy.ndarray - Boolean array containing ``True`` or ``False`` depending on - whether the corresponding secondary source is "active" or not. - secondary_source_function : callable - A function that can be used to create the sound field of a - single secondary source. See `sfs.mono.synthesize()`. - - Notes - ----- - One leg of the secondary sources has to be located on the x-axis (y0=0), - the edge at the origin. - - Derived from :cite:`Spors2016` - - """ - 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] - - r = np.linalg.norm(x0, axis=1) - phi = np.arctan2(x0[:, 1], x0[:, 0]) - phi = np.where(phi < 0, phi+2*np.pi, phi) - - if Nc is None: - Nc = np.ceil(2 * k * np.max(r) * alpha/np.pi) - - epsilon = np.ones(Nc) # weights for series expansion - epsilon[0] = 2 - - d = np.zeros(L, dtype=complex) - for m in np.arange(Nc): - nu = m*np.pi/alpha - d = d + 1/epsilon[m] * np.exp(1j*nu*np.pi/2) * np.sin(nu*phi_s) \ - * np.cos(nu*phi) * nu/r * jn(nu, k*r) - - d[phi > 0] = -d[phi > 0] - - selection = util.source_selection_all(len(x0)) - return 4*np.pi/alpha * d, selection, secondary_source_line(omega, c) - - -def esa_edge_dipole_2d_plane(omega, x0, n=[0, 1, 0], alpha=3/2*np.pi, Nc=None, - c=None): - r"""Driving function for 2-dimensional plane wave with edge dipole ESA. - - Driving function for a virtual plane wave using the 2-dimensional ESA - for an edge-shaped secondary source distribution consisting of - dipole line sources. - - Parameters - ---------- - omega : float - Angular frequency. - x0 : int(N, 3) array_like - Sequence of secondary source positions. - n : (3,) array_like, optional - Normal vector of synthesized plane wave. - alpha : float, optional - Outer angle of edge. - Nc : int, optional - Number of elements for series expansion of driving function. Estimated - if not given. - c : float, optional - Speed of sound - - Returns - ------- - d : (N,) numpy.ndarray - Complex weights of secondary sources. - selection : (N,) numpy.ndarray - Boolean array containing ``True`` or ``False`` depending on - whether the corresponding secondary source is "active" or not. - secondary_source_function : callable - A function that can be used to create the sound field of a - single secondary source. See `sfs.mono.synthesize()`. - - Notes - ----- - One leg of the secondary sources has to be located on the x-axis (y0=0), - the edge at the origin. - - Derived from :cite:`Spors2016` - - """ - 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] - - r = np.linalg.norm(x0, axis=1) - phi = np.arctan2(x0[:, 1], x0[:, 0]) - phi = np.where(phi < 0, phi+2*np.pi, phi) - - if Nc is None: - Nc = np.ceil(2 * k * np.max(r) * alpha/np.pi) - - epsilon = np.ones(Nc) # weights for series expansion - epsilon[0] = 2 - - d = np.zeros(L, dtype=complex) - for m in np.arange(Nc): - nu = m*np.pi/alpha - d = d + 1/epsilon[m] * np.exp(1j*nu*np.pi/2) * np.cos(nu*phi_s) \ - * np.cos(nu*phi) * jn(nu, k*r) - - return 4*np.pi/alpha * d - - -def esa_edge_2d_line(omega, x0, xs, alpha=3/2*np.pi, Nc=None, c=None): - r"""Driving function for 2-dimensional line source with edge ESA. - - Driving function for a virtual line source using the 2-dimensional ESA - for an edge-shaped secondary source distribution consisting of line - sources. - - Parameters - ---------- - omega : float - Angular frequency. - x0 : int(N, 3) array_like - Sequence of secondary source positions. - xs : (3,) array_like - Position of synthesized line source. - alpha : float, optional - Outer angle of edge. - Nc : int, optional - Number of elements for series expansion of driving function. Estimated - if not given. - c : float, optional - Speed of sound - - Returns - ------- - d : (N,) numpy.ndarray - Complex weights of secondary sources. - selection : (N,) numpy.ndarray - Boolean array containing ``True`` or ``False`` depending on - whether the corresponding secondary source is "active" or not. - secondary_source_function : callable - A function that can be used to create the sound field of a - single secondary source. See `sfs.mono.synthesize()`. - - Notes - ----- - One leg of the secondary sources has to be located on the x-axis (y0=0), - the edge at the origin. - - Derived from :cite:`Spors2016` - - """ - x0 = np.asarray(x0) - k = util.wavenumber(omega, c) - phi_s = np.arctan2(xs[1], xs[0]) - if phi_s < 0: - phi_s = phi_s + 2*np.pi - r_s = np.linalg.norm(xs) - L = x0.shape[0] - - r = np.linalg.norm(x0, axis=1) - phi = np.arctan2(x0[:, 1], x0[:, 0]) - phi = np.where(phi < 0, phi+2*np.pi, phi) - - if Nc is None: - Nc = np.ceil(2 * k * np.max(r) * alpha/np.pi) - - epsilon = np.ones(Nc) # weights for series expansion - epsilon[0] = 2 - - d = np.zeros(L, dtype=complex) - idx = (r <= r_s) - for m in np.arange(Nc): - nu = m*np.pi/alpha - f = 1/epsilon[m] * np.sin(nu*phi_s) * np.cos(nu*phi) * nu/r - d[idx] = d[idx] + f[idx] * jn(nu, k*r[idx]) * hankel2(nu, k*r_s) - d[~idx] = d[~idx] + f[~idx] * jn(nu, k*r_s) * hankel2(nu, k*r[~idx]) - - d[phi > 0] = -d[phi > 0] - - selection = util.source_selection_all(len(x0)) - return -1j*np.pi/alpha * d, selection, secondary_source_line(omega, c) - - -def esa_edge_dipole_2d_line(omega, x0, xs, alpha=3/2*np.pi, Nc=None, c=None): - r"""Driving function for 2-dimensional line source with edge dipole ESA. - - Driving function for a virtual line source using the 2-dimensional ESA - for an edge-shaped secondary source distribution consisting of dipole line - sources. - - Parameters - ---------- - omega : float - Angular frequency. - x0 : (N, 3) array_like - Sequence of secondary source positions. - xs : (3,) array_like - Position of synthesized line source. - alpha : float, optional - Outer angle of edge. - Nc : int, optional - Number of elements for series expansion of driving function. Estimated - if not given. - c : float, optional - Speed of sound - - Returns - ------- - d : (N,) numpy.ndarray - Complex weights of secondary sources. - selection : (N,) numpy.ndarray - Boolean array containing ``True`` or ``False`` depending on - whether the corresponding secondary source is "active" or not. - secondary_source_function : callable - A function that can be used to create the sound field of a - single secondary source. See `sfs.mono.synthesize()`. - - Notes - ----- - One leg of the secondary sources has to be located on the x-axis (y0=0), - the edge at the origin. - - Derived from :cite:`Spors2016` - - """ - x0 = np.asarray(x0) - k = util.wavenumber(omega, c) - phi_s = np.arctan2(xs[1], xs[0]) - if phi_s < 0: - phi_s = phi_s + 2*np.pi - r_s = np.linalg.norm(xs) - L = x0.shape[0] - - r = np.linalg.norm(x0, axis=1) - phi = np.arctan2(x0[:, 1], x0[:, 0]) - phi = np.where(phi < 0, phi+2*np.pi, phi) - - if Nc is None: - Nc = np.ceil(2 * k * np.max(r) * alpha/np.pi) - - epsilon = np.ones(Nc) # weights for series expansion - epsilon[0] = 2 - - d = np.zeros(L, dtype=complex) - idx = (r <= r_s) - for m in np.arange(Nc): - nu = m*np.pi/alpha - f = 1/epsilon[m] * np.cos(nu*phi_s) * np.cos(nu*phi) - d[idx] = d[idx] + f[idx] * jn(nu, k*r[idx]) * hankel2(nu, k*r_s) - d[~idx] = d[~idx] + f[~idx] * jn(nu, k*r_s) * hankel2(nu, k*r[~idx]) - - return -1j*np.pi/alpha * d - - -def esa_edge_25d_point(omega, x0, xs, xref=[2, -2, 0], alpha=3/2*np.pi, - Nc=None, c=None): - r"""Driving function for 2.5-dimensional point source with edge ESA. - - Driving function for a virtual point source using the 2.5-dimensional - ESA for an edge-shaped secondary source distribution consisting of point - sources. - - Parameters - ---------- - omega : float - Angular frequency. - x0 : int(N, 3) array_like - Sequence of secondary source positions. - xs : (3,) array_like - Position of synthesized line source. - xref: (3,) array_like or float - Reference position or reference distance - alpha : float, optional - Outer angle of edge. - Nc : int, optional - Number of elements for series expansion of driving function. Estimated - if not given. - c : float, optional - Speed of sound - - Returns - ------- - d : (N,) numpy.ndarray - Complex weights of secondary sources. - selection : (N,) numpy.ndarray - Boolean array containing ``True`` or ``False`` depending on - whether the corresponding secondary source is "active" or not. - secondary_source_function : callable - A function that can be used to create the sound field of a - single secondary source. See `sfs.mono.synthesize()`. - - Notes - ----- - One leg of the secondary sources has to be located on the x-axis (y0=0), - the edge at the origin. - - Derived from :cite:`Spors2016` - - """ - x0 = np.asarray(x0) - xs = np.asarray(xs) - xref = np.asarray(xref) - - if np.isscalar(xref): - a = np.linalg.norm(xref)/np.linalg.norm(xref-xs) - else: - a = np.linalg.norm(xref-x0, axis=1)/np.linalg.norm(xref-xs) - - d, selection, _ = esa_edge_2d_line(omega, x0, xs, alpha=alpha, Nc=Nc, c=c) - return 1j*np.sqrt(a) * d, selection, secondary_source_point(omega, c) - - -def secondary_source_point(omega, c): - """Create a point source for use in `sfs.mono.synthesize()`.""" - - def secondary_source(position, _, grid): - return _source.point(omega, position, grid, c) - - return secondary_source - - -def secondary_source_line(omega, c): - """Create a line source for use in `sfs.mono.synthesize()`.""" - - def secondary_source(position, _, grid): - return _source.line(omega, position, grid, c) - - return secondary_source diff --git a/sfs/mono/esaedge.py b/sfs/mono/esaedge.py index 7c84618..09544d3 100644 --- a/sfs/mono/esaedge.py +++ b/sfs/mono/esaedge.py @@ -37,7 +37,7 @@ def plot(d, selection, secondary_source): from scipy.special import jn, hankel2 from .. import util from .. import default -from . import source as _source +from . import secondary_source_line, secondary_source_point def plane_2d(omega, x0, n=[0, 1, 0], alpha=3/2*np.pi, Nc=None, @@ -380,21 +380,3 @@ def point_25d(omega, x0, xs, xref=[2, -2, 0], alpha=3/2*np.pi, d, selection, _ = line_2d(omega, x0, xs, alpha=alpha, Nc=Nc, c=c) return 1j*np.sqrt(a) * d, selection, secondary_source_point(omega, c) - - -def secondary_source_point(omega, c): - """Create a point source for use in `sfs.mono.synthesize()`.""" - - def secondary_source(position, _, grid): - return _source.point(omega, position, grid, c) - - return secondary_source - - -def secondary_source_line(omega, c): - """Create a line source for use in `sfs.mono.synthesize()`.""" - - def secondary_source(position, _, grid): - return _source.line(omega, position, grid, c) - - return secondary_source diff --git a/sfs/mono/nfchoa.py b/sfs/mono/nfchoa.py index 950eae9..ab65cfd 100644 --- a/sfs/mono/nfchoa.py +++ b/sfs/mono/nfchoa.py @@ -37,7 +37,7 @@ def plot(d, selection, secondary_source): from scipy.special import jn, hankel2 from .. import util from .. import default -from . import source as _source +from . import secondary_source_point def plane_2d(omega, x0, r0, n=[0, 1, 0], max_order=None, c=None): @@ -239,21 +239,3 @@ def plane_25d(omega, x0, r0, n=[0, 1, 0], max_order=None, c=None): d += (-1j)**abs(m) / (k * hn2[abs(m)]) * np.exp(1j * m * (phi0 - phi)) selection = util.source_selection_all(len(x0)) return 2*1j / r0 * d, selection, secondary_source_point(omega, c) - - -def secondary_source_point(omega, c): - """Create a point source for use in `sfs.mono.synthesize()`.""" - - def secondary_source(position, _, grid): - return _source.point(omega, position, grid, c) - - return secondary_source - - -def secondary_source_line(omega, c): - """Create a line source for use in `sfs.mono.synthesize()`.""" - - def secondary_source(position, _, grid): - return _source.line(omega, position, grid, c) - - return secondary_source diff --git a/sfs/mono/sdm.py b/sfs/mono/sdm.py index 50e3879..9c3869e 100644 --- a/sfs/mono/sdm.py +++ b/sfs/mono/sdm.py @@ -33,11 +33,9 @@ def plot(d, selection, secondary_source): """ import numpy as np -from numpy.core.umath_tests import inner1d # element-wise inner product -from scipy.special import jn, hankel2 +from scipy.special import hankel2 from .. import util -from .. import default -from . import source as _source +from . import secondary_source_line, secondary_source_point def line_2d(omega, x0, n0, xs, c=None): @@ -190,9 +188,8 @@ def plane_25d(omega, x0, n0, n=[0, 1, 0], xref=[0, 0, 0], c=None): :context: close-figs array = sfs.array.linear(32, 0.2, orientation=[0, -1, 0]) - d, selection, secondary_source = ( - sfs.mono.sdm._25d( - omega, array.x, array.n, npw, [0, -1, 0])) + d, selection, secondary_source = sfs.mono.sdm.plane_25d( + omega, array.x, array.n, npw, [0, -1, 0]) plot(d, selection, secondary_source) """ @@ -264,21 +261,3 @@ def point_25d(omega, x0, n0, xs, xref=[0, 0, 0], c=None): xs[1] / r * hankel2(1, k * r) selection = util.source_selection_all(len(x0)) return d, selection, secondary_source_point(omega, c) - - -def secondary_source_point(omega, c): - """Create a point source for use in `sfs.mono.synthesize()`.""" - - def secondary_source(position, _, grid): - return _source.point(omega, position, grid, c) - - return secondary_source - - -def secondary_source_line(omega, c): - """Create a line source for use in `sfs.mono.synthesize()`.""" - - def secondary_source(position, _, grid): - return _source.line(omega, position, grid, c) - - return secondary_source diff --git a/sfs/mono/soundfigure.py b/sfs/mono/soundfigure.py deleted file mode 100644 index 0c0cfb9..0000000 --- a/sfs/mono/soundfigure.py +++ /dev/null @@ -1,46 +0,0 @@ -"""Compute driving functions for sound figures.""" - -import numpy as np -from .. import util -from .wfs import plane_3d - - -def wfs_3d_pw(omega, x0, n0, figure, npw=[0, 0, 1], c=None): - """Compute driving function for a 2D sound figure. - - Based on - [Helwani et al., The Synthesis of Sound Figures, MSSP, 2013] - - """ - x0 = np.asarray(x0) - n0 = np.asarray(n0) - k = util.wavenumber(omega, c) - nx, ny = figure.shape - - # 2D spatial DFT of image - figure = np.fft.fftshift(figure, axes=(0, 1)) # sign of spatial DFT - figure = np.fft.fft2(figure) - # wavenumbers - kx = np.fft.fftfreq(nx, 1./nx) - ky = np.fft.fftfreq(ny, 1./ny) - # shift spectrum due to desired plane wave - figure = np.roll(figure, int(k*npw[0]), axis=0) - figure = np.roll(figure, int(k*npw[1]), axis=1) - # search and iterate over propagating plane wave components - kxx, kyy = np.meshgrid(kx, ky, sparse=True) - rho = np.sqrt((kxx) ** 2 + (kyy) ** 2) - d = 0 - for n in range(nx): - for m in range(ny): - if(rho[n, m] < k): - # dispertion relation - kz = np.sqrt(k**2 - rho[n, m]**2) - # normal vector of plane wave - npw = 1/k * np.asarray([kx[n], ky[m], kz]) - npw = npw / np.linalg.norm(npw) - # driving function of plane wave with positive kz - d_component, selection, secondary_source = wfs_3d_plane( - omega, x0, n0, npw, c) - d += selection * figure[n, m] * d_component - - return d, util.source_selection_all(len(d)), secondary_source diff --git a/sfs/mono/wfs.py b/sfs/mono/wfs.py index d54f157..341e2ef 100644 --- a/sfs/mono/wfs.py +++ b/sfs/mono/wfs.py @@ -34,10 +34,9 @@ def plot(d, selection, secondary_source): import numpy as np from numpy.core.umath_tests import inner1d # element-wise inner product -from scipy.special import jn, hankel2 +from scipy.special import hankel2 from .. import util -from .. import default -from . import source as _source +from . import secondary_source_line, secondary_source_point def line_2d(omega, x0, n0, xs, c=None): @@ -564,19 +563,42 @@ def plane_3d_delay(omega, x0, n0, n=[0, 1, 0], c=None): return d, selection, secondary_source_point(omega, c) -def secondary_source_point(omega, c): - """Create a point source for use in `sfs.mono.synthesize()`.""" +def soundfigure_2d(omega, x0, n0, figure, npw=[0, 0, 1], c=None): + """Compute driving function for a 2D sound figure. - def secondary_source(position, _, grid): - return _source.point(omega, position, grid, c) + Based on + [Helwani et al., The Synthesis of Sound Figures, MSSP, 2013] - return secondary_source - - -def secondary_source_line(omega, c): - """Create a line source for use in `sfs.mono.synthesize()`.""" - - def secondary_source(position, _, grid): - return _source.line(omega, position, grid, c) - - return secondary_source + """ + x0 = np.asarray(x0) + n0 = np.asarray(n0) + k = util.wavenumber(omega, c) + nx, ny = figure.shape + + # 2D spatial DFT of image + figure = np.fft.fftshift(figure, axes=(0, 1)) # sign of spatial DFT + figure = np.fft.fft2(figure) + # wavenumbers + kx = np.fft.fftfreq(nx, 1./nx) + ky = np.fft.fftfreq(ny, 1./ny) + # shift spectrum due to desired plane wave + figure = np.roll(figure, int(k*npw[0]), axis=0) + figure = np.roll(figure, int(k*npw[1]), axis=1) + # search and iterate over propagating plane wave components + kxx, kyy = np.meshgrid(kx, ky, sparse=True) + rho = np.sqrt((kxx) ** 2 + (kyy) ** 2) + d = 0 + for n in range(nx): + for m in range(ny): + if(rho[n, m] < k): + # dispertion relation + kz = np.sqrt(k**2 - rho[n, m]**2) + # normal vector of plane wave + npw = 1/k * np.asarray([kx[n], ky[m], kz]) + npw = npw / np.linalg.norm(npw) + # driving function of plane wave with positive kz + d_component, selection, secondary_source = plane_3d( + omega, x0, n0, npw, c) + d += selection * figure[n, m] * d_component + + return d, util.source_selection_all(len(d)), secondary_source From 9a7b09164a9a357b72a0cafa49ce86ab9edd317f Mon Sep 17 00:00:00 2001 From: Frank Schultz Date: Wed, 13 Mar 2019 11:40:05 +0100 Subject: [PATCH 03/15] Update __init__.py --- sfs/mono/__init__.py | 24 +++++++++++------------- 1 file changed, 11 insertions(+), 13 deletions(-) diff --git a/sfs/mono/__init__.py b/sfs/mono/__init__.py index dd328f4..2e5c4ac 100644 --- a/sfs/mono/__init__.py +++ b/sfs/mono/__init__.py @@ -3,22 +3,26 @@ .. autosummary:: :toctree: - esaedge + source + + wfs nfchoa sdm - wfs - - source + esaedge """ -from . import source as _source +from . import source + + +import numpy as _np +from .. import array as _array def secondary_source_point(omega, c): """Create a point source for use in `sfs.mono.synthesize()`.""" def secondary_source(position, _, grid): - return _source.point(omega, position, grid, c) + return source.point(omega, position, grid, c) return secondary_source @@ -27,22 +31,16 @@ def secondary_source_line(omega, c): """Create a line source for use in `sfs.mono.synthesize()`.""" def secondary_source(position, _, grid): - return _source.line(omega, position, grid, c) + return source.line(omega, position, grid, c) return secondary_source -import numpy as _np - -from . import source - from . import esaedge from . import nfchoa from . import sdm from . import wfs -from .. import array as _array - def shiftphase(p, phase): """Shift phase of a sound field.""" From da37636be625e9dfbb5fa98ceee3b2e336037006 Mon Sep 17 00:00:00 2001 From: Frank Schultz Date: Wed, 13 Mar 2019 11:44:50 +0100 Subject: [PATCH 04/15] Update esaedge.py --- sfs/mono/esaedge.py | 34 +--------------------------------- 1 file changed, 1 insertion(+), 33 deletions(-) diff --git a/sfs/mono/esaedge.py b/sfs/mono/esaedge.py index 09544d3..a6600d5 100644 --- a/sfs/mono/esaedge.py +++ b/sfs/mono/esaedge.py @@ -1,36 +1,4 @@ -"""Compute driving functions for various systems. - -.. include:: math-definitions.rst - -.. plot:: - :context: reset - - import matplotlib.pyplot as plt - import numpy as np - import sfs - - plt.rcParams['figure.figsize'] = 6, 6 - - xs = -1.5, 1.5, 0 - xs_focused = -0.5, 0.5, 0 - # normal vector for plane wave: - npw = sfs.util.direction_vector(np.radians(-45)) - # normal vector for focused source: - ns_focused = sfs.util.direction_vector(np.radians(-45)) - f = 300 # Hz - omega = 2 * np.pi * f - R = 1.5 # Radius of circular loudspeaker array - - grid = sfs.util.xyz_grid([-2, 2], [-2, 2], 0, spacing=0.02) - - array = sfs.array.circular(N=32, R=R) - - def plot(d, selection, secondary_source): - p = sfs.mono.synthesize(d, selection, array, secondary_source, grid=grid) - sfs.plot.soundfield(p, grid) - sfs.plot.loudspeaker_2d(array.x, array.n, selection * array.a, size=0.15) - -""" +"""Compute ESA EDGE driving functions for various systems.""" import numpy as np from numpy.core.umath_tests import inner1d # element-wise inner product From 23c4f8e575f3d7a63283412ce7de1e43d55d347d Mon Sep 17 00:00:00 2001 From: Frank Schultz Date: Wed, 13 Mar 2019 11:53:53 +0100 Subject: [PATCH 05/15] Update nfchoa.py --- sfs/mono/nfchoa.py | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/sfs/mono/nfchoa.py b/sfs/mono/nfchoa.py index ab65cfd..496725e 100644 --- a/sfs/mono/nfchoa.py +++ b/sfs/mono/nfchoa.py @@ -1,4 +1,4 @@ -"""Compute driving functions for various systems. +"""Compute NFC-HOA driving functions. .. include:: math-definitions.rst @@ -12,11 +12,8 @@ plt.rcParams['figure.figsize'] = 6, 6 xs = -1.5, 1.5, 0 - xs_focused = -0.5, 0.5, 0 # normal vector for plane wave: npw = sfs.util.direction_vector(np.radians(-45)) - # normal vector for focused source: - ns_focused = sfs.util.direction_vector(np.radians(-45)) f = 300 # Hz omega = 2 * np.pi * f R = 1.5 # Radius of circular loudspeaker array From 3779e598b626bb148d4c267c7bc34973583ce772 Mon Sep 17 00:00:00 2001 From: Frank Schultz Date: Wed, 13 Mar 2019 11:55:34 +0100 Subject: [PATCH 06/15] Update esaedge.py --- sfs/mono/esaedge.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/sfs/mono/esaedge.py b/sfs/mono/esaedge.py index a6600d5..671fd30 100644 --- a/sfs/mono/esaedge.py +++ b/sfs/mono/esaedge.py @@ -1,10 +1,8 @@ """Compute ESA EDGE driving functions for various systems.""" import numpy as np -from numpy.core.umath_tests import inner1d # element-wise inner product from scipy.special import jn, hankel2 from .. import util -from .. import default from . import secondary_source_line, secondary_source_point From 0615628109c14ffd53dce6506b1daa9af6e59e71 Mon Sep 17 00:00:00 2001 From: Frank Schultz Date: Wed, 13 Mar 2019 11:57:15 +0100 Subject: [PATCH 07/15] Update nfchoa.py --- sfs/mono/nfchoa.py | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/sfs/mono/nfchoa.py b/sfs/mono/nfchoa.py index 496725e..6677d22 100644 --- a/sfs/mono/nfchoa.py +++ b/sfs/mono/nfchoa.py @@ -30,10 +30,8 @@ def plot(d, selection, secondary_source): """ import numpy as np -from numpy.core.umath_tests import inner1d # element-wise inner product -from scipy.special import jn, hankel2 +from scipy.special import hankel2 from .. import util -from .. import default from . import secondary_source_point From 0c86e4d1bfc1a50c12e8f1f0ab18b7715121ac70 Mon Sep 17 00:00:00 2001 From: Frank Schultz Date: Wed, 13 Mar 2019 12:02:15 +0100 Subject: [PATCH 08/15] Update sdm.py --- sfs/mono/sdm.py | 12 ++---------- 1 file changed, 2 insertions(+), 10 deletions(-) diff --git a/sfs/mono/sdm.py b/sfs/mono/sdm.py index 9c3869e..e1a26cb 100644 --- a/sfs/mono/sdm.py +++ b/sfs/mono/sdm.py @@ -1,4 +1,4 @@ -"""Compute driving functions for various systems. +"""Compute SDM driving functions. .. include:: math-definitions.rst @@ -12,18 +12,14 @@ plt.rcParams['figure.figsize'] = 6, 6 xs = -1.5, 1.5, 0 - xs_focused = -0.5, 0.5, 0 # normal vector for plane wave: npw = sfs.util.direction_vector(np.radians(-45)) - # normal vector for focused source: - ns_focused = sfs.util.direction_vector(np.radians(-45)) f = 300 # Hz omega = 2 * np.pi * f - R = 1.5 # Radius of circular loudspeaker array grid = sfs.util.xyz_grid([-2, 2], [-2, 2], 0, spacing=0.02) - array = sfs.array.circular(N=32, R=R) + array = sfs.array.linear(32, 0.2, orientation=[0, -1, 0]) def plot(d, selection, secondary_source): p = sfs.mono.synthesize(d, selection, array, secondary_source, grid=grid) @@ -75,7 +71,6 @@ def line_2d(omega, x0, n0, xs, c=None): .. plot:: :context: close-figs - array = sfs.array.linear(32, 0.2, orientation=[0, -1, 0]) d, selection, secondary_source = sfs.mono.sdm.line_2d( omega, array.x, array.n, xs) plot(d, selection, secondary_source) @@ -133,7 +128,6 @@ def plane_2d(omega, x0, n0, n=[0, 1, 0], c=None): .. plot:: :context: close-figs - array = sfs.array.linear(32, 0.2, orientation=[0, -1, 0]) d, selection, secondary_source = sfs.mono.sdm.plane_2d( omega, array.x, array.n, npw) plot(d, selection, secondary_source) @@ -187,7 +181,6 @@ def plane_25d(omega, x0, n0, n=[0, 1, 0], xref=[0, 0, 0], c=None): .. plot:: :context: close-figs - array = sfs.array.linear(32, 0.2, orientation=[0, -1, 0]) d, selection, secondary_source = sfs.mono.sdm.plane_25d( omega, array.x, array.n, npw, [0, -1, 0]) plot(d, selection, secondary_source) @@ -243,7 +236,6 @@ def point_25d(omega, x0, n0, xs, xref=[0, 0, 0], c=None): .. plot:: :context: close-figs - array = sfs.array.linear(32, 0.2, orientation=[0, -1, 0]) d, selection, secondary_source = ( sfs.mono.sdm.point_25d( omega, array.x, array.n, xs, [0, -1, 0])) From aed7ab787c18b904e3d39f08a1b5f8725ee68596 Mon Sep 17 00:00:00 2001 From: Frank Schultz Date: Wed, 13 Mar 2019 12:04:52 +0100 Subject: [PATCH 09/15] Update sdm.py --- sfs/mono/sdm.py | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/sfs/mono/sdm.py b/sfs/mono/sdm.py index e1a26cb..2d775ee 100644 --- a/sfs/mono/sdm.py +++ b/sfs/mono/sdm.py @@ -236,9 +236,8 @@ def point_25d(omega, x0, n0, xs, xref=[0, 0, 0], c=None): .. plot:: :context: close-figs - d, selection, secondary_source = ( - sfs.mono.sdm.point_25d( - omega, array.x, array.n, xs, [0, -1, 0])) + d, selection, secondary_source = sfs.mono.sdm.point_25d( + omega, array.x, array.n, xs, [0, -1, 0]) plot(d, selection, secondary_source) """ From b227cbb85459101ae490d026d97cc8f36173a87f Mon Sep 17 00:00:00 2001 From: Frank Schultz Date: Wed, 13 Mar 2019 12:08:29 +0100 Subject: [PATCH 10/15] Update wfs.py --- sfs/mono/wfs.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/sfs/mono/wfs.py b/sfs/mono/wfs.py index 341e2ef..30e6e88 100644 --- a/sfs/mono/wfs.py +++ b/sfs/mono/wfs.py @@ -1,4 +1,4 @@ -"""Compute driving functions for various systems. +"""Compute WFS driving functions. .. include:: math-definitions.rst From 1a726e22f0844735f676131439533f8958d859f2 Mon Sep 17 00:00:00 2001 From: Frank Schultz Date: Wed, 13 Mar 2019 12:19:55 +0100 Subject: [PATCH 11/15] Update horizontal_plane_arrays.py --- doc/examples/horizontal_plane_arrays.py | 36 ++++++++++++------------- 1 file changed, 18 insertions(+), 18 deletions(-) diff --git a/doc/examples/horizontal_plane_arrays.py b/doc/examples/horizontal_plane_arrays.py index 75b07f9..c2abeab 100644 --- a/doc/examples/horizontal_plane_arrays.py +++ b/doc/examples/horizontal_plane_arrays.py @@ -46,34 +46,34 @@ def compute_and_plot_soundfield(title): # linear array, secondary point sources, virtual monopole array = sfs.array.linear(N, dx, center=acenter, orientation=anormal) -d, selection, secondary_source = sfs.mono.drivingfunction.wfs_3d_point( +d, selection, secondary_source = sfs.mono.wfs.point_3d( omega, array.x, array.n, xs) compute_and_plot_soundfield('linear_ps_wfs_3d_point') -d, selection, secondary_source = sfs.mono.drivingfunction.wfs_25d_point( +d, selection, secondary_source = sfs.mono.wfs.point_25d( omega, array.x, array.n, xs, xref=xnorm) compute_and_plot_soundfield('linear_ps_wfs_25d_point') -d, selection, secondary_source = sfs.mono.drivingfunction.wfs_2d_point( +d, selection, secondary_source = sfs.mono.wfs.point_2d( omega, array.x, array.n, xs) compute_and_plot_soundfield('linear_ps_wfs_2d_point') # linear array, secondary line sources, virtual line source -d, selection, secondary_source = sfs.mono.drivingfunction.wfs_2d_line( +d, selection, secondary_source = sfs.mono.wfs.line_2d( omega, array.x, array.n, xs) compute_and_plot_soundfield('linear_ls_wfs_2d_line') # linear array, secondary point sources, virtual plane wave -d, selection, secondary_source = sfs.mono.drivingfunction.wfs_3d_plane( +d, selection, secondary_source = sfs.mono.wfs.plane_3d( omega, array.x, array.n, npw) compute_and_plot_soundfield('linear_ps_wfs_3d_plane') -d, selection, secondary_source = sfs.mono.drivingfunction.wfs_25d_plane( +d, selection, secondary_source = sfs.mono.wfs.plane_25d( omega, array.x, array.n, npw, xref=xnorm) compute_and_plot_soundfield('linear_ps_wfs_25d_plane') -d, selection, secondary_source = sfs.mono.drivingfunction.wfs_2d_plane( +d, selection, secondary_source = sfs.mono.wfs.plane_2d( omega, array.x, array.n, npw) compute_and_plot_soundfield('linear_ps_wfs_2d_plane') @@ -82,11 +82,11 @@ def compute_and_plot_soundfield(title): array = sfs.array.linear_diff(N//3 * [dx] + N//3 * [dx/2] + N//3 * [dx], center=acenter, orientation=anormal) -d, selection, secondary_source = sfs.mono.drivingfunction.wfs_25d_point( +d, selection, secondary_source = sfs.mono.wfs.point_25d( omega, array.x, array.n, xs, xref=xnorm) compute_and_plot_soundfield('linear_nested_ps_wfs_25d_point') -d, selection, secondary_source = sfs.mono.drivingfunction.wfs_25d_plane( +d, selection, secondary_source = sfs.mono.wfs.plane_25d( omega, array.x, array.n, npw, xref=xnorm) compute_and_plot_soundfield('linear_nested_ps_wfs_25d_plane') @@ -95,22 +95,22 @@ def compute_and_plot_soundfield(title): array = sfs.array.linear_random(N, dx/2, 1.5*dx, center=acenter, orientation=anormal) -d, selection, secondary_source = sfs.mono.drivingfunction.wfs_25d_point( +d, selection, secondary_source = sfs.mono.wfs.point_25d( omega, array.x, array.n, xs, xref=xnorm) compute_and_plot_soundfield('linear_random_ps_wfs_25d_point') -d, selection, secondary_source = sfs.mono.drivingfunction.wfs_25d_plane( +d, selection, secondary_source = sfs.mono.wfs.plane_25d( omega, array.x, array.n, npw, xref=xnorm) compute_and_plot_soundfield('linear_random_ps_wfs_25d_plane') # rectangular array, secondary point sources array = sfs.array.rectangular((N, N//2), dx, center=acenter, orientation=anormal) -d, selection, secondary_source = sfs.mono.drivingfunction.wfs_25d_point( +d, selection, secondary_source = sfs.mono.wfs.point_25d( omega, array.x, array.n, xs, xref=xnorm) compute_and_plot_soundfield('rectangular_ps_wfs_25d_point') -d, selection, secondary_source = sfs.mono.drivingfunction.wfs_25d_plane( +d, selection, secondary_source = sfs.mono.wfs.plane_25d( omega, array.x, array.n, npw, xref=xnorm) compute_and_plot_soundfield('rectangular_ps_wfs_25d_plane') @@ -118,11 +118,11 @@ def compute_and_plot_soundfield(title): # circular array, secondary point sources N = 60 array = sfs.array.circular(N, 1, center=acenter) -d, selection, secondary_source = sfs.mono.drivingfunction.wfs_25d_point( +d, selection, secondary_source = sfs.mono.wfs.point_25d( omega, array.x, array.n, xs, xref=xnorm) compute_and_plot_soundfield('circular_ps_wfs_25d_point') -d, selection, secondary_source = sfs.mono.drivingfunction.wfs_25d_plane( +d, selection, secondary_source = sfs.mono.wfs.plane_25d( omega, array.x, array.n, npw, xref=xnorm) compute_and_plot_soundfield('circular_ps_wfs_25d_plane') @@ -132,7 +132,7 @@ def compute_and_plot_soundfield(title): xnorm = [0, 0, 0] talpha = 0 # switches off tapering -d, selection, secondary_source = sfs.mono.drivingfunction.nfchoa_2d_plane( +d, selection, secondary_source = sfs.mono.nfchoa.plane_2d( omega, array.x, 1, npw) compute_and_plot_soundfield('circular_ls_nfchoa_2d_plane') @@ -142,10 +142,10 @@ def compute_and_plot_soundfield(title): xnorm = [0, 0, 0] talpha = 0 # switches off tapering -d, selection, secondary_source = sfs.mono.drivingfunction.nfchoa_25d_point( +d, selection, secondary_source = sfs.mono.nfchoa.point_25d( omega, array.x, 1, xs) compute_and_plot_soundfield('circular_ps_nfchoa_25d_point') -d, selection, secondary_source = sfs.mono.drivingfunction.nfchoa_25d_plane( +d, selection, secondary_source = sfs.mono.nfchoa.plane_25d( omega, array.x, 1, npw) compute_and_plot_soundfield('circular_ps_nfchoa_25d_plane') From 74687a7fefdcbaae4edbfb3f6b751be808730444 Mon Sep 17 00:00:00 2001 From: Frank Schultz Date: Wed, 13 Mar 2019 12:44:42 +0100 Subject: [PATCH 12/15] soundfigure_3d update --- doc/examples/soundfigures.py | 2 +- sfs/mono/wfs.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/doc/examples/soundfigures.py b/doc/examples/soundfigures.py index e6a76a6..9405182 100644 --- a/doc/examples/soundfigures.py +++ b/doc/examples/soundfigures.py @@ -29,7 +29,7 @@ # driving function for sound figure figure = np.array(Image.open('figures/tree.png')) # read image from file figure = np.rot90(figure) # turn 0deg to the top -d, selection, secondary_source = sfs.mono.soundfigure.wfs_3d_pw( +d, selection, secondary_source = sfs.mono.wfs.soundfigure_3d( omega, array.x, array.n, figure, npw=npw) # compute synthesized sound field diff --git a/sfs/mono/wfs.py b/sfs/mono/wfs.py index 30e6e88..a05aa0f 100644 --- a/sfs/mono/wfs.py +++ b/sfs/mono/wfs.py @@ -563,7 +563,7 @@ def plane_3d_delay(omega, x0, n0, n=[0, 1, 0], c=None): return d, selection, secondary_source_point(omega, c) -def soundfigure_2d(omega, x0, n0, figure, npw=[0, 0, 1], c=None): +def soundfigure_3d(omega, x0, n0, figure, npw=[0, 0, 1], c=None): """Compute driving function for a 2D sound figure. Based on From 7424338f93fbd01cdf1aa6f23f97895b06c6cada Mon Sep 17 00:00:00 2001 From: Frank Schultz Date: Wed, 13 Mar 2019 13:15:53 +0100 Subject: [PATCH 13/15] esaedge to esa --- sfs/mono/__init__.py | 4 ++-- sfs/mono/{esaedge.py => esa.py} | 24 +++++++++++++++++------- 2 files changed, 19 insertions(+), 9 deletions(-) rename sfs/mono/{esaedge.py => esa.py} (92%) diff --git a/sfs/mono/__init__.py b/sfs/mono/__init__.py index 2e5c4ac..13fee05 100644 --- a/sfs/mono/__init__.py +++ b/sfs/mono/__init__.py @@ -8,7 +8,7 @@ wfs nfchoa sdm - esaedge + esa """ from . import source @@ -36,7 +36,7 @@ def secondary_source(position, _, grid): return secondary_source -from . import esaedge +from . import esa from . import nfchoa from . import sdm from . import wfs diff --git a/sfs/mono/esaedge.py b/sfs/mono/esa.py similarity index 92% rename from sfs/mono/esaedge.py rename to sfs/mono/esa.py index 671fd30..a6e5427 100644 --- a/sfs/mono/esaedge.py +++ b/sfs/mono/esa.py @@ -1,4 +1,14 @@ -"""Compute ESA EDGE driving functions for various systems.""" +"""Compute ESA driving functions for various systems. + + ESA is abbreviation for equivalent scattering approach. + + ESA driving functions for an edge-shaped SSD are provided below. + Further ESA for different geometries might be added here. + + Note that mode-matching (such as NFC-HOA, SDM) are equivalent + to ESA in their specific geometries (spherical/circular, planar/linear). + +""" import numpy as np from scipy.special import jn, hankel2 @@ -6,7 +16,7 @@ from . import secondary_source_line, secondary_source_point -def plane_2d(omega, x0, n=[0, 1, 0], alpha=3/2*np.pi, Nc=None, +def plane_2d_edge(omega, x0, n=[0, 1, 0], alpha=3/2*np.pi, Nc=None, c=None): r"""Driving function for 2-dimensional plane wave with edge ESA. @@ -77,7 +87,7 @@ def plane_2d(omega, x0, n=[0, 1, 0], alpha=3/2*np.pi, Nc=None, return 4*np.pi/alpha * d, selection, secondary_source_line(omega, c) -def plane_2d_with_dipole_ssd(omega, x0, n=[0, 1, 0], alpha=3/2*np.pi, Nc=None, +def plane_2d_edge_dipole_ssd(omega, x0, n=[0, 1, 0], alpha=3/2*np.pi, Nc=None, c=None): r"""Driving function for 2-dimensional plane wave with edge dipole ESA. @@ -145,7 +155,7 @@ def plane_2d_with_dipole_ssd(omega, x0, n=[0, 1, 0], alpha=3/2*np.pi, Nc=None, return 4*np.pi/alpha * d -def line_2d(omega, x0, xs, alpha=3/2*np.pi, Nc=None, c=None): +def line_2d_edge(omega, x0, xs, alpha=3/2*np.pi, Nc=None, c=None): r"""Driving function for 2-dimensional line source with edge ESA. Driving function for a virtual line source using the 2-dimensional ESA @@ -219,7 +229,7 @@ def line_2d(omega, x0, xs, alpha=3/2*np.pi, Nc=None, c=None): return -1j*np.pi/alpha * d, selection, secondary_source_line(omega, c) -def line_2d_with_dipole_ssd(omega, x0, xs, alpha=3/2*np.pi, Nc=None, c=None): +def line_2d_edge_dipole_ssd(omega, x0, xs, alpha=3/2*np.pi, Nc=None, c=None): r"""Driving function for 2-dimensional line source with edge dipole ESA. Driving function for a virtual line source using the 2-dimensional ESA @@ -290,7 +300,7 @@ def line_2d_with_dipole_ssd(omega, x0, xs, alpha=3/2*np.pi, Nc=None, c=None): return -1j*np.pi/alpha * d -def point_25d(omega, x0, xs, xref=[2, -2, 0], alpha=3/2*np.pi, +def point_25d_edge(omega, x0, xs, xref=[2, -2, 0], alpha=3/2*np.pi, Nc=None, c=None): r"""Driving function for 2.5-dimensional point source with edge ESA. @@ -344,5 +354,5 @@ def point_25d(omega, x0, xs, xref=[2, -2, 0], alpha=3/2*np.pi, else: a = np.linalg.norm(xref-x0, axis=1)/np.linalg.norm(xref-xs) - d, selection, _ = line_2d(omega, x0, xs, alpha=alpha, Nc=Nc, c=c) + d, selection, _ = line_2d_edge(omega, x0, xs, alpha=alpha, Nc=Nc, c=c) return 1j*np.sqrt(a) * d, selection, secondary_source_point(omega, c) From fad204d33db934baedec8f433f055581f89ce9ae Mon Sep 17 00:00:00 2001 From: Frank Schultz Date: Wed, 13 Mar 2019 13:58:47 +0100 Subject: [PATCH 14/15] Update __init__.py --- sfs/mono/__init__.py | 52 +++++++++++++++++++++----------------------- 1 file changed, 25 insertions(+), 27 deletions(-) diff --git a/sfs/mono/__init__.py b/sfs/mono/__init__.py index 13fee05..6961969 100644 --- a/sfs/mono/__init__.py +++ b/sfs/mono/__init__.py @@ -12,34 +12,8 @@ """ from . import source - - -import numpy as _np from .. import array as _array - - -def secondary_source_point(omega, c): - """Create a point source for use in `sfs.mono.synthesize()`.""" - - def secondary_source(position, _, grid): - return source.point(omega, position, grid, c) - - return secondary_source - - -def secondary_source_line(omega, c): - """Create a line source for use in `sfs.mono.synthesize()`.""" - - def secondary_source(position, _, grid): - return source.line(omega, position, grid, c) - - return secondary_source - - -from . import esa -from . import nfchoa -from . import sdm -from . import wfs +import numpy as _np def shiftphase(p, phase): @@ -83,3 +57,27 @@ def synthesize(d, weights, ssd, secondary_source_function, **kwargs): if weight != 0: p += a * weight * d * secondary_source_function(x, n, **kwargs) return p + + +def secondary_source_point(omega, c): + """Create a point source for use in `sfs.mono.synthesize()`.""" + + def secondary_source(position, _, grid): + return source.point(omega, position, grid, c) + + return secondary_source + + +def secondary_source_line(omega, c): + """Create a line source for use in `sfs.mono.synthesize()`.""" + + def secondary_source(position, _, grid): + return source.line(omega, position, grid, c) + + return secondary_source + + +from . import esa +from . import nfchoa +from . import sdm +from . import wfs From 5eeaecc767e6aae927979970ac196b3ef80ca708 Mon Sep 17 00:00:00 2001 From: Frank Schultz Date: Wed, 13 Mar 2019 14:13:23 +0100 Subject: [PATCH 15/15] Update esa.py --- sfs/mono/esa.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/sfs/mono/esa.py b/sfs/mono/esa.py index a6e5427..5584761 100644 --- a/sfs/mono/esa.py +++ b/sfs/mono/esa.py @@ -1,12 +1,12 @@ """Compute ESA driving functions for various systems. - ESA is abbreviation for equivalent scattering approach. +ESA is abbreviation for equivalent scattering approach. - ESA driving functions for an edge-shaped SSD are provided below. - Further ESA for different geometries might be added here. +ESA driving functions for an edge-shaped SSD are provided below. +Further ESA for different geometries might be added here. - Note that mode-matching (such as NFC-HOA, SDM) are equivalent - to ESA in their specific geometries (spherical/circular, planar/linear). +Note that mode-matching (such as NFC-HOA, SDM) are equivalent +to ESA in their specific geometries (spherical/circular, planar/linear). """