Méthode de correspondance de mode Simulation_Python

Méthode de correspondance de mode

Quelle est la méthode de correspondance de mode?

La méthode de correspondance de mode est une méthode visant à faire correspondre les modes du champ sonore souhaité synthétisé à un point de contrôle spécifique. Puisqu'il s'agit d'une solution numérique, il n'y a aucune restriction sur le placement des enceintes, mais elle peut être numériquement instable. Dans cet article, je voudrais écrire une simulation lorsque la disposition des enceintes est rendue sphérique par la méthode de correspondance de mode. Simulation Ambisonique supérieure Le champ sonore est reproduit en minimisant l'erreur quadratique du coefficient de dilatation décrit dans l'article.

Les fonctions utilisées dans la méthode de correspondance de mode sont les suivantes. De plus, la fonction de green_function est décrite dans l'article ci-dessus de Higher Ambisonics.

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

Dans cette simulation également, nous envisagerons d'utiliser des ondes sphériques.

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

La pression acoustique sur la sphère peut être exprimée comme suit. $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) $

De ces deux formules $A_{nm} = \frac{k h_n(kr_s) Y_n^m(\theta_s,\varphi_s)^* }{4\pi i^{n-1}}$

Le coefficient de dilatation peut être obtenu. Le coefficient de dilatation du champ sonore souhaité et le coefficient de dilatation du champ sonore reproduit peuvent être écrits respectivement.

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 $ et $ s2 $ représentent les sources sonores primaires et secondaires. Le problème est de trouver un signal d'entraînement qui minimise cette erreur quadratique.

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

Si vous développez cette formule $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$

Plus de signaux d'entraînement peuvent être obtenus.

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


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

#Conditions du graphe
Plot_range = 2.5
space = 0.05

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

#Position du haut-parleur
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

#Détermination optimale des commandes
N_q = int(np.floor(np.sqrt(num_sp))-1)

#Source sonore principale
ps_xyz = np.array([1.5,1.5,0]) #Placez le rayon à l'extérieur où l'enceinte est placée
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

#Calcul de la méthode de correspondance de mode
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)

#Champ sonore souhaité(Source sonore principale)Et recherchez le champ sonore de reproduction
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)

#Graphique
#Champ sonore souhaité
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()

#Champ sonore reproduit
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()

#Erreur de normalisation
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()

résultat

Champ sonore souhaité image.png Champ sonore reproduit image.png Erreur de normalisation image.png

Résumé

C'était un arrangement sphérique de la méthode d'appariement de mode. J'ai le code sur Github. https://github.com/koseki2580/mode_matching_method

Recommended Posts

Méthode de correspondance de mode Simulation_Python
Méthode binaire
Méthode spéciale
Méthode spéciale