diff --git a/docs/pyttadoc/index.rst b/docs/pyttadoc/index.rst index ca3da8d..24ebca2 100644 --- a/docs/pyttadoc/index.rst +++ b/docs/pyttadoc/index.rst @@ -57,6 +57,9 @@ Modules: * pytta.iso3741: calculations according to the standard; + * pytta.sofa: + read and write SOFA files (AES69-2022) with impulse responses; + .. toctree:: classes/index utils/index @@ -64,4 +67,5 @@ Modules: functions generate rooms + sofa properties diff --git a/docs/pyttadoc/sofa.rst b/docs/pyttadoc/sofa.rst new file mode 100644 index 0000000..bc4f8e5 --- /dev/null +++ b/docs/pyttadoc/sofa.rst @@ -0,0 +1,7 @@ +.. currentmodule:: pytta + +SOFA +==== + +.. automodule:: pytta.sofa + :members: diff --git a/pytta/__init__.py b/pytta/__init__.py index d8740aa..a9bf7eb 100644 --- a/pytta/__init__.py +++ b/pytta/__init__.py @@ -72,6 +72,9 @@ * pytta.iso3741: calculations according to the standard; + * pytta.sofa: + read and write SOFA files (AES69-2022) with impulse responses; + You can find usage examples in the toolbox's examples folder. @@ -108,6 +111,7 @@ from . import utils from . import rooms from . import iso3741 +from . import sofa from .functions import read_wav, write_wav, merge, split, \ list_devices, print_devices, get_device_from_user,\ @@ -130,6 +134,7 @@ 'generate', 'utils', 'iso3741', + 'sofa', 'rooms', # Functions diff --git a/pytta/sofa.py b/pytta/sofa.py new file mode 100644 index 0000000..4a463b6 --- /dev/null +++ b/pytta/sofa.py @@ -0,0 +1,249 @@ +""" +PyTTa SOFA input and output: +---------------------------- + +Read and write impulse responses in the Spatially Oriented Format for +Acoustics (SOFA), standardized as AES69-2022 (SOFA version 2.1). + +A SOFA file is a netCDF-4 file, which is an HDF5 file whose dimensions are +HDF5 dimension scales. This module uses h5py only and writes the netCDF +dimension metadata itself, so the files open in netCDF tools (MATLAB ncread +and ncinfo, the SOFA Toolbox, netCDF4-python). + +Supported conventions for writing (both at convention version 1.0): + + * SimpleFreeFieldHRIR: head-related impulse responses, source positions + in spherical coordinates (azimuth in degree, elevation in degree, + distance in metre); + * GeneralFIR: impulse responses for general purposes, source positions + in cartesian coordinates (metre). + +Any convention that stores Data.IR (GeneralFIR, SimpleFreeFieldHRIR, +SingleRoomSRIR, ...) can be read. + +Data mapping: + + SOFA stores Data.IR with dimensions M x R x N (measurements, receivers, + samples). Here each measurement m is one SignalObj whose R channels are + the receivers, so a SOFA file maps to a list of M SignalObj, in the + order of SourcePosition. SignalObj keeps time data in float32; the + float64 Data.IR is also returned by read as 'ir'. + +Available functions: + + >>> pytta.sofa.read(path) + >>> pytta.sofa.write(path, signal, source_positions) + +For further information, check the function specific documentation. + +""" + +from datetime import datetime + +import h5py +import numpy as np + +from .classes import SignalObj + +_DIMENSION_NAME = 'This is a netCDF dimension but not a netCDF variable.%10d' + +_CONVENTIONS = { + 'SimpleFreeFieldHRIR': { + 'globals': {'DatabaseName': '', 'ListenerShortName': ''}, + 'source': ('spherical', 'degree, degree, metre'), + 'ears': True, + }, + 'GeneralFIR': { + 'globals': {}, + 'source': ('cartesian', 'metre'), + 'ears': False, + }, +} + +_POSITIONS = {'source_position': 'SourcePosition', + 'listener_position': 'ListenerPosition', + 'receiver_position': 'ReceiverPosition', + 'emitter_position': 'EmitterPosition'} + + +def _text(value): + """Decode an HDF5 attribute value into str when it is text.""" + if isinstance(value, h5py.Empty): + return '' + if isinstance(value, bytes): + return value.decode('utf-8') + if isinstance(value, np.ndarray) and value.dtype.kind in 'SO': + return ''.join(_text(v) for v in value.ravel()) + return value + + +def read(path): + """ + Read a SOFA file with impulse responses (Data.IR). + + Parameters + ---------- + path : str + Path to the .sofa file. + + Returns + ------- + dict + * 'signal': list of M SignalObj, one per measurement, with the R + receivers as channels (float32, as SignalObj stores it; SignalObj + transposes arrays with more channels than samples, so use 'ir' + when R > N); + * 'ir': Data.IR as a float64 array of shape (M, R, N); + * 'fs': sampling rate in Hz (float); + * 'delay': Data.Delay in samples, shape (I, R) or (M, R). The delay + is returned and not applied to the signals; + * 'source_position', 'listener_position', 'receiver_position', + 'emitter_position': dicts with 'values' (the array as stored, + e.g. (M, 3) for SourcePosition and (R, 3, I) for + ReceiverPosition), 'Type' and 'Units'; + * 'attributes': dict of the global attributes, without the + 'GLOBAL:' prefix. + + Raises + ------ + ValueError + If the file has no Data.IR, or a Data.SamplingRate that changes + between measurements, or units other than hertz. + + """ + with h5py.File(path, 'r') as f: + if 'Data.IR' not in f: + raise ValueError("pytta.sofa.read only reads files with Data.IR.") + ir = np.asarray(f['Data.IR'][()], dtype='float64') + rates = np.unique(f['Data.SamplingRate'][()]) + units = _text(f['Data.SamplingRate'].attrs.get('Units', 'hertz')) + if len(rates) != 1 or units.lower() != 'hertz': + raise ValueError("Data.SamplingRate must be one value in hertz.") + fs = float(rates[0]) + delay = np.asarray(f['Data.Delay'][()]) if 'Data.Delay' in f else None + result = {'ir': ir, 'fs': fs, 'delay': delay, + 'attributes': {k: _text(v) for k, v in f.attrs.items() + if not k.startswith('_')}} + for key, name in _POSITIONS.items(): + if name in f: + result[key] = {'values': f[name][()], + 'Type': _text(f[name].attrs.get('Type', '')), + 'Units': _text(f[name].attrs.get('Units', ''))} + signals = [] + for m in range(ir.shape[0]): + signals.append(SignalObj(signalArray=ir[m].T, domain='time', + samplingRate=fs)) + result['signal'] = signals + return result + + +def write(path, signal, source_positions, convention='SimpleFreeFieldHRIR', + **attributes): + """ + Write impulse responses to a SOFA file (AES69-2022, SOFA 2.1). + + All mandatory global attributes and variables of the convention are + written. Positions other than SourcePosition get the convention defaults: + listener at the origin looking along +x with up +z, one emitter at the + origin of the source, receivers at the origin of the listener (for + SimpleFreeFieldHRIR with two receivers, the ears at +-0.09 m on y), and + Data.Delay zero. + + Parameters + ---------- + path : str + Path of the .sofa file to create (overwritten if it exists). + signal : SignalObj or list of SignalObj + One SignalObj per measurement, all with the same sampling rate, + number of samples N and number of channels R (the receivers). + source_positions : array_like, shape (M, 3) + Source position of each measurement. Spherical (azimuth deg, + elevation deg, distance m) for SimpleFreeFieldHRIR, cartesian (m) + for GeneralFIR. + convention : {'SimpleFreeFieldHRIR', 'GeneralFIR'} + SOFA convention, version 1.0. + **attributes : str + Global attributes, without the 'GLOBAL:' prefix (e.g. Title, + AuthorContact, Organization, License, DatabaseName). They replace + the defaults; mandatory ones not given are written empty. + + Raises + ------ + ValueError + For an unknown convention, measurements that differ in sampling + rate, length or number of channels, or a number of source positions + different from the number of measurements. + + """ + if convention not in _CONVENTIONS: + raise ValueError(f"convention must be one of {list(_CONVENTIONS)}.") + conv = _CONVENTIONS[convention] + signals = [signal] if isinstance(signal, SignalObj) else list(signal) + src = np.atleast_2d(np.asarray(source_positions, dtype='float64')) + if src.shape != (len(signals), 3): + raise ValueError("source_positions must have shape (M, 3), with M " + "the number of measurements.") + if len({(s.samplingRate, s.numSamples, s.numChannels) + for s in signals}) != 1: + raise ValueError("All measurements must have the same sampling rate," + " number of samples and number of channels.") + ir = np.stack([np.asarray(s.timeSignal, dtype='float64').T + for s in signals]) + M, R, N = ir.shape + receivers = np.zeros((R, 3, 1)) + if conv['ears'] and R == 2: + receivers[:, 1, 0] = [0.09, -0.09] + + from . import __version__ + now = datetime.now().strftime('%Y-%m-%d %H:%M:%S') + glob = {'Conventions': 'SOFA', 'Version': '2.1', + 'SOFAConventions': convention, 'SOFAConventionsVersion': '1.0', + 'APIName': 'pytta', 'APIVersion': __version__, + 'AuthorContact': '', 'DataType': 'FIR', + 'License': 'No license provided, ask the author for permission', + 'Organization': '', 'RoomType': 'free field', + 'DateCreated': now, 'DateModified': now, 'Title': ''} + glob.update(conv['globals']) + glob.update(attributes) + + cart = ('cartesian', 'metre') + variables = [ # name, dimensions, values, (Type, Units) + ('ListenerPosition', 'IC', [[0., 0., 0.]], cart), + ('ReceiverPosition', 'RCI', receivers, cart), + ('SourcePosition', 'MC', src, conv['source']), + ('EmitterPosition', 'ECI', np.zeros((1, 3, 1)), cart), + ('ListenerUp', 'IC', [[0., 0., 1.]], None), + ('ListenerView', 'IC', [[1., 0., 0.]], cart), + ('Data.IR', 'MRN', ir, None), + ('Data.SamplingRate', 'I', [float(signals[0].samplingRate)], None), + ('Data.Delay', 'IR', np.zeros((1, R)), None), + ] + sizes = {'C': 3, 'I': 1, 'M': M, 'R': R, 'E': 1, 'N': N} + + with h5py.File(path, 'w', track_order=True) as f: + for key, value in glob.items(): + _set_text(f.attrs, key, value) + dimid = {} + for k, (name, size) in enumerate(sizes.items()): + # netCDF-4 dimension without a coordinate variable. + dim = f.create_dataset(name, shape=(size,), dtype='>f4') + dim.make_scale(_DIMENSION_NAME % size) + dim.attrs['_Netcdf4Dimid'] = np.int32(k) + dimid[name] = k + for name, dims, values, typeunits in variables: + var = f.create_dataset(name, data=np.asarray(values, 'float64'), + track_order=True) + for axis, d in enumerate(dims): + var.dims[axis].attach_scale(f[d]) + var.attrs['_Netcdf4Coordinates'] = np.array( + [dimid[d] for d in dims], dtype='int32') + if typeunits is not None: + _set_text(var.attrs, 'Type', typeunits[0]) + _set_text(var.attrs, 'Units', typeunits[1]) + _set_text(f['Data.SamplingRate'].attrs, 'Units', 'hertz') + + +def _set_text(attrs, key, value): + """Store text as a netCDF NC_CHAR attribute (fixed-length string).""" + raw = str(value).encode('utf-8') + attrs[key] = np.bytes_(raw) if raw else h5py.Empty(np.dtype('S1')) diff --git a/tests/golden/make_sofa.m b/tests/golden/make_sofa.m new file mode 100644 index 0000000..ecebb8e --- /dev/null +++ b/tests/golden/make_sofa.m @@ -0,0 +1,83 @@ +% Generate the SOFA golden files used by tests/test_sofa.py. +% +% sofa_golden.sofa is a tiny SimpleFreeFieldHRIR 1.0 file (AES69-2022) with +% M = 3 measurements, R = 2 receivers and N = 8 samples, written with the +% netCDF library built into MATLAB (nccreate, ncwrite, ncwriteatt), so it +% does not depend on pytta. sofa_golden.mat holds what MATLAB ncread and +% ncreadatt return for that file, permuted to the C order used by netCDF +% and by pytta.sofa.read (Data.IR is M x R x N). +% +% Run from tests/golden: +% matlab -batch "make_sofa" + +f = 'sofa_golden.sofa'; +if exist(f, 'file') + delete(f); +end + +M = 3; R = 2; N = 8; +[n, r, m] = ndgrid(0:N-1, 0:R-1, 0:M-1); +ir = sin(1 + n + 10*r + 100*m) ./ (1 + n); % N x R x M in MATLAB order +fs = 44100; +src = [0 0 1.2; 90 0 1.2; 180 45 1.2]; % M x C, spherical +rec = [0 0.09 0; 0 -0.09 0]; % R x C, cartesian + +% MATLAB lists dimensions fastest first, the reverse of netCDF C order: +% {'N','R','M'} gives the netCDF variable Data.IR(M, R, N). +fmt = {'Format', 'netcdf4', 'Datatype', 'double'}; +nccreate(f, 'ListenerPosition', 'Dimensions', {'C', 3, 'I', 1}, fmt{:}); +nccreate(f, 'ReceiverPosition', 'Dimensions', {'I', 1, 'C', 3, 'R', R}, fmt{:}); +nccreate(f, 'SourcePosition', 'Dimensions', {'C', 3, 'M', M}, fmt{:}); +nccreate(f, 'EmitterPosition', 'Dimensions', {'I', 1, 'C', 3, 'E', 1}, fmt{:}); +nccreate(f, 'ListenerUp', 'Dimensions', {'C', 3, 'I', 1}, fmt{:}); +nccreate(f, 'ListenerView', 'Dimensions', {'C', 3, 'I', 1}, fmt{:}); +nccreate(f, 'Data.IR', 'Dimensions', {'N', N, 'R', R, 'M', M}, fmt{:}); +nccreate(f, 'Data.SamplingRate', 'Dimensions', {'I', 1}, fmt{:}); +nccreate(f, 'Data.Delay', 'Dimensions', {'R', R, 'I', 1}, fmt{:}); + +ncwrite(f, 'ListenerPosition', [0; 0; 0]); +ncwrite(f, 'ReceiverPosition', reshape(rec.', [1 3 R])); +ncwrite(f, 'SourcePosition', src.'); +ncwrite(f, 'EmitterPosition', [0 0 0]); +ncwrite(f, 'ListenerUp', [0; 0; 1]); +ncwrite(f, 'ListenerView', [1; 0; 0]); +ncwrite(f, 'Data.IR', ir); +ncwrite(f, 'Data.SamplingRate', fs); +ncwrite(f, 'Data.Delay', zeros(R, 1)); + +glob = {'Conventions', 'SOFA'; 'Version', '2.1'; ... + 'SOFAConventions', 'SimpleFreeFieldHRIR'; ... + 'SOFAConventionsVersion', '1.0'; 'APIName', 'MATLAB netCDF'; ... + 'APIVersion', version('-release'); 'DataType', 'FIR'; ... + 'AuthorContact', 'pytta developers'; 'Organization', 'PyTTa'; ... + 'License', 'MIT'; 'RoomType', 'free field'; ... + 'DateCreated', '2026-10-31 12:00:00'; ... + 'DateModified', '2026-10-31 12:00:00'; ... + 'Title', 'pytta.sofa golden'; 'DatabaseName', 'pytta tests'; ... + 'ListenerShortName', 'none'}; +for k = 1:size(glob, 1) + ncwriteatt(f, '/', glob{k, 1}, glob{k, 2}); +end +pos = {'ListenerPosition', 'cartesian', 'metre'; ... + 'ReceiverPosition', 'cartesian', 'metre'; ... + 'SourcePosition', 'spherical', 'degree, degree, metre'; ... + 'EmitterPosition', 'cartesian', 'metre'; ... + 'ListenerView', 'cartesian', 'metre'}; +for k = 1:size(pos, 1) + ncwriteatt(f, pos{k, 1}, 'Type', pos{k, 2}); + ncwriteatt(f, pos{k, 1}, 'Units', pos{k, 3}); +end +ncwriteatt(f, 'Data.SamplingRate', 'Units', 'hertz'); + +% What MATLAB reads back, in netCDF C order. +IR = permute(ncread(f, 'Data.IR'), [3 2 1]); % M x R x N +SamplingRate = ncread(f, 'Data.SamplingRate'); +Delay = ncread(f, 'Data.Delay').'; % I x R +SourcePosition = ncread(f, 'SourcePosition').'; % M x C +ReceiverPosition = permute(ncread(f, 'ReceiverPosition'), [3 2 1]); % R x C x I +SourcePositionType = ncreadatt(f, 'SourcePosition', 'Type'); +SourcePositionUnits = ncreadatt(f, 'SourcePosition', 'Units'); +Title = ncreadatt(f, '/', 'Title'); +save('sofa_golden.mat', 'IR', 'SamplingRate', 'Delay', 'SourcePosition', ... + 'ReceiverPosition', 'SourcePositionType', 'SourcePositionUnits', ... + 'Title', '-v7'); diff --git a/tests/golden/sofa_golden.mat b/tests/golden/sofa_golden.mat new file mode 100644 index 0000000..d25b391 Binary files /dev/null and b/tests/golden/sofa_golden.mat differ diff --git a/tests/golden/sofa_golden.sofa b/tests/golden/sofa_golden.sofa new file mode 100644 index 0000000..536d4e8 Binary files /dev/null and b/tests/golden/sofa_golden.sofa differ diff --git a/tests/test_sofa.py b/tests/test_sofa.py new file mode 100644 index 0000000..eddaea4 --- /dev/null +++ b/tests/test_sofa.py @@ -0,0 +1,127 @@ +"""Tests for pytta.sofa (AES69-2022, SOFA 2.1).""" + +import os +import tempfile +import unittest + +import h5py +import numpy as np +from scipy.io import loadmat + +import pytta +from pytta import sofa + +GOLDEN = os.path.join(os.path.dirname(__file__), 'golden') + +# Items flagged 'm' or 'rm' in SimpleFreeFieldHRIR_1.0.csv of the SOFA +# conventions (sofaconventions.org, SOFAtoolbox/conventions). +HRIR_GLOBALS = ['Conventions', 'Version', 'SOFAConventions', + 'SOFAConventionsVersion', 'APIName', 'APIVersion', + 'AuthorContact', 'DataType', 'License', 'Organization', + 'RoomType', 'DateCreated', 'DateModified', 'Title', + 'DatabaseName', 'ListenerShortName'] +HRIR_VARIABLES = {'ListenerPosition': ('I', 'C'), + 'ReceiverPosition': ('R', 'C', 'I'), + 'SourcePosition': ('M', 'C'), + 'EmitterPosition': ('E', 'C', 'I'), + 'ListenerUp': ('I', 'C'), + 'ListenerView': ('I', 'C'), + 'Data.IR': ('M', 'R', 'N'), + 'Data.SamplingRate': ('I',), + 'Data.Delay': ('I', 'R')} +HRIR_VARIABLE_ATTRS = ['ListenerPosition:Type', 'ListenerPosition:Units', + 'ReceiverPosition:Type', 'ReceiverPosition:Units', + 'SourcePosition:Type', 'SourcePosition:Units', + 'EmitterPosition:Type', 'EmitterPosition:Units', + 'ListenerView:Type', 'ListenerView:Units', + 'Data.SamplingRate:Units'] + + +class TestSofa(unittest.TestCase): + + def setUp(self): + self.tmp = tempfile.TemporaryDirectory() + + def tearDown(self): + self.tmp.cleanup() + + def test_round_trip_keeps_data_fs_and_positions(self): + """GeneralFIR with 4 receivers: write then read is exact.""" + rng = np.random.default_rng(0) + ir = rng.standard_normal((3, 4, 64)).astype('float32') + signals = [] + for m in range(3): + signals.append(pytta.SignalObj(signalArray=ir[m].T, + domain='time', + samplingRate=32000)) + src = np.array([[1.0, 0.5, 0.0], [0.0, 2.25, 1.5], [-1.0, 0.0, 0.3]]) + path = os.path.join(self.tmp.name, 'rt.sofa') + sofa.write(path, signals, src, convention='GeneralFIR', + Title='round trip') + out = sofa.read(path) + d_ir = np.max(np.abs(out['ir'] - ir)) + d_sig = max(np.max(np.abs(s.timeSignal - ir[m].T)) + for m, s in enumerate(out['signal'])) + d_src = np.max(np.abs(out['source_position']['values'] - src)) + print(f'\nround trip: max |dIR| = {d_ir}, max |dSignal| = {d_sig}, ' + f'max |dSourcePosition| = {d_src}, fs = {out["fs"]}') + self.assertEqual(d_ir, 0.0) + self.assertEqual(d_sig, 0.0) + self.assertEqual(d_src, 0.0) + self.assertEqual(out['fs'], 32000) + self.assertEqual(out['source_position']['Type'], 'cartesian') + self.assertEqual(out['attributes']['Title'], 'round trip') + self.assertEqual(out['attributes']['SOFAConventions'], 'GeneralFIR') + + def test_hrir_has_every_mandatory_item(self): + """The written file has all 'm'/'rm' items of SimpleFreeFieldHRIR.""" + signals = [] + for az in (0, 90): + signals.append(pytta.SignalObj(signalArray=np.ones((16, 2)), + domain='time', + samplingRate=48000)) + path = os.path.join(self.tmp.name, 'hrir.sofa') + sofa.write(path, signals, [[0, 0, 1.2], [90, 0, 1.2]]) + sizes = {'M': 2, 'R': 2, 'N': 16, 'C': 3, 'I': 1, 'E': 1} + with h5py.File(path, 'r') as f: + missing = [a for a in HRIR_GLOBALS if a not in f.attrs] + missing += [v for v in HRIR_VARIABLES if v not in f] + missing += [a for a in HRIR_VARIABLE_ATTRS + if a.split(':')[0] not in f + or a.split(':')[1] not in f[a.split(':')[0]].attrs] + wrong = [v for v, dims in HRIR_VARIABLES.items() if v in f + and f[v].shape != tuple(sizes[d] for d in dims)] + convention = f.attrs['SOFAConventions'] + print(f'\nmissing mandatory items: {missing}, wrong shapes: {wrong}') + self.assertEqual(missing, []) + self.assertEqual(wrong, []) + self.assertEqual(convention, b'SimpleFreeFieldHRIR') + + def test_read_matches_matlab_ncread(self): + """A SOFA file written by MATLAB reads as MATLAB ncread reads it.""" + ref = loadmat(os.path.join(GOLDEN, 'sofa_golden.mat')) + out = sofa.read(os.path.join(GOLDEN, 'sofa_golden.sofa')) + d_ir = np.max(np.abs(out['ir'] - ref['IR'])) + d_src = np.max(np.abs(out['source_position']['values'] + - ref['SourcePosition'])) + d_rec = np.max(np.abs(out['receiver_position']['values'] + - ref['ReceiverPosition'].reshape(2, 3, 1))) + d_sig = max(np.max(np.abs(s.timeSignal - ref['IR'][m].T)) + for m, s in enumerate(out['signal'])) + print(f'\nvs ncread: max |dIR| = {d_ir}, max |dSource| = {d_src}, ' + f'max |dReceiver| = {d_rec}, max |dSignal| (float32) = {d_sig}') + self.assertEqual(d_ir, 0.0) + self.assertEqual(d_src, 0.0) + self.assertEqual(d_rec, 0.0) + self.assertLess(d_sig, 1e-7) + self.assertEqual(out['fs'], ref['SamplingRate'].item()) + np.testing.assert_array_equal(out['delay'], ref['Delay']) + self.assertEqual(out['source_position']['Type'], + str(ref['SourcePositionType'][0])) + self.assertEqual(out['source_position']['Units'], + str(ref['SourcePositionUnits'][0])) + self.assertEqual(out['attributes']['Title'], str(ref['Title'][0])) + + +if __name__ == '__main__': + unittest.main()