Beamformer response by python
Delay and Sum beam pattern drawing
It's a redevelopment of the wheel, but I don't want to use it in Don Pisha, so I wrote it. The part of the array manifold vector is suspicious, but I was able to draw a figure like that.
Even so, I gave up because I could not adjust the aspect ratio in the 3D plot of plot_surface.
ax.get_It changes with proj, but the drawing area does not change, so I'm not happy.
It doesn't make much sense, but the L2 norm is defined by numba.
```python
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import numba
@numba.vectorize([numba.float64(numba.complex128),numba.float32(numba.complex64)])
def abs2(x):
return x.real**2 + x.imag**2
Linear, evenly spaced array. Line up on the X axis.
def linear_mic_array(
num = 4, # Number of array elements
spacing = 0.2, # Element separation in metres
Is3D = True,
):
PX = np.arange(num) * spacing #Microphone position
if Is3D:
return np.array([PX, np.zeros(len(PX)), np.zeros(len(PX))])
else:
return np.array([PX, np.zeros(len(PX))])
Calculation of array manifold vectors. I think plane waves are okay, but I'm not sure about spherical waves.
def calc_array_manifold_vector_from_point(mic_layout, source_from_mic, nfft=512, sampling_rate = 16000,
c=343, attn=False, ff=True):
p = np.array(source_from_mic)
if p.ndim == 1:
p = source[:, np.newaxis]
frequencies = np.linspace(0, sampling_rate // 2, nfft // 2 + 1)
omega = 2 * np.pi * frequencies
if ff:
p /= np.linalg.norm(p)
D = -1 * np.einsum("im,i->m", mic_layout, p.squeeze())
D -= D.min()
else:
D = np.sqrt(((mic_layout - source_from_mic) ** 2).sum(0))
phase = np.exp(np.einsum("k,m->km", -1j * omega / c, D))
if attn:
return 1.0 / (4 * np.pi) / D * phase
else:
return phase
Calculate the position of the direction of the sound source. Just changing the polar coordinates of the unit vector to orthogonal coordinates. From the direction, for argument generation of source_from_mic.
def get_look_direction_loc(deg_phi, deg_theta):
phi = np.deg2rad(deg_phi)
theta = np.deg2rad(deg_theta)
return np.array([[
np.cos(phi) * np.cos(theta),
np.cos(phi) * np.sin(theta),
np.sin(phi),
]]).T
Generate a Delay and Sum filter from an array manifold vector.
def calc_delay_and_sum_weights(mic_layout, source_from_mic, nfft=512, sampling_rate = 16000,
c=343, attn=False, ff=True):
a = calc_array_manifold_vector_from_point(mic_layout, source_from_mic, nfft, sampling_rate, c, attn, ff)
return a / np.einsum("km,km->k", a.conj(), a)[:, None]
Draw a beam pattern for the beamformer filter. It should be usable for MVDR and other beam formers.
def plot_freq_resp(
w,
nfft = 512, # Number of fft
angle_resolution = 500, # Number of angle points to calculate
mic_layout = linear_mic_array(),
sampling_rate = 16000, # Hz
speedSound = 343.0, # m/s
plt3D = False,
vmin = -50,
vmax = None,
):
n_freq, n_ch = w.shape
if n_freq != nfft // 2 + 1:
raise ValueError("invalid n_freq")
if n_ch != mic_layout.shape[1]:
raise ValueError("invalid n_ch")
freqs = np.linspace(0, sampling_rate // 2, nfft // 2 + 1)
angles = np.linspace(0, 180, angle_resolution)
angleRads = np.deg2rad(angles)
loc = np.array([[np.cos(theta), np.sin(theta), 0] for theta in angleRads]).T
phases = np.einsum('k,ai,am->kim', 2j * np.pi * freqs / speedSound, loc, mic_layout)
psi = np.einsum('km,kim->ki', w.conj(), np.exp(phases))
outputs = np.sqrt(abs2(psi))
logOutputs = 20 * np.log10(outputs)
if vmin is not None:
logOutputs = np.maximum(logOutputs, vmin)
if vmax is not None:
logOutputs = np.minimum(logOutputs, vmax)
plt.figure(figsize=(10,8))
if plt3D:
ax = plt.subplot(projection='3d')
surf = ax.plot_surface(*np.meshgrid(angles, freqs), logOutputs, rstride=1, cstride=1, cmap='hsv', linewidth=0.3)
#ax.get_proj = lambda: np.dot(Axes3D.get_proj(ax), np.diag([1.5, 1.5, 0.3, 1]))
ax.view_init(elev=30, azim=-45)
else:
ax = plt.subplot()
surf = ax.pcolormesh(*np.meshgrid(angles, freqs), logOutputs, cmap="inferno", shading='gouraud')
plt.colorbar(surf)
plt.xlabel("Arrival Angle[degrees]")
plt.ylabel("Frequency[Hz]")
plt.title("Frequency Response")
plt.show()
Delay and Sum Beam pattern generation dedicated to beamformers. It is assumed that the directivity is oriented in the Y-axis direction (90 degrees).
def plot_ds_freq_resp(
nfft = 512, # Number of fft
angle_resolution = 500, # Number of angle points to calculate
mic_layout = linear_mic_array(),
sampling_rate = 16000, # Hz
speedSound = 343.0, # m/s
plt3D = False,
):
freqs = np.linspace(0, sampling_rate // 2, nfft // 2 + 1)
angles = np.linspace(0, 180, angle_resolution)
angleRads = np.deg2rad(angles)
loc = np.array([[np.cos(theta), np.sin(theta), 0] for theta in angleRads]).T
phases = np.einsum('k,ai,am->kim', 2j * np.pi * freqs / speedSound, loc, mic_layout)
outputs = np.sqrt(abs2(np.exp(phases).sum(-1))) / mic_layout.shape[1]
logOutputs = np.maximum(20 * np.log10(outputs), -50)
plt.figure(figsize=(10,8))
if plt3D:
ax = plt.subplot(projection='3d')
surf = ax.plot_surface(*np.meshgrid(angles, freqs), logOutputs, rstride=1, cstride=1, cmap='hsv', linewidth=0.3)
#ax.get_proj = lambda: np.dot(Axes3D.get_proj(ax), np.diag([1.5, 1.5, 0.3, 1]))
ax.view_init(elev=30, azim=-45)
else:
ax = plt.subplot()
surf = ax.pcolormesh(*np.meshgrid(angles, freqs), logOutputs, cmap="inferno", shading='gouraud')
plt.colorbar(surf)
plt.xlabel("Arrival Angle[degrees]")
plt.ylabel("Frequency[Hz]")
plt.title("Frequency Response")
plt.show()
Plot the beam pattern (2D) of the Delay and Sum beamformer.
plot_ds_freq_resp()
Plot the beam pattern (3D) of the Delay and Sum beamformer.
plot_ds_freq_resp(plt3D=True)
Plot by changing the microphone spacing. It can be seen that the thickness of the main robe has changed and that folding distortion has occurred.
plot_ds_freq_resp(plt3D=True, mic_layout=linear_mic_array(num=8, spacing=0.06))
plot_ds_freq_resp(plt3D=True, mic_layout=linear_mic_array(num=8, spacing=0.03))
plot_ds_freq_resp(plt3D=True, mic_layout=linear_mic_array(num=8, spacing=0.015))
Plot the beam pattern by changing the direction of directivity.
mic_layout = linear_mic_array()
source_from_mic = get_look_direction_loc(0, 120)
Wds = calc_delay_and_sum_weights(mic_layout, source_from_mic)
plot_freq_resp(Wds, mic_layout=mic_layout)
Recommended Posts