Main Content

Read, Analyze and Process SOFA Files

SOFA (Spatially Oriented Format for Acoustics) [1] is a file format for storing spatially oriented acoustic data like head-related transfer functions (HRTF) and binaural or spatial room impulse responses. SOFA has been standardized by the Audio Engineering Society (AES) as AES69-2015.

In this example, you load a SOFA file containing HRTF measurements for a single subject in MATLAB. You then analyze the HRTF measurements in the time domain and the frequency domain. Finally, you use the HRTF impulse responses to spatialize an audio signal in real time by modeling a moving source based on desired azimuth and elevation values.

Explore a SOFA File in MATLAB

You use a SOFA file from the SADIE II database [2]. The file corresponds to spatially discrete free-field in-the-ear HRTF measurements for a single subject. The measurements characterize how each ear receives a sound from a point in space.

Download the SOFA file.

downloadFolder = matlab.internal.examples.downloadSupportFile("audio","SOFA/SOFA.zip");
dataFolder = tempdir;
unzip(downloadFolder,dataFolder)
netFolder = fullfile(dataFolder,"SOFA");
addpath(netFolder)
filename = "H10_48K_24bit_256tap_FIR_SOFA.sofa";

Display SOFA File Contents

SOFA files consist of binary data stored in the netCDF-4 format. You can use MATLAB to read and write netCDF files.

Display the contents of the SOFA file using ncdisp (execute ncdisp(filename)).

The file contents consist of multiple fields corresponding to different aspects of the measurements, such as the (fixed) listener position, the varying source position, the coordinate system used to capture the data, general metadata related to the measurement, as well as the measured impulse responses.

NetCDF is a "self-describing" file format, where data is stored along with attributes that can be used to assist in its interpretation. Consider the display snippet corresponding to the source position for example:

SourcePosition contains the coordinates for the varying source position used in the measurements (here, there are 2114 separate positions). The file also contains attributes (Type, Units) describing the coordinate system used to store the positions (here, spherical), as well as information about the dimensions of the data (C,M). The dimensions are defined in the AES69 standard [3]:

For the file in this example:

  • M = 2114 (the total number of measurements, each corresponding to a unique source position).

  • R = 2 (corresponding to the two ears).

  • E = 1 (one emitter or sound source per measurement).

  • N = 256 (the length of each recorded impulse response).

Read SOFA File Information

Use ncinfo to get information about the SOFA file.

SOFAInfo = ncinfo(filename);

The fields of the structure SOFAInfo hold information related to the file's dimensions, variables and attributes.

Read a Variable from the SOFA File

Use ncread to read the value of a variable in the SOFA file.

Read the variable corresponding to the measured impulse responses.

ir = ncread(filename,"Data.IR");
size(ir)
ans = 1×3

         256           2        2114

This variable holds impulse responses for the left and right ear for 2114 independent measurements. Each impulse response is of length 256.

Read the sampling rate of the measurements.

fs = ncread(filename,"Data.SamplingRate")
fs = 48000

Plot the first measured impulse response.

figure;
t = (0:size(ir,1)-1)/fs;
plot(t,ir(:,1,1))
grid on
xlabel("Time (s)")
ylabel("Impulse response")

Figure contains an axes object. The axes object with xlabel Time (s), ylabel Impulse response contains an object of type line.

Read SOFA Files with sofaread

It is possible to read and analyze the contents of the SOFA file using a combination of ncinfo and ncread. However, the process can be cumbersome and time consuming.

The function sofaread allows you to read the entire contents of a SOFA file in one line of code.

Use sofaread to get the contents of the SOFA file.

s = sofaread(filename)
s = 
  audio.sofa.SimpleFreeFieldHRIR handle with properties:

            Numerator: [2114×2×256 double]
                Delay: [0 0]
         SamplingRate: 48000
    SamplingRateUnits: "hertz"

  Show all properties

Click "Show all properties" in the display above to see the rest of the properties from the SOFA file.

You can easily access the measured impulses responses.

ir = s.Numerator;
size(ir)
ans = 1×3

        2114           2         256

You access other variables in a similar fashion. For example, read the source positions along with the coordinate system used to express them:

srcPositions = s.SourcePosition;
size(srcPositions)
ans = 1×2

        2114           3

srcPositions(1,:)
ans = 1×3

         0  -90.0000    1.2000

s.SourcePositionType
ans = 
'spherical'
s.SourcePositionUnits
ans = 
"degree, degree, meter"

View the Measurement Geometry

Call plotGeometry to view the general geometric setup of the measurements.

figure;
plotGeometry(s)

Figure contains an axes object. The axes object with xlabel X (meters), ylabel Y (meters) contains 2 objects of type line. One or more of the lines displays its values using only markers These objects represent Source position, Receiver position.

Alternatively, specify input indices to restrict the plot to desired source locations. For example, plot the source positions located in the median plane.

figure;
idx = findMeasurements(s,Plane="median");
plotGeometry(s,MeasurementIndex=idx)

Figure contains an axes object. The axes object with xlabel X (meters), ylabel Y (meters) contains 2 objects of type line. One or more of the lines displays its values using only markers These objects represent Source position, Receiver position.

Plot Measurement Impulse Responses

The file in this example uses the SimpleFreeFieldHRIR convention, which stores impulse response measurements in the time domain as FIR filters.

s.SOFAConventions
ans = 
"SimpleFreeFieldHRIR"
s.DataType
ans = 
"FIR"

Plot the impulse response of the first 3 measurements corresponding to the second receiver using impz.

figure;
impz(s, MeasurementIndex=1:3,Receiver=2)

Figure contains an axes object. The axes object with xlabel Time (s) contains 3 objects of type line. These objects represent Measurement 1. Receiver 2, Measurement 2. Receiver 2, Measurement 3. Receiver 2.

Compute HRTF Frequency Responses

It is straightforward to compute and plot the frequency response of the impulse responses using freqz.

Compute the frequency response of the first measurement for both ears. Use a frequency response length of 512.

[H,F] = freqz(s, Receiver=1:2, NPoints=512);

H is the complex frequency response array. F is the vector of corresponding frequencies (in Hertz).

size(H)
ans = 1×3

   512     1     2

size(F)
ans = 1×2

   512     1

You can also use freqz to plot the frequency response by calling the function with no output arguments.

Plot the frequency response of the first 3 measurements in the horizontal plane at zero elevation for the first receiver.

figure;
idx = findMeasurements(s, Plane="horizontal");
freqz(s, MeasurementIndex=idx(1:3), Receiver=1)

Figure contains an axes object. The axes object with xlabel Frequency (Hz), ylabel Magnitude (dB) contains 3 objects of type line. These objects represent Measurement 1097. Receiver 1, Measurement 1122. Receiver 1, Measurement 1185. Receiver 1.

Compute HRTF Spectra

It is often useful to compute and visualize the magnitude spectra of HRTF data in a specific plane in space.

For example, compute the spectrum in the horizontal plane (corresponding to an elevation angle equal to zero) for the first receiver. Use a frequency response length of 2048.

[S,F,azi] = spectrum(s, Plane="horizontal",Receiver=1, NPoints=2048); %#ok
  • S is the horizontal plane spectrum of size 2048-by-L, where L is the number of measurements in the specified plane.

  • F is the frequency vector of length 2048.

  • azi is a vector of length L containing the azimuth angles (in degrees) corresponding to the horizontal plane measurements.

You can plot the spectrum instead of returning it.

figure
spectrum(s, Plane="horizontal",Receiver=1, NPoints=2048);

Figure contains an axes object. The axes object with title Horizontal Plane Spectrum. Receiver 1., xlabel Frequency (Hz), ylabel Azimuth (degrees) contains an object of type surface.

Compute Energy-Time Curve

It is often useful to visualize the decay of the HRTF responses over time using an energy-time curve (ETC).

Measure the ETC in the horizontal plane.

ETC = energyTimeCurve(s);
size(ETC)
ans = 1×2

   256    96

The first dimension of ETC is equal to the length of the impulse response. The second dimension of ETC is equal to the number of points in the specified plane.

Plot the ETC.

figure;
energyTimeCurve(s)

Figure contains an axes object. The axes object with title Horizontal Plane Energy-Time Curve. Receiver 1., xlabel Time (s), ylabel Azimuth (degrees) contains an object of type surface.

Compute Interaural Time Difference

Interaural time difference (ITD) is the difference in arrival time of a sound between two ears. It is an important binaural cue for sound source localization.

Compute the ITD for all measurements.

ITD =  interauralTimeDifference(s, Specification="measurements");

Plot the ITD for all measurements.

figure;
interauralTimeDifference(s, Specification="measurements");

Figure contains an axes object. The axes object with title Interaural Time Difference, xlabel X (meters), ylabel Y (meters) contains an object of type patch.

You can also compute and plot the ITD in a specified plane. For example, plot the ITD in the horizontal plane at an elevation of 35 degrees.

figure
interauralTimeDifference(s, Plane="horizontal",PlaneOffsetAngle=35)

Figure contains an axes object. The axes object with title Horizontal Plane Interaural Time Difference, xlabel Azimuth (degrees), ylabel ITD (s) contains an object of type line.

Compute Interaural Level Difference

Interaural level difference (ILD) measures the difference, in decibels, of the received signals at the two ears in the desired plane.

Compute the ILD in the horizontal plane.

[ILD,F] = interauralLevelDifference(s,NPoints=1024);
size(ILD)
ans = 1×2

        1024          96

ILD is 1024-by-L, where L is the number of points in the desired plane.

Plot the ILD in the horizontal plane.

figure
interauralLevelDifference(s)

Figure contains an axes object. The axes object with title Horizontal Plane Interaural Level Difference, xlabel Frequency (Hz), ylabel Azimuth (degrees) contains an object of type surface.

Compute Frequency Directivity

HRTF measurements are highly frequency selective. It is often useful to compare the received power level of different frequencies in space.

Call directivity to plot the received level at 200 Hz for all measurements.

figure;
directivity(s, 200, Specification="measurements")

Figure contains an axes object. The axes object with title Frequency directivity at 200 Hz. Receiver 1., xlabel X (meters), ylabel Y (meters) contains an object of type patch.

Generate a similar plot at 4000 Hz.

figure;
CurrentModifier does not exist
directivity(s, 4000, Specification="measurements")

Figure contains an axes object. The axes object with title Frequency directivity at 4000 Hz. Receiver 1., xlabel X (meters), ylabel Y (meters) contains an object of type patch.

You can also compute and plot the directivity for selected planes. For example, plot the directivity at 700 Hz and 1200 Hz for the horizontal plane at an elevation of 35 degrees.

figure;
CurrentModifier does not exist
directivity(s, [700 1200], PlaneOffsetAngle=35, Receiver=1)

Figure contains an axes object with type polaraxes. The polaraxes object contains 2 objects of type line. These objects represent 700.000000 Hz, 1200.000000 Hz.

Interpolate HRTF Measurements

The HRTF measurements in the SOFA file correspond to a finite number of azimuth/elevation angle combinations. It is possible to interpolate the data to any desired spatial location using 3-D HRTF interpolation with interpolateHRTF.

Specify the desired source position (in degrees).

desiredAz = [-120;-60;0;60;120;0;-120;120];
CurrentModifier does not exist
desiredEl = [-90;90;45;0;-45;0;45;45];
desiredPosition = [desiredAz, desiredEl];

Calculate the head-related impulse response (HRIR) using the vector base amplitude panning interpolation (VBAP) algorithm at a desired source position.

interpolatedIR = interpolateHRTF(s, desiredPosition);

You can also visualize the interpolated frequency responses.

figure;
interpolateHRTF(s, desiredPosition);

Figure contains an axes object. The axes object with xlabel Frequency (Hz), ylabel Magnitude (dB) contains 16 objects of type line. These objects represent Measurement 1. Receiver 1, Measurement 1. Receiver 2, Measurement 2. Receiver 1, Measurement 2. Receiver 2, Measurement 3. Receiver 1, Measurement 3. Receiver 2, Measurement 4. Receiver 1, Measurement 4. Receiver 2, Measurement 5. Receiver 1, Measurement 5. Receiver 2, Measurement 6. Receiver 1, Measurement 6. Receiver 2, Measurement 7. Receiver 1, Measurement 7. Receiver 2, Measurement 8. Receiver 1, Measurement 8. Receiver 2.

Model Moving Source Using HRIR Filtering

Filter a mono input through the interpolated impulse responses to model a moving source.

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

desiredFs = s.SamplingRate;
CurrentModifier does not exist
[x,fs] = audioread("Counting-16-44p1-mono-15secs.wav");
y = audioresample(x,InputRate=fs,OutputRate=desiredFs);
y = y./max(abs(y));
audiowrite("Counting-16-48-mono-15secs.wav",y,desiredFs);

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 the 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)
release(spatialFilter);

Spatialize Audio in Real Time

Simulate a sound source moving in the horizontal plane, with an initial azimuth of -90 degrees. Gradually increase the azimuth as the simulation is running.

Compute the starting impulse responses based on the initial source position.

index = 1;
loc = [-90 0];

interpolatedIR = interpolateHRTF(s, ...
    loc);
spatialFilter.Numerator = squeeze(interpolatedIR);

Execute the simulation loop. Shift the source elevation by 30 degrees every 100 loop iterations. Use interpolateHRTF to estimate the new desired impulse responses.

while ~isDone(fileReader)
    index=index+1;
    frame = fileReader();
    frame = frame(:,1);
    
    audioOut = spatialFilter(frame);

    deviceWriter(audioOut);

    if mod(index,100)==0
        loc(1)=loc(1)+30;
        interpolatedIR = interpolateHRTF(s, ...
            loc);
        spatialFilter.Numerator = squeeze(interpolatedIR);
    end
end
release(deviceWriter)
release(fileReader)
release(spatialFilter);

References

[1] Majdak, P., Zotter, F., Brinkmann, F., De Muynke, J., Mihocic, M., and Noisternig, M. (2022). “Spatially Oriented Format for Acoustics 21: Introduction and Recent Advances,” J Audio Eng Soc 70, 565–584. DOI: 10.17743/jaes.2022.0026.

[2] https://www.york.ac.uk/sadie-project/database.html

[3] Majdak, P., De Muynke, J., Zotter, F., Brinkmann, F., Mihocic, M., and Ackermann, D. (2022). AES standard for file exchange - Spatial acoustic data file format (AES69-2022), Standard of the Audio Engineering Society. https://www.aes.org/publications/standards/search.cfm?docID=99

[4] Andreopoulou A, Katz BFG. Identification of perceptually relevant methods of inter-aural time difference estimation. J Acoust Soc Am. 2017 Aug;142(2):588. doi: 10.1121/1.4996457. PMID: 28863557.

See Also

| |

Related Topics