Extrapolatoin

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')

examples/sound_field_extrapolation_open_circular_array_mono.py

@@ -0,0 +1,73 @@
+"""
+    Modal analysis and extrapolation of a monochromatic sound field
+    in the cricular harmonics domain using an open circular array
+"""
+import numpy as np
+import micarray
+from micarray.util import db
+import scipy.special as special
+import matplotlib.pyplot as plt
+import matplotlib
+matplotlib.rcParams['contour.negative_linestyle'] = 'solid'
+
+# Constants
+c = 343  # speed of sound [m/s]
+
+# 2-dimensional grid
+spacing = 0.01
+x = np.expand_dims(np.arange(-1, 1, spacing), axis=0)
+y = np.expand_dims(np.arange(-1, 1, spacing), axis=1)
+r = np.sqrt(x**2+y**2).astype(complex)
+phi = np.arctan2(y, x)
+
+# Incident plane wave
+phi_pw = 0.5*np.pi  # incoming direction
+f = 1500  # temporal frequency
+k = micarray.util.asarray_1d(2*np.pi*f/c)  # corresponding wave number
+s0 = np.exp(1j*k*r*np.cos(phi-phi_pw))  # incident sound field
+
+# Microphone array and modal analysis
+N = 15  # maximum order
+order = np.roll(np.arange(-N, N+1), -N)
+threshold = 1e5  # regulaization parameter
+R = 0.5  # radius
+Phi, weights = micarray.modal.angular.grid_equal_polar_angle(N)  # array
+p = np.exp(1j*k*R*np.cos(Phi-phi_pw))  # captured signal
+bn = micarray.modal.radial.circ_radial_weights(N, k*R, setup='open')
+dn, _ = micarray.modal.radial.regularize(1/bn, threshold, 'softclip')
+pm = dn * np.fft.ifft(p)
+
+# Sound field extrapolation
+basis = special.jn(order[:, np.newaxis, np.newaxis], k * r[np.newaxis, :, :]) \
+        * np.exp(-1j*order[:, np.newaxis, np.newaxis] * phi[np.newaxis, :, :])
+s = np.tensordot(pm, basis, axes=[0, 0])
+
+
+# Plots
+plt.figure(figsize=(4, 4))
+plt.pcolormesh(x, y, np.real(s), cmap='coolwarm')
+plt.plot(R*np.cos(Phi), R*np.sin(Phi), 'k.')
+plt.axis('scaled')
+plt.axis([-1, 1, -1, 1])
+cb = plt.colorbar(fraction=0.046, pad=0.04)
+plt.clim(-1, 1)
+plt.xlabel('$x$ / m')
+plt.ylabel('$y$ / m')
+plt.title('Extrapolated Sound Field')
+plt.savefig('extrapolation_open_circ_mono.png')
+
+plt.figure(figsize=(4, 4))
+plt.pcolormesh(x, y, db(s0-s), cmap='Blues', vmin=-60)
+plt.plot(R*np.cos(Phi), R*np.sin(Phi), 'k.')
+plt.axis('scaled')
+plt.axis([-1, 1, -1, 1])
+cb = plt.colorbar(fraction=0.046, pad=0.04)
+cb.set_label('dB')
+plt.clim(-60, 0)
+plt.xlabel('$x$ / m')
+plt.ylabel('$y$ / m')
+xx, yy = np.meshgrid(x, y)
+cs = plt.contour(xx, yy, db(s0-s), np.arange(-60, 20, 20), colors='orange')
+plt.clabel(cs, fontsize=9, inline=1, fmt='%1.0f')
+plt.title('Extrapolation Error')
+plt.savefig('extrapolation_error_open_circ_mono.png')

micarray/modal/radial.py

@@ -176,7 +176,7 @@ def regularize(dn, a0, method):
     else:
         raise ValueError('method must be either: none, ' +
                          'discard, hardclip, softclip, Tikh or wng')
-    dn[0, 1:] = dn[1, 1:]
+#    dn[0, 1:] = dn[1, 1:]
     dn = dn * hn
     if not np.isfinite(dn).all():
         raise UserWarning("Filter not finite")
@@ -258,7 +258,7 @@ def circular_pw(N, k, r, setup):
         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
@@ -294,7 +294,7 @@ def circular_ls(N, k, r, rs, setup):
     """
     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:
@@ -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,7 +367,6 @@ 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