Enable flexible choice of the modal bandwidth in the examples

* The modal bandwidth of the incident sound field `Nsf` may differ from that of the modal beamforming `N`
* For open circular arrays, the incident sound can have either a finite or an infinite modal bandwidth

Author: narahahn

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

@@ -6,6 +6,7 @@
 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
@@ -18,15 +19,15 @@
 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')