Starting with hrtfpykit.hrtf¶
hrtfpykit.hrtf is the main interface for HRTF and HRIR data in hrtfpykit. It loads SOFA files that follow SimpleFreeFieldHRIR or SimpleFreeFieldHRTF and represents them as an HRTF object. That object keeps the IR samples, TF values, source positions, transforms, metrics, spherical harmonic tools, and backing SOFA object in one place.
Use Set Up first if hrtfpykit is not available in your notebook environment. The opening cells prepare one measured SONICOM SOFA file with the SONICOM dataset class, then start the HRTF workflow with load_hrtf. From there, the object is inspected, selected, transformed, reset from its backing SOFA object, compared with metric functions, represented with spherical harmonics, passed to plot workflows, synchronized back to SOFA, saved, and reloaded.
By the end, you should understand how hrtfpykit keeps HRIR and HRTF representations together in one object, how selections and transforms create new HRTF objects instead of modifying the original in place, how the backing SOFA object remains available, how metric and spherical-harmonic workflows use the current HRTF state, and how the HRTF object becomes the main working surface for acoustic processing and inspection.
Download one SONICOM HRTF file¶
The file preparation step mirrors the SOFA notebook: a SONICOM dataset object prepares one measured SimpleFreeFieldHRIR SOFA file for P0001. After the file exists locally, the workflow switches to the hrtfpykit.hrtf layer.
from pathlib import Path
from hrtfpykit.datasets import SONICOM
# Define the local dataset root.
root = Path("datasets/sonicom")
# Keep this tutorial small: use only P0001.
selected_subject_ids = ("P0001",)
# Download only the measured 44.1 kHz FreeFieldComp HRTF resource for P0001.
SONICOM(
root=root,
download=True,
download_resources="hrtf",
download_hrtf_variant={
"type": "measured",
"sample_rate": 44100,
"version": "FreeFieldComp",
},
download_server="imperial",
download_subject_ids=selected_subject_ids,
subject_ids=selected_subject_ids,
verify_checksum=True,
)
# Build the path to the downloaded SOFA file.
sofa_path = root / "P0001" / "HRTF" / "HRTF" / "44kHz" / "P0001_FreeFieldComp_44kHz.sofa"
# Stop early if the expected file is not available.
if not sofa_path.exists():
raise FileNotFoundError(f"Expected SOFA file was not found: {sofa_path}")
print(sofa_path)
Load the file as an HRTF object¶
load_hrtf is the public loader for SOFA-backed HRTF workflows. It accepts files declared with the SimpleFreeFieldHRIR and SimpleFreeFieldHRTF conventions, opens them through hrtfpykit.sofa, validates that the declared convention is supported by the HRTF layer, and returns an HRTF object. The returned object is the acoustic working object used by the rest of the tutorial.
For SimpleFreeFieldHRIR, the file stores time-domain HRIR data in Data.IR with the sampling rate in Data.SamplingRate. hrtfpykit loads those samples into IR, then derives the synchronized frequency-domain TF view with a real FFT. For SimpleFreeFieldHRTF, the file stores the frequency-domain data in Data.Real, Data.Imag, and the frequency vector N. hrtfpykit combines the real and imaginary parts into a complex TF view, normalizes missing DC when needed, and reconstructs the time-domain IR view with an inverse real FFT.
This means users do not need to branch their workflow based on whether the source file was HRIR or HRTF. After loading, the object exposes both synchronized representations, the resolved fft_length, the source positions, and the backing SOFA object.
from hrtfpykit.hrtf import load_hrtf, itd
# Load the SOFA file as a synchronized HRTF object.
hrtf = load_hrtf(sofa_path)
# Inspect the loaded convention and primary array shapes.
print("Convention:", hrtf.SOFAConventions)
print("IR shape:", hrtf.IR.values.shape)
print("TF shape:", hrtf.TF.values.shape)
print("FFT length:", hrtf.fft_length)
Access the backing SOFA object¶
The loaded HRTF object keeps the original SOFA object available through Sofa. This matters because the HRTF object is the acoustic processing layer, while the SOFA object is still the file-structure layer. With the default loader settings, global attributes, variable attributes, raw SOFA variables, dimension sizes, and provenance metadata remain available from the open backing dataset.
Processing operations update the in-memory HRTF state first. The backing SOFA object is synchronized only when update_sofa or save is called. That separation lets users inspect, select, transform, compare, and plot the HRTF before deciding what should be written back to a SOFA representation.
# Access the backing SOFA object from the HRTF instance.
sofa = hrtf.Sofa
# Inspect SOFA metadata without leaving the HRTF workflow.
print("SOFA path:", sofa.path)
print("SOFA convention:", sofa.GlobalAttributes.get("SOFAConventions").value)
print("Data.IR shape in SOFA:", sofa.Variables.get("Data.IR").value.shape)
Load with a closed SOFA handle¶
load_hrtf also accepts sofa_open=False. In that mode, hrtfpykit reads the IR, TF, and source positions into the HRTF object, then closes the backing netCDF4 dataset handle. This is useful when a workflow loads many HRTFs and does not need direct SOFA metadata immediately.
The acoustic state remains available through IR, TF, and Sources. SOFA-backed accessors such as GlobalAttributes, Variables, summary, update_sofa, and save require an open dataset. Reopen a file-backed object with Sofa.open() when direct SOFA access or persistence is needed.
# Load the acoustic state, then close the backing SOFA dataset handle.
closed_sofa_hrtf = load_hrtf(sofa_path, sofa_open=False)
print("SOFA open after load:", closed_sofa_hrtf.Sofa.is_open())
print("Closed-handle IR shape:", closed_sofa_hrtf.IR.values.shape)
print("Closed-handle source grid shape:", closed_sofa_hrtf.Sources.get_positions().shape)
# SOFA-backed metadata requires reopening the backing dataset first.
try:
closed_sofa_hrtf.Sofa.summary()
except ValueError as exc:
print("Closed SOFA access:", exc)
closed_sofa_hrtf.Sofa.open()
print("SOFA open after reopen:", closed_sofa_hrtf.Sofa.is_open())
print("SOFA convention after reopen:", closed_sofa_hrtf.Sofa.GlobalAttributes.get("SOFAConventions").value)
closed_sofa_hrtf.Sofa.close()
Inspect IR and TF views¶
IR and TF are composed interfaces that expose the two synchronized acoustic representations. IR stores time-domain HRIR samples, the sample rate, duration, and ITD helpers. TF stores complex HRTF values, frequency bins, magnitude, phase, real, and imaginary views.
These are not detached copies. They are coordinated views on the same parent HRTF state. Time-domain transforms refresh the TF view, and frequency-domain transforms refresh the IR view. Inspecting both shapes early is useful because most later operations preserve the leading source and ear axes while changing the final sample or frequency axis.
# Access the composed domain interfaces.
ir = hrtf.IR
tf = hrtf.TF
# Inspect time-domain information.
print("IR shape:", ir.values.shape)
print("Sample rate (Hz):", ir.sample_rate)
print("IR duration (ms):", ir.ir_duration)
print("ITD shape:", itd(hrtf).shape)
# Inspect frequency-domain information.
print("TF shape:", tf.values.shape)
print("Frequency bins (Hz):", tf.frequency_bins[:8])
print("Magnitude shape:", tf.magnitude.shape)
print("Phase shape:", tf.phase.shape)
Inspect source positions¶
Sources reads the SOFA SourcePosition variable and provides the source-grid operations used by selection, metrics, and plots. It can return positions in the active coordinate system, resolve named directions such as front, back, left, and right, and track which original source indices are kept after a selection.
Source positions are central to HRTF workflows because arrays are indexed by measurement position along the leading axis. If two HRTFs are compared, plotted, or exported after selection, hrtfpykit uses the current source grid to keep the acoustic data and spatial metadata aligned.
# Read source positions in spherical coordinates.
source_positions = hrtf.Sources.get_positions(angle_unit="degrees")
# Inspect the source grid.
print("Source grid shape:", source_positions.shape)
print("First five positions (azimuth degrees, elevation degrees, radius meters):")
print(source_positions[:5])
# Resolve a named source position.
front_index, front_position = hrtf.Sources.get_position_index("front")
print("Front index:", front_index)
print("Front position (azimuth degrees, elevation degrees, radius meters):", front_position)
Select positions, ears, and samples¶
select returns a new HRTF object and leaves the original object unchanged. It can select source positions, one ear, and an HRIR sample interval in the same call. Source selection can be expressed with named positions, numeric positions, a measured plane, or azimuth/elevation angle filters. Those source-selection modes are mutually exclusive because each one defines a different way to choose the leading source axis.
Ear selection with left or right removes the ear axis from the returned IR and TF arrays. The sample arguments are a selection operation: start_sample and end_sample keep that interval, shorten the IR, and recompute TF with an fft_length that reflects the selected signal. If a workflow needs more frequency resolution after sample selection, use modify_fft_length on the selected object. If the workflow needs to remove samples from the middle of an IR while preserving the original signal length, use apply_crop instead.
# Select named positions, keep only the left ear, and keep the first 128 HRIR samples.
selected = hrtf.select(
positions=["front", "left", "right"],
ear="left",
start_sample=0,
end_sample=128,
)
# Inspect the selected object.
print("Selected IR shape:", selected.IR.values.shape)
print("Selected TF shape:", selected.TF.values.shape)
print("Selected FFT length:", selected.fft_length)
# Confirm the original HRTF object was not changed.
print("Original IR shape:", hrtf.IR.values.shape)
print("Original TF shape:", hrtf.TF.values.shape)
# Increase the FFT length after sample selection when higher frequency resolution is needed.
selected_high_resolution = selected.transform.modify_fft_length(256)
print("High-resolution TF shape:", selected_high_resolution.TF.values.shape)
print("High-resolution FFT length:", selected_high_resolution.fft_length)
Select by azimuth and elevation¶
select can filter the source grid with azimuth_angles and elevation_angles. Requested values are not treated as continuous interpolation targets. Instead, each requested angle resolves to the nearest available measured spherical angle in the source grid, then the returned object keeps the sources that match the resolved azimuth and elevation filters.
This is useful when the source grid is dense and users want a clear directional subset without manually looking up exact source rows. The two angle filters can be used together, as shown below, but they cannot be mixed with positions or plane in the same call.
# Select the nearest measured sources at three azimuths on the horizontal plane.
angle_subset = hrtf.select(
azimuth_angles=[0, 90, 180],
elevation_angles=0,
)
# Inspect the selected source grid and data shapes.
print("Angle subset positions (azimuth degrees, elevation degrees, radius meters):")
print(angle_subset.Sources.get_positions(angle_unit="degrees"))
print("Angle subset IR shape:", angle_subset.IR.values.shape)
print("Angle subset TF shape:", angle_subset.TF.values.shape)
Select a measured plane¶
Plane selection keeps all sources that belong to the nearest measured plane. For plane="horizontal", plane_angle is the requested spherical elevation. For plane="median", it is the requested lateral-polar lateral angle. For plane="frontal", it is the requested spherical azimuth. The selected plane is still constrained by the positions that actually exist in the measured grid.
This is the right selection mode when a workflow needs a complete plane for inspection, cue analysis, or plotting. It preserves the plane ordering used by hrtfpykit plotting and metric helpers, instead of selecting only a few named directions.
import numpy as np
# Select the measured horizontal plane nearest 0 degrees elevation.
horizontal = hrtf.select(
plane="horizontal",
plane_angle=0.0,
)
# Inspect the selected plane.
horizontal_positions = horizontal.Sources.get_positions(angle_unit="degrees")
print("Horizontal plane shape:", horizontal_positions.shape)
print("Unique elevations (degrees):", np.unique(np.round(horizontal_positions[:, 1], 2)))
Apply transforms¶
transform exposes non-mutating processing operations. Each transform clones the parent HRTF object, applies one operation, refreshes the synchronized representation, marks the returned object as transformed, and leaves the source object unchanged.
Time-domain transforms such as apply_window, apply_padding, apply_crop, upsampling, add_itd, and delete_itd operate on IR.values and rebuild TF. Frequency-domain transforms such as apply_gain, add_ild, and delete_ild operate on TF.values and rebuild IR. modify_fft_length keeps the current IR and rebuilds TF with a different FFT length.
The short example below shows transform outputs through state fields that actually change: sample rate, FFT length, and IR sample count.
# Apply independent transforms to the loaded HRTF.
upsampled = hrtf.transform.upsampling(hrtf.IR.sample_rate * 2.0)
fft_modified = hrtf.transform.modify_fft_length(hrtf.fft_length * 2)
padded = hrtf.transform.apply_padding(40)
cropped = hrtf.transform.apply_crop(
start_sample=32,
end_sample=64,
window="hann",
)
# Inspect state changes on the returned objects.
print("Original sample rate (Hz):", hrtf.IR.sample_rate)
print("Upsampled sample rate (Hz):", upsampled.IR.sample_rate)
print("Original FFT length:", hrtf.fft_length)
print("Modified FFT length:", fft_modified.fft_length)
print("Original IR shape:", hrtf.IR.values.shape)
print("Padded IR shape:", padded.IR.values.shape)
print("Cropped transformed flag:", cropped.is_transformed())
print("Original transformed flag:", hrtf.is_transformed())
Reset from the backing SOFA object¶
reset restores the in-memory HRTF state from the attached backing SOFA object. This is useful when a selected or transformed object should return to the backed file state without reopening the SOFA file manually.
Reset is not a general undo stack. It reloads the current backing SOFA content into the HRTF object, rebuilds the synchronized IR/TF views, clears source-selection state, and marks the object as not transformed. If the backing SOFA object was already synchronized with update_sofa, reset restores that synchronized backing state.
# Create a selected and transformed object.
reset_example = hrtf.select(
positions=["front", "left", "right"],
start_sample=0,
end_sample=128,
).transform.apply_gain(-3.0, scale="db")
# Inspect the selected and transformed state.
print("Before reset IR shape:", reset_example.IR.values.shape)
print("Before reset TF shape:", reset_example.TF.values.shape)
print("Before reset transformed flag:", reset_example.is_transformed())
# Restore the object from its backing SOFA state.
reset_example.reset()
# Inspect the restored state.
print("After reset IR shape:", reset_example.IR.values.shape)
print("After reset TF shape:", reset_example.TF.values.shape)
print("After reset transformed flag:", reset_example.is_transformed())
Compute differences between matching HRTFs¶
Metric functions compare HRTF objects on the same source grid. itd_difference and ild_difference compare one reference HRTF against one or more HRTFs. hrtf_difference compares HRTFs with IR-domain metrics (rmse, mae, nrmse) or TF-domain log-spectral distortion (lsd). Source positions must match because each metric compares values source by source. If an object has been selected, every compared object should expose the same selected grid in the same order.
The examples below keep the full source grid with clone, then use separate processed copies for timing, level, and spectral changes. This makes each metric easier to interpret: add_itd() creates a timing change for ITD, add_ild() creates a controlled interaural level change for ILD, and apply_gain() creates a simple spectral magnitude change for metric="lsd".
from hrtfpykit.hrtf import hrtf_difference, ild_difference, itd_difference
# Build a full-grid reference copy.
reference = hrtf.clone()
# Create a copy with a controlled timing change.
timing_processed = reference.transform.add_itd(20, unit="samples")
# Create a copy with a controlled interaural level change.
level_processed = reference.transform.add_ild(6.0)
# Create a copy with a global spectral magnitude change.
spectral_processed = reference.transform.apply_gain(-1.0, scale="db")
# Compare timing, level, and spectral differences.
itd_diff = itd_difference(reference, timing_processed, output="samples", absolute=True, reduction_axis="itds")
ild_diff = ild_difference(reference, level_processed, mode="broad-band", absolute=True, reduction_axis="ilds")
lsd_global = hrtf_difference(reference, spectral_processed, metric="lsd", reduction_axis="global")
# Print compact summaries instead of full per-position arrays.
print("ITD difference shape:", itd_diff.shape)
print("ITD difference range (samples):", f"{np.min(itd_diff):.1f}", f"{np.max(itd_diff):.1f}")
print("ILD difference shape:", ild_diff.shape)
print("ILD difference range (dB):", f"{np.min(ild_diff):.1f}", f"{np.max(ild_diff):.1f}")
print("Global LSD (dB):", f"{lsd_global:.1f}")
Build spherical-harmonic representations¶
sht builds a spherical-harmonic representation of HRTF magnitudes from the current HRTF object. This belongs in the HRTF workflow because the decomposition uses the object’s current source grid, frequency bins, ear selection, and magnitude data. If the HRTF has been selected or transformed before this step, the SH representation is built from that current state.
The SH object stores two main arrays: C, the coefficient tensor, and Y, the real spherical-harmonic basis evaluated on the source positions. For order L, the coefficient count is (L + 1) ** 2. A decomposition with ear="both" keeps a two-ear axis in C, so the coefficient tensor has shape (coefficients, ears, frequency_bins).
The code below decomposes the full loaded HRTF magnitude with a low SH order, reconstructs magnitude values on the same measured source grid with sht_inverse, and compares the reconstruction against the original magnitude with sht_error. The reconstruction is magnitude-only. It does not reconstruct phase and it does not create a new complex HRTF object.
from hrtfpykit.hrtf import sht, sht_error, sht_inverse
# Build a spherical-harmonic representation from the current HRTF magnitude.
sh_order = 4
sh = sht(reference, sh_order=sh_order, ear="both")
# Reconstruct magnitudes on the same source grid used for the decomposition.
reconstructed_magnitude = sht_inverse(sh)
original_magnitude = np.abs(reference.TF.values[:, 0:2, :])
# Evaluate global reconstruction error in the dB domain.
abs_err, rel_err, rms_err, max_err = sht_error(
original_magnitude=original_magnitude,
reconstructed_magnitude=reconstructed_magnitude,
magnitude="db",
reference="max",
)
print("SH order:", sh_order)
print("Coefficient count:", (sh_order + 1) ** 2)
print("SH coefficients shape:", sh.C.shape)
print("SH basis shape:", sh.Y.shape)
print("Original magnitude shape:", original_magnitude.shape)
print("Reconstructed magnitude shape:", reconstructed_magnitude.shape)
print("RMS reconstruction error (dB):", f"{rms_err:.2f}")
print("Maximum reconstruction error (dB):", f"{max_err:.2f}")
Synchronize the HRTF object back to SOFA¶
The HRTF object stores working IR and TF arrays separately from its backing SOFA object. This avoids writing every intermediate selection or transform back into the file structure automatically. update_sofa is the explicit synchronization step that writes the current HRTF state into the backing SOFA object.
When selections or transforms change fixed SOFA dimensions such as M for source positions or N for samples/frequency bins, pass change_sofa_dimensions=True. Without that flag, hrtfpykit refuses to silently force arrays into incompatible netCDF dimensions. The backing SOFA dataset must also be open; if the object was loaded with sofa_open=False or closed manually, call Sofa.open() before synchronization.
# Create an export object with selected positions and shortened HRIR samples.
export_hrtf = hrtf.select(
positions=["front", "left", "right"],
start_sample=0,
end_sample=128,
)
# Synchronize the selected HRTF state into the backing SOFA object.
export_hrtf.update_sofa(change_sofa_dimensions=True)
# Inspect the synchronized SOFA dimensions and variables.
print("M dimension:", export_hrtf.Sofa.Dimensions.get("M").value)
print("N dimension:", export_hrtf.Sofa.Dimensions.get("N").value)
print("Data.IR shape:", export_hrtf.Sofa.Variables.get("Data.IR").value.shape)
print("SourcePosition shape:", export_hrtf.Sofa.Variables.get("SourcePosition").value.shape)
Save and reload the processed HRTF¶
save synchronizes the current HRTF state and writes the backing SOFA object to disk. The output convention defaults to the same convention as the backed file, but save can also forward convention options to update_sofa when a workflow needs HRIR or HRTF SOFA output explicitly.
SOFA persistence is a two-ear export path. HRTFs created with select(ear="left") or select(ear="right") are useful analysis objects, but they are not valid objects for this SOFA save workflow because the exported file must contain both receiver channels. Keep ear="both" for HRTFs that will be synchronized with SOFA or saved to disk.
Reloading the saved file with load_hrtf is a practical final check. It confirms that the exported SOFA file can be opened again, converted into a synchronized HRTF object, and used by the same inspection, plotting, metric, and dataset workflows.
# Create an output folder for tutorial artifacts.
output_dir = Path("tutorial_outputs")
output_dir.mkdir(exist_ok=True)
# Save the selected and shortened HRTF as a new SOFA file.
saved_path = export_hrtf.save(
path=output_dir / "P0001_selected_hrtf.sofa",
overwrite=True,
change_sofa_dimensions=True,
)
# Reload the saved file through the HRTF API.
reloaded = load_hrtf(saved_path)
# Confirm that the saved file preserves the exported HRTF state.
print("Saved path:", saved_path)
print("Reloaded convention:", reloaded.SOFAConventions)
print("Reloaded IR shape:", reloaded.IR.values.shape)
print("Reloaded TF shape:", reloaded.TF.values.shape)
print("Reloaded FFT length:", reloaded.fft_length)