Modusanpassungsmethode Simulation_Python

Modusanpassungsmethode

Was ist die Modusanpassungsmethode?

Das Modenanpassungsverfahren ist ein Verfahren, das darauf abzielt, die Moden des gewünschten Schallfeldes anzupassen, das an einem bestimmten Kontrollpunkt synthetisiert wird. Da es sich um eine numerische Lösung handelt, gibt es keine Einschränkungen für die Platzierung von Lautsprechern, sie kann jedoch numerisch instabil sein. In diesem Artikel möchte ich eine Simulation schreiben, wenn die Lautsprecheranordnung durch die Modusanpassungsmethode sphärisch gemacht wird. Higher Ambisonics Simulation Das Schallfeld wird reproduziert, indem der quadratische Fehler des im Artikel beschriebenen Expansionskoeffizienten minimiert wird.

Die in der Modusanpassungsmethode verwendeten Funktionen sind wie folgt. Darüber hinaus wird die Funktion von green_function im obigen Artikel von Higher Ambisonics beschrieben.

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 = shankel(n,1,kr)
    return y

def shankel(n,k,z):
    if k == 1:
        sha = spherical_jn(n,z) + 1j*spherical_yn(n,z)
    else:
        sha = spherical_jn(n,z) - 1j*spherical_yn(n,z)
    return sha

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 calcCmn(N,k,r,pos):
    C_N = np.zeros(((N+1)**2,int(len(pos))),dtype = np.complex)
    if len(pos) == 1:
        hn = BnMat2(N,k,r)
        C_N = hn * np.conj(calcSH(N,pos[0,0],pos[0,1])).T
    else:
        hn = BnMat2(N,k,r)
        for i in range(int(len(pos[:,0]))):
            C_N[:,i] = hn * np.conj(calcSH(N,pos[i,0],pos[i,1])).T
            
    return C_N

Simulation

Auch in dieser Simulation werden wir die Verwendung von sphärischen Wellen in Betracht ziehen.

\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)^*

Der Schalldruck auf die Kugel kann wie folgt ausgedrückt werden. $p\bigl(\theta,\varphi,r,\omega\bigr) = \sum_{n=0}^\infty\sum_{m=-n}^n A_{nm} 4\pi i^n j_n(kr) Y_n^m(\theta,\varphi) $

Aus diesen beiden Formeln $A_{nm} = \frac{k h_n(kr_s) Y_n^m(\theta_s,\varphi_s)^* }{4\pi i^{n-1}}$

Der Expansionskoeffizient kann erhalten werden. Der Ausdehnungskoeffizient des gewünschten Schallfeldes und der Ausdehnungskoeffizient des wiedergegebenen Schallfeldes können jeweils geschrieben werden.

A_{nm}^{des} = \frac{k h_n(kr_{s1}) Y_n^m(\theta_{s1},\varphi_{s1})^*}{4\pi i^{n-1}}
A_{nm}^{rep} = \frac{k h_n(kr_{s2}) Y_n^m(\theta_{s2},\varphi_{s2})^*d}{4\pi i^{n-1}}

$ s1 $ und $ s2 $ repräsentieren die primären und sekundären Schallquellen. Das Problem besteht darin, ein Ansteuersignal zu finden, das diesen quadratischen Fehler minimiert.

min {|A_{nm}^{rep} - A_{nm}^{des}|}^2 + \lambda d^H d

Wenn Sie diese Formel erweitern $C = h_n(kr_{s2}) Y_n^m(\theta_{s2},\varphi_{s2})^* $ $b = h_n(kr_{s1}) Y_n^m(\theta_{s1},\varphi_{s1})^* $ $d = \Bigl(C^H C + \lambda I\Bigr)^{-1} + C^H b$

Es können mehr Antriebssignale erhalten werden.

mode_matching_method.py


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

pos = np.array([0,0,0])


#Anfangsbedingungen
freq = 550
cv = 344
num_sp = 64
lambda_ = 10**-8
r_sp = 1.5

#Diagrammbedingungen
Plot_range = 2.5
space = 0.05

#Berechnung
omega = 2 * np.pi * freq
k = omega/cv
pos_azi, pos_ele, pos_r = fu.cart2sph(pos[0],pos[1],pos[2])
vec_R = np.array([pos_azi, pos_ele, pos_r])

#Lautsprecherposition
sp_xyz = fu.EquiDistSphere(num_sp)
sp_xyz = sp_xyz * r_sp
sp_azi,sp_ele, sp_r = fu.cart2sph(sp_xyz[:,0],sp_xyz[:,1],sp_xyz[:,2])
sp_ele1 = np.pi/2 - sp_ele
sp_sph = np.vstack([np.vstack([sp_azi,sp_ele1]),sp_r]).T

#Optimale Bestellbestimmung
N_q = int(np.floor(np.sqrt(num_sp))-1)

#Primäre Schallquelle
ps_xyz = np.array([1.5,1.5,0]) #Stellen Sie den Radius außerhalb der Position des Lautsprechers auf
ps_azi,ps_ele, ps_r = fu.cart2sph(ps_xyz[0],ps_xyz[1],ps_xyz[2])
ps_ele1 = np.pi/2 - ps_ele
ps_sph = np.vstack([np.vstack([ps_azi,ps_ele1]),ps_r]).T

#Berechnung der Modusanpassungsmethode
C = fu.calcCmn(N_q,k,r_sp,sp_sph);     
b = fu.calcCmn(N_q,k,ps_sph[0,2],ps_sph);     

I = np.identity(num_sp)
CC = np.conj(C.T) @ C
Cb = np.conj(C.T) @ b

d_MM = 	np.linalg.solve(CC + lambda_ * I,Cb)

#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)
P_org = np.zeros((X.size,Y.size), dtype = np.complex)
MM = 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.green2(1,freq,cv,ps_xyz,pos_r)
        for l in range(num_sp):
            G_transfer = gr.green2(1,freq,cv,sp_xyz[l,:],pos_r)
            MM[i,j] += G_transfer *  d_MM[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(MM[i,j]-P_org[i,j]))**2)/((np.abs(P_org[i,j]))**2))

#Sweet Spot
rad = (N_q*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.jet,origin='lower', 
               extent=[-Plot_range, Plot_range, -Plot_range, Plot_range],
               vmax=0.1, vmin=-0.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(MM), interpolation='gaussian',cmap=cm.jet,origin='lower', 
               extent=[-Plot_range, Plot_range, -Plot_range, Plot_range],
               vmax=0.1, vmin=-0.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_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()

Ergebnis

Gewünschtes Schallfeld image.png Wiedergegebenes Schallfeld image.png Normalisierungsfehler image.png

Zusammenfassung

Es war eine sphärische Anordnung der Modenanpassungsmethode. Ich habe den Code auf Github. https://github.com/koseki2580/mode_matching_method

Recommended Posts

Modusanpassungsmethode Simulation_Python
Binäre Methode
Spezielle Methode
Spezielle Methode