Merge pull request #1 from spatialaudio/spherical_ps

Spherical ps

Author: spors

doc/index.rst

@@ -4,6 +4,10 @@
 
 .. toctree::
 
+    usage
+    modal
+    utilities
+
 .. only:: html
 
     * :ref:`genindex`

doc/modal.rst

@@ -0,0 +1,14 @@
+Modal Beamforming
+=================
+
+.. automodule:: micarray.modal
+
+Angular
+-------
+
+.. automodule:: micarray.modal.angular
+
+Radial
+------
+
+.. automodule:: micarray.modal.radial

doc/readthedocs-environment.yml

@@ -0,0 +1,9 @@
+channels:
+  - conda-forge
+dependencies:
+  - python==3.5
+  - sphinx>=1.3
+  - sphinx_rtd_theme
+  - numpy
+  - scipy
+  - matplotlib>=1.5

doc/usage.rst

@@ -0,0 +1,24 @@
+Usage
+=====
+
+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 calculations.
+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/
+
+How to Get Started
+------------------
+
+Various examples are located in the directory 

doc/utilities.rst

@@ -0,0 +1,4 @@
+Utilities
+=========
+
+.. automodule:: micarray.util

micarray/modal/radial.py

@@ -1,11 +1,111 @@
 from __future__ import division
 import numpy as np
 from scipy import special
+from .. import util
 
 
-def spherical(N, kr, setup, plane_wave):
-    if np.isscalar(kr):
-        kr = np.asarray([kr])
+def spherical_pw(N, k, r, setup):
+    r"""Radial coefficients for a plane wave
+
+    Computes the radial component of the spherical harmonics expansion of a
+    plane wave impinging on a spherical array.
+
+    .. math::
+
+        \mathring{P}_n(k) = 4 \pi i^n b_n(kr)
+
+    Parameters
+    ----------
+    N : int
+        Maximum order.
+    k : array_like
+        Wavenumber.
+    r : float
+        Radius of microphone array.
+    setup : {'open', 'card', 'rigid'}
+        Array configuration (open, cardioids, rigid).
+
+    Returns
+    -------
+    numpy.ndarray
+        Radial weights for all orders up to N and the given wavenumbers.
+    """
+    kr = util.asarray_1d(k*r)
+    n = np.arange(N+1)
+
+    bn = weights(N, kr, setup)
+    for i, x in enumerate(kr):
+        bn[i, :] = bn[i, :] * 4*np.pi * (1j)**n
+    return bn
+
+
+def spherical_ps(N, k, r, rs, setup):
+    r"""Radial coefficients for a point source
+
+    Computes the radial component of the spherical harmonics expansion of a
+    point source impinging on a spherical array.
+
+    .. math::
+
+        \mathring{P}_n(k) = 4 \pi (-i) k h_n^{(2)}(k r_s) b_n(kr)
+
+    Parameters
+    ----------
+    N : int
+        Maximum order.
+    k : array_like
+        Wavenumber.
+    r : float
+        Radius of microphone array.
+    rs : float
+        Distance of source.
+    setup : {'open', 'card', 'rigid'}
+        Array configuration (open, cardioids, rigid).
+
+    Returns
+    -------
+    numpy.ndarray
+        Radial weights for all orders up to N and the given wavenumbers.
+    """
+    k = util.asarray_1d(k)
+    krs = k*rs
+    n = np.arange(N+1)
+
+    bn = weights(N, k*r, setup)
+    for i, x in enumerate(krs):
+        hn = special.spherical_jn(n, x) - 1j * special.spherical_yn(n, x)
+        bn[i, :] = bn[i, :] * 4*np.pi * (-1j) * hn * k[i]
+
+    return bn
+
+
+def weights(N, kr, setup):
+    r"""Radial weighing functions
+
+    Computes the radial weighting functions for diferent array types
+    (cf. eq.(2.62), Rafaely 2015).
+
+    For instance for an rigid array
+
+    .. math::
+
+        b_n(kr) = j_n(kr) - \frac{j_n^\prime(kr)}{h_n^{(2)\prime}(kr)}h_n^{(2)}(kr)
+
+    Parameters
+    ----------
+    N : int
+        Maximum order.
+    kr : array_like
+        Wavenumber * radius.
+    setup : {'open', 'card', 'rigid'}
+        Array configuration (open, cardioids, rigid).
+
+    Returns
+    -------
+    numpy.ndarray
+        Radial weights for all orders up to N and the given wavenumbers.
+
+    """
     n = np.arange(N+1)
     bns = np.zeros((len(kr), N+1), dtype=complex)
     for i, x in enumerate(kr):
@@ -21,11 +121,8 @@ def spherical(N, kr, setup, plane_wave):
             bn = jn - jnd/hnd*hn
         else:
             raise ValueError('setup must be either: open, card or rigid')
-        if plane_wave:
-            bn = bn * 4*np.pi * (1j)**n
         bns[i, :] = bn
-        bns = np.squeeze(bns)
-    return bns
+    return np.squeeze(bns)
 
 
 def regularize(dn, a0, method):

readthedocs.yml

@@ -0,0 +1,4 @@
+conda:
+    file: doc/readthedocs-environment.yml
+python:
+    setup_py_install: true