Réponse de Beamformer par python

Dessin de motif de faisceau de retard et de somme

C'est un réaménagement des roues, mais je ne veux pas l'utiliser à Don Pisha, alors je l'ai écrit. La partie du vecteur manifold du tableau est suspecte, mais j'ai pu dessiner une figure comme celle-là.

Malgré tout, j'ai abandonné car je ne pouvais pas ajuster le rapport hauteur / largeur dans le tracé 3D de plot_surface.

ax.get_Cela change avec proj, mais la zone de dessin ne change pas, donc je ne suis pas content.

 Cela n'a pas beaucoup de sens, mais la norme L2 est définie par 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

Réseau linéaire et régulièrement espacé. Alignez-vous sur l'axe 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 #Position du microphone
    if Is3D:
        return np.array([PX, np.zeros(len(PX)), np.zeros(len(PX))])
    else:
        return np.array([PX, np.zeros(len(PX))])

Calcul d'un vecteur de variété de tableau. Je pense que les ondes planes sont correctes, mais je ne suis pas sûr des ondes sphériques.

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

Calculez la position de la direction de la source sonore. Changement juste des coordonnées polaires du vecteur unitaire aux coordonnées orthogonales. Depuis la direction, pour générer l'argument de 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

Générez un filtre de délai et de somme à partir d'un vecteur de variété de tableau.

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]

Dessinez un motif de faisceau pour le filtre du formateur de faisceaux. Il doit être utilisable pour MVDR et autres profileurs de poutre.

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()

Génération de motifs de retard et de somme de faisceaux dédiés aux formateurs de faisceau. On suppose que la directionnalité est orientée dans la direction de l'axe Y (90 degrés).

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()

Retard et somme Tracez le motif de faisceau (2D) du formateur de faisceau.

plot_ds_freq_resp()

Délai et somme Tracez le motif de faisceau (3D) du formateur de faisceau.

plot_ds_freq_resp(plt3D=True)

Tracez en modifiant l'espacement du microphone. On peut voir que l'épaisseur du lobe principal a changé et qu'une distorsion de pliage s'est produite.

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))

Tracez le motif du faisceau en modifiant la direction de la directionnalité.

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)