From 169e51f5e6a3c91b185cd91b707bf83ff30f855a Mon Sep 17 00:00:00 2001 From: Sergio Aguirre Date: Sat, 31 Oct 2026 21:56:42 +0100 Subject: [PATCH] [FEATURE] add pytta.sofa to read and write SOFA files (AES69) New module pytta/sofa.py reads and writes impulse responses in SOFA (AES69-2022, SOFA 2.1) with h5py only. The writer creates the netCDF-4 dimension scales and attributes itself, so the files open with MATLAB ncinfo/ncread and other netCDF tools. read(path) returns one SignalObj per measurement with the receivers as channels, the float64 Data.IR, fs, Data.Delay, the source, listener, receiver and emitter positions with Type and Units, and the global attributes. Any convention with Data.IR can be read. write(path, signal, source_positions, convention, **attributes) writes every mandatory attribute and variable of SimpleFreeFieldHRIR 1.0 or GeneralFIR 1.0. Tests: exact write/read round trip, presence and shape of all mandatory SimpleFreeFieldHRIR items, and parity with MATLAB ncread on a small golden file written by tests/golden/make_sofa.m. Co-Authored-By: Claude Fable 5.1 --- docs/pyttadoc/index.rst | 4 + docs/pyttadoc/sofa.rst | 7 + pytta/__init__.py | 5 + pytta/sofa.py | 249 ++++++++++++++++++++++++++++++++++ tests/golden/make_sofa.m | 83 ++++++++++++ tests/golden/sofa_golden.mat | Bin 0 -> 1085 bytes tests/golden/sofa_golden.sofa | Bin 0 -> 60118 bytes tests/test_sofa.py | 127 +++++++++++++++++ 8 files changed, 475 insertions(+) create mode 100644 docs/pyttadoc/sofa.rst create mode 100644 pytta/sofa.py create mode 100644 tests/golden/make_sofa.m create mode 100644 tests/golden/sofa_golden.mat create mode 100644 tests/golden/sofa_golden.sofa create mode 100644 tests/test_sofa.py 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 0000000000000000000000000000000000000000..d25b391f68257b91d8b46e70e3ac8ab01c20193e GIT binary patch literal 1085 zcmeZu4DoSvQZUssQ1EpO(M`+DN!3vZ$Vn_o%P-2c0*X01nwjV*I2WZRmZYXAkb1@gjQ0QdfZQEGE*p^M0Adgg z0<1vH1H{Zg%mlU@Z(rGW`&dl;`}U5#MrNO>$cFX% zUr*L6d;S0Y{+NC1&!n$7jq~uHLt=GwS)Nd+w(F z($wqUXK7#DKjro*{u@|Z|jykgnZE&Er6^_Gd;Uuo|e zS7XGt`MkYRx{hJL+(~<`?$nnn?6%u~d+_Xgde4dd5sPe@8g87}zyJQ!ziPgR_V3DA z`eoyP~ zR50b>{(S3225Oc&_w!9%@@v<+1^cVdYb?05V!wT(l(R?uioN!%7uDEnYB$@f%L!+Q zaqQT?=~u@WwVeI?O>Rsy($w6yUwFp(RiG$!1!71N&_GTCjE2Sn#!MR-1lX6VEl-=a zRE+hBc>1)IXU_|nj`BztFkb(9lyzIo?Y9fR)Rj~+EPl%u4KmCaZdd`bVM>W|3VeQN zbWfZ-`Bs8KAoHBz=5>J0V@MKX<49=z zB-uDa*ozy@#5j*sM=h7V3?Z?Fj&50jX+=qVS7aDEw7KttjPimT<$-2Ygoqb6Gpgg_ zGUC!a4I@Bid1X}vX07UD>N)4DtLeqs#mpe*B;X7((H3qZ2fB#@Ud>0GjJ%728n-Zn f>{+mBiAvyl3wzx_MO%AyRa5>|@(j;+a7F_F>FSk5 literal 0 HcmV?d00001 diff --git a/tests/golden/sofa_golden.sofa b/tests/golden/sofa_golden.sofa new file mode 100644 index 0000000000000000000000000000000000000000..536d4e8858b4352ef12300cf315383eed06ec529 GIT binary patch literal 60118 zcmeHQ3w%_?xu4A@Aqi^&2!f!vJT4CjF)t{#H!R6!m!u?{bT=YWuefAS$hKiOCJTyU zMY$Dwr7cL%MhX{fsgZ|o3V788Q9(;96dzbXxZe5zR1~XJdX@WqXXZP*8xq-E14a1e zNAjJSGmmfn^PMx{KQrI?SxK?;g0#MA2`MRp!U+{#KO|GK)r}vFvVT-jQtVE!Jee?1 zPcvtzTPwaVFY+cxd5w~33BrUay)F@nAk5vx)%#_RgGy!`ka2OXlAR#Bb;C+h-;!0i zbPa@0!`dkH3!$t*JoL-s_CyAQZbS6B$E=$#^Z-v+so^uCzWu=Xxu zyH6H;Ta~9CrV4_uwX}@DtY(ocQm}4&S)s=xE`%i1HUWN$OUoTDPici~ipO49<`9<2 z0(3EpBw4^!SnkmCk5l=pYWQrrR$G+J*4*f&%aG48TbPnEbBnnAno4L_(jB0x2e)N4|`!t3a`%(G}$ zpUoBak_jn2;I=TE=#i6#D}-e>8ZqjSjXJm;PKVp!vOA{eI*f)&PzRqk<1uJM#~W$9Pyq;wv11UZ4<9`3z6G_i$Yf!uLwGw=Bf^^`X+-g@TTk>xjHF4d z=qnXb?`L6A1D-D5F;iUj%7$RQ))4Z;P~*#Kpz6#_tpVRG1-Y5dW;cFsqQ3-; z6uw{t7y(9r5nu!u0Y-okU<4QeMt~7u1Q>y{1p%FiGzpy$Z;OLN+Oa({TvFjBit#E} zT|%D8*H1=%q0*3P?jZ4fx@hlQfeN;)k(&KVDJb4~=$2QV9*;xZrkF94M~Jn^hFpJO z(x7A@4}lv=B(l%s7JaCs5HGyG_D*6UBOD|#C9>~vaqLL4+&i&psWB+PsU?^r#Hi7{ zYMV%z*6I47wZjyj%oBh9zO2r#YVFkn8Bd~8Nm@!f$P1HFBuQDs6o=pg#YEnungx=! zb+~kEo=)#7>M~c}2SxOdJP3s;k;vX|ZmnAL15a5QZ1jhKt0XQKv+e2qLIsYI;7lUIZ7~ApGXkyoe|pc+>_zh&y<{n!_7y(A$3?X1p&A@QCw;lcz{C|#AHJbl-=gG9Q zkpIuq$>+ZJN(JvCm^YU;<{&dlOh56Frcz=w=ngR+rgGR+tcOLJKRrgQY8#xYzV;Cp z|BvL$mfhm}#-I#_L4i>rQl?#ri!^cZ{kT!uZ2Ip#Qt=WgMBC&vZVmoMJHbg^r0)+Y z^1fCziD^>~%1E_J6Bp2hl64O~*c3e|blYSn`YCw;rkN6x)va_-a7IJ5W~=}x+w50! z?NF=;_oRR09qD6aTWBG=m8dwbZ!Fp9`wZJ@8_=_b?I$kJG9L~P=y-J2?syt;lC&3H zjgG7y(9r5nu!u0Y-okU<4R}vkw7-YF7L|{r{h29Q*$jhZoKNmtHS>dKU8kzzBD5 zV;Ek4{@}nilH*P~ng2If*7*5A z+S?|tzw!Qm!UL7WW_MP|_P%ES|GeWLuXMUgD#6_ndQ^cz6zEg}9MaGOo!T(`t`J|P z?#H_e6_hkeN;IuN@;7@`HKKgThmRChl@tSS4*9x5vTxWoz|x)&mx}*L4CIPlkUB_sHZPkdE`!v|#6=IRxLu5e%bud=u2uJl zxq7ibBs{xH+^1ibCgraa2jFs`_+e%3MzKJztjc}I7U71=zGCCp9(zO?WEdp2=6r9z z2oLqr2ix1`qm2x zRF`29C9*!UE=oTlli`Ib1EtflyGrJR1Ox<2Mt~7u1Q-EEfDvE>7y(9r5nu!u0Y>1T zh(P<>p2_7(XKgYjGK(Z_BBst)Z2hcFr{KgdS$-|R2Y+X};3iIkByD>(G*o=85m(#N zjksDuq4Nl@lLO=^9QE6 z!Qp)~@|+M%iNrvr>3NE^I`6FJE%N4R_$}9fhU#Aki1+&d9XSJ1xK}FON9o?qst=H? z-{=eoSygA8ylDQh)D`Q*ABy8MVF@pNiq_yj<&!-dJ&is}I`g3?BAlf007@l|vkrv7 z;!`~IhGnbrhAFFfFW3L1cj?p~(9=z_vwQIV@2`7LDbIotU<4QeMj&1UTCb;hSyo|Fm3XS%V+#?o9rlSaXuZo|7rjmko%Zp=B+K z6W>qyiSMU)!8_r9rE8Xq03*N%FaqZd0?hxPH$%-GWds-jMt~7u1Q-EEfDvE>7y(9r z5jc+#koaIH{>iX89;kAAS)s=xE`&npYYCa^ z;?i=5%Trq6n&Pn+mN|rFGGfvM zj}KRS*pB++qf@6hbW?$I5|%*Ol~-~=Apbvv96x5H(3OwrpN4a#A= zJ~&I$sD&a57-KF{@V;DuvFdu!kdUM^5s*y}bRdV~xo>@$dPek62ehbN(nF1~n7zO? zj8>N-XG^RqJB?0fp+v2cBy9+_7r!W(K$@~-1Q>yH5P{bGKTQiXo6(0kXo1^ER|ooN zS(?)66eK~#ALb?`STX{P03*N%FanGKBftnS0*nA7zz8q`|2zaFt8F&;eimCP z-+adX{!MDXf061?JDhc42RUm}CR$p`Bxf8iR~=mzwZnMx6H-rkZ`tF2J*>09^E;P~ zCQ5}Q?bw%X!>72LqO}06w#3FSp1b5)nbOg{{p3`NTKaqO(gk$Qk`Z787y(A$yhflk z|4-8-@%iNchd|n$L02pp0Y-okU<4QeMt~7u1Q-EEfDvE>7y(A$+k!xQ{(oo_F!2Et zkAFQt?qdI9pF{$;Xk%z+MZnR)6#-I*B&vzmo#Oxh52&bIA;5Z@RUGU8=#JIfXooha zZ214D?U9!4LZo);UUvGV1PnfL0$60SsLi|L(ZF;@^w_BftnS0_P9{@$&yP z(bBqrpbz>eEZG15Ib;;v6-Iy&U<4QeMt~7u1Q-EEfDvE>7=iN+fwVV%Q(us_-qyO~ zDm~mfo#J_V_*9mgGjUwo=`}Kaj&od0e#(4G`LWWAWTe1HmegNXeR#;`d~)ZZ&-?7N z&AO!Yp5?t-Y!CO}yKm!nwuL`G@b0l`^S6ap&po{1gEt-wZ(2IL{+;9vwp(5qw*1d6 zVcX$FSu+oxI1*l|y>aaR5j(@z%$=A$>BG%7$91!F)1P?7_Ehc4zinK($u{Bk;aw*! zYzY?(4^FLb+8iF*|H9gfQZ|I|E*rD(%MbV35=K4u)w0+2+b;g;-*0-n>1Er4OJ5l< zbMcemNe>hx%-Ovy{Mc*X-F{VaIJ~Ie(_bCxvcz_KUwg@IEzj7lefE*<(^ef0H|=>j zWB1K_!t;EAcIm-gwp9~DnrEmpiXgU2D=anwYu0lhrw z*{~dD%gxHMjKSq2Gx!P8b74VyFfdaKgdnr56rD}!mHOeE+nTdM)0}>-&R60tbz81L z??GnpY}8iiwwlbz%C=mNE|Uz!3M)%p-szgI7=Gj%n<_7?Dl05fE~+g5BEJUTRW<_o zqbo))8hT%{b4O3L43OWDNR?UfGvHOT>Sf*3X=kas@Xn@xINc?ckP8db0vkk5s6gvB z1h?{$M81+BrwrZc;uI1FhSs_C<1<1G13B-nRBTu9)B2*Sl48&ZI)kSMO+q3jsshA7Ff?RL6&b4yS`VCGG)Fr4Fv=T{DmOzy%#Z91zT(E=Xx3}b-ekudW zmaVV;x~Qtm3EDl3SdSFoaSF^QcXn@a5$1(8NDOs(`$me3F-JbECxZC=hL5wvDr|Nr zF76589ldw7D90S5WC!l+=W>V}F-M+sjy3U(9pl9vm?KjN{Gj_4=E}+96U+gMb_KC_ z)%Y5*2Xl;&{;hsob?c2H2Xp9_BE+qh`OCx)A&v={{$BRT)nX!E4L}2d*4OzT6?5Q1 zSf%--E)Q)NyD>pGX(4PSM_w1pF+qPqg!tqqx9$@=u@0*c_-Xr+Q71kVT`)m6Xd%{m zT_;>C-dzJJ=U&jqKQ(5!V<3(}GX}UASYrT-8MK)3i7|;73yK+n7`uqEyci>kv8Nal zj#G@6kdu{@pOtThcZFU_>zOA}Re`gTF<(T#{OO$z z6WRZF3hUtHWvFq=Y0C^RE4^1=H({CKt*B>eZ>ll?E}Rjn4K@HmA#Zg^pDc0OLV)XM zg+gAdPn)UL;Q%z6VP??ltZ10#4ft>N;#97eNKz##XH`{sjh>YGtHA?g>?B!UT4ku_ z4hE+ivnam!u+bX$Gr(Gde+L+&2VKzY4I2FR>8;V<4Y@EgkIcy)nU`bD$^CA2HvBVG z#4_bUpC4+9qmc4}S>+GaNp*3eG$6=o4Ayw9(}HzAEfC#+cPe;4H2%cE8N?`J4QY)b z0IY`utjymS(gIq82apaKCCCZ`!GH#Z&(jdSb*Fic1c4SZErSqK!|AADEI8vM*45qy z0NU>jSgi7|O*N2L1^glTvLMTxu7w)3qt0VrQgde0W2&)i*`~n8;%@5e8=BXB;d*T7 z>#J<6uhkm-)!sS_s9-o^jp*J>_0p$JgNGs$oiw(EtJk%m4|o78A=TaIwMumh5zn@9Dwj{QG|al8Ql literal 0 HcmV?d00001 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()