Merge pull request #8 from spatialaudio/consistent_usage_of_transform_matrices

Fix inconsistent usage of transform matrices

Author: trettberg

examples/modal_beamforming_open_array.py

@@ -35,7 +35,7 @@ def dot_product_sph(v, u):
 # compute microphone signals for an incident broad-band plane wave
 p = np.exp(1j * k[:, np.newaxis]*r * dot_product_sph((azi, elev), pw_angle))
 # compute the plane wave dcomposition
-A_pwd = np.matmul(np.matmul(np.conj(Y_q.T), D), Y_p)
+A_pwd = np.matmul(np.matmul(Y_q, D), np.conj(Y_p.T))
 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)
 
@@ -46,11 +46,11 @@ def dot_product_sph(v, u):
 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('spherical_modal_open_beamformer_pwd_fd.png')
 
 plt.figure()
 plt.pcolormesh(range(2*len(k)-2), azi_pwd/np.pi, 20*np.log10(np.abs(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('spherical_modal_open_beamformer_pwd_td.png')

examples/modal_beamforming_open_circular_array.py

@@ -18,12 +18,12 @@
 # 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
+# pressure on the surface of an open 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.matmul(np.matmul(Psi_p, D), np.conj(Psi_pw.T))
 p = np.squeeze(p)
 
 # incident plane wave exhibiting infinite spatial bandwidth
@@ -35,7 +35,7 @@
 D = micarray.modal.radial.circ_diagonal_mode_mat(Dn)
 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)
+A_pwd = np.matmul(np.matmul(Psi_q, D), np.conj(Psi_p.T))
 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)
 

examples/modal_beamforming_rigid_array.py

@@ -22,7 +22,7 @@
 D = micarray.modal.radial.diagonal_mode_mat(bn)
 Y_p = micarray.modal.angular.sht_matrix(N, azi, elev)
 Y_pw = micarray.modal.angular.sht_matrix(N, azi_pw, np.pi/2)
-p = np.matmul(np.matmul(np.conj(Y_pw.T), D), Y_p)
+p = np.matmul(np.matmul(Y_p, D), np.conj(Y_pw.T))
 p = np.squeeze(p)
 
 # plane wave decomposition using modal beamforming
@@ -35,7 +35,7 @@
 dn, _ = micarray.modal.radial.regularize(1/bn, 100, 'softclip')
 D = micarray.modal.radial.diagonal_mode_mat(dn)
 # compute the PWD
-A_mb = np.matmul(np.matmul(np.conj(Y_q.T), D), Y_p)
+A_mb = np.matmul(np.matmul(Y_q, D), np.conj(Y_p.T))
 q_mb = np.squeeze(np.matmul(A_mb, np.expand_dims(p, 2)))
 q_mb_t = np.fft.fftshift(np.fft.irfft(q_mb, axis=0), axes=0)
 

examples/modal_beamforming_rigid_circular_array.py

@@ -23,7 +23,7 @@
 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.matmul(np.matmul(Psi_p, D), np.conj(Psi_pw.T))
 p = np.squeeze(p)
 
 # plane wave decomposition using modal beamforming
@@ -32,7 +32,7 @@
 Bn = micarray.modal.radial.circular_pw(N, k, r, setup='rigid')
 Dn, _ = micarray.modal.radial.regularize(1/Bn, 100, 'softclip')
 D = micarray.modal.radial.circ_diagonal_mode_mat(Dn)
-A_pwd = np.matmul(np.matmul(np.conj(Psi_q.T), D), Psi_p)
+A_pwd = np.matmul(np.matmul(Psi_q, D), np.conj(Psi_p.T))
 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)
 

micarray/modal/angular.py

@@ -12,7 +12,7 @@ def sht_matrix(N, azi, elev, weights=None):
 
     .. math::
 
-        \mathbf{Y} = \left[ \begin{array}{cccccc} 
+        \mathbf{Y} = \left[ \begin{array}{cccccc}
         Y_0^0(\theta[0], \phi[0]) & Y_1^{-1}(\theta[0], \phi[0]) & Y_1^0(\theta[0], \phi[0]) & Y_1^1(\theta[0], \phi[0]) & \dots & Y_N^N(\theta[0], \phi[0])  \\
         Y_0^0(\theta[1], \phi[1]) & Y_1^{-1}(\theta[1], \phi[1]) & Y_1^0(\theta[1], \phi[1]) & Y_1^1(\theta[1], \phi[1]) & \dots & Y_N^N(\theta[1], \phi[1])  \\
         \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\
@@ -25,6 +25,10 @@ def sht_matrix(N, azi, elev, weights=None):
 
         Y_n^m(\theta, \phi) = \sqrt{\frac{2n + 1}{4 \pi} \frac{(n-m)!}{(n+m)!}} P_n^m(\cos \theta) e^{i m \phi}
 
+
+    (Note: :math:`\mathbf{Y}` is interpreted as the inverse transform (or synthesis)
+    matrix in examples and documentation.)
+
     Parameters
     ----------
     N : int
@@ -38,22 +42,23 @@ def sht_matrix(N, azi, elev, weights=None):
 
     Returns
     -------
-    Ymn : ((N+1)**2, Q) numpy.ndarray
+    Ymn : (Q, (N+1)**2) numpy.ndarray
         Matrix of spherical harmonics.
+
     """
     azi = util.asarray_1d(azi)
     elev = util.asarray_1d(elev)
     if azi.ndim == 0:
-        M = 1
+        Q = 1
     else:
-        M = len(azi)
+        Q = len(azi)
     if weights is None:
-        weights = np.ones(M)
-    Ymn = np.zeros([(N+1)**2, M], dtype=complex)
+        weights = np.ones(Q)
+    Ymn = np.zeros([Q, (N+1)**2], dtype=complex)
     i = 0
     for n in range(N+1):
         for m in range(-n, n+1):
-            Ymn[i, :] = weights * special.sph_harm(m, n, azi, elev)
+            Ymn[:, i] = weights * special.sph_harm(m, n, azi, elev)
             i += 1
     return Ymn
 
@@ -105,6 +110,10 @@ def cht_matrix(N, pol, weights=None):
         1 & e^{i\varphi[Q-1]} & \cdots & e^{iN\varphi[Q-1]} & e^{-iN\varphi[Q-1]} & \cdots & e^{-i\varphi[Q-1]}
         \end{array} \right]
 
+    (Note: :math:`\Psi` is interpreted as the inverse transform (or synthesis)
+    matrix in examples and documentation.)
+
+
     Parameters
     ----------
     N : int
@@ -116,8 +125,9 @@ def cht_matrix(N, pol, weights=None):
 
     Returns
     -------
-    Psi : (2N+1, Q) numpy.ndarray
+    Psi : (Q, 2N+1) numpy.ndarray
         Matrix of circular harmonics.
+
     """
     pol = util.asarray_1d(pol)
     if pol.ndim == 0:
@@ -126,10 +136,10 @@ def cht_matrix(N, pol, weights=None):
         Q = len(pol)
     if weights is None:
         weights = np.ones(Q)
-    Psi = np.zeros([(2*N+1), Q], dtype=complex)
+    Psi = np.zeros([Q, (2*N+1)], dtype=complex)
     order = np.roll(np.arange(-N, N+1), -N)
     for i, n in enumerate(order):
-        Psi[i, :] = weights * np.exp(1j * n * pol)
+        Psi[:, i] = weights * np.exp(1j * n * pol)
     return Psi