#Delay and Sumのビームパターンの描画
車輪の再開発だけど、ドンピシャで使いたいのがないので、書いてみた。
アレイ・マニフォールド・ベクトル(array manifold vector)の部分が怪しいけど、それっぽい図形が描けた。
それにしても、plot_surface の3Dプロットで、アスペクト比が調整できなくて断念。
ax.get_projで変わるけど、描画領域が変わらず、うれしくない。
あまり意味はないけどL2ノルムをnumbaで定義。
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
直線状等間隔アレイ。X軸上に並ぶ。
def linear_mic_array(
num = 4, # Number of array elements
spacing = 0.2, # Element separation in metres
Is3D = True,
):
PX = np.arange(num) * spacing # マイク位置
if Is3D:
return np.array([PX, np.zeros(len(PX)), np.zeros(len(PX))])
else:
return np.array([PX, np.zeros(len(PX))])
アレイ・マニフォールド・ベクトル(array manifold vector)の計算。
平面波は大丈夫だと思うけど、球面波はよくわからん。
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
音源の方向の位置を計算。単位ベクトルの極座標から直行座標に変えているだけ。
向きから、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
アレイ・マニフォールド・ベクトル(array manifold vector)からDelay and Sumのフィルタを生成。
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]
ビームフォーマのフィルタに対するビームパターンを描画。MVDRとか他のビームフォーマにも使えるはず。
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ビームフォーマ専用のビームパターン生成。Y軸方向(90度)に指向性が向いている想定。
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()
Delay and Sumビームフォーマのビームパターン(2D)をプロット。
plot_ds_freq_resp()
Delay and Sumビームフォーマのビームパターン(3D)をプロット。
plot_ds_freq_resp(plt3D=True)
マイク間隔を変更してプロット。メインローブの太さが変わっていたり、折り返し歪が発生しているのがわかる。
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))
指向性の向きを変えてビームパターンをプロット。
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)
