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

本文主要是介绍波束图（beam pattern）的python和matlab实现，希望对大家解决编程问题提供一定的参考价值，需要的开发者们随着小编来一起学习吧！

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

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 differenceWc = 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 numbadef abs2(x):return x.real**2 + x.imag**2def linear_mic_array(num = 4,       # number of array elementsspacing = 0.2, # element separation in metresIs3D = True,
):PX = np.arange(num) * spacing   # microphone positionif 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 * frequenciesif 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*phaseelse:return phasedef 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),]]).Tdef 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.shapeif 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]).Tphases = 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 fftangle_resolution = 500,    # Number of angle points to calculatemic_layout       = linear_mic_array(),sampling_rate    = 16000,  # HzspeedSound       = 343.0,  # m/splt3D            = 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]).Tphases = 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)