##１．はじめに
この記事はMATLABで音響シミュレーションを行なう例として書いてみました。物理読み物のような感じで見ていただければと思います。
##２．シミュレーションの手法
例題として次のようなものをあげてみます。
この図の下のように振動板があって、図のような波形で一様に振動するとします。このとき、観測点Pにおいて観測される波形はどのようになりますか？ という問題です。
解き方は大きく分けて２通りあり、１つは有限要素法のようにメッシュを切って（伝搬領域を細かく分割するの意味）、解析するものです。演算量は多くなりますが、例えば回折の影響なども加味した結果が得られるのが特徴です。
もう１つの方法は光線法です。今回は光線法による解き方について述べてみます。
（ジュジュツ光線？ ちゃう。ちゃうで）
##３．光線法の手法とMATLABでの記述
まず、振動板がこのままだと解析できないので、細かく分割します。分割のピッチは波形の波長を考慮して決めます。ｎ個に分割したとしましょう。
細かくした微量振動板の座標を(x1,y1),・・・,(xn,yn）とし、微小振動板の間隔をｗとしますと、各振動板のｘ座標は、
x = (1:n)*w;
x = x-mean(x);%原点中心に持ってくる
こんな感じで表せて（各振動板のy座標はゼロ）、各振動板と観測点（Px,Py）の距離Lは、
L = abs(x-Px + 1i*Py); %上下方向は虚軸としてabsでユークリッド距離を求める
となり、各振動板から観測点Pまでの伝搬時間tpは、媒質の音速がVcなので、
tp = L/Vc;
波形も伝搬時間も離散時間系を用いる（サンプリング周期tsとする）と、
tc = tp/ts;
で、到達時間が下の図のようになったとすると、
この到達時間hv（hv=[6 2 2 0 2 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1]のようになる）と波形Aの畳み込みを行ない、
Ap = conv(A,hv);
到達波形Apを得ます。
※到達時間のヒストグラムについては、昨年のアドベントカレンダーで書いたMATLABで「たたみこみ」を語るを参照ください。
##４．２層の媒質による伝搬経路
ということで、ここまでは比較的順調に来ました。
ところで、伝搬媒質が１種類でなく２種類になった場合、例えば下の図のようになったとしたらどのようにして各振動板から観測点までの経路を求めるのでしょうか。
媒質１と媒質２の境界で屈折するわけですが、２つの媒質の音速が定まっているので（Vc1,Vc2)、スネルの法則により屈折角を求めることができます。
しかしながら、スネルの法則で簡単に求められるのは、振動板から出て、媒質境界のどの部分を通過したら屈折角がいくつで（結果として壁のどの部分に到達するか）ということであり、振動板と観測点の位置が定められているときに、媒質境界のどこを通るかというのはなかなか計算が難しいものです。いわんや境界面が下図のように曲面においておや、ということになります。
筆者は（いきなり論文口調）、当初、はさみうち法において境界の通過点を求めようとしたのですが、残念ながらあまり能率の良い手法とはなりませんでした。
##５．そこでマシンガン
すみません、マシンガンじゃなくて重機関銃だったみたいです。これです。
台座があって、銃が固定されてて、回転することで向きが変えられるやつです。これをそれぞれの微小振動板のところに置き、角度を変えながら撃ちます。
どの銃も同時に発砲するものとします。
こんな感じに角度を変えて撃ち、また変えて撃つ、というのを繰り返していきます。
結果として、壁面にこんな感じに弾痕が残ります。
（実際にはもっと密度が濃くならないと誤差が大きいです）。
これらの弾痕には、「どの銃（振動板）から撃たれたか」「境界面のどこを通ってきたか」の情報が含まれています。
と、いうことは、その弾痕をあけた弾が、銃から壁までにかかったフライトタイムがわかります。
次に壁を微小間隔で区切ります。
今、観測点Ｐがここにあるとすると、Ｐを含む区切りの中の弾痕に関する弾のフライトタイムを計算し、まとめます。結果は、こうなります。
なので、あとは波形と畳み込みすれば、Ｐ点（正確にはＰ点付近）で観測される音の波形が分かる、というわけです。
##６．MATLAB的な見せ場は？
すみません、意外にありませんでした。
直線と円弧の交点座標の計算は２次方程式の解の公式で解いていますし。
%直線 y=Ax+B と x軸上の円の弧（直径:R,弧とx軸の交点のx座標:L,
%T:弧が左（下）に凸が-1、右（上）に凸が+1
%媒質１，２の音速:Vc1,Vc2
%直線と弧の交点の座標:x2,y2
%媒質２における直線の方程式 y-A2 x +B2
C=L-R*T;
aa=A.^2+1;
bb=2*(A.*B-C);
cc=C^2+B.^2-R^2;
x2dum(1,:)=(-bb+sqrt(bb.^2-4*aa.*cc))./(2*aa);
x2dum(2,:)=(-bb-sqrt(bb.^2-4*aa.*cc))./(2*aa);
x2dum1=x2dum(1,:);
Ix1=find(imag(x2dum1)~=0);
if(isempty(Ix1)~=1)
x2dum(1,Ix1)=NaN;
x2dum(2,Ix1)=NaN;
end
y2dum(1,:)=A.*x2dum(1,:)+B;
y2dum(2,:)=A.*x2dum(2,:)+B;
if(T>0)
[x2 x2num]=max(x2dum);
else
[x2 x2num]=min(x2dum);
end
Ix2=find(x2num==1);
x2=x2dum(2,:);
y2=y2dum(2,:);
x2(Ix2)=x2dum(1,Ix2);
y2(Ix2)=y2dum(1,Ix2);
あとは、壁全体のデータからＰ点付近のデータを取り出して計算するあたり。
% ANGall:全弾丸の弾痕の位置（角度表示）
% Pang:観測点Pの位置（角度表示）
% ptcang:壁の分割ピッチ（角度表示）
% CNTall:全弾丸のフライト時間（サンプリング時間でカウント）
% A:振動波形
% 観測点Ｐ付近のデータを抜き出しCjに入力
Cj= ANGall(:)>=(Pang-ptcang/2) & ANGall(:)<(Pang+ptcang/2);
% 該当する弾丸のフライト時間をCNTdumに入力
CNTsub=CNTall(Cj);
%インパルス応答（フライトタイムの集計）を計算
h=histogram(CNTsub,min(CNTsub):max(CNTsub));
HISTcnt=h.Value;
%振動波形との畳み込み
AA=conv(A,HISTcnt);
ということで。
おつかれさまでした。