Über höhere Ambisonics

Ambisonics ist eine Methode, die eine dreidimensionale Schallwiedergabe realisiert, indem die Ankunftsrichtung von Schallwellen, die auf einen bestimmten Schallempfangspunkt einfallen, dh die Einfallsrichtung von Schallwellen, reproduziert wird. Die einfallenden Richtungseigenschaften werden reproduziert, indem die dreidimensionalen Richtungseigenschaften des an einem bestimmten Punkt einfallenden Schalls durch eine sphärische harmonische Funktion erweitert werden und die Ausgabe jedes Lautsprechers in der Lautsprecheranordnung bestimmt wird, um die entwickelten Richtungseigenschaften während der Wiedergabe zu reproduzieren. Diejenigen, die sich auf die 0. und 1. Ordnung konzentrieren, werden Ambisonics genannt, und diejenigen, die auch sphärische harmonische Funktionen 2. oder höherer Ordnung verwenden, werden Ambisonics höherer Ordnung genannt.

In diesem Artikel beschreibe ich den Code des Schallflusses mit einem Mikrofon, der virtuellen Anordnung der Lautsprecher und der Wiedergabe des Schallfelds.

Über die in der Simulation verwendeten Wellen

Die in diesem Artikel beschriebene Simulation behandelt nur ebene Wellen, es gibt jedoch drei Arten: ebene Wellen, sphärische Wellen und gerichtete sphärische Wellen. Die Wellengleichung für jedes kartesische Koordinatensystem kann wie folgt geschrieben werden: $e^{i\vec{k_i}\vec{r}}$ $\frac{e^{ik|r-r_s|}}{4\pi|r-r_s|}$ $\frac{e^{ik|r-r_s|}}{4\pi|r-r_s|}\Bigl(\alpha-\bigl(1-\alpha\bigr)\bigl(1+\frac{i}{k|r-r_s|}\bigr)\cos\gamma\Bigr)$ $ \ Gamma $ ist der Winkel zwischen der Empfangsposition und dem Schallquellenpositionsvektor, und $ \ alpha $ ist der Richtungskoeffizient ($ \ alpha $ = 1 entspricht der sphärischen Welle). Dann kann die Gleichung für jede Welle im Polarkoordinatensystem wie folgt geschrieben werden: $e^{i\vec{k_i}\vec{r}} = \sum_{n=0}^\infty\sum_{m=-n}^n4\pi i^n j_n(kr) Y_n^m(\theta,\varphi) Y_n^m(\theta_i,\varphi_i)^*$

\frac{e^{ik|r-r_s|}}{4\pi|r-r_s|} = \sum_{n=0}^\infty\sum_{m=-n}^n ik j_n(kr) h_n(kr_s) Y_n^m(\theta,\varphi) Y_n^m(\theta_s,\varphi_s)^*

\frac{e^{ik|r-r_s|}}{4\pi|r-r_s|}\Bigl(\alpha-\bigl(1-\alpha\bigr)\bigl(1+\frac{i}{k|r-r_s|}\bigr)\cos\gamma\Bigr) = \sum_{n=0}^\infty\sum_{m=-n}^n k j_n(kr)[i\alpha h_n(kr_s) +(1-\alpha)h_n^{'}(kr_s)]Y_n^m(\theta,\varphi) Y_n^m(\theta_s,\varphi_s)^*

$ j_n (kr) $ ist die erstklassige Kugelgefäßfunktion und $ h_n (kr) $ ist die erstklassige Hankelfunktion. $ ^ * $ Stellt eine komplexe Co-Konvertierung dar.

`green_function.py`


import numpy as np 
import function as fu


def Plane_wave(freq,cv, pos_s,vec_dn, pos_r):
    # freq: frequency
    # cv: sound speed
    # pos_s: source position (Not used in this function)
    # vec_dn: plane wave direction vector
    # pos_r: receiver position
    omega = 2*np.pi*freq
    k = omega/cv
    vec_dn_unit = np.array([0]*3,dtype = 'float')
    vec_dn_unit[0] = vec_dn[0] / np.linalg.norm(vec_dn, ord=2)
    vec_dn_unit[1] = vec_dn[1] / np.linalg.norm(vec_dn, ord=2)
    vec_dn_unit[2] = vec_dn[2] / np.linalg.norm(vec_dn, ord=2)
    azi,ele,_ = fu.cart2sph(-vec_dn_unit[0],-vec_dn_unit[1],-vec_dn_unit[2])
    ele_1 = np.pi / 2 - ele
    kx = k * np.cos(azi) * np.sin(ele_1)
    ky = k * np.sin(azi) * np.sin(ele_1)
    kz = k * np.cos(ele_1)
    x = pos_r[0] - pos_s[0]
    y = pos_r[1] - pos_s[1]
    z = pos_r[2] - pos_s[2]
    ret = np.exp(1j*(kx*x + ky*y + kz*z))
    return ret

def green(freq, cv, vec1, vec2):
    vec = np.array([0]*3,dtype = 'float')
    vec[0] = vec1[0] - vec2[0] 
    vec[1] = vec1[1] - vec2[1]
    vec[2] = vec1[2] - vec2[2]
    rr = np.linalg.norm(vec, ord=2)
    omega = 2 * np.pi * freq
    k = omega / cv
    jkr = 1j * k * rr
    ret = np.exp(jkr)/(4 * np.pi * rr)
    return ret

def green2(alpha, freq, cv, vecrs, vecr):
    # alpha:directivity coefficient 0~1 
    #                               0   :dipole 
    #                               0.5 :cardiod 
    #                               1   :monopole
    # CV: sound speed
    # freq: frequency [Hz](can be vector)
    # vecrs: sound source position
    # vecr: receiver position
    vecrd = np.array([0]*3,dtype = 'float')
    vecrd[0] = vecrs[0] - vecr[0]
    vecrd[1] = vecrs[1] - vecr[1]
    vecrd[2] = vecrs[2] - vecr[2]
    rd = np.linalg.norm(vecrd, ord=2)
    omega = 2 * np.pi * freq
    k = omega/cv
    cosgamma = np.dot(vecrd,-vecrs) / (np.linalg.norm(vecrd, ord=2) * np.linalg.norm(vecrs, ord=2))
    ret = green(freq, cv, vecrs, vecr) * (alpha - (1 - alpha)*(1 + (1j/(k*rd)))*cosgamma)
    return ret

Die hauptsächlich in Ambisonics verwendeten Funktionen sind nachstehend zusammengefasst.

`function.py`


import scipy as sp
from scipy.special import sph_harm
import numpy as np
from scipy.special import spherical_jn
from scipy.special import spherical_yn
from numpy import arctan2, sqrt
import numexpr as ne

def cart2sph(x,y,z):
    """ x, y, z :  ndarray coordinates
        ceval: backend to use: 
              - eval :  pure Numpy
              - numexpr.evaluate:  Numexpr """
              
    azimuth = ne.evaluate('arctan2(y,x)')
    xy2 = ne.evaluate('x**2 + y**2')
    elevation = ne.evaluate('arctan2(z, sqrt(xy2))')
    r = eval('sqrt(xy2 + z**2)')
    return azimuth, elevation, r

def EquiDistSphere(n):
    golden_angle = np.pi * (3 - np.sqrt(5))
    theta = golden_angle * (np.arange(n)+1)
    z = np.linspace(1 - 1.0 / n, 1.0 / n - 1, n)
    radius = np.sqrt(1 - z * z)
    
    points = np.zeros((n, 3))
    points[:,0] = radius * np.cos(theta)
    points[:,1] = radius * np.sin(theta)
    points[:,2] = z
    return points

def BnMat2(N,k,r):
    B = np.zeros((N+1)**2, dtype = np.complex).reshape((N+1)**2,)
    dx = 0
    for i in range(N+1):
        A = Bn2(i,k,r)
        for j in np.arange(-i,i+1,1):
            B[dx] = A
            dx +=1
    return B

def Bn2(n,k,r):
    kr = k*r
    y = 4 * np.pi * (1j**(n)) * spherical_jn(n,kr)
    return y
def calcSH(N,azi,ele):
    y_n = np.array([],dtype = np.complex)
    for n in range(N+1):
        for m in np.arange(-n,n+1,1):
            y_n = np.append(y_n,sph_harm(m,n,azi,ele))
    return y_n

def encode(ambi_order, azi, ele,k,r, mic_sig):
    Y_N = np.zeros((len(azi),(ambi_order+1)**2),dtype = np.complex)
    for i in range(azi.size):
        Y_N[i,:] = calcSH(ambi_order, azi[i],ele[i])
    invY = np.linalg.pinv(Y_N)             
    Wnkr = BnMat2(ambi_order,k,r)
    ambi_signal = np.diag(1/Wnkr) @ invY @ mic_sig
    return ambi_signal

def decode(ambi_order, azi, ele,r,lambda_,ambi_signal):
    Y_N = np.zeros((len(azi),(ambi_order+1)**2),dtype = np.complex)
    for i in range(azi.size):
        Y_N[i,:] = calcSH(ambi_order, azi[i],ele[i])
    invY = np.linalg.pinv(np.conj(Y_N.T)) 
    r_res = r % lambda_
    spkr_output = np.exp(1j*(2*np.pi*r_res/lambda_)) * (invY @ ambi_signal)
    return spkr_output

Ambisonics (internes Problem ver)

Der Schalldruck an jeder Position auf der Kugel kann wie folgt ausgedrückt werden.

p\bigl(\theta,\varphi,r,\omega\bigr) = \sum_{n=0}^\infty\sum_{m=-n}^n A_{nm} W_{nkr} Y_n^m(\theta,\varphi)

$ A_ {nm} $ ist der Expansionskoeffizient (Ambisonics-Koeffizient) und $ W_ {nkr} $ ist die Form des Mikrofons (in dieser Simulation wird ein hohles Array verwendet). Ambisonics macht zwei Dinge: Kodieren und Dekodieren. Zunächst wird der Expansionskoeffizient aus der Codierung erhalten, und das vom Lautsprecher ausgegebene Ansteuersignal wird durch Decodierung erhalten.

Kodieren

Finden Sie den Expansionskoeffizienten $ A_ {nm} . Die folgenden Formeln werden durch Matrizen dargestellt. $A_{nm}= diag(\frac{1}{W_{nkr}}) * Y_n^m(\theta,\varphi)^+ * p\bigl(\theta,\varphi,r,\omega\bigr) $$

Sie können den Ausdehnungskoeffizienten mehr finden. $ ^ + $ Repräsentiert eine verallgemeinerte inverse Matrix. Das Problem hierbei ist, dass die Welle nicht mit der Frequenz reproduziert wird, bei der $ W_ {nkr} $ 0 wird, weil sie durch $ W_ {nkr} $ geteilt wird. Um dieses verbotene Frequenzproblem zu vermeiden, werden Mikrofonarrays wie Hardball-Arrays und hohle Nieren-Arrays in Betracht gezogen.

Der folgende Code beschreibt zunächst den Fluss der primären Schallquelle (Schallfeld, das Sie abspielen möchten) aus der Mikrofonplatzierung, der Tonsammlung und der Codierung. Ab dem folgenden Artikel sind die Mikrofone gleichmäßig auf der Kugeloberfläche angeordnet.

blog.marmakoide.org/?p=1

`~encode.py`


import green_function as gr
import  function as fu
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.patches as patches
import matplotlib.cm as cm

#Anfangsbedingungen
ambi_order = 7              #Spaltungsauftrag
freq = 550                  #Frequenz
pos_s = np.array([2,0,0])   #Primäre Schallquellenkoordinaten x,y,z
vec_dn = np.array([1,0,0])  #Richtung der ebenen Welle
cv = 344                    #Schallgeschwindigkeit
r_sp = 1.5                  #Radius zum Aufstellen des Lautsprechers
r_mic = 0.05                #Mikrofonradius
#Diagrammbedingungen
Plot_range = 2.0
space = 0.05

#Berechnung
lambda_ = cv/freq
omega = 2*np.pi*freq
k = omega/cv


#Legen Sie die Lautsprecherposition fest
num_sp = (ambi_order+1)**2
xyz = fu.EquiDistSphere((ambi_order+1)**2)
xyz = np.array([xyz]).reshape(num_sp,3)

pos_sp = xyz * r_sp

azi_sp, ele_sp, rr_sp = fu.cart2sph(pos_sp[:,0],pos_sp[:,1],pos_sp[:,2])
ele_sp1 = np.pi /2 - ele_sp
#Entscheiden Sie sich für die Mikrofonanordnung
num_mic = (ambi_order+1)**2
pos_mic = xyz * r_mic
azi_mic, ele_mic, rr_mic = fu.cart2sph(pos_mic[:,0],pos_mic[:,1],pos_mic[:,2])
ele_mic1 = np.pi /2 - ele_mic
#Tonsammlung
mic_sig = np.array([])
for i in range(num_mic):
    mic_sig_ = gr.Plane_wave(freq,cv,pos_s,vec_dn,pos_mic[i])
    mic_sig = np.append(mic_sig, mic_sig_)

`encode.py`


#Kodieren
ambi_signal = fu.encode(ambi_order,azi_mic,ele_mic1,k,r_mic,mic_sig)
"""
def encode(ambi_order, azi, ele,k,r, mic_sig):
    Y_N = np.zeros((len(azi),(ambi_order+1)**2),dtype = np.complex)
    for i in range(azi.size):
        Y_N[i,:] = calcSH(ambi_order, azi[i],ele[i])
    invY = np.linalg.pinv(Y_N)             
    Wnkr = BnMat2(ambi_order,k,r)
    ambi_signal = np.diag(1/Wnkr) @ invY @ mic_sig
    return ambi_signal
"""

Dekodieren

Führt eine Verarbeitung durch, um den Expansionskoeffizienten des Lautsprechers mit dem durch Codierung erhaltenen Expansionskoeffizienten abzugleichen. Insbesondere ist es die folgende Formel.

A_{nm}= Y_n^m(\theta,\varphi)^* * w_l

$ w_l $ ist das Antriebssignal. Aus dieser Formel

w_l= Y_n^m(\theta,\varphi)^{*^{+}} * A_{nm}

Sie können mehr Antriebssignale erhalten. Bei dieser Decodierung stimmen die primäre Schallquelle und die sekundäre Schallquelle nur überein, wenn es sich um ebene Wellen handelt. Wenn Sie jedoch die Phase berücksichtigen, stimmen sie in dieser Formel möglicherweise nicht überein. Die Decodierungsformel, die auch die Phasenverschiebung berücksichtigt, lautet wie folgt. $r_{res} \equiv mod(r_{sp},\lambda)$ $w_l= e^{i\bigl(\frac{2\pi r_{res}}{\lambda}\bigr)}\times Y_n^m(\theta,\varphi)^{*^{+}} * A_{nm} $

Verwenden Sie die folgende Formel, um andere Wellen zu reproduzieren.

w_l= Y_n^m(\theta,\varphi)^{*^{+}} * diag(\frac{1}{F_{nkr}}) * A_{nm}

$ F_ {nkr} $ ist ein Wert, der von der Welle bestimmt wird.

In dem folgenden Code wird das Ansteuersignal aus dem durch Codieren erhaltenen Expansionskoeffizienten erhalten, und das Ansteuersignal wird erhalten.

`decode.py`


#Dekodieren
spkr_output = fu.decode(ambi_order,azi_sp,ele_sp1,r_sp,lambda_,ambi_signal)
"""
def decode(ambi_order, azi, ele,r,lambda_,ambi_signal):
    Y_N = np.zeros((len(azi),(ambi_order+1)**2),dtype = np.complex)
    for i in range(azi.size):
        Y_N[i,:] = calcSH(ambi_order, azi[i],ele[i])
    invY = np.linalg.pinv(np.conj(Y_N.T)) 
    r_res = r % lambda_
    spkr_output = np.exp(1j*(2*np.pi*r_res/lambda_)) * (invY @ ambi_signal)
    return spkr_output
"""

Über das Diagramm

Der Graph wird als Normalisierungsfehler bezeichnet. Verwenden Sie die folgende Formel. $ N_E(r,\omega) = 10\log_{10}\frac{{|p_{rep}(r,\omega)-p_{des}(r,\omega)|}^2}{{|p_{des}(r,\omega)|}^2} $

$ p_ {rep} (r, \ omega) $ ist der Schalldruck im wiedergegebenen Schallfeld, und $ p_ {des} (r, \ omega) $ ist der Schalldruck im primären Schallfeld, dh das gewünschte Schallfeld. Zusätzlich kann der Radius, der den Sweet-Spot-Bereich des Bereichs darstellt, in dem der Wiedergabefehler 4% oder weniger beträgt, aus der folgenden Formel berechnet werden.

r = \frac{c}{2 \pi f}N

$ c $ ist die Schallgeschwindigkeit, $ f $ ist die Frequenz und $ N $ ist die Cutoff-Reihenfolge. Die Cutoff-Reihenfolge beträgt in der Regel $ floor (\ sqrt {Anzahl der Lautsprecher oder Anzahl der Mikrofone}) - 1 $. Es kann $ floor (\ sqrt {Anzahl der Lautsprecher oder Anzahl der Mikrofone}) - 2 $ sein. Es ist möglich, unterschiedliche Aufträge für die Tonsammlung und -wiedergabe zu verwenden. In dieser Simulation stimmen sie überein.

`plot.py`


X = np.arange(-Plot_range,Plot_range+space,space)
Y = np.arange(-Plot_range,Plot_range+space,space)
XX,YY = np.meshgrid(X,Y)
P_org = np.zeros((X.size,Y.size), dtype = np.complex)
P_HOA = np.zeros((X.size,Y.size), dtype = np.complex)
for j in range(X.size):
    for i in range(Y.size):
        pos_r = np.array([X[j],Y[i],0])
        P_org[i,j] = gr.Plane_wave(freq,cv,pos_s,vec_dn,pos_r)
        for l in range(num_sp):
            G_transfer = gr.Plane_wave(freq,cv,pos_sp[l,:],-pos_sp[l,:],pos_r)
            P_HOA[i,j] += G_transfer *  spkr_output[l]

NE = np.zeros((X.size,Y.size), dtype = np.complex)
for i in range(X.size):
    for j in range(Y.size):
        NE[i,j] = 10*np.log10(((np.abs(P_HOA[i,j]-P_org[i,j]))**2)/((np.abs(P_org[i,j]))**2))

#Sweet Spot
rad = (ambi_order*cv)/(2*np.pi*freq)

#Graph
#Gewünschtes Schallfeld
sweet_spot = patches.Circle(xy=(0, 0), radius=rad, ec='r',fill=False)
r_sp_ = patches.Circle(xy=(0, 0), radius=r_sp, ec='b',fill=False)
plt.figure()
ax = plt.axes()
im = ax.imshow(np.real(P_org), interpolation='gaussian',cmap=cm.hot,origin='lower', 
               extent=[-Plot_range, Plot_range, -Plot_range, Plot_range],
               vmax=1, vmin=-1)
ax.grid(False)
ax.add_patch(sweet_spot)
ax.add_patch(r_sp_)
plt.axis('scaled')
ax.set_aspect('equal')
plt.colorbar(im)
plt.savefig('original_sound.pdf',bbox_inches="tight",dpi = 64,
            facecolor = "lightgray", tight_layout = True) 
plt.show()

#Wiedergegebenes Schallfeld
sweet_spot = patches.Circle(xy=(0, 0), radius=rad, ec='r',fill=False)
r_sp_ = patches.Circle(xy=(0, 0), radius=r_sp, ec='b',fill=False)
plt.figure()
ax2 = plt.axes()
im2 = ax2.imshow(np.real(P_HOA), interpolation='gaussian',cmap=cm.hot,origin='lower', 
               extent=[-Plot_range, Plot_range, -Plot_range, Plot_range],
               vmax=1, vmin=-1)
ax2.grid(False)
ax2.add_patch(sweet_spot)
ax2.add_patch(r_sp_)
plt.axis('scaled')
ax2.set_aspect('equal')
plt.colorbar(im2)
plt.savefig('reproduced_sound.pdf',bbox_inches="tight",dpi = 64,
            facecolor = "lightgray", tight_layout = True)
plt.show()

#Normalisierungsfehler
sweet_spot = patches.Circle(xy=(0, 0), radius=rad, ec='r',fill=False)
r_sp_ = patches.Circle(xy=(0, 0), radius=r_sp, ec='b',fill=False)
plt.figure()
ax3 = plt.axes()
im3 = ax3.imshow(np.real(NE), interpolation='gaussian',cmap=cm.pink,origin='lower', 
               extent=[-Plot_range, Plot_range, -Plot_range, Plot_range],
               vmax=0, vmin=-60)
ax3.grid(False)
ax3.add_patch(sweet_spot)
ax3.add_patch(r_sp_)
plt.axis('scaled')
ax3.set_aspect('equal')
plt.colorbar(im3)
plt.savefig('normalization_error.pdf',bbox_inches="tight",dpi = 64,
            facecolor = "lightgray", tight_layout = True)
plt.show()

Ergebnis

Es ist das Ergebnis jedes Schallfeldes. Gewünschtes Schallfeld

Wiedergegebenes Schallfeld

Normalisierungsfehler

`ambisonics.py`


import green_function as gr
import  function as fu
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.patches as patches
import matplotlib.cm as cm

#Anfangsbedingungen
ambi_order = 7              #Spaltungsauftrag
freq = 550                  #Frequenz
pos_s = np.array([2,0,0])   #Primäre Schallquellenkoordinaten x,y,z
vec_dn = np.array([1,0,0])  #Richtung der ebenen Welle
cv = 344                    #Schallgeschwindigkeit
r_sp = 1.5                  #Radius zum Aufstellen des Lautsprechers
r_mic = 0.05                #Mikrofonradius
#Diagrammbedingungen
Plot_range = 2.0
space = 0.05

#Berechnung
lambda_ = cv/freq
omega = 2*np.pi*freq
k = omega/cv


#Legen Sie die Lautsprecherposition fest
num_sp = (ambi_order+1)**2
xyz = fu.EquiDistSphere((ambi_order+1)**2)
xyz = np.array([xyz]).reshape(num_sp,3)

pos_sp = xyz * r_sp

azi_sp, ele_sp, rr_sp = fu.cart2sph(pos_sp[:,0],pos_sp[:,1],pos_sp[:,2])
ele_sp1 = np.pi /2 - ele_sp
#Entscheiden Sie sich für die Mikrofonanordnung
num_mic = (ambi_order+1)**2
pos_mic = xyz * r_mic
azi_mic, ele_mic, rr_mic = fu.cart2sph(pos_mic[:,0],pos_mic[:,1],pos_mic[:,2])
ele_mic1 = np.pi /2 - ele_mic
mic_sig = np.array([])
for i in range(num_mic):
    mic_sig_ = gr.Plane_wave(freq,cv,pos_s,vec_dn,pos_mic[i])
    mic_sig = np.append(mic_sig, mic_sig_)
#Kodieren
ambi_signal = fu.encode(ambi_order,azi_mic,ele_mic1,k,r_mic,mic_sig)

#Dekodieren
spkr_output = fu.decode(ambi_order,azi_sp,ele_sp1,r_sp,lambda_,ambi_signal)

#Gewünschtes Schallfeld(Primäre Schallquelle)Und suchen Sie das Wiedergabeschallfeld
X = np.arange(-Plot_range,Plot_range+space,space)
Y = np.arange(-Plot_range,Plot_range+space,space)
XX,YY = np.meshgrid(X,Y)
P_org = np.zeros((X.size,Y.size), dtype = np.complex)
P_HOA = np.zeros((X.size,Y.size), dtype = np.complex)
for j in range(X.size):
    for i in range(Y.size):
        pos_r = np.array([X[j],Y[i],0])
        P_org[i,j] = gr.Plane_wave(freq,cv,pos_s,vec_dn,pos_r)
        for l in range(num_sp):
            G_transfer = gr.Plane_wave(freq,cv,pos_sp[l,:],-pos_sp[l,:],pos_r)
            P_HOA[i,j] += G_transfer *  spkr_output[l]
#Normalisierungsfehler
NE = np.zeros((X.size,Y.size), dtype = np.complex)
for i in range(X.size):
    for j in range(Y.size):
        NE[i,j] = 10*np.log10(((np.abs(P_HOA[i,j]-P_org[i,j]))**2)/((np.abs(P_org[i,j]))**2))

#Sweet Spot
rad = (ambi_order*cv)/(2*np.pi*freq)

#Graph
#Gewünschtes Schallfeld
sweet_spot = patches.Circle(xy=(0, 0), radius=rad, ec='r',fill=False)
r_sp_ = patches.Circle(xy=(0, 0), radius=r_sp, ec='b',fill=False)
plt.figure()
ax = plt.axes()
im = ax.imshow(np.real(P_org), interpolation='gaussian',cmap=cm.hot,origin='lower', 
               extent=[-Plot_range, Plot_range, -Plot_range, Plot_range],
               vmax=1, vmin=-1)
ax.grid(False)
ax.add_patch(sweet_spot)
ax.add_patch(r_sp_)
plt.axis('scaled')
ax.set_aspect('equal')
plt.colorbar(im)
plt.savefig('original_sound.pdf',bbox_inches="tight",dpi = 64,
            facecolor = "lightgray", tight_layout = True) 
plt.show()

#Wiedergegebenes Schallfeld
sweet_spot = patches.Circle(xy=(0, 0), radius=rad, ec='r',fill=False)
r_sp_ = patches.Circle(xy=(0, 0), radius=r_sp, ec='b',fill=False)
plt.figure()
ax2 = plt.axes()
im2 = ax2.imshow(np.real(P_HOA), interpolation='gaussian',cmap=cm.jet,origin='lower', 
               extent=[-Plot_range, Plot_range, -Plot_range, Plot_range],
               vmax=1, vmin=-1)
ax2.grid(False)
ax2.add_patch(sweet_spot)
ax2.add_patch(r_sp_)
plt.axis('scaled')
ax2.set_aspect('equal')
plt.colorbar(im2)
plt.savefig('reproduced_sound.pdf',bbox_inches="tight",dpi = 64,
            facecolor = "lightgray", tight_layout = True)
plt.show()

#Finden Sie den Normalisierungsfehler
sweet_spot = patches.Circle(xy=(0, 0), radius=rad, ec='r',fill=False)
r_sp_ = patches.Circle(xy=(0, 0), radius=r_sp, ec='b',fill=False)
plt.figure()
ax3 = plt.axes()
im3 = ax3.imshow(np.real(NE), interpolation='gaussian',cmap=cm.pink_r,origin='lower', 
               extent=[-Plot_range, Plot_range, -Plot_range, Plot_range],
               vmax=0, vmin=-60)
ax3.grid(False)
ax3.add_patch(sweet_spot)
ax3.add_patch(r_sp_)
plt.axis('scaled')
ax3.set_aspect('equal')
plt.colorbar(im3)
plt.savefig('normalization_error.pdf',bbox_inches="tight",dpi = 64,
            facecolor = "lightgray", tight_layout = True)
plt.show()

Zusammenfassung

Da ich das, was ursprünglich in MATLAB in Python gemacht wurde, umgeschrieben habe, denke ich, dass es Stellen gibt, an denen es seltsam geschrieben ist, aber wenn es etwas gibt, das besser ist, kommentieren Sie es bitte.

Ich gebe den Code an Github https://github.com/koseki2580/ambisonics

Ambisonics Simulation Python

Über höhere Ambisonics

Über die in der Simulation verwendeten Wellen

green_function.py

function.py