From cb99ad746103f172c18a40e6329538d79715b294 Mon Sep 17 00:00:00 2001 From: Nara Hahn Date: Sun, 3 Dec 2017 23:03:43 +0100 Subject: [PATCH 1/4] Generate generalized spiral points on a sphere --- examples/generalized_sprial_points.py | 20 ++++++++++++++++++ micarray/modal/angular.py | 30 +++++++++++++++++++++++++++ 2 files changed, 50 insertions(+) create mode 100644 examples/generalized_sprial_points.py diff --git a/examples/generalized_sprial_points.py b/examples/generalized_sprial_points.py new file mode 100644 index 0000000..c9b83a7 --- /dev/null +++ b/examples/generalized_sprial_points.py @@ -0,0 +1,20 @@ +"""Compute the generalized sprial points on a sphere.""" +import numpy as np +import matplotlib.pyplot as plt +from mpl_toolkits.mplot3d import Axes3D +import micarray + +M = 700 +azi, elev, _ = micarray.modal.angular.grid_generalized_spiral(M, C=3.6) + +x = np.sin(elev) * np.cos(azi) +y = np.sin(elev) * np.sin(azi) +z = np.cos(elev) + +fig = plt.figure(figsize=(8, 8)) +ax = fig.gca(projection='3d') +ax.scatter(x, y, z) +ax.set_xlabel('$\hat{x}$') +ax.set_ylabel('$\hat{y}$') +ax.set_zlabel('$\hat{z}$') +ax.set_title('Generalized Spiral Points ($M={}$)'.format(M)) diff --git a/micarray/modal/angular.py b/micarray/modal/angular.py index 81d9c55..b30ed07 100644 --- a/micarray/modal/angular.py +++ b/micarray/modal/angular.py @@ -197,6 +197,36 @@ def grid_gauss(n): return azi, elev, weights +def grid_generalized_spiral(M, C=3.6): + """Generalized spiral points on sphere. + + According to (cf. Rakhmanov 1994, sec. 5). + + Parameters + ---------- + M : int + Number of point on the sphere. + + Returns + ------- + azi : array_like + Azimuth. + elev : array_like + Elevation. + weights : array_like + Quadrature weights. + + """ + k = np.arange(1, M+1) + h = -1 + 2*(k-1)/(M-1) + elev = np.arccos(h) + azi = np.zeros(M) + for n in k[:-1]: + azi[n] = azi[n-1] + C/np.sqrt(M)*1/np.sqrt(1-h[n]**2) + weights = np.ones(M) + return azi, elev, weights + + def grid_equal_polar_angle(n): """Equi-angular sampling points on a circle. From 967773070be9d0a85f99625b2ac2629d833d31ff Mon Sep 17 00:00:00 2001 From: Nara Hahn Date: Mon, 4 Dec 2017 09:43:54 +0100 Subject: [PATCH 2/4] Fix computation for k=M --- micarray/modal/angular.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/micarray/modal/angular.py b/micarray/modal/angular.py index b30ed07..b4bdecf 100644 --- a/micarray/modal/angular.py +++ b/micarray/modal/angular.py @@ -221,7 +221,7 @@ def grid_generalized_spiral(M, C=3.6): h = -1 + 2*(k-1)/(M-1) elev = np.arccos(h) azi = np.zeros(M) - for n in k[:-1]: + for n in k[:-2]: azi[n] = azi[n-1] + C/np.sqrt(M)*1/np.sqrt(1-h[n]**2) weights = np.ones(M) return azi, elev, weights From fbe4748e5530c3a99990579251621ab5dd8ee05d Mon Sep 17 00:00:00 2001 From: Nara Hahn Date: Mon, 4 Dec 2017 11:48:04 +0100 Subject: [PATCH 3/4] Don't return quadrature weight --- micarray/modal/angular.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/micarray/modal/angular.py b/micarray/modal/angular.py index b4bdecf..8ae226e 100644 --- a/micarray/modal/angular.py +++ b/micarray/modal/angular.py @@ -223,7 +223,7 @@ def grid_generalized_spiral(M, C=3.6): azi = np.zeros(M) for n in k[:-2]: azi[n] = azi[n-1] + C/np.sqrt(M)*1/np.sqrt(1-h[n]**2) - weights = np.ones(M) + weights = None return azi, elev, weights From 6afeee82f7b172afc9198d9a98e24c65b0947608 Mon Sep 17 00:00:00 2001 From: Nara Hahn Date: Mon, 4 Dec 2017 11:49:50 +0100 Subject: [PATCH 4/4] Robustness analysis of the generalized spiral points --- examples/generalized_sprial_points.py | 42 ++++++++++++++++++++++++--- 1 file changed, 38 insertions(+), 4 deletions(-) diff --git a/examples/generalized_sprial_points.py b/examples/generalized_sprial_points.py index c9b83a7..9d13f12 100644 --- a/examples/generalized_sprial_points.py +++ b/examples/generalized_sprial_points.py @@ -3,18 +3,52 @@ import matplotlib.pyplot as plt from mpl_toolkits.mplot3d import Axes3D import micarray +from micarray.util import db -M = 700 + +def sph2cart(alpha, beta, r): + """Spherical to cartesian coordinates.""" + x = r * np.cos(alpha) * np.sin(beta) + y = r * np.sin(alpha) * np.sin(beta) + z = r * np.cos(beta) + return x, y, z + + +# Microphone array +R = 0.5 # radius +N = 6 # modal bandwidth +M = 1*(N+1)**2 # number of microphones azi, elev, _ = micarray.modal.angular.grid_generalized_spiral(M, C=3.6) +x, y, z = sph2cart(azi, elev, R) +Y = micarray.modal.angular.sht_matrix(N, azi, elev) # synthesis matrix +Y_inv = np.linalg.pinv(Y) # analysis matrix +k = np.linspace(0.1, 40, 100) # wavenumber +bn = micarray.modal.radial.spherical_pw(N, k, R, setup='open') +D = micarray.modal.radial.diagonal_mode_mat(bn) +B = np.matmul(D, Y) +condnum = np.linalg.cond(B) # Condition number -x = np.sin(elev) * np.cos(azi) -y = np.sin(elev) * np.sin(azi) -z = np.cos(elev) +# Fig. Microphone array fig = plt.figure(figsize=(8, 8)) ax = fig.gca(projection='3d') +ax.plot(x, y, z, c='lightgray') ax.scatter(x, y, z) ax.set_xlabel('$\hat{x}$') ax.set_ylabel('$\hat{y}$') ax.set_zlabel('$\hat{z}$') ax.set_title('Generalized Spiral Points ($M={}$)'.format(M)) + +# Fig. Pseudo inverse matrix +plt.figure() +plt.pcolormesh(db(np.matmul(Y_inv, Y))) +plt.axis('scaled') +plt.colorbar(label='dB') +plt.title(r'$E = Y^\dagger Y$') + +# Fig. Condition number +plt.figure() +plt.semilogy(k*R, condnum) +plt.xlabel('$kr$') +plt.ylabel('Condition number') +plt.ylim(top=10e4)