How to use the spaudiopy.sph function in spaudiopy

w_rE = sph.max_rE_weights(N)
# Choose here:
w_taper = w_Hann


# %% Spatial dirac in SH domain
dirac_azi = np.deg2rad(90)
dirac_colat = np.deg2rad(90)

# cross section
azi = np.linspace(0, 2 * np.pi, 720, endpoint=True)
colat = np.pi / 2 * np.ones_like(azi)

# Bandlimited Dirac pulse
dirac_untapered = 4 * np.pi / (N + 1) ** 2 * \
                  sph.bandlimited_dirac(N, azi - dirac_azi)
dirac_tapered = 4 * np.pi / (N + 1) ** 2 * \
                sph.bandlimited_dirac(N, azi - dirac_azi, w_n=w_taper)

# Coloration compensation of windowing
compensation_untapered = sph.binaural_coloration_compensation(N, f)
compensation_tapered = sph.binaural_coloration_compensation(N, f,
                                                            w_taper=w_taper)

# Get an FIR filter
ntaps = 128 + 1
assert (ntaps % 2), "Does not produce uneven number of filter taps."
filter_taps_untapered = firwin2(ntaps, f / (fs // 2), compensation_untapered)
filter_taps_tapered = firwin2(ntaps, f / (fs // 2), compensation_tapered)

# %% --- PLOTS ---
plots.polar(azi, dirac_untapered, title='Dirac untapered')

'../data/FABIAN_HRTF_DATABASE_V2/1 HRIRs/SphericalHarmonics/FABIAN_HRIR_measured_HATO_0.mat')

# Extracting the data is a bit ugly here...
SamplingRate = int(file['SamplingRate'])
SH_l = file['SH'][0][0][0]
SH_r = file['SH'][0][0][1]
f = np.squeeze(file['SH'][0][0][5])
print("Positive frequncy bins:", len(f))
print("SH shape:", SH_l.shape)
plt.plot(np.abs(SH_l[:, 8]))
plt.xlabel('SH coefficient')
plt.ylabel('Amplitude')
plt.title('Left HRTF SH at f={:.3} Hz'.format(f[8]))

# %% Inverse SHT
HRTF_l = sph.inverse_sht(SH_l, azi, colat, 'complex')
HRTF_r = sph.inverse_sht(SH_r, azi, colat, 'complex')

assert HRTF_l.shape == HRTF_r.shape
print("HRTF shape:", HRTF_l.shape)
plt_idx = int(HRTF_l.shape[0] / 2)
plots.f_amp(f, [HRTF_l[plt_idx, :],
                HRTF_r[plt_idx, :]],
            title=r"HRTF for $\phi={:.2}, \theta={:.2}$".format(
                azi[plt_idx], colat[plt_idx]),
            labels=['left', 'right'])

# %% [markdown]
# The inverse spherical harmonics transform renders from the spherical harmonics representation to the defined grid.
# Since the original input to the SHT are HRTFs we obtain again the time-frequency transfer functions.
#
# Now, from time-frequency domain to time domain HRIRs by the inverse Fourier transform.

el_steps = np.deg2rad(el_steps)
    phi_plot, theta_plot = np.meshgrid(np.arange(0., 2 * np.pi + azi_steps,
                                                 azi_steps),
                                       np.arange(10e-3, np.pi + el_steps,
                                                 el_steps))

    fig = plt.figure(figsize=plt.figaspect(1 / N_plots),
                     constrained_layout=True)
    ax_l = []
    for i_p, ff in enumerate(F_l):
        F_nm = utils.asarray_1d(ff)
        F_nm = F_nm[:, np.newaxis]
        if SH_type is None:
            SH_type = 'complex' if np.iscomplexobj(F_nm) else 'real'

        f_plot = sph.inverse_sht(F_nm, phi_plot.ravel(), theta_plot.ravel(),
                                 SH_type)
        f_r = np.abs(f_plot)
        f_ang = np.angle(f_plot)

        x_plot, y_plot, z_plot = utils.sph2cart(phi_plot.ravel(),
                                                theta_plot.ravel(),
                                                f_r.ravel())

        ax = fig.add_subplot(1, N_plots, i_p + 1, projection='3d')

        m = cm.ScalarMappable(cmap=cm.hsv,
                              norm=colors.Normalize(vmin=-np.pi, vmax=np.pi))
        m.set_array(f_ang)
        c = m.to_rgba(f_ang.reshape(phi_plot.shape))

        ax.plot_surface(x_plot.reshape(phi_plot.shape),

def sh_coeffs(F_nm, SH_type=None, azi_steps=5, el_steps=3, title=None):
    """Plot spherical harmonics coefficients as function on the sphere."""
    F_nm = utils.asarray_1d(F_nm)
    F_nm = F_nm[:, np.newaxis]
    if SH_type is None:
        SH_type = 'complex' if np.iscomplexobj(F_nm) else 'real'

    azi_steps = np.deg2rad(azi_steps)
    el_steps = np.deg2rad(el_steps)
    phi_plot, theta_plot = np.meshgrid(np.arange(0., 2 * np.pi + azi_steps,
                                                 azi_steps),
                                       np.arange(10e-3, np.pi + el_steps,
                                                 el_steps))

    f_plot = sph.inverse_sht(F_nm, phi_plot.ravel(), theta_plot.ravel(),
                             SH_type)
    f_r = np.abs(f_plot)
    f_ang = np.angle(f_plot)

    x_plot, y_plot, z_plot = utils.sph2cart(phi_plot.ravel(),
                                            theta_plot.ravel(),
                                            f_r.ravel())

    fig = plt.figure()
    ax = fig.gca(projection='3d')

    m = cm.ScalarMappable(cmap=cm.hsv,
                          norm=colors.Normalize(vmin=-np.pi, vmax=np.pi))
    m.set_array(f_ang)
    c = m.to_rgba(f_ang.reshape(phi_plot.shape))

raise ValueError('Run hull.setup_for_ambisonic() first!')
    if N_sph is None:
        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)
    # sqrt(E) normalization (eq.6)
    a_w = np.sqrt(np.sum((2 * (np.arange(N_sph + 1) + 1)) / (4 * np.pi) *
                         a_n**2))
    a_n /= a_w
    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')

    # tapered dirac
    Y_k = Y_td @ np.diag(a_n) @ F_nm
    # ALLRAD2 Decoder
    ls_sig = np.sqrt(4 * np.pi / J * np.square(G.T) @ np.square(Y_k))

    # remove imaginary loudspeakers
    ls_sig = np.delete(ls_sig, ambisonics_hull.imaginary_speaker, axis=0)
    return ls_sig

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import firwin2
from spaudiopy import sph, utils, plots
from spaudiopy import process as pcs

# sampling rate (Hz)
fs = 48000
# evaluated frequencies, 1000 points
f = np.linspace(0, fs / 2, 1000)
# target spherical harmonics order N (>= 3)
N = 5

# tapering windows
w_Hann = pcs.half_sided_Hann(N)
w_rE = sph.max_rE_weights(N)
# Choose here:
w_taper = w_Hann


# %% Spatial dirac in SH domain
dirac_azi = np.deg2rad(90)
dirac_colat = np.deg2rad(90)

# cross section
azi = np.linspace(0, 2 * np.pi, 720, endpoint=True)
colat = np.pi / 2 * np.ones_like(azi)

# Bandlimited Dirac pulse
dirac_untapered = 4 * np.pi / (N + 1) ** 2 * \
                  sph.bandlimited_dirac(N, azi - dirac_azi)
dirac_tapered = 4 * np.pi / (N + 1) ** 2 * \

def characteristic_ambisonic_order(hull):
    """Zotter, F., & Frank, M. (2012). All-Round Ambisonic Panning and
    Decoding. Journal of Audio Engineering Society, Sec. 7.
    """
    _hull = copy.copy(hull)
    # projection for loudspeakers not on unit sphere
    _xp, _yp, _zp = sph.project_on_sphere(_hull.points[:, 0],
                                          _hull.points[:, 1],
                                          _hull.points[:, 2],)
    _hull.points = np.c_[_xp, _yp, _zp]
    # project centroids (for non-uniform setups)
    src = np.asarray(sph.project_on_sphere(hull.centroids[:, 0],
                                           hull.centroids[:, 1],
                                           hull.centroids[:, 2])).T
    # all loudspeaker triplets enclosing listener
    gains = vbap(src, _hull, valid_simplices=_hull._encloses_listener)
    # Energy vector of each center
    rE, rE_mag = sph.r_E(_hull.points, gains)
    # eq. (16)
    spread = 2 * np.arccos(rE_mag) * (180 / np.pi)
    N_e = 2 * 137.9 / np.average(spread) - 1.51
    # ceil might be optimistic
    return int(np.ceil(N_e))

def sh_to_b(cls, multisig):
        """Alternative constructor, convert from sig.Multisignal.
        Assumes ACN channel order.
        """
        assert isinstance(multisig, MultiSignal)
        _B = sph.sh_to_b(multisig.copy().get_signals())
        return cls([*_B], fs=multisig.fs)

# %% Spherical Harmonics Order
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')