波束图（beam pattern）的python和matlab实现

【代码】波束图（beam pattern）的python和matlab实现。

chinabing

8929人浏览 · 2022-08-21 22:30:39

chinabing · 2022-08-21 22:30:39 发布

关注、点赞、收藏是对我最大的支持，谢谢！

1、什么是波束图

2、波束图的原理

3、波束图的实现

1、什么是波束图

通过波束图可以知晓哪个方向的信号被增强，哪个方向的信号被抑制。

2、波束图的原理

声源到各麦克风的时间是不一样的，存在时间差，以mic1为参考点，mic2和micM均会提前，提前的时间为 $(\Delta t, 2\Delta t,,...,(M-1)\Delta t,)$ ,其中 $\Delta=dcos\theta/c$ 。

假设声波波长 $\lambda=d/2$ ，频率为 $f$ ，相位差为 $(\psi, 2\psi,...,(M-1)\psi)$ ，其中

$\psi = 2\pi \Delta t/T= 2\pi f \Delta t= 2\pi f dcos\theta/c =2\pi dcos\theta/\lambda = 4 \pi cos\theta$

设定期望阵列流形矢量为

$w=(e^{j\psi }, e^{j2\psi },...,e^{j(M-1)\psi })$

其它方向阵列流形矢量 $p(\theta) = (e^{j\psi^{'}}, e^{j2\psi^{'}},...,e^{j(M-1)\psi^{'}})$ ，各方向的波束响应 $B(\theta)$ 可以用波束图来描述，

$B(\theta) = w^Tp(\theta)$

3、波束图的实现

clear; close all; clc;

c=340;
f=1000;
lambda=c/f;  % wave length
k=2*pi/lambda;
d=lambda/2; % mic distance 
M=10;
theta_d = 80*pi/180;  % angle of incidence
theta_angle = 0:0.1:180;
theta = theta_angle*pi/180;

psi = 2*pi*d*cos(theta_d)/lambda;  % phase difference

Wc = exp(-1i*psi*[0:M-1])/M;
B = zeros(size(theta));

%% compute beam pattern
for i=1:length(theta)
    psi2 = 2*pi*d*cos(theta(i))/lambda; 
    p = exp(1i * psi2 * [0:M-1]);
    B(i) = conj(Wc)*p';
end
B_db = 20*log10(B);
limit_dB = -50;
index = B_db <limit_dB;
B_db(index) = limit_dB;

figure;
plot(theta_angle, B_db, 'linewidth', 1.5);
grid on;
xlabel('azimuth(^°)'); ylabel('20lg(B)/dB');
title('beam pattern');

波束图


import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import numba

def abs2(x):
    return x.real**2 + x.imag**2

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

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_from_mic[:, 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

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

def cal_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]

def plot_freq_resp(
        w,
        nfft   =  512,
        angle_resolution = 500,
        mic_layout = linear_mic_array(),
        sampling_rate = 16000,
        speedSound = 434,
        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()

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

if __name__ == '__main__':
    plot_ds_freq_resp()
    plot_ds_freq_resp(plt3D=True)