simulation ambisonique Python

À propos de Higher Ambisonics

L'ambisonique est une méthode qui réalise une reproduction sonore tridimensionnelle en reproduisant la direction d'arrivée des ondes sonores incidentes sur un certain point de réception du son, c'est-à-dire la directionnalité incidente des ondes sonores. Les caractéristiques directionnelles incidentes sont reproduites en développant les caractéristiques directionnelles tridimensionnelles du son incident à un certain point par une fonction harmonique sphérique et en déterminant la sortie de chaque haut-parleur dans le réseau de haut-parleurs afin de reproduire les caractéristiques directionnelles développées pendant la reproduction. Ceux qui se concentrent sur les ordres 0 et 1er sont appelés ambisoniques, et ceux qui utilisent également des fonctions harmoniques sphériques d'ordre 2 ou supérieur sont appelés ambisoniques d'ordre supérieur.

Dans cet article, je décrirai le code du flux de collecte du son avec un microphone, puis de disposition virtuelle des enceintes, et de reproduction du champ sonore.

À propos des ondes utilisées dans la simulation

La simulation décrite dans cet article ne traite que des ondes planes, mais il en existe trois types: les ondes planes, les ondes sphériques et les ondes sphériques directionnelles. L'équation d'onde pour chaque système de coordonnées cartésiennes peut s'écrire: $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 $ est l'angle entre la position de réception et le vecteur de position de la source sonore, et $ \ alpha $ est le coefficient directionnel ($ \ alpha $ = 1 correspond à l'onde sphérique). Ensuite, l'équation pour chaque onde dans le système de coordonnées polaires peut être écrite comme suit: $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)^*$

$ j_n (kr) $ est la fonction Vessel de sphère de première classe, et $ h_n (kr) $ est la fonction Hankel de première classe. $ ^ * $ Représente une co-conversion complexe.

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

Les fonctions principalement utilisées dans Ambisonics sont résumées ci-dessous.

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 (problème interne ver)

La pression acoustique à n'importe quelle position sur la sphère peut être exprimée comme suit.

$ A_ {nm} $ est le coefficient d'expansion (coefficient ambisonique), et $ W_ {nkr} $ est la forme du microphone (un tableau creux est utilisé dans cette simulation). Ambisonics fait deux choses: l'encodage et le décodage. Tout d'abord, le coefficient de dilatation est obtenu à partir du codage, et le signal de commande émis par le haut-parleur est obtenu par décodage.

Encoder

Trouvez le coefficient d'expansion $ A_ {nm} . Les formules suivantes sont représentées par des matrices. $A_{nm}= diag(\frac{1}{W_{nkr}}) * Y_n^m(\theta,\varphi)^+ * p\bigl(\theta,\varphi,r,\omega\bigr) $$

Vous pouvez trouver le coefficient de dilatation plus. $ ^ + $ Représente une matrice inverse généralisée. Le problème ici est que l'onde n'est pas reproduite à la fréquence où $ W_ {nkr} $ devient 0 car elle est divisée par $ W_ {nkr} $. Des matrices de microphones telles que des matrices de hardball et des matrices cardioïdes creuses sont envisagées pour éviter ce problème de fréquence interdit.

Le code suivant décrit d'abord le flux de la source sonore principale (champ sonore que vous souhaitez lire) à partir du placement du microphone, de la collection de sons et du codage. D'après l'article suivant, les microphones sont disposés uniformément sur la surface sphérique.

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

#Conditions initiales
ambi_order = 7              #Ordre de clivage
freq = 550                  #la fréquence
pos_s = np.array([2,0,0])   #Coordonnées de la source sonore principale x,y,z
vec_dn = np.array([1,0,0])  #Direction de l'onde plane
cv = 344                    #Vitesse du son
r_sp = 1.5                  #Rayon pour placer l'enceinte
r_mic = 0.05                #Rayon du microphone
#Conditions du graphe
Plot_range = 2.0
space = 0.05

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


#Déterminer l'emplacement des enceintes
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
#Décidez de la disposition du microphone
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
#Collection sonore
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

#Encoder
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
"""

Décoder

Effectue un traitement pour faire correspondre le coefficient de dilatation du haut-parleur avec le coefficient de dilatation obtenu par codage. Plus précisément, c'est la formule suivante.

$ w_l $ est le signal du lecteur. De cette formule

Vous pouvez obtenir plus de signaux d'entraînement. Dans ce décodage, la source sonore principale et la source sonore secondaire ne correspondent que lorsqu'il s'agit d'ondes planes, mais lorsque l'on considère la phase, elles peuvent ne pas correspondre dans cette formule. La formule de décodage qui prend également en compte le déphasage est la suivante. $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} $

Pour reproduire d'autres ondes, utilisez la formule suivante.

$ F_ {nkr} $ est une valeur déterminée par la vague.

Dans le code suivant, le signal d'attaque est obtenu à partir du coefficient d'expansion obtenu par codage et le signal d'attaque est obtenu.

decode.py

#Décoder
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
"""

À propos du graphique

Le graphique est appelé l'erreur de normalisation. Utilisez la formule suivante. $ 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) $ est la pression acoustique dans le champ sonore reproduit, et $ p_ {des} (r, \ omega) $ est la pression acoustique dans le champ sonore primaire, c'est-à-dire le champ sonore souhaité. De plus, le rayon représentant la zone du point idéal de la zone où l'erreur de reproduction est de 4% ou moins peut être calculé à partir de la formule suivante.

$ c $ est la vitesse du son, $ f $ la fréquence et $ N $ est l'ordre de coupure. L'ordre de coupure est généralement $ floor (\ sqrt {nombre de haut-parleurs ou nombre de microphones}) - 1 $. Il peut s'agir de $ floor (\ sqrt {nombre de haut-parleurs ou nombre de microphones}) - 2 $. Il est possible d'utiliser différents ordres pour la collecte et la lecture du son. Dans cette simulation, ils correspondent.

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)

#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.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()

#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(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()

#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,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

C'est le résultat de chaque champ sonore. Champ sonore souhaité

Champ sonore reproduit

Erreur de normalisation

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

#Conditions initiales
ambi_order = 7              #Ordre de clivage
freq = 550                  #la fréquence
pos_s = np.array([2,0,0])   #Coordonnées de la source sonore principale x,y,z
vec_dn = np.array([1,0,0])  #Direction de l'onde plane
cv = 344                    #Vitesse du son
r_sp = 1.5                  #Rayon pour placer l'enceinte
r_mic = 0.05                #Rayon du microphone
#Conditions du graphe
Plot_range = 2.0
space = 0.05

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


#Déterminer l'emplacement des enceintes
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
#Décidez de la disposition du microphone
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_)
#Encoder
ambi_signal = fu.encode(ambi_order,azi_mic,ele_mic1,k,r_mic,mic_sig)

#Décoder
spkr_output = fu.decode(ambi_order,azi_sp,ele_sp1,r_sp,lambda_,ambi_signal)

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

#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.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()

#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(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()

#Trouvez l'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ésumé

Fabriqué à l'origine avec MATLAB, il a été réécrit en Python, donc je pense qu'il y a des endroits où il est écrit étrangement, mais s'il y a quelque chose que vous aimez, veuillez commenter.

Je donne le code à github https://github.com/koseki2580/ambisonics