前回の復習
前回の記事の内容を確認してから本記事を読みことをお勧めいたいします。
スポンサーリンク
最小二乗法を用いた音場最適化方法（Least Squares Optimization）
“最小二乗法を用いた音場最適化方法”は目標応答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)
\)
ここで、学生の方は
学生
どーゆーことー？
そもそもq_targetとqは同じでしょ？
という感じですよね。
パイセン
q_targetとqは違うんですよ
順を追ってわかりやすく説明します
例えば、図1のような音場を構成するために、目標応答p_targetをダークゾーンの音圧は0でブライトゾーンの音圧は1となるように設定したとします。ただ、現実問題として、そのような解が存在しないことは明白です。なぜなら、音は伝搬（進行）するので、図のようにブライトゾーンの周りにダークゾーンが作ることはできないからです。これを行列計算で行っても同じことが言えます。目標応答p_targetを再現するqは算出できません。しかし、p_targetを再現”できそうな”qは得られます。そのqに伝達関数行列Gを乗じた応答pは、p_targetっぽいけど違います。p_targetっぽいpを更にそれっぽくなるようにする方法が”最小二乗法を用いた音場最適化方法”です。
図1 目標応答音場の例
最適化理論
先ほど、目標応答p_targetと最終的な計算結果のpは異なることを説明しました。両者の誤差を下式の誤差ベクトルeで表します。
$$
e = p_{Target} – p
$$
この誤差ベクトルeが最小になる入力qを求めればokっという考えが”最小二乗法を用いた音場最適化方法”です。最適化を計算するためには、評価関数というものを定義します。その評価関数は下式になります。ちなみに評価関数J_LSはスカラーです
$$
J_{LS} = e^{H} e + λ_{M}(q^H q – E_m)
$$
上記最適化式のλ_M….の項は入力パワーを制約するための条件です。スピーカーから超大音量を出して目的の音場を構成するような解が生じないように入力パワーq^H*qも小さくしましょうね、という意味があります。実際の計算ではλ_M=0のとき誤差e^H*eを最小化して、入力パワーも考慮させるためにλ_M≠0を設定します。この値はどうすればいいかは後日記載したいと思いますが、自分でパラスタ（パラメータースタディ）する方が簡単に理解できるかと思います。Emはシステムが必要とする入力パワーを制約するための定数です。
上記のような最適化（線形最適化）では「評価関数を設計変数の偏微分 ＝ 0」を解けばよいので、下式を解きます。
$$
\frac{\partial J_{LS}}{\partial q} = 0
$$
ブログ上で数式書くのはしんどいので、詳細は割愛しますが最終的に下記になります。
$$
q = (G^H Z + λ_{M} I)^{-1} G^H p_{Target}
$$
パイセン
最小二乗法の最適化の式は最終的にこの形に行きつくことが多いです。
まぁ、線形システムとして考えて、入力で応答を最適化しようと思って、逆行列計算しようとおもったらこの形になるのは想像できるかと思いますが…
MATLABでの検証
では、プログラムを組んで理論を検証しましょう。
ここでは入力－応答間の伝達関数を球面波で求めることにします。過去に球面波のプログラム方法について説明しているので、よかったらそちらも確認ください。
スポンサーリンク
先にプログラムでの検証結果を図2に示します。
- 音源の位置は、青丸で示したy軸±1.5m
- ブライトゾーンは、x軸0.1～0.7mでy軸0.4～1.0mのゾーン
- ダークゾーンは、x軸-0.7～-0.1mでy軸-1.0～-0.4mのゾーン
図2 音場制御の検証結果（最小二乗法を用いた音場最適化方法）
どうっすか？
ちなみに、最適化ではなく、単純に逆行列計算した結果を図3に示します。
図3 単純に逆行列計算した結果
どちらの結果も確かに2つのゾーンに分かれており、ダークゾーンは音圧が低くなっているので、音場制御自体はうまくいっています。ただ、正直なところ、両者の差がよくわかりませんね。
“最小二乗法を用いた音場最適化方法”は目標応答p_targetをうまいこと設定すると、効力をはっきしてきます。この辺については次回説明したいと思います。
今回はこのへんで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<0.1
r_bz=0.1; %r=0の時にP(ii1,ii2)が無限大になるのを防ぐ目的
end
if r_dz<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=[ones(size(Gb,1),1);zeros(size(Gd,1),1)];
% %-----------前記事からの追加分-------------% %
q=inv(G'*G+1*eye(size(G'*G)))*G'*p;
% %-----------前記事からの追加分-------------% %
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<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
コメント