interpolateHRTF

3-D head-related transfer function (HRTF) interpolation

Syntax

interpolatedHRTF = interpolateHRTF(HRTF,sourcePositions,desiredSourcePositions)

interpolatedHRTF = interpolateHRTF(sofaObj,desiredSourcePositions)

interpolateHRTF(sofaObj,desiredSourcePositions)

interpolatedHRTF = interpolateHRTF(___,Name=Value)

Description

interpolatedHRTF = interpolateHRTF(HRTF,sourcePositions,desiredSourcePositions) returns the interpolated head-related transfer function (HRTF) at the desired position.

interpolatedHRTF = interpolateHRTF(sofaObj,desiredSourcePositions) returns the interpolated HRTF using the SOFA file data in sofaObj.

interpolateHRTF(sofaObj,desiredSourcePositions) with no output arguments plots the frequency response of the HRTF at the desired source positions.

interpolatedHRTF = interpolateHRTF(___,Name=Value) specifies options using one or more name-value arguments.

Examples

Render 3-D Audio on Headphones

Modify the 3-D audio image of a sound file by filtering it through a head-related transfer function (HRTF). Set the location of the sound source by specifying the desired azimuth and elevation.

load 'ReferenceHRTF.mat' hrtfData sourcePosition

hrtfData = permute(double(hrtfData),[2,3,1]);

sourcePosition = sourcePosition(:,[1,2]);

Calculate the head-related impulse response (HRIR) using the VBAP algorithm at a desired source position.

desiredAz = 110;
desiredEl = -45;
desiredPosition = [desiredAz desiredEl];

interpolatedIR  = interpolateHRTF(hrtfData,sourcePosition,desiredPosition, ...
                                  Algorithm="vbap");

Create a dsp.AudioFileReader object to read in a file frame by frame. Create an audioDeviceWriter object to play audio to your sound card frame by frame. Create a dsp.FrequencyDomainFIRFilter with the combined left-ear and right-ear filter pair. Since you want to keep the received signals at each ear separate, set SumFilteredOutputs to false.

fileReader = dsp.AudioFileReader("RockDrums-48-stereo-11secs.mp3");
deviceWriter = audioDeviceWriter(SampleRate=fileReader.SampleRate);

spatialFilter = dsp.FrequencyDomainFIRFilter(squeeze(interpolatedIR),SumFilteredOutputs=false);

In an audio stream loop:

1. Read in a frame of audio data.

2. Feed the stereo audio data through the filter.

3. Write the audio to your output device.

while ~isDone(fileReader)
    audioIn = fileReader();

    audioOut = spatialFilter(audioIn);
    deviceWriter(squeeze(audioOut));
end

As a best practice, release your System objects when complete.

release(deviceWriter)
release(fileReader)

Model Moving Source Using HRIR Filtering

Load a SOFA file containing HRTF measurements.

s = sofaread("ReferenceHRTF.sofa");

Specify the desired source positions and then calculate the HRTF at these locations using the interpolateHRTF function.

desiredAz = [-120;-60;0;60;120;0;-120;120];
desiredEl = [-90;90;45;0;-45;0;45;45];
desiredPosition = [desiredAz desiredEl];

interpolatedIR  = interpolateHRTF(s,desiredPosition);

Create an audio file sampled at 48 kHz for compatibility with the SOFA file.

[audio,fs] = audioread("Counting-16-44p1-mono-15secs.wav");
audio = audioresample(audio,InputRate=fs,OutputRate=s.SamplingRate);
audio = audio./max(abs(audio));
audiowrite("Counting-16-48-mono-15secs.wav",audio,s.SamplingRate);

Create a dsp.AudioFileReader object to read in a file frame by frame. Create an audioDeviceWriter object to play audio to your sound card frame by frame. Create a dsp.FrequencyDomainFIRFilter object. Set the Numerator property to the combined left-ear and right-ear filter pair. Since you want to keep the received signals at each ear separate, set SumFilteredOutputs to false.

fileReader = dsp.AudioFileReader("Counting-16-48-mono-15secs.wav");
deviceWriter = audioDeviceWriter(SampleRate=fileReader.SampleRate);

spatialFilter = dsp.FrequencyDomainFIRFilter(squeeze(interpolatedIR(1,:,:)),SumFilteredOutputs=false);

In an audio stream loop:

1. Read in a frame of audio data.

2. Feed the audio data through filter.

3. Write the audio to your output device. If you have a stereo output hardware, such as headphones, you can hear the source shifting position over time.

4. Modify the desired source position in 2-second intervals by updating the filter coefficients.

durationPerPosition = 2;
samplesPerPosition = durationPerPosition*fileReader.SampleRate;
samplesPerPosition = samplesPerPosition - rem(samplesPerPosition,fileReader.SamplesPerFrame);

sourcePositionIndex = 1;
samplesRead = 0;
while ~isDone(fileReader)
    audioIn = fileReader();
    samplesRead = samplesRead + fileReader.SamplesPerFrame;
    
    audioOut = spatialFilter(audioIn);
    deviceWriter(audioOut);
    
    if mod(samplesRead,samplesPerPosition) == 0
        sourcePositionIndex = sourcePositionIndex + 1;
        spatialFilter.Numerator = squeeze(interpolatedIR(sourcePositionIndex,:,:));
    end
end

As a best practice, release your System objects when complete.

release(deviceWriter)
release(fileReader)

Interpolate SOFA HRTF and Visualize Frequency Response

Load a SOFA file containing HRTF measurements.

s = sofaread("ReferenceHRTF.sofa");

Specify the desired source positions and then calculate the HRTF at these locations.

desiredAz = 0:40:180;
desiredEl = -30:15:30;
desiredPosition = [desiredAz', desiredEl'];

interpolatedIR = interpolateHRTF(s, desiredPosition);

Call interpolateHRTF with no output arguments to plot the frequency response of the HRTF.

interpolateHRTF(s, desiredPosition)
legend(Location="southeast")

Input Arguments

HRTF — HRTF values measured at source positions

N-by-2-by-M array

HRTF values measured at the source positions, specified as a N-by-2-by-M array.

• N –– Number of known HRTF pairs

• M –– Number of samples in each known HRTF

If you specify HRTF with real numbers, the function assumes that the input represents an impulse response, and M corresponds to the length of the impulse response. If you specify HRTF with complex numbers, the function assumes that the input represents a transfer function, and M corresponds to the number of bins in the frequency response. The output of the interpolateHRTF function has the same complexity and interpretation as the input.

Data Types: single | double
Complex Number Support: Yes

sourcePositions — Source positions corresponding to measured HRTF values

N-by-2 matrix

Source positions corresponding to measured HRTF values, specified as a N-by-2 matrix. N is the number of known HRTF pairs. The two columns correspond to the azimuth and elevation of the source in degrees, respectively.

Azimuth must be in the range [−180,360]. You can use the −180 to 180 convention or the 0 to 360 convention.

Elevation must be in the range [−90,180]. You can use the −90 to 90 convention or the 0 to 180 convention.

Data Types: single | double

sofaObj — SOFA object

SimpleFreeFieldHRIR object | SimpleFreeFieldHRTF object

Object containing HRTF data from a SOFA file, specified as a SimpleFreeFieldHRIR or SimpleFreeFieldHRTF object. Use the sofaread function to read in data from SOFA files.

desiredSourcePositions — Desired source positions for HRTF interpolation

P-by-2 matrix

Desired source position for HRTF interpolation, specified as a P-by-2 matrix. P is the number of desired source positions. The columns correspond to the desired azimuth and elevation of the source in degrees, respectively.

Azimuth must be in the range [−180,360]. You can use the −180 to 180 convention or the 0 to 360 convention.

Elevation must be in the range [−90,180]. You can use the −90 to 90 convention or the 0 to 180 convention.

Data Types: single | double

Name-Value Arguments

Specify optional pairs of arguments as Name1=Value1,...,NameN=ValueN, where Name is the argument name and Value is the corresponding value. Name-value arguments must appear after other arguments, but the order of the pairs does not matter.

Before R2021a, use commas to separate each name and value, and enclose Name in quotes.

Example: Algorithm="vbap"

Algorithm — Interpolation algorithm

"bilinear" (default) | "vbap" | "nearest"

Interpolation algorithm, specified as "bilinear", "vbap", or "nearest".

• "bilinear" –– 3-D bilinear interpolation, as specified by [1].

• "vbap" –– Vector base amplitude panning interpolation, as specified by [2].

• "nearest" — Nearest neighbor interpolation, using the Euclidean distance to find the nearest neighbor.

Data Types: char | string

Output Arguments

interpolatedHRTF — Interpolated HRTF

P-by-2-by-M

Interpolated HRTF, returned as a P-by-2-by-M array.

• P –– Number of desired source positions, specified by the number of rows in the desiredSourcePositions input argument.

• M –– Number of samples in each known HRTF, specified by the number of pages in the HRTF input argument.

interpolatedHRTF has the same complexity and interpretation as the input. If you specify the input, HRTF, with real numbers, the function assumes that the input represents an impulse response. If you specify the input with complex numbers, the function assumes that the input represents a transfer function.

Data Types: single | double
Complex Number Support: Yes

More About

SOFA Files

Spatially oriented format for acoustics (SOFA) is a file format for storing spatially oriented acoustic data like HRTFs and binaural or spatial room impulse responses. The AES69 standard [5] defines the SOFA file format.

References

[1] F.P. Freeland, L.W.P. Biscainho and P.S.R. Diniz, "Interpolation of Head-Related Transfer Functions (HRTFS): A multi-source approach." 2004 12th European Signal Processing Conference. Vienna, 2004, pp. 1761–1764.

[2] Pulkki, Ville. "Virtual Sound Source Positioning Using Vector Based Amplitude Panning." Journal of Audio Engineering Society. Vol. 45. Issue 6, pp. 456–466.

[3] Majdak, P., Zotter, F. Brinkmann, F., De Muynke, J., Mihocic, M., and Noisternig, M., “Spatially Oriented Format for Acoustics 2.1: Introduction and Recent Advances.” Journal of the Audio Engineering Society 70, no. 7/8 (July 25, 2022): 565–84. https://doi.org/10.17743/jaes.2022.0026.

[4] Majdak, P., Carpentier, T., Nicol, R., Roginska, A., Suzuki, Y., Watanabe, K., Wierstorf, H., et al., “Spatially Oriented Format for Acoustics: A Data Exchange Format Representing Head-Related Transfer Functions.” In Proceedings of the 134th Convention of the Audio Engineering Society (AES), Paper 8880. Roma, Italy, 2013.

[5] AES69-2022. "AES standard for file exchange - Spatial acoustic data file format." Standard of the Audio Engineering Society. https://www.aes.org/publications/standards/search.cfm?docID=99

Extended Capabilities

C/C++ Code Generation
Generate C and C++ code using MATLAB® Coder™.

Version History

Introduced in R2018b

R2024a: Visualize frequency response of HRTF

Call interpolateHRTF with no output arguments to visualize the frequency response of data for desired azimuth and elevation values. This syntax works with the SimpleFreeFieldHRIR or SimpleFreeFieldHRTF SOFA conventions.

R2023b: Interpolate HRTFs from data stored in SOFA files

Interpolate an HRTF from a SOFA file containing HRTF data. Read in a SOFA file using sofaread and pass the resulting SOFA object to interpolateHRTF.

R2023b: Interpolate HRTFs using the nearest neighbor algorithm

Use nearest neighbor interpolation by setting the Algorithm name-value argument to "nearest". Nearest neighbor interpolation is faster than the 3-D bilinear and vector-based amplitude panning (VBAP) algorithms at the cost of being less accurate.