Es ist eine Neuentwicklung der Räder, aber ich möchte es nicht in Don Pisha verwenden, also habe ich es geschrieben. Der Teil des Array-Mannigfaltigkeitsvektors ist verdächtig, aber ich konnte eine solche Figur zeichnen.
Trotzdem gab ich auf, weil ich das Seitenverhältnis im 3D-Plot von plot_surface nicht anpassen konnte.
ax.get_Es ändert sich mit proj, aber der Zeichenbereich ändert sich nicht, also bin ich nicht glücklich.
Es macht nicht viel Sinn, aber die L2-Norm wird durch numba definiert.
```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
Lineares, gleichmäßig verteiltes Array. Richten Sie sich auf der X-Achse aus.
def linear_mic_array(
num = 4, # Number of array elements
spacing = 0.2, # Element separation in metres
Is3D = True,
):
PX = np.arange(num) * spacing #Mikrofonposition
if Is3D:
return np.array([PX, np.zeros(len(PX)), np.zeros(len(PX))])
else:
return np.array([PX, np.zeros(len(PX))])
Berechnung eines Array-Mannigfaltigkeitsvektors. Ich denke, ebene Wellen sind in Ordnung, aber ich bin mir bei sphärischen Wellen nicht sicher.
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
Berechnen Sie die Position der Schallquellenrichtung. Einfach von den Polarkoordinaten des Einheitsvektors zu den orthogonalen Koordinaten wechseln. Aus der Richtung, um das Argument von source_from_mic zu generieren.
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
Generieren Sie einen Verzögerungs- und Summenfilter aus einem Array-Mannigfaltigkeitsvektor.
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]
Zeichnen Sie ein Strahlmuster für den Strahlformerfilter. Es sollte für MVDR und andere Strahlformer verwendbar sein.
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()
Erzeugung von Verzögerungs- und Summenstrahlmustern für Strahlformer. Es wird angenommen, dass die Richtwirkung in Richtung der Y-Achse (90 Grad) ausgerichtet ist.
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()
Verzögerung und Summe Zeichnen Sie das Strahlmuster (2D) des Strahlformers.
plot_ds_freq_resp()
Verzögerung und Summe Zeichnen Sie das Strahlmuster (3D) des Strahlformers.
plot_ds_freq_resp(plt3D=True)
Zeichnen Sie durch Ändern des Mikrofonabstands. Es ist ersichtlich, dass sich die Dicke des Hauptlappens geändert hat und dass eine Faltverzerrung aufgetreten ist.
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))
Zeichnen Sie das Strahlmuster, indem Sie die Richtung der Richtung ändern.
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)
Recommended Posts