モードマッチング法

##　モードマッチング法とは
モードマッチング法は特定の制御点で合成された所望の音場のモードを一致させることを目的とした手法です。数値的な解法であるため、スピーカーの配置に制限がありませんが数値的に不安定になることもあります。
モードマッチング法でスピーカー配置を球状にした場合のシミュレーションを今回の記事では書きたいと思います。
高次アンビソニックスシミュレーション
の記事で書いてある展開係数の二乗誤差を最小化させる事で音場を再現します。

モードマッチング法で用いる関数は以下です。また、green_functionの関数に関しましては上記の高次アンビソニックスの記事に記載しております。

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

シミュレーション

今回のシミュレーションでも球面波を用いて考えます。

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

球面上の音圧は次のように表すことが可能です。
$$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) $$

この２式より
$$A_{nm} = \frac{k h_n(kr_s) Y_n^m(\theta_s,\varphi_s)^* }{4\pi i^{n-1}}$$

展開係数が求まります。
所望音場の展開係数と再現音場での展開係数はそれぞれ書くことができます。

$$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$と$s2$は1次音源、２次音源を表します。この二乗誤差を最小化するような駆動信号を求める問題となります。

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

この式を展開していくと
$$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$$

より駆動信号が求まる。

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


#初期条件
freq = 550
cv = 344
num_sp = 64
lambda_ = 10**-8
r_sp = 1.5

#グラフ条件
Plot_range = 2.5
space = 0.05

#計算
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])

#スピーカー位置
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

#最適な次数決定 
N_q = int(np.floor(np.sqrt(num_sp))-1)

#一次音源
ps_xyz = np.array([1.5,1.5,0]) # スピーカーを配置してる半径をよりも外に配置
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

#モードマッチング法の計算
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)

#所望音場(一次音源)と再現音場を求める
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))

#スイートスポット
rad = (N_q*cv)/(2*np.pi*freq)

#グラフ
#所望音場
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()

#再現音場
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()

#正規化誤差
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()

##　結果
所望音場

再現音場

正規化誤差

##まとめ
モードマッチング法の球状配置verでした。
Githubにコード載せてます。
https://github.com/koseki2580/mode_matching_method