1010N = 20 # order of modal beamformer/microphone array
1111pw_angle = (np .pi , np .pi / 2 ) # incidence angle of plane wave
1212azi_pwd = np .linspace (0 , 2 * np .pi , 91 , endpoint = False ) # angles for plane wave decomposition
13- kr = np .linspace (0.1 , 20 , 100 ) # wavenumber-radius vector
13+ k = np .linspace (0.1 , 20 , 100 ) # wavenumber vector
14+ r = 1 # radius of array
1415
1516
1617def dot_product_sph (v , u ):
1718 # evaluate dot-product between u and v in spherical coordinates
18- return (np .cos (v [0 ])* np .sin (v [1 ])* np .cos (u [0 ])* np .sin (u [1 ])
19- + np .sin (v [0 ])* np .sin (v [1 ])* np .sin (u [0 ])* np .sin (u [1 ])
20- + np .cos (v [1 ])* np .cos (u [1 ]))
19+ return (np .cos (v [0 ])* np .sin (v [1 ])* np .cos (u [0 ])* np .sin (u [1 ]) +
20+ np .sin (v [0 ])* np .sin (v [1 ])* np .sin (u [0 ])* np .sin (u [1 ]) +
21+ np .cos (v [1 ])* np .cos (u [1 ]))
2122
2223
2324# get quadrature grid (microphone positions) of order N
@@ -27,28 +28,28 @@ def dot_product_sph(v, u):
2728# get spherical harmonics matrix for a source ensemble of azimuthal plane waves
2829Y_q = micarray .modal .angular .sht_matrix (N , azi_pwd , np .pi / 2 )
2930# get radial filters
30- bn = micarray .modal .radial .spherical (N , kr , setup = 'open' , plane_wave = True )
31+ bn = micarray .modal .radial .spherical_pw (N , k , r , setup = 'open' )
3132dn , _ = micarray .modal .radial .regularize (1 / bn , 100 , 'softclip' )
3233D = micarray .modal .radial .diagonal_mode_mat (dn )
3334
3435# compute microphone signals for an incident broad-band plane wave
35- p = np .exp (1j * kr [:, np .newaxis ] * dot_product_sph ((azi , elev ), pw_angle ))
36+ p = np .exp (1j * k [:, np .newaxis ]* r * dot_product_sph ((azi , elev ), pw_angle ))
3637# compute the plane wave dcomposition
3738A_pwd = np .matmul (np .matmul (np .conj (Y_q .T ), D ), Y_p )
3839q_pwd = np .squeeze (np .matmul (A_pwd , np .expand_dims (p , 2 )))
3940q_pwd_t = np .fft .fftshift (np .fft .irfft (q_pwd , axis = 0 ), axes = 0 )
4041
4142# visualize plane wave decomposition (aka beampattern)
4243plt .figure ()
43- plt .pcolormesh (kr , azi_pwd / np .pi , 20 * np .log10 (np .abs (q_pwd .T )), vmin = - 10 )
44+ plt .pcolormesh (k , azi_pwd / np .pi , 20 * np .log10 (np .abs (q_pwd .T )), vmin = - 40 )
4445plt .colorbar ()
4546plt .xlabel (r'$kr$' )
4647plt .ylabel (r'$\phi / \pi$' )
4748plt .title ('Plane wave docomposition by modal beamformer (frequency domain)' )
4849plt .savefig ('modal_open_beamformer_pwd_fd.png' )
4950
5051plt .figure ()
51- plt .pcolormesh (range (2 * len (kr )- 2 ), azi_pwd / np .pi , 20 * np .log10 (np .abs (q_pwd_t .T )), vmin = - 10 )
52+ plt .pcolormesh (range (2 * len (k )- 2 ), azi_pwd / np .pi , 20 * np .log10 (np .abs (q_pwd_t .T )), vmin = - 40 )
5253plt .colorbar ()
5354plt .ylabel (r'$\phi / \pi$' )
5455plt .title ('Plane wave docomposition by modal beamformer (time domain)' )
0 commit comments