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

"""
from collections import namedtuple as _namedtuple

import numpy as _np

from . import util as _util

[docs]class SecondarySourceDistribution(_namedtuple('SecondarySourceDistribution',
'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 'SecondarySourceDistribution(\n' + ',\n'.join(
'    {}={}'.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 SecondarySourceDistribution(
self.x[indices], self.n[indices], self.a[indices])

[docs]def as_secondary_source_distribution(arg, **kwargs):
r"""Create a `SecondarySourceDistribution`.

Parameters
----------
arg : sequence of between 1 and 3 array_like objects
All elements are converted to NumPy arrays.
If only 1 element is given, all normal vectors are set to *NaN*.
If only 1 or 2 elements are given, all weights are set to ``1.0``.
**kwargs
All keyword arguments are forwarded to :func:`numpy.asarray`.

Returns
-------
`SecondarySourceDistribution`
A named tuple consisting of three `numpy.ndarray`\s containing
positions, normal vectors and weights.

"""
if len(arg) == 3:
x, n, a = arg
elif len(arg) == 2:
x, n = arg
a = 1.0
elif len(arg) == 1:
x, = arg
n = _np.nan, _np.nan, _np.nan
a = 1.0
else:
raise TypeError('Between 1 and 3 elements are required')
x = _util.asarray_of_rows(x, **kwargs)
n = _util.asarray_of_rows(n, **kwargs)
if len(n) == 1:
n = _np.tile(n, (len(x), 1))
a = _util.asarray_1d(a, **kwargs)
if len(a) == 1:
a = _np.tile(a, len(x))
return SecondarySourceDistribution(x, n, a)

[docs]def linear(N, spacing, *, center=[0, 0, 0], orientation=[1, 0, 0]):
"""Return linear, equidistantly sampled 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
-------
`SecondarySourceDistribution`
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.plot2d.loudspeakers(x0, n0, a0)
plt.axis('equal')
plt.xlabel('x / m')
plt.ylabel('y / m')

"""
return _linear_helper(_np.arange(N) * spacing, center, orientation)

[docs]def linear_diff(distances, *, center=[0, 0, 0], orientation=[1, 0, 0]):
"""Return 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
-------
`SecondarySourceDistribution`
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.plot2d.loudspeakers(x0, n0, a0)
plt.axis('equal')
plt.xlabel('x / m')
plt.ylabel('y / m')

"""
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):
"""Return 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
-------
`SecondarySourceDistribution`
Positions, orientations and weights of secondary sources.

Examples
--------
.. plot::
:context: close-figs

x0, n0, a0 = sfs.array.linear_random(
N=12,
min_spacing=0.15, max_spacing=0.4,
orientation=[0, -1, 0])
sfs.plot2d.loudspeakers(x0, n0, a0)
plt.axis('equal')
plt.xlabel('x / m')
plt.ylabel('y / m')

"""
r = _np.random.RandomState(seed)
distances = r.uniform(min_spacing, max_spacing, size=N-1)
return linear_diff(distances, center=center, orientation=orientation)

[docs]def circular(N, R, *, center=[0, 0, 0]):
"""Return circular secondary source distribution parallel to the xy-plane.

Parameters
----------
N : int
Number of secondary sources.
R : float
center
See `linear()`.

Returns
-------
`SecondarySourceDistribution`
Positions, orientations and weights of secondary sources.

Examples
--------
.. plot::
:context: close-figs

x0, n0, a0 = sfs.array.circular(16, 1)
sfs.plot2d.loudspeakers(x0, n0, a0, size=0.2, show_numbers=True)
plt.axis('equal')
plt.xlabel('x / m')
plt.ylabel('y / m')

"""
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 SecondarySourceDistribution(positions, normals, weights)

[docs]def rectangular(N, spacing, *, center=[0, 0, 0], orientation=[1, 0, 0]):
"""Return 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
-------
`SecondarySourceDistribution`
Positions, orientations and weights of secondary sources.

Examples
--------
.. plot::
:context: close-figs

x0, n0, a0 = sfs.array.rectangular((4, 8), 0.2)
sfs.plot2d.loudspeakers(x0, n0, a0, show_numbers=True)
plt.axis('equal')
plt.xlabel('x / m')
plt.ylabel('y / m')

"""
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(
# left
linear(N1, spacing, center=[-offset1, 0, 0], orientation=[1, 0, 0]),
# upper
linear(N2, spacing, center=[0, offset2, 0], orientation=[0, -1, 0]),
# right
linear(N1, spacing, center=[offset1, 0, 0], orientation=[-1, 0, 0]),
# lower
linear(N2, spacing, center=[0, -offset2, 0], orientation=[0, 1, 0]),
)
positions, normals = _rotate_array(positions, normals,
[1, 0, 0], orientation)
positions += center
return SecondarySourceDistribution(positions, normals, weights)

[docs]def rounded_edge(Nxy, Nr, spacing, *, center=[0, 0, 0], orientation=[1, 0, 0]):
"""Return SSD 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.
spacing : float
Distance (in metres) between secondary sources.
center : (3,) array_like, optional
Position of edge.
orientation : (3,) array_like, optional
Normal vector of array.  Default orientation is along xy-axis.

Returns
-------
`SecondarySourceDistribution`
Positions, orientations and weights of secondary sources.

Examples
--------
.. plot::
:context: close-figs

x0, n0, a0 = sfs.array.rounded_edge(8, 5, 0.2)
sfs.plot2d.loudspeakers(x0, n0, a0)
plt.axis('equal')
plt.xlabel('x / m')
plt.ylabel('y / m')

"""
Nr += 1
R = 2/_np.pi * Nr * spacing

# array along y-axis
x00, n00, a00 = linear(Nxy, spacing,
center=[0, Nxy//2*spacing+spacing/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] = spacing
positions = _np.concatenate((positions, x00))
directions = _np.concatenate((directions, n00))
weights = _np.concatenate((weights, a00))

# array along x-axis
x00, n00, a00 = linear(Nxy, spacing,
center=[Nxy//2*spacing-spacing/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 SecondarySourceDistribution(positions, directions, weights)

[docs]def edge(Nxy, spacing, *, center=[0, 0, 0], orientation=[1, 0, 0]):
"""Return SSD along the xy-axis with sharp edge at the origin.

Parameters
----------
Nxy : int
Number of secondary sources along x- and y-axis.
spacing : float
Distance (in metres) between secondary sources.
center : (3,) array_like, optional
Position of edge.
orientation : (3,) array_like, optional
Normal vector of array.  Default orientation is along xy-axis.

Returns
-------
`SecondarySourceDistribution`
Positions, orientations and weights of secondary sources.

Examples
--------
.. plot::
:context: close-figs

x0, n0, a0 = sfs.array.edge(8, 0.2)
sfs.plot2d.loudspeakers(x0, n0, a0)
plt.axis('equal')
plt.xlabel('x / m')
plt.ylabel('y / m')

"""
# array along y-axis
x00, n00, a00 = linear(Nxy, spacing,
center=[0, Nxy//2*spacing+spacing/2, 0])
x00 = _np.flipud(x00)
positions = x00
directions = n00
weights = a00

# array along x-axis
x00, n00, a00 = linear(Nxy, spacing,
center=[Nxy//2*spacing-spacing/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 SecondarySourceDistribution(positions, directions, weights)

[docs]def planar(N, spacing, *, center=[0, 0, 0], orientation=[1, 0, 0]):
"""Return 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
-------
`SecondarySourceDistribution`
Positions, orientations and weights of secondary sources.

Examples
--------
.. plot::
:context: close-figs

x0, n0, a0 = sfs.array.planar(
(4,3), 0.5, orientation=[0, 0, 1])  # 4 sources along y, 3 sources along x
x0, n0, a0 = sfs.array.planar(
(4,3), 0.5, orientation=[1, 0, 0])  # 4 sources along y, 3 sources along z

x0, n0, a0 = sfs.array.planar(
(4,3), 0.5, orientation=[0, 1, 0])  # 4 sources along x, 3 sources along z
sfs.plot2d.loudspeakers(x0, n0, a0)  # plot the last ssd in 2D
plt.axis('equal')
plt.xlabel('x / m')
plt.ylabel('y / m')

"""
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 SecondarySourceDistribution(positions, normals, weights)

[docs]def cube(N, spacing, *, center=[0, 0, 0], orientation=[1, 0, 0]):
"""Return 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
-------
`SecondarySourceDistribution`
Positions, orientations and weights of secondary sources.

Examples
--------
.. plot::
:context: close-figs

x0, n0, a0 = sfs.array.cube(
N=2, spacing=0.5,
center=[0, 0, 0], orientation=[1, 0, 0])
sfs.plot2d.loudspeakers(x0, n0, a0)
plt.axis('equal')
plt.xlabel('x / m')
plt.ylabel('y / m')
plt.title('view onto xy-plane')

"""
N1, N2, N3 = (N, N, N) if _np.isscalar(N) else N
d = spacing
offset1 = d * (N2 - 1) / 2 + d / _np.sqrt(2)
offset2 = d * (N1 - 1) / 2 + d / _np.sqrt(2)
offset3 = d * (N3 - 1) / 2 + d / _np.sqrt(2)
positions, directions, weights = concatenate(
# west
planar((N1, N3), d, center=[-offset1, 0, 0], orientation=[1, 0, 0]),
# north
planar((N2, N3), d, center=[0, offset2, 0], orientation=[0, -1, 0]),
# east
planar((N1, N3), d, center=[offset1, 0, 0], orientation=[-1, 0, 0]),
# south
planar((N2, N3), d, center=[0, -offset2, 0], orientation=[0, 1, 0]),
# bottom
planar((N2, N1), d, center=[0, 0, -offset3], orientation=[0, 0, 1]),
# top
planar((N2, N1), d, center=[0, 0, offset3], orientation=[0, 0, -1]),
)
positions, directions = _rotate_array(positions, directions,
[1, 0, 0], orientation)
positions += center
return SecondarySourceDistribution(positions, directions, weights)

"""Load spherical secondary source distribution from file.

ASCII Format (see MATLAB SFS Toolbox) with 4 numbers (3 for the cartesian
position vector, 1 for the integration weight) per secondary source located
on the unit circle which is resized by the given radius and shifted to the
given center.

Returns
-------
`SecondarySourceDistribution`
Positions, orientations and weights of secondary sources.

Examples
--------
content of ``example_array_6LS_3D.txt``::

1 0 0 1
-1 0 0 1
0 1 0 1
0 -1 0 1
0 0 1 1
0 0 -1 1

corresponds to the `3-dimensional 6-point spherical 3-design
<http://neilsloane.com/sphdesigns/dim3/des.3.6.3.txt>`_.

.. plot::
:context: close-figs

'../data/arrays/example_array_6LS_3D.txt',
center=[0, 0, 0])
sfs.plot2d.loudspeakers(x0, n0, a0, size=0.25)
plt.axis('equal')
plt.xlabel('x / m')
plt.ylabel('y / m')
plt.title('view onto xy-plane')

"""
positions, weights = data[:, :3], data[:, 3]
normals = -positions
positions += center
return SecondarySourceDistribution(positions, normals, weights)

[docs]def load(file, *, center=[0, 0, 0], orientation=[1, 0, 0]):
"""Load secondary source distribution from file.

Comma Separated Values (CSV) format with 7 values
(3 for the cartesian position vector, 3 for the cartesian inward normal
vector, 1 for the integration weight) per secondary source.

Returns
-------
`SecondarySourceDistribution`
Positions, orientations and weights of secondary sources.

Examples
--------
content of ``example_array_4LS_2D.csv``::

1,0,0,-1,0,0,1
0,1,0,0,-1,0,1
-1,0,0,1,0,0,1
0,-1,0,0,1,0,1

corresponds to 4 sources at 1, j, -1, -j in the complex plane. This setup
is typically used for Quadraphonic audio reproduction.

.. plot::
:context: close-figs

sfs.plot2d.loudspeakers(x0, n0, a0)
plt.axis('equal')
plt.xlabel('x / m')
plt.ylabel('y / m')

.. plot::
:context: close-figs

'../data/arrays/wfs_university_rostock_2018.csv')
sfs.plot2d.loudspeakers(x0, n0, a0)
plt.axis('equal')
plt.xlabel('x / m')
plt.ylabel('y / m')
plt.title('top view of 64 channel WFS system at university of Rostock')

"""
positions, normals, weights = data[:, :3], data[:, 3:6], data[:, 6]
positions, normals = _rotate_array(positions, normals,
[1, 0, 0], orientation)
positions += center
return SecondarySourceDistribution(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 along the
contour.

closed : bool
``True`` if the loudspeaker contour is closed.

Returns
-------
(N,) numpy.ndarray
Weights of secondary sources.

Examples
--------
>>> import sfs
>>> x0, n0, a0 = sfs.array.circular(2**5, 1)
>>> a = sfs.array.weights_midpoint(x0, closed=True)
>>> max(abs(a0-a))
0.0003152601902411123

"""
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 SecondarySourceDistribution(positions, normals, weights)

[docs]def concatenate(*arrays):
"""Concatenate `SecondarySourceDistribution` objects.

Returns
-------
`SecondarySourceDistribution`
Positions, orientations and weights
of the concatenated secondary sources.

Examples
--------
.. plot::
:context: close-figs

ssd1 = sfs.array.edge(10, 0.2)
ssd2 = sfs.array.edge(20, 0.1, center=[2, 2, 0], orientation=[-1, 0, 0])
x0, n0, a0 = sfs.array.concatenate(ssd1, ssd2)
sfs.plot2d.loudspeakers(x0, n0, a0)
plt.axis('equal')
plt.xlabel('x / m')
plt.ylabel('y / m')

"""
return SecondarySourceDistribution._make(_np.concatenate(i)
for i in zip(*arrays))
```