pythonによるビームフォーマの応答

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)をプロット。
python
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)