# Source code for sfs.array

"""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

[docs]class ArrayData(namedtuple('ArrayData', 'x n a')): """Named tuple returned by array functions. See 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)) + ')'
[docs] def take(self, indices): """Return a sub-array given by *indices*.""" return ArrayData(self.x[indices], self.n[indices], self.a[indices])
[docs]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. Examples -------- .. 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)
[docs]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 linear(). Returns ------- ArrayData Positions, orientations and weights of secondary sources. Examples -------- .. 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)
[docs]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 linear(). seed : {None, int, array_like}, optional Random seed. See numpy.random.RandomState. Returns ------- ArrayData Positions, orientations and weights of secondary sources. Examples -------- .. 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)
[docs]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 linear(). Returns ------- ArrayData Positions, orientations and weights of secondary sources. Examples -------- .. 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)
[docs]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 linear(). The *orientation* corresponds to the first linear segment. Returns ------- ArrayData Positions, orientations and weights of secondary sources. Examples -------- .. 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)
[docs]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. Examples -------- .. 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)
[docs]def edge(Nxy, dx, center=[0, 0, 0], orientation=[1, 0, 0]): """Array along the xy-axis with edge at the origin. Parameters ---------- Nxy : int Number of secondary sources along x- and y-axis. 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. Examples -------- .. plot:: :context: close-figs x0, n0, a0 = sfs.array.edge(8, 0.2) sfs.plot.loudspeaker_2d(x0, n0, a0) plt.axis('equal') """ # array along y-axis x00, n00, a00 = linear(Nxy, dx, center=[0, Nxy//2*dx+dx/2, 0]) x00 = np.flipud(x00) positions = x00 directions = n00 weights = a00 # array along x-axis x00, n00, a00 = linear(Nxy, dx, center=[Nxy//2*dx-dx/2, 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)
[docs]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 linear(). Returns ------- ArrayData Positions, orientations and weights of secondary sources. """ 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)
[docs]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 rectangular(), the last one specifies the vertical edge. spacing : float Distance (in metres) between secondary sources. center, orientation See linear(). The *orientation* corresponds to the first planar segment. Returns ------- ArrayData Positions, orientations and weights of secondary sources. """ 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)
[docs]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. """ data = np.loadtxt(fname) positions, weights = data[:, :3], data[:, 3] normals = -positions positions *= radius positions += center return ArrayData(positions, normals, weights)
[docs]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. """ 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)
[docs]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.""" center = util.asarray_1d(center) 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)
[docs]def concatenate(*arrays): """Concatenate ArrayData objects.""" return ArrayData._make(np.concatenate(i) for i in zip(*arrays))