How to use the spaudiopy.process function in spaudiopy

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_resample_hrirs(test_jobs):
    hrirs = spa.IO.load_hrirs(fs=44100, filename='dummy')
    hrir_l_rsmp_r, hrir_r_rsmp_r, _ = spa.process.resample_hrirs(hrirs.left,
                                                                 hrirs.right,
                                                                 44100, 48000,
                                                                 jobs_count=1)
    hrir_l_rsmp_t, hrir_r_rsmp_t, _ = spa.process.resample_hrirs(hrirs.left,
                                                                 hrirs.right,
                                                                 44100, 48000,
                                                                 jobs_count=
                                                                 test_jobs)
    assert_allclose([hrir_l_rsmp_t, hrir_r_rsmp_t],
                    [hrir_l_rsmp_r, hrir_r_rsmp_r])

def test_resample_hrirs(test_jobs):
    hrirs = spa.IO.load_hrirs(fs=44100, filename='dummy')
    hrir_l_rsmp_r, hrir_r_rsmp_r, _ = spa.process.resample_hrirs(hrirs.left,
                                                                 hrirs.right,
                                                                 44100, 48000,
                                                                 jobs_count=1)
    hrir_l_rsmp_t, hrir_r_rsmp_t, _ = spa.process.resample_hrirs(hrirs.left,
                                                                 hrirs.right,
                                                                 44100, 48000,
                                                                 jobs_count=
                                                                 test_jobs)
    assert_allclose([hrir_l_rsmp_t, hrir_r_rsmp_t],
                    [hrir_l_rsmp_r, hrir_r_rsmp_r])

# default grid:
    if (grid_azi is None) and (grid_colat is None):
        grid_azi, grid_colat, _ = grids.gauss(35)  # grid positions

    # %% Inverse SHT
    HRTF_l = sph.inverse_sht(SH_l, grid_azi, grid_colat, 'complex')
    HRTF_r = sph.inverse_sht(SH_r, grid_azi, grid_colat, 'complex')
    assert HRTF_l.shape == HRTF_r.shape
    # %%
    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

    # %% Resample
    fs_target = 48000
    hrir_l_48k, hrir_r_48k, _ = process.resample_hrirs(hrir_l, hrir_r,
                                                       SamplingRate,
                                                       fs_target)
    fs_target = 96000
    hrir_l_96k, hrir_r_96k, _ = process.resample_hrirs(hrir_l, hrir_r,
                                                       SamplingRate,
                                                       fs_target)

    savemat(os.path.join(current_file_dir, '../data/HRTF_default_44100.mat'),
            {'hrir_l': hrir_l,
             'hrir_r': hrir_r,
             'azi': grid_azi, 'colat': grid_colat,
             'fs': 44100})
    savemat(os.path.join(current_file_dir, '../data/HRTF_default_48000.mat'),
            {'hrir_l': hrir_l_48k,
             'hrir_r': hrir_r_48k,
             'azi': grid_azi, 'colat': grid_colat,

# %% Inverse SHT
    HRTF_l = sph.inverse_sht(SH_l, grid_azi, grid_colat, 'complex')
    HRTF_r = sph.inverse_sht(SH_r, grid_azi, grid_colat, 'complex')
    assert HRTF_l.shape == HRTF_r.shape
    # %%
    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

    # %% Resample
    fs_target = 48000
    hrir_l_48k, hrir_r_48k, _ = process.resample_hrirs(hrir_l, hrir_r,
                                                       SamplingRate,
                                                       fs_target)
    fs_target = 96000
    hrir_l_96k, hrir_r_96k, _ = process.resample_hrirs(hrir_l, hrir_r,
                                                       SamplingRate,
                                                       fs_target)

    savemat(os.path.join(current_file_dir, '../data/HRTF_default_44100.mat'),
            {'hrir_l': hrir_l,
             'hrir_r': hrir_r,
             'azi': grid_azi, 'colat': grid_colat,
             'fs': 44100})
    savemat(os.path.join(current_file_dir, '../data/HRTF_default_48000.mat'),
            {'hrir_l': hrir_l_48k,
             'hrir_r': hrir_r_48k,
             'azi': grid_azi, 'colat': grid_colat,
             'fs': 48000})
    savemat(os.path.join(current_file_dir, '../data/HRTF_default_96000.mat'),
            {'hrir_l': hrir_l_96k,
             'hrir_r': hrir_r_96k,

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)

axis=1, arr=hrir_l)
hrir_r_hp = np.apply_along_axis(lambda m:
                                scysignal.convolve(m, h_headphone),
                                axis=1, arr=hrir_r)

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,