Frequency data phase unwrapping issue

[stealthBanana]

Hi,

I'm writing to you because I encountered a problem trying to extrapolate the phase angle of the frequency data provided by MEsh2HRTF, I'm trying to plot the unwrapped phase at each azimuth of the evaluation grid, however, in some cases, when unwrapping the phase the cyclicity of the signal is lost, meaning that the phase calculated at 0 and 360 degrees differ to a large amount. This issue appears when I plot the values at 0 deg elevation at a 3500Hz frequency.

The problem can be reproduced with the following code:

import pyfar as pf
import numpy as np
import matplotlib.pyplot as plt

signal, sources, _ = pf.io.read_sofa('./tutorial.sofa')

# get angles for 0 deg elevation
elevation = 0
_, elevation_mask = sources.find_slice(coordinate="elevation", unit="deg", value=elevation, tol=0.1)
azimuth = sources.get_sph('top_elev', 'deg')[elevation_mask, elevation]

# get id for 3500Hz frequency
freq_id = signal.find_nearest_frequency(value=3500)

# get frequency data
channel_id = 0
frequency_data = signal[elevation_mask, channel_id][np.argsort(np.argsort(azimuth))].freq[:,freq_id].flatten()

# calculate phase angle
phase = np.arctan2(frequency_data.imag, frequency_data.real)

# unwrap phase
unwrapped_phase = np.unwrap(phase)

# plot phase
fig, ax = plt.subplots(nrows=1, ncols=2, sharey=True, figsize=(14, 5))
ax[0].plot(azimuth, phase, color='g', label='Wrapped phase')
ax[1].plot(azimuth, unwrapped_phase, color='r', label='Unwrapped phase')

ax[0].set_ylabel('Phase [rad]')
ax[0].set_xlabel('Azimuth[deg]')
ax[1].set_xlabel('Azimuth[deg]')
fig.legend()

plt.show()

This is the resulting plot that I'm getting, the expected result would be that the phase at 0 and 360 degrees are the same value:

The SOFA file has been generated using the Mesh2HRTF tutorial head mesh and the export settings are:

• BEM method: ML-FMM BEM
• Source Type: Point source
• Pictures: False
• Reference: True
• Compute HRIRs: True
• c: 343.18
• Evaluation Grids name: Rings_1_0m
• Min. frequency: 3500
• Max. frequency: 3500

The evaluation grid is a custom grid created with the write_evaluation_grid.py script, it has a node for each azimuth degree at 5 different elevations, -30, -15, 0, 15 and 30 degrees with 1m radius:

The head mesh has been remeshed with the hrtf_mesh_grading algorithm and the point source has been placed 1cm off the ear.

Here I have attached the SOFA, blender file and evaluation grid. tutorial.zip

[f-brinkmann]

Hi, phase unwrapping is tricky. Without having thought about the details it might be that you ran into the problem described by Johannes Zaar who proposed spatial unwrapping. For what reason are you looking at the phase?