add 'db' in 'util.py'

Author: narahahn

examples/modal_beamforming_open_circular_array.py

@@ -12,7 +12,7 @@
 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
-M = 181 # number of microphones
+M = 181  # number of microphones
 
 # get uniform grid (microphone positions) of number M
 pol, weights = micarray.modal.angular.grid_equal_polar_angle(M)
@@ -26,23 +26,23 @@
 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))
+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, 20*np.log10(np.abs(q_pwd.T)), vmin=-40)
+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, 20*np.log10(np.abs(q_pwd_t.T)), vmin=-40)
+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)')

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