PKVG2܁sfs/__init__.py"""Sound Field Synthesis Toolbox. http://sfs.rtfd.org/ """ __version__ = "0.2.0" from . import tapering from . import array from . import util from . import defs try: from . import plot except ImportError: pass from . import mono PK+j~G4%:: sfs/array.py"""Compute positions of various secondary source distributions. .. plot:: :context: reset import sfs import matplotlib.pyplot as plt plt.rcParams['figure.figsize'] = 8, 4.5 # inch plt.rcParams['axes.grid'] = True .. autoclass:: ArrayData :members: take """ from __future__ import division # for Python 2.x from collections import namedtuple import numpy as np from . import util class ArrayData(namedtuple('ArrayData', 'x n a')): """Named tuple returned by array functions. See :obj:`collections.namedtuple`. Attributes ---------- x : (N, 3) numpy.ndarray Positions of secondary sources n : (N, 3) numpy.ndarray Orientations (normal vectors) of secondary sources a : (N,) numpy.ndarray Weights of secondary sources """ __slots__ = () def __repr__(self): return 'ArrayData(\n' + ',\n'.join( ' {0}={1}'.format(name, repr(data).replace('\n', '\n ')) for name, data in zip('xna', self)) + ')' def take(self, indices): """Return a sub-array given by `indices`.""" return ArrayData(self.x[indices], self.n[indices], self.a[indices]) def linear(N, spacing, center=[0, 0, 0], orientation=[1, 0, 0]): """Linear secondary source distribution. Parameters ---------- N : int Number of secondary sources. spacing : float Distance (in metres) between secondary sources. center : (3,) array_like, optional Coordinates of array center. orientation : (3,) array_like, optional Orientation of the array. By default, the loudspeakers have their main axis pointing into positive x-direction. Returns ------- ArrayData Positions, orientations and weights of secondary sources. See :class:`ArrayData`. Example ------- .. plot:: :context: close-figs x0, n0, a0 = sfs.array.linear(16, 0.2, orientation=[0, -1, 0]) sfs.plot.loudspeaker_2d(x0, n0, a0) plt.axis('equal') """ return _linear_helper(np.arange(N) * spacing, center, orientation) def linear_diff(distances, center=[0, 0, 0], orientation=[1, 0, 0]): """Linear secondary source distribution from a list of distances. Parameters ---------- distances : (N-1,) array_like Sequence of secondary sources distances in metres. center, orientation See :func:`linear` Returns ------- ArrayData Positions, orientations and weights of secondary sources. See :class:`ArrayData`. Example ------- .. plot:: :context: close-figs x0, n0, a0 = sfs.array.linear_diff(4 * [0.3] + 6 * [0.15] + 4 * [0.3], orientation=[0, -1, 0]) sfs.plot.loudspeaker_2d(x0, n0, a0) plt.axis('equal') """ distances = util.asarray_1d(distances) ycoordinates = np.concatenate(([0], np.cumsum(distances))) return _linear_helper(ycoordinates, center, orientation) def linear_random(N, min_spacing, max_spacing, center=[0, 0, 0], orientation=[1, 0, 0], seed=None): """Randomly sampled linear array. Parameters ---------- N : int Number of secondary sources. min_spacing, max_spacing : float Minimal and maximal distance (in metres) between secondary sources. center, orientation See :func:`linear` seed : {None, int, array_like}, optional Random seed. See :class:`numpy.random.RandomState`. Returns ------- ArrayData Positions, orientations and weights of secondary sources. See :class:`ArrayData`. Example ------- .. plot:: :context: close-figs x0, n0, a0 = sfs.array.linear_random(12, 0.15, 0.4, orientation=[0, -1, 0]) sfs.plot.loudspeaker_2d(x0, n0, a0) plt.axis('equal') """ r = np.random.RandomState(seed) distances = r.uniform(min_spacing, max_spacing, size=N-1) return linear_diff(distances, center, orientation) def circular(N, R, center=[0, 0, 0]): """Circular secondary source distribution parallel to the xy-plane. Parameters ---------- N : int Number of secondary sources. R : float Radius in metres. center See :func:`linear`. Returns ------- ArrayData Positions, orientations and weights of secondary sources. See :class:`ArrayData`. Example ------- .. plot:: :context: close-figs x0, n0, a0 = sfs.array.circular(16, 1) sfs.plot.loudspeaker_2d(x0, n0, a0, size=0.2, show_numbers=True) plt.axis('equal') """ center = util.asarray_1d(center) alpha = np.linspace(0, 2 * np.pi, N, endpoint=False) positions = np.zeros((N, len(center))) positions[:, 0] = R * np.cos(alpha) positions[:, 1] = R * np.sin(alpha) positions += center normals = np.zeros_like(positions) normals[:, 0] = np.cos(alpha + np.pi) normals[:, 1] = np.sin(alpha + np.pi) weights = np.ones(N) * 2 * np.pi * R / N return ArrayData(positions, normals, weights) def rectangular(N, spacing, center=[0, 0, 0], orientation=[1, 0, 0]): """Rectangular secondary source distribution. Parameters ---------- N : int or pair of int Number of secondary sources on each side of the rectangle. If a pair of numbers is given, the first one specifies the first and third segment, the second number specifies the second and fourth segment. spacing : float Distance (in metres) between secondary sources. center, orientation See :func:`linear`. The `orientation` corresponds to the first linear segment. Returns ------- ArrayData Positions, orientations and weights of secondary sources. See :class:`ArrayData`. Example ------- .. plot:: :context: close-figs x0, n0, a0 = sfs.array.rectangular((4, 8), 0.2) sfs.plot.loudspeaker_2d(x0, n0, a0, show_numbers=True) plt.axis('equal') """ N1, N2 = (N, N) if np.isscalar(N) else N offset1 = spacing * (N2 - 1) / 2 + spacing / np.sqrt(2) offset2 = spacing * (N1 - 1) / 2 + spacing / np.sqrt(2) positions, normals, weights = concatenate( linear(N1, spacing, [-offset1, 0, 0], [1, 0, 0]), # left linear(N2, spacing, [0, offset2, 0], [0, -1, 0]), # upper linear(N1, spacing, [offset1, 0, 0], [-1, 0, 0]), # right linear(N2, spacing, [0, -offset2, 0], [0, 1, 0]), # lower ) positions, normals = _rotate_array(positions, normals, [1, 0, 0], orientation) positions += center return ArrayData(positions, normals, weights) def rounded_edge(Nxy, Nr, dx, center=[0, 0, 0], orientation=[1, 0, 0]): """Array along the xy-axis with rounded edge at the origin. Parameters ---------- Nxy : int Number of secondary sources along x- and y-axis. Nr : int Number of secondary sources in rounded edge. Radius of edge is adjusted to equdistant sampling along entire array. center : (3,) array_like, optional Position of edge. orientation : (3,) array_like, optional Normal vector of array. Default orientation is along xy-axis. Returns ------- ArrayData Positions, orientations and weights of secondary sources. See :class:`ArrayData`. Example ------- .. plot:: :context: close-figs x0, n0, a0 = sfs.array.rounded_edge(8, 5, 0.2) sfs.plot.loudspeaker_2d(x0, n0, a0) plt.axis('equal') """ # radius of rounded edge Nr += 1 R = 2/np.pi * Nr * dx # array along y-axis x00, n00, a00 = linear(Nxy, dx, center=[0, Nxy//2*dx+dx/2+R, 0]) x00 = np.flipud(x00) positions = x00 directions = n00 weights = a00 # round part x00 = np.zeros((Nr, 3)) n00 = np.zeros((Nr, 3)) a00 = np.zeros(Nr) for n in range(0, Nr): alpha = np.pi/2 * n/Nr x00[n, 0] = R * (1-np.cos(alpha)) x00[n, 1] = R * (1-np.sin(alpha)) n00[n, 0] = np.cos(alpha) n00[n, 1] = np.sin(alpha) a00[n] = dx positions = np.concatenate((positions, x00)) directions = np.concatenate((directions, n00)) weights = np.concatenate((weights, a00)) # array along x-axis x00, n00, a00 = linear(Nxy, dx, center=[Nxy//2*dx-dx/2+R, 0, 0], orientation=[0, 1, 0]) x00 = np.flipud(x00) positions = np.concatenate((positions, x00)) directions = np.concatenate((directions, n00)) weights = np.concatenate((weights, a00)) # rotate array positions, directions = _rotate_array(positions, directions, [1, 0, 0], orientation) # shift array to desired position positions += center return ArrayData(positions, directions, weights) def planar(N, spacing, center=[0, 0, 0], orientation=[1, 0, 0]): """Planar secondary source distribtion. Parameters ---------- N : int or pair of int Number of secondary sources along each edge. If a pair of numbers is given, the first one specifies the number on the horizontal edge, the second one specifies the number on the vertical edge. spacing : float Distance (in metres) between secondary sources. center, orientation See :func:`linear`. Returns ------- ArrayData Positions, orientations and weights of secondary sources. See :class:`ArrayData`. """ N1, N2 = (N, N) if np.isscalar(N) else N zcoordinates = np.arange(N2) * spacing zcoordinates -= np.mean(zcoordinates[[0, -1]]) # move center to origin subarrays = [linear(N1, spacing, center=[0, 0, z]) for z in zcoordinates] positions, normals, weights = concatenate(*subarrays) weights *= spacing positions, normals = _rotate_array(positions, normals, [1, 0, 0], orientation) positions += center return ArrayData(positions, normals, weights) def cube(N, spacing, center=[0, 0, 0], orientation=[1, 0, 0]): """Cube-shaped secondary source distribtion. Parameters ---------- N : int or triple of int Number of secondary sources along each edge. If a triple of numbers is given, the first two specify the edges like in :func:`rectangular`, the last one specifies the vertical edge. spacing : float Distance (in metres) between secondary sources. center, orientation See :func:`linear`. The `orientation` corresponds to the first planar segment. Returns ------- ArrayData Positions, orientations and weights of secondary sources. See :class:`ArrayData`. """ N1, N2, N3 = (N, N, N) if np.isscalar(N) else N offset1 = spacing * (N2 - 1) / 2 + spacing / np.sqrt(2) offset2 = spacing * (N1 - 1) / 2 + spacing / np.sqrt(2) offset3 = spacing * (N3 - 1) / 2 + spacing / np.sqrt(2) positions, directions, weights = concatenate( planar((N1, N3), spacing, [-offset1, 0, 0], [1, 0, 0]), # west planar((N2, N3), spacing, [0, offset2, 0], [0, -1, 0]), # north planar((N1, N3), spacing, [offset1, 0, 0], [-1, 0, 0]), # east planar((N2, N3), spacing, [0, -offset2, 0], [0, 1, 0]), # south planar((N2, N1), spacing, [0, 0, -offset3], [0, 0, 1]), # bottom planar((N2, N1), spacing, [0, 0, offset3], [0, 0, -1]), # top ) positions, directions = _rotate_array(positions, directions, [1, 0, 0], orientation) positions += center return ArrayData(positions, directions, weights) def sphere_load(fname, radius, center=[0, 0, 0]): """Spherical secondary source distribution loaded from datafile. ASCII Format (see MATLAB SFS Toolbox) with 4 numbers (3 position, 1 weight) per secondary source located on the unit circle. Returns ------- ArrayData Positions, orientations and weights of secondary sources. See :class:`ArrayData`. """ data = np.loadtxt(fname) positions, weights = data[:, :3], data[:, 3] normals = -positions positions *= radius positions += center return ArrayData(positions, normals, weights) def load(fname, center=[0, 0, 0], orientation=[1, 0, 0]): """Load secondary source positions from datafile. Comma Seperated Values (CSV) format with 7 values (3 positions, 3 normal vectors, 1 weight) per secondary source. Returns ------- ArrayData Positions, orientations and weights of secondary sources. See :class:`ArrayData`. """ data = np.loadtxt(fname, delimiter=',') positions, normals, weights = data[:, :3], data[:, 3:6], data[:, 6] positions, normals = _rotate_array(positions, normals, [1, 0, 0], orientation) positions += center return ArrayData(positions, normals, weights) def weights_midpoint(positions, closed): """Calculate loudspeaker weights for a simply connected array. The weights are calculated according to the midpoint rule. Parameters ---------- positions : (N, 3) array_like Sequence of secondary source positions. .. note:: The loudspeaker positions have to be ordered on the contour! closed : bool ``True`` if the loudspeaker contour is closed. Returns ------- (N,) numpy.ndarray Weights of secondary sources. """ positions = util.asarray_of_rows(positions) if closed: before, after = -1, 0 # cyclic else: before, after = 1, -2 # mirrored positions = np.row_stack((positions[before], positions, positions[after])) distances = np.linalg.norm(np.diff(positions, axis=0), axis=1) return (distances[:-1] + distances[1:]) / 2 def _rotate_array(positions, normals, n1, n2): """Rotate secondary sources from n1 to n2.""" R = util.rotation_matrix(n1, n2) positions = np.inner(positions, R) normals = np.inner(normals, R) return positions, normals def _linear_helper(ycoordinates, center, orientation): """Create a full linear array from an array of y-coordinates.""" N = len(ycoordinates) positions = np.zeros((N, 3)) positions[:, 1] = ycoordinates - np.mean(ycoordinates[[0, -1]]) positions, normals = _rotate_array(positions, [1, 0, 0], [1, 0, 0], orientation) positions += center normals = np.tile(normals, (N, 1)) weights = weights_midpoint(positions, closed=False) return ArrayData(positions, normals, weights) def concatenate(*arrays): """Concatenate :class:`ArrayData` objects.""" return ArrayData._make(np.concatenate(i) for i in zip(*arrays)) PKG)) sfs/util.py"""Various utility functions.""" import numpy as np from . import defs def rotation_matrix(n1, n2): """Compute rotation matrix for rotation from `n1` to `n2`. Parameters ---------- n1, n2 : (3,) array_like Two vectors. They don't have to be normalized. Returns ------- (3, 3) numpy.ndarray Rotation matrix. """ n1 = normal_vector(n1) n2 = normal_vector(n2) I = np.identity(3) if np.all(n1 == n2): return I # no rotation elif np.all(n1 == -n2): return -I # flip # TODO: check for *very close to* parallel vectors # Algorithm from http://math.stackexchange.com/a/476311 v = v0, v1, v2 = np.cross(n1, n2) s = np.linalg.norm(v) # sine c = np.inner(n1, n2) # cosine vx = np.matrix([[0, -v2, v1], [v2, 0, -v0], [-v1, v0, 0]]) # skew-symmetric cross-product matrix return I + vx + vx**2 * (1 - c) / s**2 def wavenumber(omega, c=None): """Compute the wavenumber for a given radial frequency.""" if c is None: c = defs.c return omega / c def direction_vector(alpha, beta=np.pi/2): """Compute normal vector from azimuth, colatitude.""" return sph2cart(alpha, beta, 1) def sph2cart(alpha, beta, r): """Spherical to cartesian coordinates.""" x = r * np.cos(alpha) * np.sin(beta) y = r * np.sin(alpha) * np.sin(beta) z = r * np.cos(beta) return x, y, z def cart2sph(x, y, z): """Cartesian to spherical coordinates.""" alpha = np.arctan2(y, x) beta = np.arccos(z / np.sqrt(x**2 + y**2)) r = np.sqrt(x**2 + y**2 + z**2) return alpha, beta, r def asarray_1d(a, **kwargs): """Squeeze the input and check if the result is one-dimensional. Returns `a` converted to a :class:`numpy.ndarray` and stripped of all singleton dimensions. Scalars are "upgraded" to 1D arrays. The result must have exactly one dimension. If not, an error is raised. """ result = np.squeeze(np.asarray(a, **kwargs)) if result.ndim == 0: result = result.reshape((1,)) elif result.ndim > 1: raise ValueError("array must be one-dimensional") return result def asarray_of_rows(a, **kwargs): """Convert to 2D array, turn column vector into row vector. Returns `a` converted to a :class:`numpy.ndarray` and stripped of all singleton dimensions. If the result has exactly one dimension, it is re-shaped into a 2D row vector. """ result = np.squeeze(np.asarray(a, **kwargs)) if result.ndim == 1: result = result.reshape(1, -1) return result def strict_arange(start, stop, step=1, endpoint=False, dtype=None, **kwargs): """Like :func:`numpy.arange`, but compensating numeric errors. Unlike :func:`numpy.arange`, but similar to :func:`numpy.linspace`, providing `endpoint=True` includes both endpoints. Parameters ---------- start, stop, step, dtype See :func:`numpy.arange`. endpoint See :func:`numpy.linspace`. .. note:: With `endpoint=True`, the difference between `start` and `end` value must be an integer multiple of the corresponding `spacing` value! **kwargs All further arguments are forwarded to :func:`numpy.isclose`. Returns ------- numpy.ndarray Array of evenly spaced values. See :func:`numpy.arange`. """ remainder = (stop - start) % step if np.any(np.isclose(remainder, (0.0, step), **kwargs)): if endpoint: stop += step * 0.5 else: stop -= step * 0.5 elif endpoint: raise ValueError("Invalid stop value for endpoint=True") return np.arange(start, stop, step, dtype) def xyz_grid(x, y, z, spacing, endpoint=True, **kwargs): """Create a grid with given range and spacing. Parameters ---------- x, y, z : float or pair of float Inclusive range of the respective coordinate or a single value if only a slice along this dimension is needed. spacing : float or triple of float Grid spacing. If a single value is specified, it is used for all dimensions, if multiple values are given, one value is used per dimension. If a dimension (`x`, `y` or `z`) has only a single value, the corresponding spacing is ignored. endpoint : bool, optional If ``True`` (the default), the endpoint of each range is included in the grid. Use ``False`` to get a result similar to :func:`numpy.arange`. See :func:`sfs.util.strict_arange`. **kwargs All further arguments are forwarded to :func:`sfs.util.strict_arange`. Returns ------- XyzComponents A grid that can be used for sound field calculations. See :class:`sfs.util.XyzComponents`. See Also -------- strict_arange, numpy.meshgrid """ if np.isscalar(spacing): spacing = [spacing] * 3 ranges = [] scalars = [] for i, coord in enumerate([x, y, z]): if np.isscalar(coord): scalars.append((i, coord)) else: start, stop = coord ranges.append(strict_arange(start, stop, spacing[i], endpoint=endpoint, **kwargs)) grid = np.meshgrid(*ranges, sparse=True, copy=False) for i, s in scalars: grid.insert(i, s) return XyzComponents(grid) def normalize(p, grid, xnorm): """Normalize sound field wrt position `xnorm`.""" return p / np.abs(probe(p, grid, xnorm)) def probe(p, grid, x): """Determine the value at position `x` in the sound field `p`.""" grid = XyzComponents(grid) x = asarray_1d(x) r = np.linalg.norm(grid - x) idx = np.unravel_index(r.argmin(), r.shape) return p[idx] def broadcast_zip(*args): """Broadcast arguments to the same shape and then use :func:`zip`.""" return zip(*np.broadcast_arrays(*args)) def normal_vector(x): """Normalize a 1D vector.""" x = asarray_1d(x) return x / np.linalg.norm(x) def displacement(v, omega): """Particle displacement .. math:: d(x, t) = \int_0^t v(x, t) dt """ v = XyzComponents(v) return v / (1j * omega) def db(x, power=False): """Convert `x` to decibel. Parameters ---------- x : array_like Input data. Values of 0 lead to negative infinity. power : bool, optional If `power=False` (the default), `x` is squared before conversion. """ with np.errstate(divide='ignore'): return 10 if power else 20 * np.log10(np.abs(x)) class XyzComponents(np.ndarray): """See __init__().""" def __init__(self, components, **kwargs): """Triple (or pair) of arrays: x, y, and optionally z. Instances of this class can be used to store coordinate grids (either regular grids like in :func:`sfs.util.xyz_grid` or arbitrary point clouds) or vector fields (e.g. particle velocity). This class is a subclass of :class:`numpy.ndarray`. It is one-dimensional (like a plain :class:`list`) and has a length of 3 (or 2, if no z-component is available). It uses `dtype=object` in order to be able to store other `numpy.ndarray`\s of arbitrary shapes but also scalars, if needed. Because it is a NumPy array subclass, it can be used in operations with scalars and "normal" NumPy arrays, as long as they have a compatible shape. Like any NumPy array, instances of this class are iterable and can be used, e.g., in for-loops and tuple unpacking. If slicing or broadcasting leads to an incompatible shape, a plain `numpy.ndarray` with `dtype=object` is returned. Parameters ---------- components : triple or pair of array_like The values to be used as X, Y and Z arrays. Z is optional. **kwargs All further arguments are forwarded to :func:`numpy.asarray`, which is applied to the elements of `components`. """ # This method does nothing, it's only here for the documentation! def __new__(cls, components, **kwargs): # object arrays cannot be created and populated in a single step: obj = np.ndarray.__new__(cls, len(components), dtype=object) for i, component in enumerate(components): obj[i] = np.asarray(component, **kwargs) return obj def __array_finalize__(self, obj): if self.ndim == 0: pass # this is allowed, e.g. for np.inner() elif self.ndim > 1 or len(self) not in (2, 3): raise ValueError("XyzComponents can only have 2 or 3 components") def __array_prepare__(self, obj, context=None): if obj.ndim == 1 and len(obj) in (2, 3): return obj.view(XyzComponents) return obj def __array_wrap__(self, obj, context=None): if obj.ndim != 1 or len(obj) not in (2, 3): return obj.view(np.ndarray) return obj def __getitem__(self, index): if isinstance(index, slice): start, stop, step = index.indices(len(self)) if start == 0 and stop in (2, 3) and step == 1: return np.ndarray.__getitem__(self, index) # Slices other than xy and xyz are "downgraded" to ndarray return np.ndarray.__getitem__(self.view(np.ndarray), index) def __repr__(self): return 'XyzComponents(\n' + ',\n'.join( ' {0}={1}'.format(name, repr(data).replace('\n', '\n ')) for name, data in zip('xyz', self)) + ')' def make_property(index, doc): def getter(self): return self[index] def setter(self, value): self[index] = value return property(getter, setter, doc=doc) x = make_property(0, doc='x-component.') y = make_property(1, doc='y-component.') z = make_property(2, doc='z-component (optional).') del make_property def apply(self, func, *args, **kwargs): """Apply a function to each component. The function `func` will be called once for each component, passing the current component as first argument. All further arguments are passed after that. The results are returned as a new :class:`XyzComponents` object. """ return XyzComponents([func(i, *args, **kwargs) for i in self]) PKGd711 sfs/plot.py"""Plot sound fields etc.""" from __future__ import division import matplotlib.pyplot as plt from matplotlib.patches import PathPatch from matplotlib.path import Path from matplotlib.collections import PatchCollection from mpl_toolkits import axes_grid1 from mpl_toolkits.mplot3d import Axes3D import numpy as np from . import util def _register_coolwarm_clip(alpha): """Create color map with "over" and "under" values. The 'coolwarm' colormap is based on the paper "Diverging Color Maps for Scientific Visualization" by Kenneth Moreland http://www.sandia.gov/~kmorel/documents/ColorMaps/ """ from matplotlib.colors import LinearSegmentedColormap cdict = plt.cm.datad['coolwarm'] cmap = LinearSegmentedColormap('coolwarm_clip', cdict) cmap.set_over([alpha * c + 1 - alpha for c in cmap(1.0)[:3]]) cmap.set_under([alpha * c + 1 - alpha for c in cmap(0.0)[:3]]) plt.cm.register_cmap(cmap=cmap) _register_coolwarm_clip(0.7) del _register_coolwarm_clip def virtualsource_2d(xs, ns=None, type='point', ax=None): """Draw position/orientation of virtual source.""" xs = np.asarray(xs) ns = np.asarray(ns) if ax is None: ax = plt.gca() if type == 'point': vps = plt.Circle(xs, .05, edgecolor='k', facecolor='k') ax.add_artist(vps) for n in range(1, 3): vps = plt.Circle(xs, .05+n*0.05, edgecolor='k', fill=False) ax.add_artist(vps) elif type == 'plane': ns = 0.2 * ns ax.arrow(xs[0], xs[1], ns[0], ns[1], head_width=0.05, head_length=0.1, fc='k', ec='k') def reference_2d(xref, size=0.1, ax=None): """Draw reference/normalization point.""" xref = np.asarray(xref) if ax is None: ax = plt.gca() ax.plot((xref[0]-size, xref[0]+size), (xref[1]-size, xref[1]+size), 'k-') ax.plot((xref[0]-size, xref[0]+size), (xref[1]+size, xref[1]-size), 'k-') def secondarysource_2d(x0, n0, grid=None): """Simple plot of secondary source locations.""" x0 = np.asarray(x0) n0 = np.asarray(n0) ax = plt.gca() # plot only secondary sources inside simulated area if grid is not None: x0, n0 = _visible_secondarysources_2d(x0, n0, grid) # plot symbols for x00 in x0: ss = plt.Circle(x00[0:2], .05, edgecolor='k', facecolor='k') ax.add_artist(ss) def loudspeaker_2d(x0, n0, a0=0.5, size=0.08, show_numbers=False, grid=None, ax=None): """Draw loudspeaker symbols at given locations and angles. Parameters ---------- x0 : (N, 3) array_like Loudspeaker positions. n0 : (N, 3) or (3,) array_like Normal vector(s) of loudspeakers. a0 : float or (N,) array_like, optional Weighting factor(s) of loudspeakers. size : float, optional Size of loudspeakers in metres. show_numbers : bool, optional If ``True``, loudspeaker numbers are shown. grid : triple of numpy.ndarray, optional If specified, only loudspeakers within the `grid` are shown. ax : Axes object, optional The loudspeakers are plotted into this :class:`~matplotlib.axes.Axes` object or -- if not specified -- into the current axes. """ x0 = util.asarray_of_rows(x0) n0 = util.asarray_of_rows(n0) a0 = util.asarray_1d(a0).reshape(-1, 1) # plot only secondary sources inside simulated area if grid is not None: x0, n0 = _visible_secondarysources_2d(x0, n0, grid) # normalized coordinates of loudspeaker symbol (see IEC 60617-9) codes, coordinates = zip(*( (Path.MOVETO, [-0.62, 0.21]), (Path.LINETO, [-0.31, 0.21]), (Path.LINETO, [0, 0.5]), (Path.LINETO, [0, -0.5]), (Path.LINETO, [-0.31, -0.21]), (Path.LINETO, [-0.62, -0.21]), (Path.CLOSEPOLY, [0, 0]), (Path.MOVETO, [-0.31, 0.21]), (Path.LINETO, [-0.31, -0.21]), )) coordinates = np.column_stack([coordinates, np.zeros(len(coordinates))]) coordinates *= size patches = [] for x00, n00 in util.broadcast_zip(x0, n0): # rotate and translate coordinates R = util.rotation_matrix([1, 0, 0], n00) transformed_coordinates = np.inner(coordinates, R) + x00 patches.append(PathPatch(Path(transformed_coordinates[:, :2], codes))) # add collection of patches to current axis p = PatchCollection(patches, edgecolor='0', facecolor=np.tile(1 - a0, 3)) if ax is None: ax = plt.gca() ax.add_collection(p) if show_numbers: for idx, (x00, n00) in enumerate(util.broadcast_zip(x0, n0)): x, y = x00[:2] - 1.2 * size * n00[:2] ax.text(x, y, idx + 1, horizontalalignment='center', verticalalignment='center') def _visible_secondarysources_2d(x0, n0, grid): """Determine secondary sources which lie within `grid`.""" grid = util.XyzComponents(grid) x, y = grid[:2] idx = np.where((x0[:, 0] > x.min()) & (x0[:, 0] < x.max()) & (x0[:, 1] > y.min()) & (x0[:, 1] < x.max())) idx = np.squeeze(idx) return x0[idx, :], n0[idx, :] def loudspeaker_3d(x0, n0, a0=None, w=0.08, h=0.08): """Plot positions and normals of a 3D secondary source distribution.""" fig = plt.figure(figsize=(15, 15)) ax = fig.add_subplot(111, projection='3d') ax.quiver(x0[:, 0], x0[:, 1], x0[:, 2], n0[:, 0], n0[:, 1], n0[:, 2], length=0.1) plt.xlabel('x (m)') plt.ylabel('y (m)') plt.title('Secondary Sources') fig.show() def soundfield(p, grid, xnorm=None, cmap='coolwarm_clip', vmin=-2.0, vmax=2.0, xlabel=None, ylabel=None, colorbar=True, colorbar_kwargs={}, ax=None, **kwargs): """Two-dimensional plot of sound field. Parameters ---------- p : array_like Sound pressure values (or any other scalar quantity if you like). If the values are complex, the imaginary part is ignored. Typically, `p` is two-dimensional with a shape of `(Ny, Nx)`, `(Nz, Nx)` or `(Nz, Ny)`. This is the case if :func:`sfs.util.xyz_grid` was used with a single number for `z`, `y` or `x`, respectively. However, `p` can also be three-dimensional with a shape of `(Ny, Nx, 1)`, `(1, Nx, Nz)` or `(Ny, 1, Nz)`. This is the case if :func:`numpy.meshgrid` was used with a scalar for `z`, `y` or `x`, respectively (and of course with the default ``indexing='xy'``). .. note:: If you want to plot a single slice of a pre-computed "full" 3D sound field, make sure that the slice still has three dimensions (including one singleton dimension). This way, you can use the original `grid` of the full volume without changes. This works because the grid component corresponding to the singleton dimension is simply ignored. grid : triple or pair of numpy.ndarray The grid that was used to calculate `p`, see :func:`sfs.util.xyz_grid`. If `p` is two-dimensional, but `grid` has 3 components, one of them must be scalar. xnorm : array_like, optional Coordinates of a point to which the sound field should be normalized before plotting. If not specified, no normalization is used. See :func:`sfs.util.normalize`. Returns ------- AxesImage See :func:`matplotlib.pyplot.imshow`. Other Parameters ---------------- xlabel, ylabel : str Overwrite default x/y labels. Use ``xlabel=''`` and ``ylabel=''`` to remove x/y labels. The labels can be changed afterwards with :func:`matplotlib.pyplot.xlabel` and :func:`matplotlib.pyplot.ylabel`. colorbar : bool, optional If ``False``, no colorbar is created. colorbar_kwargs : dict, optional Further colorbar arguments, see :func:`add_colorbar`. ax : Axes, optional If given, the plot is created on `ax` instead of the current axis (see :func:`matplotlib.pyplot.gca`). cmap, vmin, vmax, **kwargs All further parameters are forwarded to :func:`matplotlib.pyplot.imshow`. See Also -------- sfs.plot.level """ p = np.asarray(p) grid = util.XyzComponents(grid) # normalize sound field wrt xnorm if xnorm is not None: p = util.normalize(p, grid, xnorm) if p.ndim == 3: if p.shape[2] == 1: p = p[:, :, 0] # first axis: y; second axis: x plotting_plane = 'xy' elif p.shape[1] == 1: p = p[:, 0, :].T # first axis: z; second axis: y plotting_plane = 'yz' elif p.shape[0] == 1: p = p[0, :, :].T # first axis: z; second axis: x plotting_plane = 'xz' else: raise ValueError("If p is 3D, one dimension must have length 1") elif len(grid) == 3: if grid[2].ndim == 0: plotting_plane = 'xy' elif grid[1].ndim == 0: plotting_plane = 'xz' elif grid[0].ndim == 0: plotting_plane = 'yz' else: raise ValueError( "If p is 2D and grid is 3D, one grid component must be scalar") else: # 2-dimensional case plotting_plane = 'xy' if plotting_plane == 'xy': x, y = grid[[0, 1]] elif plotting_plane == 'xz': x, y = grid[[0, 2]] elif plotting_plane == 'yz': x, y = grid[[1, 2]] if ax is None: ax = plt.gca() im = ax.imshow(np.real(p), cmap=cmap, origin='lower', extent=[x.min(), x.max(), y.min(), y.max()], vmax=vmax, vmin=vmin, **kwargs) if xlabel is None: xlabel = plotting_plane[0] + ' / m' if ylabel is None: ylabel = plotting_plane[1] + ' / m' ax.set_xlabel(xlabel) ax.set_ylabel(ylabel) if colorbar: add_colorbar(im, **colorbar_kwargs) return im def level(p, grid, xnorm=None, power=False, cmap=None, vmax=3, vmin=-50, **kwargs): """Two-dimensional plot of level (dB) of sound field. Takes the same parameters as :func:`sfs.plot.soundfield`. Other Parameters ---------------- power : bool, optional See :func:`sfs.util.db`. """ # normalize before converting to dB! if xnorm is not None: p = util.normalize(p, grid, xnorm) L = util.db(p, power=power) return soundfield(L, grid=grid, xnorm=None, cmap=cmap, vmax=vmax, vmin=vmin, **kwargs) def particles(x, trim=None, ax=None, xlabel='x (m)', ylabel='y (m)', edgecolor='', **kwargs): """Plot particle positions as scatter plot""" XX, YY = [np.real(c) for c in x[:2]] if trim is not None: xmin, xmax, ymin, ymax = trim idx = np.where((XX > xmin) & (XX < xmax) & (YY > ymin) & (YY < ymax)) XX = XX[idx] YY = YY[idx] if ax is None: ax = plt.gca() ax.scatter(XX, YY, edgecolor=edgecolor, **kwargs) if xlabel: ax.set_xlabel(xlabel) if ylabel: ax.set_ylabel(ylabel) def add_colorbar(im, aspect=20, pad=0.5, **kwargs): """Add a vertical color bar to a plot. Parameters ---------- im : ScalarMappable The output of :func:`sfs.plot.soundfield`, :func:`sfs.plot.level` or any other :class:`matplotlib.cm.ScalarMappable`. aspect : float, optional Aspect ratio of the colorbar. Strictly speaking, since the colorbar is vertical, it's actually the inverse of the aspect ratio. pad : float, optional Space between image plot and colorbar, as a fraction of the width of the colorbar. .. note:: The `pad` argument of :meth:`matplotlib.figure.Figure.colorbar` has a slightly different meaning ("fraction of original axes")! \**kwargs All further arguments are forwarded to :meth:`matplotlib.figure.Figure.colorbar`. See Also -------- matplotlib.pyplot.colorbar """ ax = im.axes divider = axes_grid1.make_axes_locatable(ax) width = axes_grid1.axes_size.AxesY(ax, aspect=1/aspect) pad = axes_grid1.axes_size.Fraction(pad, width) current_ax = plt.gca() cax = divider.append_axes("right", size=width, pad=pad) plt.sca(current_ax) return ax.figure.colorbar(im, cax=cax, orientation='vertical', **kwargs) PK+j~GAAi sfs/defs.py"""Definition of constants.""" # speed of sound c = 343 # static density of air rho0 = 1.2250 # tolerance used for secondary source selection selection_tolerance = 1e-6 PKhEVrsfs/tapering.py"""Weights (tapering) for the driving function.""" # from scipy import signal import numpy as np def none(active): """No tapering window.""" return np.asarray(active, dtype=np.float64) def kaiser(active): """Kaiser tapering window.""" idx = _windowidx(active) # compute coefficients window = np.zeros(active.shape) window[idx] = np.kaiser(len(idx), 2) return window def tukey(active, alpha): """Tukey tapering window.""" idx = _windowidx(active) # alpha out of limits if alpha <= 0 or alpha >= 1: return none(active) # design Tukey window x = np.linspace(0, 1, len(idx)+2) w = np.ones(x.shape) first_condition = x < alpha/2 w[first_condition] = 0.5 * (1 + np.cos(2*np.pi/alpha * (x[first_condition] - alpha/2))) third_condition = x >= (1 - alpha/2) w[third_condition] = 0.5 * (1 + np.cos(2*np.pi/alpha * (x[third_condition] - 1 + alpha/2))) # fit window into tapering function window = np.zeros(active.shape) window[idx] = w[1:-1] return window def _windowidx(active): """Returns list of connected indices for window function.""" active = np.asarray(active, dtype=np.float64) # find index were active loudspeakers begin (works for connected contours) if (active[0] == 1 and active[-1] == 0) or np.all(active): a0 = 0 else: a0 = np.argmax(np.diff(active)) + 1 # shift generic index vector to get a connected list of indices idx = np.roll(np.arange(len(active)), -a0) # remove indices of inactive secondary sources idx = idx[0:len(np.squeeze(np.where(active == 1)))] return idx PKz[^GVWMggsfs/mono/__init__.pyfrom . import drivingfunction from . import source from . import synthesized from . import soundfigure PK+j~G=ݺ " "sfs/mono/drivingfunction.py"""Compute driving functions for various systems.""" import numpy as np from numpy.core.umath_tests import inner1d # element-wise inner product from scipy.special import hankel2 from scipy.special import sph_jn, sph_yn from .. import util from .. import defs def wfs_2d_line(omega, x0, n0, xs, c=None): """Line source by 2-dimensional WFS. :: D(x0,k) = j k (x0-xs) n0 / |x0-xs| * H1(k |x0-xs|) """ x0 = np.asarray(x0) n0 = np.asarray(n0) xs = np.squeeze(np.asarray(xs)) k = util.wavenumber(omega, c) ds = x0 - xs r = np.linalg.norm(ds, axis=1) return 1j * k * inner1d(ds, n0) / r * hankel2(1, k * r) def _wfs_point(omega, x0, n0, xs, c=None): """Point source by two- or three-dimensional WFS. :: (x0-xs) n0 D(x0,k) = j k ------------- e^(-j k |x0-xs|) |x0-xs|^(3/2) """ x0 = np.asarray(x0) n0 = np.asarray(n0) xs = np.squeeze(np.asarray(xs)) k = util.wavenumber(omega, c) ds = x0 - xs r = np.linalg.norm(ds, axis=1) return 1j * k * inner1d(ds, n0) / r ** (3 / 2) * np.exp(-1j * k * r) wfs_2d_point = _wfs_point def wfs_25d_point(omega, x0, n0, xs, xref=[0, 0, 0], c=None, omalias=None): """Point source by 2.5-dimensional WFS. :: ____________ (x0-xs) n0 D(x0,k) = \|j k |xref-x0| ------------- e^(-j k |x0-xs|) |x0-xs|^(3/2) """ x0 = np.asarray(x0) n0 = np.asarray(n0) xs = np.squeeze(np.asarray(xs)) xref = np.squeeze(np.asarray(xref)) k = util.wavenumber(omega, c) ds = x0 - xs r = np.linalg.norm(ds, axis=1) return wfs_25d_preeq(omega, omalias, c) * \ np.sqrt(np.linalg.norm(xref - x0)) * inner1d(ds, n0) / \ r ** (3 / 2) * np.exp(-1j * k * r) wfs_3d_point = _wfs_point def _wfs_plane(omega, x0, n0, n=[0, 1, 0], c=None): """Plane wave by two- or three-dimensional WFS. Eq.(17) from [Spors et al, 2008]:: D(x0,k) = j k n n0 e^(-j k n x0) """ x0 = np.asarray(x0) n0 = np.asarray(n0) n = np.squeeze(np.asarray(n)) k = util.wavenumber(omega, c) return 2j * k * np.inner(n, n0) * np.exp(-1j * k * np.inner(n, x0)) wfs_2d_plane = _wfs_plane def wfs_25d_plane(omega, x0, n0, n=[0, 1, 0], xref=[0, 0, 0], c=None, omalias=None): """Plane wave by 2.5-dimensional WFS. :: ____________ D_2.5D(x0,w) = \|j k |xref-x0| n n0 e^(-j k n x0) """ x0 = np.asarray(x0) n0 = np.asarray(n0) n = np.squeeze(np.asarray(n)) xref = np.squeeze(np.asarray(xref)) k = util.wavenumber(omega, c) return wfs_25d_preeq(omega, omalias, c) * \ np.sqrt(2*np.pi * np.linalg.norm(xref - x0)) * \ np.inner(n, n0) * np.exp(-1j * k * np.inner(n, x0)) wfs_3d_plane = _wfs_plane def wfs_25d_preeq(omega, omalias, c): """Preqeualization for 2.5D WFS.""" 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): """Plane wave by simple delay of secondary sources.""" x0 = np.asarray(x0) n = np.squeeze(np.asarray(n)) k = util.wavenumber(omega, c) return np.exp(-1j * k * np.inner(n, x0)) def source_selection_plane(n0, n): """Secondary source selection for a plane wave. Eq.(13) from [Spors et al, 2008] """ n0 = np.asarray(n0) n = np.squeeze(np.asarray(n)) return np.inner(n, n0) >= defs.selection_tolerance def source_selection_point(n0, x0, xs): """Secondary source selection for a point source. Eq.(15) from [Spors et al, 2008] """ n0 = np.asarray(n0) x0 = np.asarray(x0) xs = np.squeeze(np.asarray(xs)) ds = x0 - xs return inner1d(ds, n0) >= defs.selection_tolerance def source_selection_all(N): """Select all secondary sources.""" return np.ones(N) >= 0 def nfchoa_2d_plane(omega, x0, r0, n=[0, 1, 0], c=None): """Point source by 2.5-dimensional WFS.""" x0 = np.asarray(x0) k = util.wavenumber(omega, c) alpha, beta, r = util.cart2sph(n[0], n[1], n[2]) alpha0, beta0, tmp = util.cart2sph(x0[:, 0], x0[:, 1], x0[:, 2]) # determine max order of circular harmonics M = _hoa_order_2d(len(x0)) # compute driving function d = 0 for m in np.arange(-M, M): d = d + 1j**(-m) / hankel2(m, k * r0) * \ np.exp(1j * m * (alpha0 - alpha)) return - 2j / (np.pi*r0) * d def nfchoa_25d_point(omega, x0, r0, xs, c=None): """Point source by 2.5-dimensional WFS. :: __ (2) 1 \ h|m| (w/c r) D(phi0,w) = ----- /__ ------------- e^(i m (phi0-phi)) 2pi r0 m=-N..N (2) h|m| (w/c r0) """ x0 = np.asarray(x0) k = util.wavenumber(omega, c) alpha, beta, r = util.cart2sph(xs[0], xs[1], xs[2]) alpha0, beta0, tmp = util.cart2sph(x0[:, 0], x0[:, 1], x0[:, 2]) # determine max order of circular harmonics M = _hoa_order_2d(len(x0)) # compute driving function d = 0 a = _sph_hn2(M, k * r) / _sph_hn2(M, k * r0) for m in np.arange(-M, M): d += a[0, abs(m)] * np.exp(1j * m * (alpha0 - alpha)) return 1 / (2 * np.pi * r0) * d def nfchoa_25d_plane(omega, x0, r0, n=[0, 1, 0], c=None): """Plane wave by 2.5-dimensional WFS. :: __ 2i \ i^|m| D_25D(phi0,w) = -- /__ ------------------ e^(i m (phi0-phi_pw) ) r0 m=-N..N (2) w/c h|m| (w/c r0) """ x0 = np.asarray(x0) k = util.wavenumber(omega, c) alpha, beta, r = util.cart2sph(n[0], n[1], n[2]) alpha0, beta0, tmp = util.cart2sph(x0[:, 0], x0[:, 1], x0[:, 2]) # determine max order of circular harmonics M = _hoa_order_2d(len(x0)) # compute driving function d = 0 a = 1 / _sph_hn2(M, k * r0) for m in np.arange(-M, M): d += (1j)**(-abs(m)) * a[0, abs(m)] * \ np.exp(1j * m * (alpha0 - alpha)) return -2 / r0 * d def sdm_2d_line(omega, x0, n0, xs, c=None): """Line source by two-dimensional SDM. The secondary sources have to be located on the x-axis (y0=0). Derived from [Spors 2009, 126th AES Convention], Eq.(9), Eq.(4):: D(x0,k) = """ x0 = np.asarray(x0) n0 = np.asarray(n0) xs = np.squeeze(np.asarray(xs)) k = util.wavenumber(omega, c) ds = x0 - xs r = np.linalg.norm(ds, axis=1) return - 1j/2 * k * xs[1] / r * hankel2(1, k * r) def sdm_2d_plane(omega, x0, n0, n=[0, 1, 0], c=None): """Plane wave by two-dimensional SDM. The secondary sources have to be located on the x-axis (y0=0). Derived from [Ahrens 2011, Springer], Eq.(3.73), Eq.(C.5), Eq.(C.11):: D(x0,k) = kpw,y * e^(-j*kpw,x*x) """ x0 = np.asarray(x0) n0 = np.asarray(n0) n = np.squeeze(np.asarray(n)) k = util.wavenumber(omega, c) return k * n[1] * np.exp(-1j * k * n[0] * x0[:, 0]) def sdm_25d_plane(omega, x0, n0, n=[0, 1, 0], xref=[0, 0, 0], c=None): """Plane wave by 2.5-dimensional SDM. The secondary sources have to be located on the x-axis (y0=0). Eq.(3.79) from [Ahrens 2011, Springer]:: D_2.5D(x0,w) = """ x0 = np.asarray(x0) n0 = np.asarray(n0) n = np.squeeze(np.asarray(n)) xref = np.squeeze(np.asarray(xref)) k = util.wavenumber(omega, c) return 4j * np.exp(-1j*k*n[1]*xref[1]) / hankel2(0, k*n[1]*xref[1]) * \ np.exp(-1j*k*n[0]*x0[:, 0]) def sdm_25d_point(omega, x0, n0, xs, xref=[0, 0, 0], c=None): """Point source by 2.5-dimensional SDM. The secondary sources have to be located on the x-axis (y0=0). Driving funcnction from [Spors 2010, 128th AES Covention], Eq.(24):: D(x0,k) = """ x0 = np.asarray(x0) n0 = np.asarray(n0) xs = np.squeeze(np.asarray(xs)) xref = np.squeeze(np.asarray(xref)) k = util.wavenumber(omega, c) ds = x0 - xs r = np.linalg.norm(ds, axis=1) return 1/2 * 1j * k * np.sqrt(xref[1] / (xref[1] - xs[1])) * \ xs[1] / r * hankel2(1, k * r) def _sph_hn2(n, z): """Spherical Hankel function of 2nd kind.""" return np.asarray(sph_jn(n, z)) - 1j * np.asarray(sph_yn(n, z)) def _hoa_order_2d(N): """Computes order of HOA.""" if N % 2 == 0: return N//2 else: return (N-1)//2 PK+j~Gƅ[sfs/mono/soundfigure.py"""Compute driving functions for sound figures.""" import numpy as np from .. import util from . import drivingfunction 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 a = drivingfunction.source_selection_plane(n0, npw) a = a * figure[n, m] d += a * drivingfunction.wfs_3d_plane(omega, x0, n0, npw, c) return d PK+j~Gxwahhsfs/mono/synthesized.py"""Computation of synthesized sound fields.""" import numpy as np from .source import point def generic(omega, x0, n0, d, grid, c=None, source=point): """Compute sound field for a generic driving function.""" d = np.squeeze(np.asarray(d)) if len(d) != len(x0): raise ValueError("length mismatch") p = 0 for weight, position, direction in zip(d, x0, n0): if weight != 0: p += weight * source(omega, position, direction, grid, c) return p def shiftphase(p, phase): """Shift phase of a sound field.""" p = np.asarray(p) return p * np.exp(1j * phase) PK+j~Gb;(;(sfs/mono/source.py"""Compute the sound field generated by a sound source. .. plot:: :context: reset import sfs import numpy as np import matplotlib.pyplot as plt plt.rcParams['figure.figsize'] = 8, 4.5 # inch x0 = 1.5, 1, 0 f = 500 # Hz omega = 2 * np.pi * f grid = sfs.util.xyz_grid([-2, 3], [-1, 2], 0, spacing=0.02) """ import itertools import numpy as np from scipy import special from .. import util from .. import defs def point(omega, x0, n0, grid, c=None): """Point source. :: 1 e^(-j w/c |x-x0|) G(x-x0, w) = --- ----------------- 4pi |x-x0| Examples -------- .. plot:: :context: close-figs p = sfs.mono.source.point(omega, x0, None, grid) sfs.plot.soundfield(p, grid) plt.title("Point Source at {} m".format(x0)) Normalization ... multiply by :math:`4\pi` ... .. plot:: :context: close-figs p *= 4 * np.pi sfs.plot.soundfield(p, grid) plt.title("Point Source at {} m (normalized)".format(x0)) """ k = util.wavenumber(omega, c) x0 = util.asarray_1d(x0) grid = util.XyzComponents(grid) r = np.linalg.norm(grid - x0) return 1 / (4*np.pi) * np.exp(-1j * k * r) / r def point_velocity(omega, x0, n0, grid, c=None): """Velocity of a point source. Returns ------- XyzComponents Particle velocity at positions given by `grid`. See :class:`sfs.util.XyzComponents`. """ k = util.wavenumber(omega, c) x0 = util.asarray_1d(x0) grid = util.XyzComponents(grid) offset = grid - x0 r = np.linalg.norm(offset) v = point(omega, x0, n0, grid, c=c) v *= (1+1j*k*r) / (defs.rho0 * defs.c * 1j*k*r) return util.XyzComponents([v * o / r for o in offset]) def point_modal(omega, x0, n0, grid, L, N=None, deltan=0, c=None): """Point source in a rectangular room using a modal room model. Parameters ---------- omega : float Frequency of source. x0 : (3,) array_like Position of source. n0 : (3,) array_like Normal vector (direction) of source (only required for compatibility). grid : triple of numpy.ndarray The grid that is used for the sound field calculations. See :func:`sfs.util.xyz_grid`. L : (3,) array_like Dimensionons of the rectangular room. N : (3,) array_like or int, optional Combination of modal orders in the three-spatial dimensions to calculate the sound field for or maximum order for all dimensions. If not given, the maximum modal order is approximately determined and the sound field is computed up to this maximum order. deltan : float, optional Absorption coefficient of the walls. c : float, optional Speed of sound. Returns ------- numpy.ndarray Sound pressure at positions given by `grid`. """ k = util.wavenumber(omega, c) x0 = util.asarray_1d(x0) x, y, z = util.XyzComponents(grid) if N is None: # determine maximum modal order per dimension Nx = int(np.ceil(L[0]/np.pi * k)) Ny = int(np.ceil(L[1]/np.pi * k)) Nz = int(np.ceil(L[2]/np.pi * k)) mm = range(Nx) nn = range(Ny) ll = range(Nz) elif np.isscalar(N): # compute up to a given order mm = range(N) nn = range(N) ll = range(N) else: # compute field for one order combination only mm = [N[0]] nn = [N[1]] ll = [N[2]] kmp0 = [((kx + 1j * deltan)**2, np.cos(kx * x) * np.cos(kx * x0[0])) for kx in [m * np.pi / L[0] for m in mm]] kmp1 = [((ky + 1j * deltan)**2, np.cos(ky * y) * np.cos(ky * x0[1])) for ky in [n * np.pi / L[1] for n in nn]] kmp2 = [((kz + 1j * deltan)**2, np.cos(kz * z) * np.cos(kz * x0[2])) for kz in [l * np.pi / L[2] for l in ll]] ksquared = k**2 p = 0 for (km0, p0), (km1, p1), (km2, p2) in itertools.product(kmp0, kmp1, kmp2): km = km0 + km1 + km2 p = p + 8 / (ksquared - km) * p0 * p1 * p2 return p def point_modal_velocity(omega, x0, n0, grid, L, N=None, deltan=0, c=None): """Velocity of point source in a rectangular room using a modal room model. Parameters ---------- omega : float Frequency of source. x0 : (3,) array_like Position of source. n0 : (3,) array_like Normal vector (direction) of source (only required for compatibility). grid : triple of numpy.ndarray The grid that is used for the sound field calculations. See :func:`sfs.util.xyz_grid`. L : (3,) array_like Dimensionons of the rectangular room. N : (3,) array_like or int, optional Combination of modal orders in the three-spatial dimensions to calculate the sound field for or maximum order for all dimensions. If not given, the maximum modal order is approximately determined and the sound field is computed up to this maximum order. deltan : float, optional Absorption coefficient of the walls. c : float, optional Speed of sound. Returns ------- XyzComponents Particle velocity at positions given by `grid`. See :class:`sfs.util.XyzComponents`. """ k = util.wavenumber(omega, c) x0 = util.asarray_1d(x0) x, y, z = util.XyzComponents(grid) if N is None: # determine maximum modal order per dimension Nx = int(np.ceil(L[0]/np.pi * k)) Ny = int(np.ceil(L[1]/np.pi * k)) Nz = int(np.ceil(L[2]/np.pi * k)) mm = range(Nx) nn = range(Ny) ll = range(Nz) elif np.isscalar(N): # compute up to a given order mm = range(N) nn = range(N) ll = range(N) else: # compute field for one order combination only mm = [N[0]] nn = [N[1]] ll = [N[2]] kmp0 = [((kx + 1j * deltan)**2, np.sin(kx * x) * np.cos(kx * x0[0])) for kx in [m * np.pi / L[0] for m in mm]] kmp1 = [((ky + 1j * deltan)**2, np.sin(ky * y) * np.cos(ky * x0[1])) for ky in [n * np.pi / L[1] for n in nn]] kmp2 = [((kz + 1j * deltan)**2, np.sin(kz * z) * np.cos(kz * x0[2])) for kz in [l * np.pi / L[2] for l in ll]] ksquared = k**2 vx = 0+0j vy = 0+0j vz = 0+0j for (km0, p0), (km1, p1), (km2, p2) in itertools.product(kmp0, kmp1, kmp2): km = km0 + km1 + km2 vx = vx - 8*1j / (ksquared - km) * p0 vy = vy - 8*1j / (ksquared - km) * p1 vz = vz - 8*1j / (ksquared - km) * p2 return util.XyzComponents([vx, vy, vz]) def line(omega, x0, n0, grid, c=None): """Line source parallel to the z-axis. Note: third component of x0 is ignored. :: (2) G(x-x0, w) = -j/4 H0 (w/c |x-x0|) Examples -------- .. plot:: :context: close-figs p = sfs.mono.source.line(omega, x0, None, grid) sfs.plot.soundfield(p, grid) plt.title("Line Source at {} m".format(x0[:2])) Normalization ... .. plot:: :context: close-figs p *= np.sqrt(8 * np.pi * omega / sfs.defs.c) * np.exp(1j * np.pi / 4) sfs.plot.soundfield(p, grid) plt.title("Line Source at {} m (normalized)".format(x0[:2])) """ k = util.wavenumber(omega, c) x0 = util.asarray_1d(x0) x0 = x0[:2] # ignore z-component grid = util.XyzComponents(grid) r = np.linalg.norm(grid[:2] - x0) p = -1j/4 * special.hankel2(0, k * r) return _duplicate_zdirection(p, grid) def line_velocity(omega, x0, n0, grid, c=None): """Velocity of line source parallel to the z-axis. Returns ------- XyzComponents Particle velocity at positions given by `grid`. See :class:`sfs.util.XyzComponents`. """ k = util.wavenumber(omega, c) x0 = util.asarray_1d(x0) x0 = x0[:2] # ignore z-component grid = util.XyzComponents(grid) offset = grid[:2] - x0 r = np.linalg.norm(offset) v = -1/(4*defs.c*defs.rho0) * special.hankel2(1, k * r) v = [v * o / r for o in offset] assert v[0].shape == v[1].shape if len(grid) > 2: v.append(np.zeros_like(v[0])) return util.XyzComponents([_duplicate_zdirection(vi, grid) for vi in v]) def line_dipole(omega, x0, n0, grid, c=None): """Line source with dipole characteristics parallel to the z-axis. Note: third component of x0 is ignored. :: (2) G(x-x0, w) = jk/4 H1 (w/c |x-x0|) cos(phi) """ k = util.wavenumber(omega, c) x0 = util.asarray_1d(x0) x0 = x0[:2] # ignore z-component n0 = n0[:2] grid = util.XyzComponents(grid) dx = grid[:2] - x0 r = np.linalg.norm(dx) p = 1j*k/4 * special.hankel2(1, k * r) * np.inner(dx, n0) / r return _duplicate_zdirection(p, grid) def plane(omega, x0, n0, grid, c=None): """Plane wave. :: G(x, w) = e^(-i w/c n x) Example ------- .. plot:: :context: close-figs direction = 45 # degree n0 = sfs.util.direction_vector(np.radians(direction)) p_plane = sfs.mono.source.plane(omega, x0, n0, grid) sfs.plot.soundfield(p_plane, grid); plt.title("Plane wave with direction {} degree".format(direction)) """ k = util.wavenumber(omega, c) x0 = util.asarray_1d(x0) n0 = util.asarray_1d(n0) grid = util.XyzComponents(grid) return np.exp(-1j * k * np.inner(grid - x0, n0)) def plane_velocity(omega, x0, n0, grid, c=None): """Velocity of a plane wave. :: V(x, w) = 1/(rho c) e^(-i w/c n x) n Returns ------- XyzComponents Particle velocity at positions given by `grid`. See :class:`sfs.util.XyzComponents`. """ v = plane(omega, x0, n0, grid, c=c) / (defs.rho0 * defs.c) return util.XyzComponents([v * n for n in n0]) def _duplicate_zdirection(p, grid): """If necessary, duplicate field in z-direction.""" gridshape = np.broadcast(*grid).shape if len(gridshape) > 2: return np.tile(p, [1, 1, gridshape[2]]) else: return p PKVGxx#sfs-0.2.0.dist-info/DESCRIPTION.rstSound Field Synthesis Toolbox for Python ======================================== Python implementation of the `Sound Field Synthesis Toolbox`_. .. _Sound Field Synthesis Toolbox: http://github.com/sfstoolbox/sfs/ Documentation: http://sfs.rtfd.org/ Code: http://github.com/sfstoolbox/sfs-python/ Python Package Index: http://pypi.python.org/pypi/sfs/ License: MIT -- see the file ``LICENSE`` for details. Requirements ------------ Obviously, you'll need Python_. We normally use Python 3.x, but it *should* also work with Python 2.x. NumPy_ and SciPy_ are needed for the calcuations. If you also want to plot the resulting sound fields, you'll need matplotlib_. Instead of installing all of them separately, you should probably get a Python distribution that already includes everything, e.g. Anaconda_. .. _Python: http://www.python.org/ .. _NumPy: http://www.numpy.org/ .. _SciPy: http://www.scipy.org/scipylib/ .. _matplotlib: http://matplotlib.org/ .. _Anaconda: http://docs.continuum.io/anaconda/ Installation ------------ Once you have installed the above-mentioned dependencies, you can use pip_ to download and install the latest release of the Sound Field Synthesis Toolbox with a single command:: pip install sfs --user If you want to install it system-wide for all users (assuming you have the necessary rights), you can just drop the ``--user`` option. To un-install, use:: pip uninstall sfs .. _pip: http://www.pip-installer.org/en/latest/installing.html How to Get Started ------------------ Various examples are located in the directory examples/ * sound_field_synthesis.py: Illustrates the general usage of the toolbox * horizontal_plane_arrays.py: Computes the sound fields for various techniques, virtual sources and loudspeaker array configurations * soundfigures.py: Illustrates the synthesis of sound figures with Wave Field Synthesis PKVGu!sfs-0.2.0.dist-info/metadata.json{"classifiers": ["Development Status :: 3 - Alpha", "License :: OSI Approved :: MIT License", "Operating System :: OS Independent", "Programming Language :: Python", "Topic :: Scientific/Engineering"], "extensions": {"python.details": {"contacts": [{"email": "sfstoolbox@gmail.com", "name": "SFS Toolbox Developers", "role": "author"}], "document_names": {"description": "DESCRIPTION.rst"}, "project_urls": {"Home": "http://github.com/sfstoolbox/"}}}, "extras": [], "generator": "bdist_wheel (0.26.0)", "keywords": ["audio", "SFS", "WFS", "Ambisonics"], "license": "MIT", "metadata_version": "2.0", "name": "sfs", "platform": "any", "run_requires": [{"requires": ["numpy", "scipy"]}], "summary": "Sound Field Synthesis Toolbox", "test_requires": [{"requires": ["pytest"]}], "version": "0.2.0"}PKVGפS!sfs-0.2.0.dist-info/top_level.txtsfs PKVGndnnsfs-0.2.0.dist-info/WHEELWheel-Version: 1.0 Generator: bdist_wheel (0.26.0) Root-Is-Purelib: true Tag: py2-none-any Tag: py3-none-any PKVG+&[ sfs-0.2.0.dist-info/METADATAMetadata-Version: 2.0 Name: sfs Version: 0.2.0 Summary: Sound Field Synthesis Toolbox Home-page: http://github.com/sfstoolbox/ Author: SFS Toolbox Developers Author-email: sfstoolbox@gmail.com License: MIT Keywords: audio,SFS,WFS,Ambisonics Platform: any Classifier: Development Status :: 3 - Alpha Classifier: License :: OSI Approved :: MIT License Classifier: Operating System :: OS Independent Classifier: Programming Language :: Python Classifier: Topic :: Scientific/Engineering Requires-Dist: numpy Requires-Dist: scipy Sound Field Synthesis Toolbox for Python ======================================== Python implementation of the `Sound Field Synthesis Toolbox`_. .. _Sound Field Synthesis Toolbox: http://github.com/sfstoolbox/sfs/ Documentation: http://sfs.rtfd.org/ Code: http://github.com/sfstoolbox/sfs-python/ Python Package Index: http://pypi.python.org/pypi/sfs/ License: MIT -- see the file ``LICENSE`` for details. Requirements ------------ Obviously, you'll need Python_. We normally use Python 3.x, but it *should* also work with Python 2.x. NumPy_ and SciPy_ are needed for the calcuations. If you also want to plot the resulting sound fields, you'll need matplotlib_. Instead of installing all of them separately, you should probably get a Python distribution that already includes everything, e.g. Anaconda_. .. _Python: http://www.python.org/ .. _NumPy: http://www.numpy.org/ .. _SciPy: http://www.scipy.org/scipylib/ .. _matplotlib: http://matplotlib.org/ .. _Anaconda: http://docs.continuum.io/anaconda/ Installation ------------ Once you have installed the above-mentioned dependencies, you can use pip_ to download and install the latest release of the Sound Field Synthesis Toolbox with a single command:: pip install sfs --user If you want to install it system-wide for all users (assuming you have the necessary rights), you can just drop the ``--user`` option. To un-install, use:: pip uninstall sfs .. _pip: http://www.pip-installer.org/en/latest/installing.html How to Get Started ------------------ Various examples are located in the directory examples/ * sound_field_synthesis.py: Illustrates the general usage of the toolbox * horizontal_plane_arrays.py: Computes the sound fields for various techniques, virtual sources and loudspeaker array configurations * soundfigures.py: Illustrates the synthesis of sound figures with Wave Field Synthesis PKVGWز  sfs-0.2.0.dist-info/RECORDsfs/__init__.py,sha256=Q6SD5iWBd7n62IsBJeM6-pD7T_WVTkros2LOu8MYLpQ,242 sfs/array.py,sha256=SpbtoQSpMLn59jx_BxuZ5T00bL6ElqmU2hI3mc-f_Cs,15003 sfs/defs.py,sha256=TPT6wZz4Ax2bKe5TR12py88jwnXPa1_lR1h8z5BCoRI,172 sfs/plot.py,sha256=l6uD0GuSpDmfaASXttm-AlqvUlBVac2b3WQcjYYA3mo,12572 sfs/tapering.py,sha256=rK4pghz8hOO_WgDvk-665v9nLHf6-7sLDzrdt0LBxIA,1708 sfs/util.py,sha256=WnE_0kCTV0b6zEN-ORQw1aqHUhQSLaXSvlZKak99Fzk,10512 sfs/mono/__init__.py,sha256=41iYheUFsZdwqZCrsd7E195o6418XugcNidySA3RVDA,103 sfs/mono/drivingfunction.py,sha256=6meIkE5TnlKPHGUVUtxiamNsN_RIKWMMlBdDyK600Vo,8717 sfs/mono/soundfigure.py,sha256=JgkmzFYOwAF99Jx6_dwN0WWpd1RN47g_0ko-VG1mrIY,1563 sfs/mono/source.py,sha256=_BkK-mHh-XA8dLpE0E5NslXFPCtt_GAAJ6heKQCD6n0,10299 sfs/mono/synthesized.py,sha256=B_GWVbN86R25STZuCXrvpTKTaFBqfiot9pV_CEfO4Ek,616 sfs-0.2.0.dist-info/DESCRIPTION.rst,sha256=EXRv8DfoxD6h2NPgwAtm8PKDC-YyxbruVc8VrIhhx2g,1912 sfs-0.2.0.dist-info/METADATA,sha256=6CdiTyr3TQX7LAr70Ss49IAS2LoSHIj9usuWg9OvYIg,2439 sfs-0.2.0.dist-info/RECORD,, sfs-0.2.0.dist-info/WHEEL,sha256=GrqQvamwgBV4nLoJe0vhYRSWzWsx7xjlt74FT0SWYfE,110 sfs-0.2.0.dist-info/metadata.json,sha256=ZOe7YF8VAKj73eXn0HvpjwQhvCo-uzBK3dCShrZjxfk,793 sfs-0.2.0.dist-info/top_level.txt,sha256=6mBnHRnczkH-3MoFY3vtAo_cer661XzjlrnrFOWAW28,4 PKVG2܁sfs/__init__.pyPK+j~G4%:: sfs/array.pyPKG)) ;sfs/util.pyPKGd711 esfs/plot.pyPK+j~GAAi bsfs/defs.pyPKhEVr7sfs/tapering.pyPKz[^GVWMggsfs/mono/__init__.pyPK+j~G=ݺ " "sfs/mono/drivingfunction.pyPK+j~Gƅ[sfs/mono/soundfigure.pyPK+j~Gxwahh?sfs/mono/synthesized.pyPK+j~Gb;(;(sfs/mono/source.pyPKVGxx#Gsfs-0.2.0.dist-info/DESCRIPTION.rstPKVGu!sfs-0.2.0.dist-info/metadata.jsonPKVGפS!Xsfs-0.2.0.dist-info/top_level.txtPKVGndnnsfs-0.2.0.dist-info/WHEELPKVG+&[ @sfs-0.2.0.dist-info/METADATAPKVGWز  sfs-0.2.0.dist-info/RECORDPK|F