How to use the spaudiopy.plots function in spaudiopy

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')
plots.polar(azi, dirac_tapered, title='Dirac tapered')
plots.spectrum([filter_taps_untapered, filter_taps_tapered], fs, scale_mag=True,
               title='Coloration Equalization', labels=['untapered', 'tapered'])

plt.show()

# 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')
plots.polar(azi, dirac_tapered, title='Dirac tapered')
plots.spectrum([filter_taps_untapered, filter_taps_tapered], fs, scale_mag=True,
               title='Coloration Equalization', labels=['untapered', 'tapered'])

plt.show()

print("Compensated HRIR:", hrir_l_hp.shape)

freq = np.fft.rfftfreq(hrir_l_hp.shape[1], d=1. / SamplingRate)
plots.f_amp(freq, [np.fft.rfft(hrir_l_hp[plt_idx, :]),
                   np.fft.rfft(hrir_r_hp[plt_idx, :])],
            labels=['HRTF left', 'HRTF right'],
            title='Compensated HRTF')

# %% Resample to 48k
fs_target = 48000
hrir_l_hp48k, hrir_r_hp48k, _ = process.resample_HRIRs(hrir_l_hp, hrir_r_hp,
                                                       SamplingRate,
                                                       fs_target)
print("Resampled HRIR:", hrir_l_hp48k.shape)
freq = np.fft.rfftfreq(hrir_l_hp48k.shape[1], d=1. / SamplingRate)
plots.f_amp(freq, [np.fft.rfft(hrir_l_hp48k[plt_idx, :]),
                   np.fft.rfft(hrir_r_hp48k[plt_idx, :])],
            labels=['HRTF left', 'HRTF right'],
            title='Resampled HRTF')

# %% Save to .mat
savemat('../data/HRTF_default', {'hrir_l': hrir_l_hp,
                                 'hrir_r': hrir_r_hp,
                                 'azi': azi, 'elev': colat,
                                 'SamplingRate': SamplingRate})
savemat('../data/HRTF_default48k', {'hrir_l': hrir_l_hp48k,
                                   'hrir_r': hrir_r_hp48k,
                                   'azi': azi, 'elev': colat,
                                   'SamplingRate': fs_target})

# %%
plt.show()

def show(self):
        """Plot hull object."""
        plots.hull(self, title='Loudspeaker Setup')

print("Complex valued SHT of time signal:")
utils.test_diff(sig_t, f_c_t)
print("Real valued least-squares SHT of time signal:")
utils.test_diff(sig_t, f_lst_t)

# %%
# Check B format conversion
B_sig = np.array([1, 1, 0, 0])  # W, X, Y, Z
F_B = sph.b_to_sh(B_sig)
B_sig_re = sph.sh_to_b(F_B)
print("B format to SH conversion:")
utils.test_diff(B_sig, B_sig_re)

# %%
# Some plots
plots.sh_coeffs(F_nm, title="Ambeo: all channels max")
plots.sh_coeffs(F_B, title="b_to_sh: W+X")

# %%
plots.sh_coeffs_subplot([np.array([1, 0, 0, 0]),
                         np.array([0, 1, 0, 0]),
                         np.array([0, 0, 1, 0]),
                         np.array([0, 0, 0, 1])],
                        titles=["0", "1, -1", "1, 0", "1, 1"])

# %%
plots.sh_coeffs(np.sqrt(2) * np.array([1, 0, 0, -1]), 'complex',
                title="sqrt(2) * [1, 0, 0, -1] complex coeffs")

# %%
# Look at simple B format generator
sig2 = np.ones(8)

# %% 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)

# %% Look at some performance measures
plots.decoder_performance(ls_setup, 'NLS')
plots.decoder_performance(ls_setup, 'VBAP')
plots.decoder_performance(ls_setup, 'VBAP', retain_outside=True)
plt.suptitle('VBAP with imaginary loudspeaker')
plots.decoder_performance(ls_setup, 'VBIP', retain_outside=True)
plt.suptitle('VBIP with imaginary loudspeaker')
plots.decoder_performance(ls_setup, 'ALLRAP')
plots.decoder_performance(ls_setup, 'ALLRAP2')

# %% Binauralize
fs = 44100
hrirs = IO.load_hrirs(fs)

l_vbap_ir, r_vbap_ir = ls_setup.binauralize(ls_setup.loudspeaker_signals(
                                            gains_vbap), fs)

l_allrap_ir, r_allrap_ir = ls_setup.binauralize(ls_setup.loudspeaker_signals(
                                                gains_allrap), fs)
l_allrap2_ir, r_allrap2_ir = ls_setup.binauralize(ls_setup.loudspeaker_signals(
                                                  gains_allrap2), fs)

l_nls_ir, r_nls_ir = ls_setup.binauralize(ls_setup.loudspeaker_signals(
                                          gains_nls), fs)

# %%
# Check B format conversion
B_sig = np.array([1, 1, 0, 0])  # W, X, Y, Z
F_B = sph.b_to_sh(B_sig)
B_sig_re = sph.sh_to_b(F_B)
print("B format to SH conversion:")
utils.test_diff(B_sig, B_sig_re)

# %%
# Some plots
plots.sh_coeffs(F_nm, title="Ambeo: all channels max")
plots.sh_coeffs(F_B, title="b_to_sh: W+X")

# %%
plots.sh_coeffs_subplot([np.array([1, 0, 0, 0]),
                         np.array([0, 1, 0, 0]),
                         np.array([0, 0, 1, 0]),
                         np.array([0, 0, 0, 1])],
                        titles=["0", "1, -1", "1, 0", "1, 1"])

# %%
plots.sh_coeffs(np.sqrt(2) * np.array([1, 0, 0, -1]), 'complex',
                title="sqrt(2) * [1, 0, 0, -1] complex coeffs")

# %%
# Look at simple B format generator
sig2 = np.ones(8)
B = sph.src_to_B(sig2, np.pi / 4, np.pi / 4)
B_nm = sph.b_to_sh(B)
plots.sh_coeffs(B_nm[:, 0], title="Sig 2 B")

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.

# %%
hrir_l = np.fft.irfft(HRTF_l)  # creates 256 samples(t)
hrir_r = np.fft.irfft(HRTF_r)  # creates 256 samples(t)

assert hrir_l.shape == hrir_r.shape

ls_setup = decoder.LoudspeakerSetup(ls_x, ls_y, ls_z)
ls_setup.show()

# Load SH impulse response
ambi_ir = sig.MultiSignal.from_file('../data/IR_Gewandhaus_SH1.wav')
# convert to B-format
ambi_ir = sig.AmbiBSignal.sh_to_b(ambi_ir)

fs = ambi_ir.fs

# - SDM Encoding:
sdm_p = ambi_ir.W
sdm_azi, sdm_colat, _ = sdm.pseudo_intensity(ambi_ir, f_bp=(100, 5000))

# Show first 10000 samples DOA
plots.doa(sdm_azi[:10000], sdm_colat[:10000], fs, p=sdm_p[:10000])


# - SDM Decoding:
# very quick stereo SDM decoding. This is only for testing!
ir_st_l, ir_st_r = sdm.render_stereo_sdm(sdm_p, sdm_azi, sdm_colat)

# Loudspeaker decoding
s_pos = np.array(utils.sph2cart(sdm_azi, sdm_colat)).T
ls_gains = decoder.nearest_loudspeaker(s_pos, ls_setup)
assert len(ls_gains) == len(sdm_p)
ir_ls_l, ir_ls_r = sdm.render_binaural_loudspeaker_sdm(sdm_p, ls_gains,
                                                       ls_setup, fs)

# Render some examples
s_in = sig.MonoSignal.from_file('../data/piano_mono.flac', fs)
s_in.trim(2.6, 6)

normal_limit = 85
    aperture_limit = 90
    opening_limit = 135
    blacklist = None
    ls_setup = IO.load_layout("../data/ls_layouts/Graz.json",
                              listener_position=listener_position)
    ls_setup.pop_triangles(normal_limit, aperture_limit, opening_limit,
                           blacklist)

else:
    raise ValueError


# %% Show setup
ls_setup.show()
plots.hull_normals(ls_setup)

# Test source location
src = np.array([1, 0.5, 2.5])
src_azi, src_colat, _ = utils.cart2sph(*src.tolist())

# %% VBAP
gains_vbap = decoder.vbap(src, ls_setup)


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