Read .sofa HRIR

SOFA (Spatially Oriented Format for Acoustics) files (https://www.sofaconventions.org) are HDF5 containers that hold data like Head-Related Impulse Responses (HRIRs), Binaural Room Impulse Responses (BRIRs) or Directional Impulse Responses (DRIRs, array impulse responses).

Additionally, they hold important meta information and are a great way to archive such recording data. Unfortunately, there currently is no Python API available (You may find at Matlab/Octave and C++ implementation here: https://www.sofaconventions.org/mediawiki/index.php/Software_and_APIs).

Luckily, sofa files can be read using the netCDF4 package, the extracted impulse responses can then be saved into a format that is (currently) easier to work with like numpy's .npy. Note that this is much less optimized for filesize!

Dependencies

This example mainly relies on the netCDF4 package to read in the sofa format. sound_field_analysis is only used for some format definitions.

To be able to install netCDF4, having the packages hdf5 and netcdf installed was necessary. On OSX this worked by:

brew install hdf5
brew install netcdf
pip install netCDF4

In [2]:
import os
from netCDF4 import Dataset
import numpy as np
import sys
sys.path.insert(0, '../')
from sound_field_analysis import io


---------------------------------------------------------------------------
ImportError                               Traceback (most recent call last)
<ipython-input-2-e8f1a70cff48> in <module>()
      1 import os
----> 2 from netCDF4 import Dataset
      3 import numpy as np
      4 import sys
      5 sys.path.insert(0, '../../sound_field_analysis-py')

/anaconda3/lib/python3.6/site-packages/netCDF4/__init__.py in <module>()
      1 # init for netCDF4. package
      2 # Docstring comes from extension module _netCDF4.
----> 3 from ._netCDF4 import *
      4 # Need explicit imports for names beginning with underscores
      5 from ._netCDF4 import __doc__, __pdoc__

ImportError: dlopen(/anaconda3/lib/python3.6/site-packages/netCDF4/_netCDF4.cpython-36m-darwin.so, 2): Symbol not found: _clock_gettime
  Referenced from: /anaconda3/lib/python3.6/site-packages/netCDF4/.dylibs/libcurl.4.dylib (which was built for Mac OS X 10.13)
  Expected in: /usr/lib/libSystem.B.dylib
 in /anaconda3/lib/python3.6/site-packages/netCDF4/.dylibs/libcurl.4.dylib

Load .sofa file

Plenty of sofa files are listed at https://www.sofaconventions.org/mediawiki/index.php/Files In this example, the "mit_kemar_large_pinna.sofa" was used.


In [2]:
sofa_file_name = 'sofa/mit_kemar_large_pinna.sofa'
sofa_file = Dataset(sofa_file_name, 'r', format='NETCDF4')
print('sofa_file: ' + str(sofa_file))


Out[2]:
<class 'netCDF4._netCDF4.Dataset'>
root group (NETCDF4 data model, file format HDF5):
    Conventions: SOFA
    Version: 0.6
    SOFAConventions: SimpleFreeFieldHRIR
    SOFAConventionsVersion: 0.4
    APIName: ARI SOFA API for Matlab/Octave
    APIVersion: 0.4.3
    ApplicationName: ARI converter
    ApplicationVersion: 
    AuthorContact: piotr@majdak.com
    Comment: Converted by Piotr Majdak, Acoustics Research Institute, OeAW
    DataType: FIR
    History: Converted from the MIT database
    License: This HRTF data is provided without any usage restrictions. It is requested that you cite Gardner and Martin (1995) when using this data for research or commercial applications.
    Organization: MIT Media Lab
    References: Gardner, W. G., and Martin, K. D. (1995). "HRTF measurements of a KEMAR," J Acoust Soc Am 97, 3907-3908.
    RoomType: free field
    Origin: http://sound.media.mit.edu/resources/KEMAR.html
    DateCreated: 1999-11-16 20:01:54
    DateModified: 2014-05-22 09:14:59
    Title: HRTF
    DatabaseName: MIT
    ListenerShortName: KEMAR, large pinna
    dimensions(sizes): I(1), C(3), R(2), E(1), N(512), M(710)
    variables(dimensions): float64 ListenerPosition(I,C), float64 ReceiverPosition(R,C,I), float64 SourcePosition(M,C), float64 EmitterPosition(E,C,I), float64 ListenerUp(I,C), float64 ListenerView(I,C), float64 Data.IR(M,R,N), float64 Data.SamplingRate(I), float64 Data.Delay(I,R)
    groups: 

SOFA content

If everything went correctly, you should see a description of the selection .sofa file above.

Generally, HRIR / BRIR sets will have R=2 receiver positions (left and right ear) at I=1 listener position (usually 0,0,0), with impulse responses of M source positions (with E=1 emitter position, usually 0,0,0) of length N.

Other data such as array recordings can have several receiver positions for a single source position and the script below needs to be adjusted accordingly.

As an example, the ear distance is calculated as the difference between the y coordinates of the left and the right ear.


In [3]:
print('SourcePosition: ' + str(sofa_file['SourcePosition']))
print('Data.IR: ' + str(sofa_file['Data.IR']))
print('Ear distance: ' + str(sofa_file['ReceiverPosition'][1, 1] - sofa_file['ReceiverPosition'][0, 1]) + ' m')


SourcePosition: <class 'netCDF4._netCDF4.Variable'>
float64 SourcePosition(M, C)
    Type: spherical
    Units: degree, degree, meter
unlimited dimensions: 
current shape = (710, 3)
filling on, default _FillValue of 9.969209968386869e+36 used

Data.IR: <class 'netCDF4._netCDF4.Variable'>
float64 Data.IR(M, R, N)
unlimited dimensions: 
current shape = (710, 2, 512)
filling on, default _FillValue of 9.969209968386869e+36 used

Ear distance: [ 0.18] m

In [4]:
# extract IRs
HRIRs_l = np.squeeze(sofa_file['Data.IR'][:, 0, :])
HRIRs_r = np.squeeze(sofa_file['Data.IR'][:, 1, :])
Az, El, R = np.squeeze(np.hsplit(sofa_file['SourcePosition'][:], 3))
fs = sofa_file['Data.SamplingRate'][:][0]

In [5]:
sofa_file.close()

Save as npy file

To now write the data into the handy binary .npy format, we could simply call np.save().

Specifically for the sound_field_analysis toolbox, we convert the angles into radiants (and elevation to colatitude), create an io.ArraySignal and finally save that.


In [6]:
azimuth = Az / 180 * np.pi
colatitude = np.pi / 2 - El / 180 * np.pi

hrir_full_l = io.ArraySignal(io.TimeSignal(HRIRs_l, fs), io.SphericalGrid(azimuth, colatitude, R))
hrir_full_r = io.ArraySignal(io.TimeSignal(HRIRs_r, fs), io.SphericalGrid(azimuth, colatitude, R))

In [7]:
np.save(os.path.splitext(os.path.basename(sofa_file_name))[0] + '_L', hrir_full_l)
np.save(os.path.splitext(os.path.basename(sofa_file_name))[0] + '_R', hrir_full_r)