Starting with hrtfpykit.hrtf

hrtfpykit.hrtf is the HRTF/HRIR object layer built on top of SOFA files that follow SimpleFreeFieldHRIR or SimpleFreeFieldHRTF. The goal is to show how a measured SOFA file becomes an HRTF object with synchronized time-domain HRIR data, frequency-domain HRTF data, source positions, immutable transforms, built-in plots, metrics, spherical-harmonic representations, and access to the backing SOFA object.

Use Set Up first if hrtfpykit is not available in your notebook environment. The opening cells prepare one measured HUTUBS SOFA file with the HUTUBS 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, plotted, 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 HUTUBS HRTF file

The file preparation step mirrors the SOFA notebook: a HUTUBS dataset object prepares one measured SimpleFreeFieldHRIR SOFA file for pp1. After the file exists locally, the workflow switches to the hrtfpykit.hrtf layer.

from pathlib import Path

from hrtfpykit.datasets import HUTUBS

# Define the local dataset root.
root = Path("datasets/hutubs")

# Exclude all HUTUBS subjects except pp1.
exclude_subject_ids = tuple(f"pp{i}" for i in range(2, 97))

# Download only the measured HRTF resource for pp1.
HUTUBS(
    root=root,
    download=True,
    download_resources="hrtf",
    download_hrtf_variant="measured",
    exclude_subject_ids=exclude_subject_ids,
    verify_checksum=True,
)

# Build the path to the downloaded SOFA file.
sofa_path = root / "pp1_HRIRs_measured.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

# 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. If a workflow needs global attributes, variable attributes, raw SOFA variables, dimension sizes, or provenance metadata, those remain available without reopening the file separately.

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 is intentional: users can 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)

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 (s):", ir.ir_duration)
print("ITD shape:", ir.get_itd().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 range 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. HRIR sample cropping is applied along the final time axis. After cropping, the returned object recomputes TF using the cropped IR length as the FFT length, so the TF bin count and fft_length reflect the selected signal. If a workflow needs more frequency resolution after cropping, use modify_fft_length on the selected object.

# Select named positions, keep only the left ear, and crop 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 cropping 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 elevation. For plane="median" and plane="frontal", plane_angle is the requested 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 available for comparison, reset, or another processing branch.

Time-domain transforms operate on IR.values and rebuild TF. Frequency-domain transforms operate on TF.values and rebuild IR. This keeps later selection, plotting, metric, and export workflows based on a coherent HRTF state.

# Apply independent transforms to the same selected HRTF state.
windowed = angle_subset.transform.apply_window("hann")
padded = angle_subset.transform.apply_padding(40)
gained = angle_subset.transform.apply_gain(-3.0, scale="db")

# Compare the derived objects with the source object.
print("Source IR shape:", angle_subset.IR.values.shape)
print("Windowed IR shape:", windowed.IR.values.shape)
print("Padded IR shape:", padded.IR.values.shape)
print("Gained transformed flag:", gained.is_transformed())
print("Source transformed flag:", angle_subset.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 such as itd_difference, ild_difference, and lsd compare two HRTF objects on the same source grid. The source positions must match because each metric compares values position by position. If either object has been selected, both objects 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, asymmetric ear gain creates a level change for ILD, and global gain creates a simple spectral magnitude change for LSD.

from hrtfpykit.hrtf import ild_difference, itd_difference, lsd

# 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 an asymmetric ear-level change.
ear_gain_db = np.array([0.0, -6.0], dtype=float).reshape(1, 2, 1)
level_processed = reference.transform.apply_gain(ear_gain_db, scale="db")

# 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")
ild_diff = ild_difference(reference, level_processed, mode="broad-band")
lsd_global = lsd(reference, spectral_processed, reduction="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 or from its TF view. 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}")

Plot source grids and spatial planes

Built-in HRTF plotting methods are called from the loaded HRTF object. plot_source_grid visualizes the current source grid, and plot_plane_grid highlights canonical spatial planes. These methods read the current Sources state, so selected HRTF objects plot their selected grid rather than the original full grid.

These plots are often the first sanity check after loading or selecting data. They make it clear whether the file uses the expected spatial coverage, whether a subset contains the intended directions, and which source positions belong to the horizontal, median, or frontal planes.

import matplotlib.pyplot as plt

# Plot the full source grid attached to the loaded HRTF.
hrtf.plot_source_grid(show=False)

# Plot the source grid with canonical spatial planes highlighted.
hrtf.plot_plane_grid(
    plane=["horizontal", "median", "frontal"],
    show=False,
)

# Plot the source grid after azimuth/elevation selection.
angle_subset.plot_source_grid(show=False)

# Display the figures in an interactive notebook session.
plt.show()

Plot amplitude and magnitude responses

plot_amplitude shows time-domain HRIR waveforms, plot_magnitude shows frequency-domain HRTF magnitudes, and plot_amplitude_and_magnitude places both views for one source direction in the same figure. These plots use the current IR and TF state.

Amplitude plots are useful for checking delays, padding, windowing, and time-domain artifacts. Magnitude plots are useful for checking spectral cues, filtering, gain changes, and frequency resolution. Looking at both views for the same source direction helps confirm that a transform had the intended effect in both domains.

# Plot HRIR amplitude for one source direction.
reference.plot_amplitude(
    positions=["front"],
    ear="both",
    x_axis="samples",
    show=False,
)

# Plot HRTF magnitude for the same source direction.
reference.plot_magnitude(
    positions=["front"],
    ear="both",
    x_axis="log",
    show=False,
)

# Plot the HRIR amplitude and HRTF magnitude together.
reference.plot_amplitude_and_magnitude(
    position="front",
    ear="both",
    amplitude_x_axis="samples",
    magnitude_x_axis="log",
    show=False,
)

# Display the figures in an interactive notebook session.
plt.show()

Plot spectral planes

plot_spectrum_plane shows HRTF magnitude over a measured horizontal or median plane. plot_elevation_spectrum shows how the magnitude response changes across elevation for a fixed azimuth slice.

These heatmaps are useful when the question is spatial rather than one-directional. Instead of inspecting a single front or left response, the plot shows how spectral structure changes across an entire measured plane or vertical slice.

# Plot a horizontal-plane magnitude spectrum for the left ear.
hrtf.plot_spectrum_plane(
    plane="horizontal",
    elevation_angle=0.0,
    ear="left",
    x_axis="linear",
    show=False,
)

# Plot the front-facing elevation spectrum for the left ear.
hrtf.plot_elevation_spectrum(
    azimuth="front",
    ear="left",
    x_axis="linear",
    show=False,
)

# Display the figures in an interactive notebook session.
plt.show()

Plot ITD and ILD cues

Cue plots expose binaural timing and level differences around the source grid. plot_itd_curve and plot_ild_curve show signed curves over a horizontal plane. plot_absolute_itd and plot_absolute_ild show absolute cue magnitudes in polar form. plot_ild_plane shows frequency-dependent ILD over a measured plane.

These plots are useful after loading, selecting, or transforming an HRTF because ITD and ILD are compact cues that can reveal obvious spatial inconsistencies. Curve plots show broad directional trends, while plane plots show how frequency-dependent level differences vary across space.

# Plot signed ITD and absolute ITD around the horizontal plane.
hrtf.plot_itd_curve(elevation_angle=0.0, show=False)
hrtf.plot_absolute_itd(elevation_angle=0.0, show=False)

# Plot signed ILD, absolute ILD, and frequency-dependent ILD.
hrtf.plot_ild_curve(elevation_angle=0.0, show=False)
hrtf.plot_absolute_ild(elevation_angle=0.0, show=False)
hrtf.plot_ild_plane(
    plane="horizontal",
    elevation_angle=0.0,
    show=False,
)

# Display the figures in an interactive notebook session.
plt.show()

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.

# Create an export object with selected positions and cropped 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.

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 cropped HRTF as a new SOFA file.
saved_path = export_hrtf.save(
    path=output_dir / "pp1_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)