Merge pull request #5 from spatialaudio/sfa_circular_harmonics

Radial weight vectors in CHT

Author: trettberg

examples/modal_beamforming_open_circular_array.py

@@ -6,43 +6,51 @@
 import numpy as np
 import matplotlib.pyplot as plt
 import micarray
+from micarray.util import db
 
-N = 90  # order of modal beamformer/microphone array
+Nsf = 50  # order of the incident sound field
+N = 30  # 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
+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)
-# 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
+
+# pressure on the surface of a rigid cylinder for an incident plane wave
+Bn = micarray.modal.radial.circular_pw(Nsf, k, r, setup='open')
+D = micarray.modal.radial.circ_diagonal_mode_mat(Bn)
+Psi_p = micarray.modal.angular.cht_matrix(Nsf, pol)
+Psi_pw = micarray.modal.angular.cht_matrix(Nsf, pw_angle)
+p = np.matmul(np.matmul(np.conj(Psi_p.T), D), Psi_pw)
+p = np.squeeze(p)
+
+# incident plane wave exhibiting infinite spatial bandwidth
+# p = np.exp(1j * k[:, np.newaxis]*r * np.cos(pol - pw_angle))
+
+# plane wave decomposition using modal beamforming
 Bn = micarray.modal.radial.circular_pw(N, k, r, setup='open')
-Dn, _ = micarray.modal.radial.regularize(1/Bn, 100, 'softclip')
+Dn, _ = micarray.modal.radial.regularize(1/Bn, 3000, '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
+Psi_p = micarray.modal.angular.cht_matrix(N, pol, weights)
+Psi_q = micarray.modal.angular.cht_matrix(N, pol_pwd)
 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.pcolormesh(k, pol_pwd/np.pi, 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.savefig('modal_circ_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.pcolormesh(range(2*len(k)-2), pol_pwd/np.pi, 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')
+plt.savefig('modal_circ_open_beamformer_pwd_td.png')

examples/modal_beamforming_rigid_circular_array.py

@@ -1,32 +1,33 @@
 """
     Compute the plane wave decomposition for an incident broadband plane wave
-    on an rigid circular array using a modal beamformer of finite order.
+    on a rigid circular array using a modal beamformer of finite order.
 """
 
 import numpy as np
 import matplotlib.pyplot as plt
 import micarray
+from micarray.util import db
 
 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
+pol_pwd = np.linspace(0, 2*np.pi, 180, endpoint=False)  # angles for PWD
 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)
+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)
 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.matmul(np.matmul(np.conj(Psi_p.T), D), Psi_pw)
 p = np.squeeze(p)
 
 # plane wave decomposition using modal beamforming
-Psi_p = micarray.modal.angular.cht_matrix(N, pol)
+Psi_p = micarray.modal.angular.cht_matrix(N, pol, weights)
 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')
@@ -37,16 +38,16 @@
 
 # visualize plane wave decomposition (aka beampattern)
 plt.figure()
-plt.pcolormesh(k, pol_pwd/np.pi, micarray.util.db(q_pwd.T), vmin=-40)
+plt.pcolormesh(k, pol_pwd/np.pi, 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.savefig('modal_circ_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.pcolormesh(range(2*len(k)-2), pol_pwd/np.pi, 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')
+plt.savefig('modal_circ_open_beamformer_pwd_td.png')

micarray/modal/radial.py

@@ -239,7 +239,7 @@ def circular_pw(N, k, r, setup):
 
     .. math::
 
-        \mathring{P}_n(k) = i^{-n} b_n(kr)
+        \mathring{P}_n(k) = i^n b_n(kr)
 
     Parameters
     ----------
@@ -254,14 +254,14 @@ def circular_pw(N, k, r, setup):
 
     Returns
     -------
-    bn : (M, N+1) numpy.ndarray
+    bn : (M, 2*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)
+    n = np.roll(np.arange(-N, N+1), -N)
 
     bn = circ_radial_weights(N, kr, setup)
-    return (1j)**(n) * bn
+    return 1j**n * bn
 
 
 def circular_ls(N, k, r, rs, setup):
@@ -289,20 +289,20 @@ def circular_ls(N, k, r, rs, setup):
 
     Returns
     -------
-    bn : (M, N+1) numpy.ndarray
+    bn : (M, 2*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)
+    n = np.roll(np.arange(-N, N+1), -N)
 
     bn = circ_radial_weights(N, k*r, setup)
     if len(k) == 1:
         bn = bn[np.newaxis, :]
     for i, x in enumerate(krs):
         Hn = special.hankel2(n, x)
-        bn[i, :] = bn[i, :] * -1j/4 * Hn
-    return np.squeeze(bn)
+        bn[i, :] = bn[i, :] * Hn
+    return -1j/4 * np.squeeze(bn)
 
 
 def circ_radial_weights(N, kr, setup):
@@ -327,7 +327,7 @@ def circ_radial_weights(N, kr, setup):
 
     Returns
     -------
-    bn : (M, N+1) numpy.ndarray
+    bn : (M, 2*N+1) numpy.ndarray
         Radial weights for all orders up to N and the given wavenumbers.
 
     """
@@ -348,6 +348,7 @@ def circ_radial_weights(N, kr, setup):
         else:
             raise ValueError('setup must be either: open, card or rigid')
         Bns[i, :] = bn
+    Bns = np.concatenate((Bns, (Bns*(-1)**np.arange(N+1))[:, :0:-1]), axis=-1)
     return np.squeeze(Bns)
 
 
@@ -366,35 +367,10 @@ def circ_diagonal_mode_mat(bk):
         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)