Merge pull request #3 from spatialaudio/sfa_circular_harmonics

SFA circular harmonics

Author: trettberg

examples/modal_beamforming_open_circular_array.py

@@ -0,0 +1,48 @@
+"""
+    Compute the plane wave decomposition for an incident broadband plane wave
+    on an open circular array using a modal beamformer of finite order.
+"""
+
+import numpy as np
+import matplotlib.pyplot as plt
+import micarray
+
+N = 90  # order of modal beamformer/microphone array
+pw_angle = 1.23 * np.pi  # incidence angle of plane wave
+pol_pwd = np.linspace(0, 2*np.pi, 91, endpoint=False)  # angles for plane wave decomposition
+k = np.linspace(0.1, 20, 100)  # wavenumber vector
+r = 1  # radius of array
+
+# get uniform grid (microphone positions) of order N
+pol, weights = micarray.modal.angular.grid_equal_polar_angle(N)
+# get circular harmonics matrix for sensors
+Psi_p = micarray.modal.angular.cht_matrix(N, pol, weights)
+# get circular harmonics matrix for a source ensemble of azimuthal plane wave
+Psi_q = micarray.modal.angular.cht_matrix(N, pol_pwd)
+# get radial filters
+Bn = micarray.modal.radial.circular_pw(N, k, r, setup='open')
+Dn, _ = micarray.modal.radial.regularize(1/Bn, 100, 'softclip')
+D = micarray.modal.radial.circ_diagonal_mode_mat(Dn)
+
+# compute microphone signals for an incident broad-band plane wave
+p = np.exp(1j * k[:, np.newaxis]*r * np.cos(pol - pw_angle))
+# compute plane wave decomposition
+A_pwd = np.matmul(np.matmul(np.conj(Psi_q.T), D), Psi_p)
+q_pwd = np.squeeze(np.matmul(A_pwd, np.expand_dims(p, 2)))
+q_pwd_t = np.fft.fftshift(np.fft.irfft(q_pwd, axis=0), axes=0)
+
+# visualize plane wave decomposition (aka beampattern)
+plt.figure()
+plt.pcolormesh(k, pol_pwd/np.pi, micarray.util.db(q_pwd.T), vmin=-40)
+plt.colorbar()
+plt.xlabel(r'$kr$')
+plt.ylabel(r'$\phi / \pi$')
+plt.title('Plane wave docomposition by modal beamformer (frequency domain)')
+plt.savefig('modal_open_beamformer_pwd_fd.png')
+
+plt.figure()
+plt.pcolormesh(range(2*len(k)-2), pol_pwd/np.pi, micarray.util.db(q_pwd_t.T), vmin=-40)
+plt.colorbar()
+plt.ylabel(r'$\phi / \pi$')
+plt.title('Plane wave docomposition by modal beamformer (time domain)')
+plt.savefig('modal_open_beamformer_pwd_td.png')

examples/modal_beamforming_rigid_circular_array.py

@@ -0,0 +1,52 @@
+"""
+    Compute the plane wave decomposition for an incident broadband plane wave
+    on an rigid circular array using a modal beamformer of finite order.
+"""
+
+import numpy as np
+import matplotlib.pyplot as plt
+import micarray
+
+Nsf = 50  # order of the incident sound field
+N = 30  # order of modal beamformer/microphone array
+pw_angle = 1 * np.pi  # incidence angle of plane wave
+pol_pwd = np.linspace(0, 2*np.pi, 180, endpoint=False)  # angles for plane wave decomposition
+k = np.linspace(0.1, 20, 100)  # wavenumber vector
+r = 1  # radius of array
+
+# get uniform grid (microphone positions) of order N
+pol, weights = micarray.modal.angular.grid_equal_polar_angle(N)
+
+# pressure on the surface of a rigid cylinder for an incident plane wave
+bn = micarray.modal.radial.circular_pw(Nsf, k, r, setup='rigid')
+D = micarray.modal.radial.circ_diagonal_mode_mat(bn)
+Psi_p = micarray.modal.angular.cht_matrix(Nsf, pol, weights)
+Psi_pw = micarray.modal.angular.cht_matrix(Nsf, pw_angle)
+p = np.matmul(np.matmul(np.conj(Psi_pw.T), D), Psi_p)
+p = np.squeeze(p)
+
+# plane wave decomposition using modal beamforming
+Psi_p = micarray.modal.angular.cht_matrix(N, pol)
+Psi_q = micarray.modal.angular.cht_matrix(N, pol_pwd)
+Bn = micarray.modal.radial.circular_pw(N, k, r, setup='rigid')
+Dn, _ = micarray.modal.radial.regularize(1/Bn, 100, 'softclip')
+D = micarray.modal.radial.circ_diagonal_mode_mat(Dn)
+A_pwd = np.matmul(np.matmul(np.conj(Psi_q.T), D), Psi_p)
+q_pwd = np.squeeze(np.matmul(A_pwd, np.expand_dims(p, 2)))
+q_pwd_t = np.fft.fftshift(np.fft.irfft(q_pwd, axis=0), axes=0)
+
+# visualize plane wave decomposition (aka beampattern)
+plt.figure()
+plt.pcolormesh(k, pol_pwd/np.pi, micarray.util.db(q_pwd.T), vmin=-40)
+plt.colorbar()
+plt.xlabel(r'$kr$')
+plt.ylabel(r'$\phi / \pi$')
+plt.title('Plane wave docomposition by modal beamformer (frequency domain)')
+plt.savefig('modal_open_beamformer_pwd_fd.png')
+
+plt.figure()
+plt.pcolormesh(range(2*len(k)-2), pol_pwd/np.pi, micarray.util.db(q_pwd_t.T), vmin=-40)
+plt.colorbar()
+plt.ylabel(r'$\phi / \pi$')
+plt.title('Plane wave docomposition by modal beamformer (time domain)')
+plt.savefig('modal_open_beamformer_pwd_td.png')

micarray/modal/angular.py

@@ -91,6 +91,48 @@ def Legendre_matrix(N, ctheta):
     return Lmn
 
 
+def cht_matrix(N, pol, weights=None):
+    r"""Matrix of circular harmonics up to order N for given angles.
+
+    Computes a matrix of circular harmonics up to order :math:`N`
+    for the given angles/grid.
+
+    .. math::
+        \Psi = \left[ \begin{array}{ccccccc}
+        1 & e^{i\varphi[0]} & \cdots & e^{iN\varphi[0]} & e^{-iN\varphi[0]} & \cdots & e^{-i\varphi[0]} \\
+        1 & e^{i\varphi[1]} & \cdots & e^{iN\varphi[1]} & e^{-iN\varphi[1]} & \cdots & e^{-i\varphi[1]} \\
+        \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\
+        1 & e^{i\varphi[Q-1]} & \cdots & e^{iN\varphi[Q-1]} & e^{-iN\varphi[Q-1]} & \cdots & e^{-i\varphi[Q-1]}
+        \end{array} \right]
+
+    Parameters
+    ----------
+    N : int
+        Maximum order.
+    pol : (Q,) array_like
+        Polar angle.
+    weights : (Q,) array_like, optional
+        Weights.
+
+    Returns
+    -------
+    Psi : (2N+1, Q) numpy.ndarray
+        Matrix of circular harmonics.
+    """
+    pol = util.asarray_1d(pol)
+    if pol.ndim == 0:
+        Q = 1
+    else:
+        Q = len(pol)
+    if weights is None:
+        weights = np.ones(Q)
+    Psi = np.zeros([(2*N+1), Q], dtype=complex)
+    order = np.roll(np.arange(-N, N+1), -N)
+    for i, n in enumerate(order):
+        Psi[i, :] = weights * np.exp(1j * n * pol)
+    return Psi
+
+
 def grid_equal_angle(n):
     """Equi-angular sampling points on a sphere.
 
@@ -153,3 +195,24 @@ def grid_gauss(n):
     weights = np.repeat(weights, 2*n+2)
     weights *= np.pi / (n+1)      # sum(weights) == 4pi
     return azi, elev, weights
+
+
+def grid_equal_polar_angle(n):
+    """Equi-angular sampling points on a circle.
+
+    Parameters
+    ----------
+    n : int
+        Maximum order
+
+    Returns
+    -------
+    pol : array_like
+        Polar angle.
+    weights : array_like
+        Weights.
+    """
+    num_mic = 2*n+1
+    pol = np.linspace(0, 2*np.pi, num=num_mic, endpoint=False)
+    weights = 1/num_mic * np.ones(num_mic)
+    return pol, weights

micarray/modal/radial.py

@@ -227,3 +227,171 @@ def repeat_n_m(v):
     """
     krlist = [np.tile(v, (2*i+1, 1)).T for i, v in enumerate(v.T.tolist())]
     return np.squeeze(np.concatenate(krlist, axis=-1))
+
+
+def circular_pw(N, k, r, setup):
+    r"""Radial coefficients for a plane wave.
+
+    Computes the radial component of the circular harmonics expansion of a
+    plane wave impinging on a circular array.
+
+    .. math::
+
+        \mathring{P}_n(k) = i^{-n} b_n(kr)
+
+    Parameters
+    ----------
+    N : int
+        Maximum order.
+    k : (M,) array_like
+        Wavenumber.
+    r : float
+        Radius of microphone array.
+    setup : {'open', 'card', 'rigid'}
+        Array configuration (open, cardioids, rigid).
+
+    Returns
+    -------
+    bn : (M, N+1) 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 = circ_radial_weights(N, kr, setup)
+    for i, x in enumerate(kr):
+        bn[i, :] = bn[i, :] * (1j)**(n)
+    return bn
+
+
+def circular_ls(N, k, r, rs, setup):
+    r"""Radial coefficients for a line source.
+
+    Computes the radial component of the circular harmonics expansion of a
+    line source impinging on a circular array.
+
+    .. math::
+
+        \mathring{P}_n(k) = \frac{-i}{4} H_n^{(2)}(k r_s) b_n(kr)
+
+    Parameters
+    ----------
+    N : int
+        Maximum order.
+    k : (M,) array_like
+        Wavenumber.
+    r : float
+        Radius of microphone array.
+    rs : float
+        Distance of source.
+    setup : {'open', 'card', 'rigid'}
+        Array configuration (open, cardioids, rigid).
+
+    Returns
+    -------
+    bn : (M, N+1) 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 = circ_radial_weights(N, k*r, setup)
+    for i, x in enumerate(krs):
+        Hn = special.hankel2(n, x)
+        bn[i, :] = bn[i, :] * -1j/4 * Hn
+    return bn
+
+
+def circ_radial_weights(N, kr, setup):
+    r"""Radial weighing functions.
+
+    Computes the radial weighting functions for diferent array types
+
+    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 : (M,) array_like
+        Wavenumber * radius.
+    setup : {'open', 'card', 'rigid'}
+        Array configuration (open, cardioids, rigid).
+
+    Returns
+    -------
+    bn : (M, N+1) 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):
+        Jn = special.jv(n, x)
+        if setup == 'open':
+            bn = Jn
+        elif setup == 'card':
+            bn = Jn - 1j * special.jvp(n, x, n=1)
+        elif setup == 'rigid':
+            Jnd = special.jvp(n, x, n=1)
+            Hn = special.hankel2(n, x)
+            Hnd = special.h2vp(n, x)
+            bn = Jn - Jnd/Hnd*Hn
+        else:
+            raise ValueError('setup must be either: open, card or rigid')
+        Bns[i, :] = bn
+    return np.squeeze(Bns)
+
+
+def circ_diagonal_mode_mat(bk):
+    """Diagonal matrix of radial coefficients for all modes/wavenumbers.
+
+    Parameters
+    ----------
+    bk : (M, N+1) numpy.ndarray
+        Vector containing values for all wavenumbers :math:`M` and modes up to
+        order :math:`N`
+
+    Returns
+    -------
+    Bk : (M, 2*N+1, 2*N+1) numpy.ndarray
+        Multidimensional array containing diagnonal matrices with input
+        vector on main diagonal.
+    """
+    bk = mirror_vec(bk)
+    if len(bk.shape) == 1:
+        bk = bk[np.newaxis, :]
+    K, N = bk.shape
+    Bk = np.zeros([K, N, N], dtype=complex)
+    for k in range(K):
+        Bk[k, :, :] = np.diag(bk[k, :])
+    return np.squeeze(Bk)
+
+
+def mirror_vec(v):
+    """Mirror elements in a vector.
+
+    Returns a vector of length *2*len(v)-1* with symmetric elements.
+    The first *len(v)* elements are the same as *v* and the last *len(v)-1*
+    elements are *v[:0:-1]*. The function can be used to order the circular
+    harmonic coefficients. If *v* is a matrix, it is treated as a stack of
+    vectors residing in the last index and broadcast accordingly.
+
+    Parameters
+    ----------
+    v : (, N+1) numpy.ndarray
+        Input vector of stack of input vectors
+
+    Returns
+    -------
+     : (, 2*N+1) numpy.ndarray
+        Vector of stack of vectors containing mirrored elements
+    """
+    if len(v.shape) == 1:
+        v = v[np.newaxis, :]
+    return np.concatenate((v, v[:, :0:-1]), axis=1)

micarray/util.py

@@ -91,3 +91,19 @@ def matdiagmul(A, b):
     for k in range(K):
         C[k, :, :] = A * b[k, :]
     return C
+
+
+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))