MATLABで学ぶ信号処理・音響工学音場制御のプログラミングについて3 （効果的な目標音場の設定）

前回の復習
効果的なブライトゾーンの目標音場の設定方法
MATLABでの検証
1. 関連

前回の復習

前回の記事の内容を確認してから本記事を読みことをお勧めいたいします。

MATLABで学ぶ信号処理・音響工学　音場制御のプログラミングについて1

本記事ではフィードフォワードの音場制御方法の理論について、数式とMATLABのコードを交えながら説明いたします。理論としてはすごくシンプルなので、高校生でも理解できると思います。

MATLABで学ぶ信号処理・音響工学　音場制御のプログラミングについて2 （最小二乗法を用いた最適化）

最小二乗法を用いた音場最適化制御方法について、理論式やMATLABプログラムを示しながら説明します。ここまで簡単に詳しく説明してくれている文献はほかにないので、ぜひ読んでみてください。

効果的なブライトゾーンの目標音場の設定方法

“最小二乗法を用いた音場最適化方法”は目標応答p_targetと下式で算出されたqに伝達関数行列Gを乗じた応答pとの誤差を最小にする方法です。

$\left( \begin{array}{c} q_1 \\ \vdots \\ q_n \end{array} \right)= \begin{pmatrix} G_{11} & \cdots & G_{1n} \\ \vdots & \ddots & \vdots \\ G_{m1} & \cdots & G_{mn} \end{pmatrix} ^{\mathrm{-1}} \left( \begin{array}{c} p_{target1} \\ \vdots \\ p_{targetm} \end{array} \right)$

上式のような逆行列計算ではうまくいかないことがわかっているので、下式のように誤差ベクトルeを定義して、その誤差が最小になる入力ベクトルqを求めるんでしたね。前回の記事に記載した式を一部だけ下に示します。

$e = p_{Target} – p$
$J_{LS} = e^{H} e + λ_{M}(q^H q – E_m)$
$q = (G^H Z + λ_{M} I)^{-1} G^H p_{Target}$

そして前回の記事の最後に、”最小二乗法を用いた音場最適化方法”は目標音場p_targetをうまいこと設定すると、効力を発揮すると述べました。本記事では、この目標音場の設定方法について説明します。

これまで、目標音場の設定ではブライトゾーンの音圧が1、ダークゾーンの音圧が0になるようにしていました。ただし、この目標音場（解）は不自然であることがわかります。実際に音場を制御したときを考えてみてください。ブライトゾーンで設定しているすべての制御点（応答点）の音圧が1にすることは不可能です。これは、スピーカーから音が出力したときに、近接するスピーカーで音の干渉が生じるので、容易に想像できると思います。

実際には不可能な目標音場を設定していたとしても、”最小二乗法を用いた音場最適化方法”では、誤差が最小になる解を導きます。ただし、効率的に音場を制御したいならば、音場制御で目標音場が実現可能な音場を設定することが望ましく、その方が最適化の効果が大きくなることがわかっています。

では、実現可能な音場をどのように設定するのでしょうか？

答えは簡単です。
入力点（スピーカー）から同位相で音を出力したときの音場を目標音場にするのです。
もっと、テクニカルな話をすると、ブライトゾーンに近い入力点（スピーカー）からだけ同位相で音を出力させ、そのときの音場を目標音場にします。
（このとき思想としては、ブライトゾーン近傍のスピーカーはブライトゾーンに音を伝搬させ、ダークゾーンにも伝搬してしまう音に対してはダークゾーン近傍のスピーカーで音を打ち消すというものです。）

MATLABでの検証

では、プログラムを組んで理論を検証しましょう。
ここでは入力－応答間の伝達関数を球面波で求めることにします。過去に球面波のプログラム方法について説明しているので、よかったらそちらも確認ください。

MATLABで学ぶ音響工学　球面波1

はじめまして　MATLABパイセンです。今回はMATLABで球面波のプログラム方法を紹介したいと思います。参考にした文献は下記です。下記の物性記号を用いて、球面波の...

MATLABで学ぶ音響工学　球面波2

2次元での球面波のプログラム方法を紹介しています。

先にプログラムでの検証結果を図1に示します。

音源の位置は、青丸で示したy軸±1.5m
ブライトゾーンは、x軸0.1～0.7mでy軸0.4～1.0mのゾーン
ダークゾーンは、x軸-0.7～-0.1mでy軸-1.0～-0.4mのゾーン

図1　音場制御の検証結果@5000Hz（効果的なBZの目標音場の設定）

どうっすか？
ちなみに、ブライトゾーンの目標音場を1としたときの計算結果を図2に示します。

図2　音場制御の検証結果@5000Hz（BZの目標音場=1）

ただ、正直なところ、両者の差がよくわかりませんね。
では、もっと低周波数で比較してみましょう。ここでは500Hzを比較します。

図3　音場制御の検証結果@500Hz（効果的なBZの目標音場の設定）

図4　音場制御の検証結果@500Hz（BZの目標音場=1）

どうでしょうか？
この解析結果ではわかりにくいですが、効果的なブライトゾーンの目標音場の設定をした方が、ブライトゾーンとダークゾーンの音圧差が大きくなっていそうです。

ここでまでの議論で、音場制御がうまくいっているかどうかを評価する指標を説明していなかったので、「差が大きくなった」だとか「うまくいっていそう」といった曖昧な表現を使っていました。次回の記事では音場制御がうまくいっているのかの評価指標について説明します。

今回はこのへんでGood luck

clear all;close all;clc
 
 
% % % % % % 初期設定
x=-2:0.1:2;
y=-2:0.1:2;
rho=1.293;
A=1;
c=340; %音速
freq=5000; %周波数
w=2*pi*freq;
k=w/c;
 
% % % % % % 音源位置の設定
x_sp=-1:0.02:1;
y_sp=ones(1,length(x_sp));
input=[x_sp x_sp
    y_sp*(-1.5) y_sp*(+1.5)].';
 
% % % % % % ブライトゾーンとダークゾーンの設定
bz_geo=[];
dz_geo=[];
bz_x=-0.3:0.1:0.3;
bz_y=-0.3:0.1:0.3;
 
for ii1=1:length(bz_x)
    for ii2=1:length(bz_y)
        bz_geo=[bz_geo;bz_x(ii1)+0.4 bz_y(ii2)+0.7];
        dz_geo=[dz_geo;bz_x(ii1)-0.4 bz_y(ii2)-0.7];
    end
end
zone=[bz_geo;dz_geo];
 
Gb=zeros(size(bz_geo,1),size(input,1));
Gd=zeros(size(dz_geo,1),size(input,1));
 
for ii1=1:size(bz_geo,1)
    for ii2=1:size(input,1)
        r_bz=sqrt( (bz_geo(ii1,1)-input(ii2,1))^2 + (bz_geo(ii1,2)-input(ii2,2))^2 );
        r_dz=sqrt( (dz_geo(ii1,1)-input(ii2,1))^2 + (dz_geo(ii1,2)-input(ii2,2))^2 );
        if r_bz&amp;amp;amp;amp;amp;lt;0.1
            r_bz=0.1; %r=0の時にP(ii1,ii2)が無限大になるのを防ぐ目的
        end
        if r_dz&amp;amp;amp;amp;amp;lt;0.1
            r_dz=0.1; %r=0の時にP(ii1,ii2)が無限大になるのを防ぐ目的
        end
        Gb(ii1,ii2)=j*w*rho*A./(4*pi*r_bz).*exp(j*(k*r_bz));
        Gd(ii1,ii2)=j*w*rho*A./(4*pi*r_dz).*exp(j*(k*r_dz));
    end
end
G=[Gb;Gd];
% %-----------前記事からの追加分-------------% %
p_target=G*ones(size(input,1),1);
p_target(size(bz_geo,1)+1:end,1)=zeros(size(dz_geo,1),1);
 
q=inv(G'*G+1*eye(size(G'*G)))*G'*p_target;
% %-----------前記事からの追加分-------------% %
 
P=zeros(length(y),length(x));
for ii1=1:length(x)
    for ii2=1:length(y)
        for ii3=1:size(input,1)
            r=sqrt( (x(ii1)-input(ii3,1))^2 + (y(ii2)-input(ii3,2))^2 );
        if r&amp;amp;amp;amp;amp;lt;0.2
            r=0.2; %r=0の時にP(ii1,ii2)が無限大になるのを防ぐ目的
        end
        P(ii2,ii1)=P(ii2,ii1)+q(ii3)*j*w*rho*A./(4*pi*r).*exp(j*(k*r));
        end
    end
end
 
% x(ii1)およびy(ii2)という位置における振幅abs(P(ii1,ii2))は一意に決まるため
% P(ii1,ii2)を複素数平面上で回転させれば、音波が伝搬するアニメーションが作成できる
 
% figure(1)
% surf(x,y,abs(P))
% shading interp
% ylabel('X [m]')
% xlabel('Y [m]')
 
t=linspace(0,2*pi,6)/w;
for ii1=1:length(t)
clf
% surf(x,y,real( P*exp(j*(-w*t(ii1))) ))
surf(x,y,log10(abs(real( P*exp(j*(-w*t(ii1))))) ))
shading interp
hold on
plot3(input(:,1),input(:,2),ones(size(input,1),1)*50,'o')
plot3(zone(:,1),zone(:,2),ones(size(zone,1),1)*50,'ro')
ylabel('Y [m]')
xlabel('X [m]')
view([0 90])
% pause(0.1)
saveas(gcf,[num2str(ii1+200) '.png'],'png')
end