How to use the spaudiopy.sph.sh_matrix function in spaudiopy

def test_cpxSH(expected_Ynm):
    vecs = spa.grids.load_t_design(2*N)
    dirs = spa.utils.vecs2dirs(vecs)
    Y_nm = spa.sph.sh_matrix(N, dirs[:, 0], dirs[:, 1], SH_type='complex')
    assert(np.allclose(Y_nm, expected_Ynm))

def test_realSH(expected_Ynm):
    vecs = spa.grids.load_t_design(2*N)
    dirs = spa.utils.vecs2dirs(vecs)
    Y_nm = spa.sph.sh_matrix(N, dirs[:, 0], dirs[:, 1], SH_type='real')
    assert(np.allclose(Y_nm, expected_Ynm))

# %% Ambisonic decoding
# Ambisonic setup
N_e = ls_setup.get_characteristic_order()
ls_setup.ambisonics_setup(update_hull=True, N_kernel=20)

# Show ALLRAP hulls
plots.hull(ls_setup.ambisonics_hull, title='Ambisonic hull')
plots.hull(ls_setup.kernel_hull, mark_invalid=False, title='Kernel hull')

# ALLRAP
gains_allrap = decoder.allrap(src, ls_setup, N_sph=N_e)
# ALLRAP2
gains_allrap2 = decoder.allrap2(src, ls_setup, N_sph=N_e)
# ALLRAD
input_F_nm = sph.sh_matrix(N_e, src_azi, src_colat, 'real').T  # SH dirac
out_allrad = decoder.allrad(input_F_nm, ls_setup, N_sph=N_e)
out_allrad2 = decoder.allrad2(input_F_nm, ls_setup, N_sph=N_e)


utils.test_diff(gains_allrap, out_allrad, msg="ALLRAD and ALLRAP:")
utils.test_diff(gains_allrap2, out_allrad2, msg="ALLRAD2 and ALLRAP2:")

# Nearest Loudspeaker
gains_nls = decoder.nearest_loudspeaker(src, ls_setup)

# %% test multiple sources
_grid, _weights = grids.load_Fliege_Maier_nodes(10)
G_vbap = decoder.vbap(_grid, ls_setup)
G_allrap = decoder.allrap(_grid, ls_setup)
G_allrap2 = decoder.allrap2(_grid, ls_setup)
G_vbip = decoder.vbip(_grid, ls_setup)

# %%
# Grids
t = grids.load_t_design(2)
t_az, t_colat, t_r = utils.cart2sph(t[:, 0], t[:, 1], t[:, 2])
azi = t_az
colat = t_colat

# %%
# First, check condition number to which SH order the SHT is stable
# Tetraeder is not suited for SHT N>1:
sph.check_cond_sht(3, t_az, t_colat, 'real')

# %% Real and Complex SHs
Y_nm_c = sph.sh_matrix(N, azi, colat, 'complex')
Y_nm_r = sph.sh_matrix(N, azi, colat, 'real')

# %%
# Look at some SHTs
sig = np.array([1, 1, 1, 1])
sig_t = np.c_[np.eye(4), np.eye(4)]  # second axis s(t)
sig_B = sph.soundfield_to_b(sig)
F_nm = sph.sht(sig, N, azi, colat, SH_type='real')
F_nm_t = sph.sht(sig_t, N, azi, colat, SH_type='real')

# %%
F_nm_c = sph.sht(sig, N, azi, colat, SH_type='complex')
F_nm_c_t = sph.sht(sig_t, N, azi, colat, SH_type='complex')

# %%
F_nm_lst = sph.sht_lstsq(sig, N, azi, colat, SH_type='complex')
F_nm_lst_t = sph.sht_lstsq(sig_t, N, azi, colat, SH_type='real')

N_sph = hull.characteristic_order

    # virtual t-design loudspeakers
    J = len(kernel_hull.points)
    # virtual speakers expressed as VBAP phantom sources
    G = vbap(src=kernel_hull.points, hull=ambisonics_hull)

    # SH tapering coefficients
    a_n = sph.max_rE_weights(N_sph)
    a_n = sph.repeat_order_coeffs(a_n)

    # virtual Ambisonic decoder
    _t_azi, _t_colat, _t_r = utils.cart2sph(kernel_hull.points[:, 0],
                                            kernel_hull.points[:, 1],
                                            kernel_hull.points[:, 2])
    Y_td = sph.sh_matrix(N_sph, _t_azi, _t_colat, SH_type='real')
    # ALLRAD Decoder
    D = 4 * np.pi / J * G.T @ Y_td
    # apply tapering to decoder matrix
    D = D @ np.diag(a_n)
    # loudspeaker output signals
    ls_sig = D @ F_nm
    # remove imaginary loudspeakers
    ls_sig = np.delete(ls_sig, ambisonics_hull.imaginary_speaker, axis=0)
    return ls_sig

N = 1

# %%
# Grids
t = grids.load_t_design(2)
t_az, t_colat, t_r = utils.cart2sph(t[:, 0], t[:, 1], t[:, 2])
azi = t_az
colat = t_colat

# %%
# First, check condition number to which SH order the SHT is stable
# Tetraeder is not suited for SHT N>1:
sph.check_cond_sht(3, t_az, t_colat, 'real')

# %% Real and Complex SHs
Y_nm_c = sph.sh_matrix(N, azi, colat, 'complex')
Y_nm_r = sph.sh_matrix(N, azi, colat, 'real')

# %%
# Look at some SHTs
sig = np.array([1, 1, 1, 1])
sig_t = np.c_[np.eye(4), np.eye(4)]  # second axis s(t)
sig_B = sph.soundfield_to_b(sig)
F_nm = sph.sht(sig, N, azi, colat, SH_type='real')
F_nm_t = sph.sht(sig_t, N, azi, colat, SH_type='real')

# %%
F_nm_c = sph.sht(sig, N, azi, colat, SH_type='complex')
F_nm_c_t = sph.sht(sig_t, N, azi, colat, SH_type='complex')

# %%
F_nm_lst = sph.sht_lstsq(sig, N, azi, colat, SH_type='complex')