diff --git a/examples/modal_beamforming_open_circular_array.py b/examples/modal_beamforming_open_circular_array.py index 8d91eda..dfc07eb 100644 --- a/examples/modal_beamforming_open_circular_array.py +++ b/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') diff --git a/examples/modal_beamforming_rigid_circular_array.py b/examples/modal_beamforming_rigid_circular_array.py index 1597bdc..1c930b7 100644 --- a/examples/modal_beamforming_rigid_circular_array.py +++ b/examples/modal_beamforming_rigid_circular_array.py @@ -1,16 +1,17 @@ """ 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 @@ -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') diff --git a/examples/sound_field_extrapolation_open_circular_array_mono.py b/examples/sound_field_extrapolation_open_circular_array_mono.py new file mode 100644 index 0000000..76cfc82 --- /dev/null +++ b/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') diff --git a/micarray/modal/radial.py b/micarray/modal/radial.py index ae1c408..75f78f9 100644 --- a/micarray/modal/radial.py +++ b/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