In [1]:
import h5py
import numpy as np
from ear.core import hoa
import matplotlib.pyplot as plt
%matplotlib inline
In [2]:
plt.rcParams["figure.figsize"] = (15, 8)
In [3]:
f = h5py.File("HRIR_L2702_NF150.sofa")
In [4]:
def fftcorrelate(a, b, axis):
assert a.shape == b.shape
a = np.moveaxis(a, axis, 0)
b = np.moveaxis(b, axis, 0)
af = np.fft.rfft(a[::-1], len(a) * 2, axis=0)
bf = np.fft.rfft(b, len(a) * 2, axis=0)
cf = af * bf
c = np.fft.irfft(cf, axis=0)
return np.moveaxis(c, 0, axis)
def fftdelay(a, b, axis):
cc = fftcorrelate(a, b, axis)
return np.argmax(cc, axis) - a.shape[axis] - 1
In [5]:
itd_samples = fftdelay(f["Data.IR"][:, 0], f["Data.IR"][:, 1], 1)
In [6]:
az, el = f["SourcePosition"][:, 0], f["SourcePosition"][:, 1]
In [7]:
cs = plt.tricontour(az, el, itd_samples, 20)
plt.clabel(cs, inline=1, fontsize=10, fmt="%d")
plt.xlabel("azimuth (deg)")
plt.xlabel("elevation (deg)")
plt.axis("equal");
In [12]:
n, m = hoa.from_acn(np.arange(13**2))
Y = hoa.sph_harm(n[np.newaxis, :], m[np.newaxis, :],
np.radians(az)[:, np.newaxis], np.radians(el)[:, np.newaxis],
norm=hoa.norm_N3D)
v = np.dot(itd_samples, np.linalg.pinv(Y.T))
In [13]:
saz, sel = np.mgrid[0:360, -90:91]
In [14]:
sY = hoa.sph_harm(n[np.newaxis, :], m[np.newaxis, :],
np.radians(saz)[..., np.newaxis], np.radians(sel)[..., np.newaxis],
norm=hoa.norm_N3D)
sitd_samples = np.dot(sY, v)
In [15]:
cs = plt.contour(saz, sel, sitd_samples, np.linspace(-40, 40, 21))
plt.clabel(cs, inline=1, fontsize=10, fmt="%d")
plt.colorbar()
plt.axis("equal")
Out[15]:
