Binaural Rendering of Array Data

This notebook imports a local version of sound_field_analysis-py to prototype quickly inside the package. It produces a binaural rendering of a spherical microphone array measurement.

To run it, you need to download the necessary impulse responses either as SOFA or MIRO files and reflect your choice in the variable is_load_sofa. After downloading from the provided links the files need to be placed in the specified directory.

In [1]:
# this automatically reloads the sfa-py package if anything is changed there
%load_ext autoreload
%autoreload 2

# Requirements
from IPython import display
import numpy as np

# Load the local sound_field_analysis-py package
import sys
sys.path.insert(0, '../')
from sound_field_analysis import io, gen, process, plot, sph, utils

Settings

In [2]:
is_load_sofa    = True  # whether SOFA source data should be loaded, MIRO source data otherwise
sh_max_order    = 8  # maximal utilized spherical harmonics rendering order
rf_nfft         = 4096  # target length of radial filter in samples
rf_amp_max_db   = 0  # soft limiting level of radial filter in dB
is_apply_rfi    = True  # whether radial filter improvement should be applied
azim_offset_deg = -37  # azimuth head roation offset in degrees (to make the listener face the sound source)

Load data

SOFA (recommended):

Variable File Description Author Source
DRIR examples/ data/ DRIR_CR1_VSA_110RS_L.sofa WDR Control Room 1, left loudspeaker
spherical Lebedev grid, 110 nodes
up to SH order 29
TH Köln WDR Impulse Response Compilation (SOFA)
HRIR examples/ data/ HRIR_L2702.sofa Neumann K100, artificial head
spherical Lebedev grid, 2702 nodes
up to SH order 8
TH Köln sofaconventions.org

MIRO:

Variable File Description Author Source
DRIR examples/ data/ CR1_VSA_110RS_L.mat ** WDR Control Room 1, left loudspeaker
spherical Lebedev grid, 110 nodes
up to SH order 29
TH Köln WDR CR1_VSA
HRIR_L,
HRIR_R
examples/ data/ HRIR_L2702.mat ** Neumann K100, artificial head
spherical Lebedev grid, 2702 nodes
up to SH order 8
TH Köln KU100 HRIR Dataset
examples/ miro.m MIRO Matlab Class Definition TH Köln miro Class Definition (Required for KU/VSA Datasets)

** These files contain data in the Matlab MIRO format, which cannot be accessed in Python. To read them, you need to download the specified MIRO Matlab Class Definition file. You will have to convert the data to a Matlab structure first by executing the prepared Matlab script examples/ miro_to_struct.m .

In [3]:
if is_load_sofa:
    # load impulse responses from SOFA file
    DRIR   = io.read_SOFA_file('data/DRIR_CR1_VSA_110RS_L.sofa')
    HRIR   = io.read_SOFA_file('data/HRIR_L2702.sofa')
    NFFT   = HRIR.l.signal.shape[-1]
else:
    # load impulse responses from MIRO struct
    DRIR   = io.read_miro_struct('data/CR1_VSA_110RS_L_struct.mat')
    HRIR_L = io.read_miro_struct('data/HRIR_L2702_struct.mat', channel='irChOne')
    HRIR_R = io.read_miro_struct('data/HRIR_L2702_struct.mat', channel='irChTwo')
    NFFT   = HRIR_L.signal.signal.shape[-1]

# automatically calculate target processing length
# by summing impulse response lengths of DRIR, HRIR and radial filters
NFFT += DRIR.signal.signal.shape[1] + rf_nfft
open SOFA file "data/DRIR_CR1_VSA_110RS_L.sofa"
 --> samplerate: 48000 Hz, receivers: 110, emitters: 1, measurements: 1, samples: 17000, format: float64
 --> convention: SingleRoomDRIR, version: 0.3
 --> listener: Earthworks M30 Omni - S/N3551E (Rigid Sphere)
 --> author: Benjamin Bernsch�tz, Philipp Stade, Max R�hl
('degree, degree, metre', 'spherical')

open SOFA file "data/HRIR_L2702.sofa"
 --> samplerate: 48000 Hz, receivers: 2, emitters: 1, measurements: 2702, samples: 128, format: float64
 --> convention: SimpleFreeFieldHRIR, version: 1.0
 --> listener: Neumann KU100
 --> author: Benjamin Bernsch�tz
('degree, degree, metre', 'spherical')

Processing

In [4]:
# FFT and spatFT
if is_load_sofa:
    # transform SOFA data
    HRTF_L = process.FFT(HRIR.l.signal, fs=HRIR.l.fs, NFFT=NFFT, calculate_freqs=False)
    HRTF_R = process.FFT(HRIR.r.signal, fs=HRIR.r.fs, NFFT=NFFT, calculate_freqs=False)
    Hnm_L = process.spatFT(HRTF_L, HRIR.grid, order_max=sh_max_order)
    Hnm_R = process.spatFT(HRTF_R, HRIR.grid, order_max=sh_max_order)
else:
    # transform MIRO data
    HRTF_L = process.FFT(HRIR_L.signal.signal, fs=HRIR_L.signal.fs, NFFT=NFFT, calculate_freqs=False)
    HRTF_R = process.FFT(HRIR_R.signal.signal, fs=HRIR_R.signal.fs, NFFT=NFFT, calculate_freqs=False)
    Hnm_L = process.spatFT(HRTF_L, HRIR_L.grid, order_max=sh_max_order)
    Hnm_R = process.spatFT(HRTF_R, HRIR_R.grid, order_max=sh_max_order)

DRTF = process.FFT(DRIR.signal.signal[:,:int(NFFT)], fs=DRIR.signal.fs, NFFT=NFFT, calculate_freqs=False)
Pnm = process.spatFT(DRTF, DRIR.grid, order_max=sh_max_order)

# compute radial filter
dn = gen.radial_filter_fullspec(max_order=sh_max_order, NFFT=rf_nfft, fs=DRIR.signal.fs,
                                array_configuration=DRIR.configuration, amp_maxdB=rf_amp_max_db)

# show frequency domain preview of radial filter


if is_apply_rfi:
    # improve radial filter (remove DC offset and make casual)
    dn, _, dn_delay_samples = process.rfi(dn, kernelSize=rf_nfft)
else:
    # make radial filter causal
    dn_delay_samples = rf_nfft / 2
    dn *= gen.delay_fd(target_length_fd=dn.shape[1], delay_samples=dn_delay_samples)

# show frequency domain preview of radial filter
dn_legend = list(f'n = {n}' for n in range(sh_max_order + 1))
plot.plot2D(dn, fs=DRIR.signal.fs, viz_type='LogFFT', title='radial filter', line_names=dn_legend)

# adjust radial filter length
dn = utils.zero_pad_fd(dn, target_length_td=NFFT)

# # show time domain preview of radial filter
# plot.plot2D(np.fft.irfft(dn), fs=DRIR.signal.fs, viz_type='Time', title='radial filter', line_names=dn_legend)

# adjust radial filter size according to SH order as grades
dn = np.repeat(dn, np.arange(sh_max_order * 2 + 1, step=2) + 1, axis=0)