How to use spaudiopy - 10 common examples

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

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()

Compensated loudspeaker signals.

    References
    ----------
    Tervo, S., et. al. (2015).
    Spatial Analysis and Synthesis of Car Audio System and Car Cabin Acoustics
    with a Compact Microphone Array. Journal of the Audio Engineering Society.

    """
    ls_distance = ls_setup.d  # ls distance
    a = ls_setup.a  # distance attenuation exponent

    CHECK_SANITY = False

    # prepare filterbank
    filter_gs, ff = pcs.frac_octave_filterbank(n=1, N_out=2**16,
                                               fs=fs, f_low=62.5, f_high=16000,
                                               mode='amplitude')

    # band dependent block size
    band_blocksizes = np.zeros(ff.shape[0])
    # proposed by Tervo
    band_blocksizes[1:] = np.round(7 / ff[1:, 0] * fs)
    band_blocksizes[0] = np.round(7 / ff[0, 1] * fs)
    # make sure they are even
    band_blocksizes = (np.ceil(band_blocksizes / 2) * 2).astype(int)

    padsize = band_blocksizes.max()

    ntaps = padsize // 2 - 1
    assert(ntaps % 2), "N does not produce uneven number of filter taps."
    irs = np.zeros([filter_gs.shape[0], ntaps])

mag_diff = np.abs(H_p) / \
                            np.clip(sdm_mag_coherent, 10e-10, None)
            elif band_idx == 1:
                mag_diff = np.abs(H_p) / \
                           (0.5 * np.clip(sdm_mag_coherent, 10e-10, None) +
                            0.5 * np.clip(sdm_mag_incoherent, 10e-10, None))
            elif band_idx == 2:
                mag_diff = np.abs(H_p) / \
                           (0.25 * np.clip(sdm_mag_coherent, 10e-10, None) +
                            0.75 * np.clip(sdm_mag_incoherent, 10e-10, None))
            else:
                mag_diff = np.abs(H_p) / np.clip(sdm_mag_incoherent, 10e-10,
                                                 None)
            # soft clip gain
            if soft_clip:
                mag_diff = pcs.gain_clipping(mag_diff, 1)

            # apply to ls input
            Y = H_sdm * mag_diff[np.newaxis, :]

            # inverse STFT
            X = np.real(np.fft.ifft(Y, axis=1))
            # Zero Phase
            assert(np.mod(X.shape[1], 2))
            # delay
            zp_delay = X.shape[1] // 2
            X = np.roll(X, zp_delay, axis=1)

            # overlap add
            ls_sigs_band[:, padsize + start_idx - zp_delay:
                            padsize + start_idx - zp_delay + nfft,
                         band_idx] += X

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_allrap2(test_jobs):
    vecs = spa.grids.load_t_design(degree=5)
    hull = spa.decoder.LoudspeakerSetup(*vecs.T)
    hull.ambisonics_setup(update_hull=False)
    src = np.random.randn(1000, 3)
    gains_r = spa.decoder.allrap2(src, hull, jobs_count=1)
    gains_t = spa.decoder.allrap2(src, hull, jobs_count=test_jobs)
    assert_allclose(gains_t, gains_r)

def test_vbap(test_jobs):
    vecs = spa.grids.load_t_design(degree=5)
    hull = spa.decoder.LoudspeakerSetup(*vecs.T)
    src = np.random.randn(1000, 3)
    gains_r = spa.decoder.vbap(src, hull, jobs_count=1)
    gains_t = spa.decoder.vbap(src, hull, jobs_count=test_jobs)
    assert_allclose(gains_t, gains_r)

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

def test_allrap(test_jobs):
    vecs = spa.grids.load_t_design(degree=5)
    hull = spa.decoder.LoudspeakerSetup(*vecs.T)
    hull.ambisonics_setup(update_hull=False)
    src = np.random.randn(1000, 3)
    gains_r = spa.decoder.allrap(src, hull, jobs_count=1)
    gains_t = spa.decoder.allrap(src, hull, jobs_count=test_jobs)
    assert_allclose(gains_t, gains_r)