💨
MATLABで学ぶ音響工学 球面波1
Zennのお試しということで、今回はMATLABで球面波のプログラム方法を紹介したいと思います。
下記の物性記号を用いて、球面波の圧力pは下式で表せます。
ρ：空気の密度
j：√-1
ω：角周波数
k：音波の波数（2π／波長）
今回は簡単にrをxに置き換えて、3次元から1次元に次元を下げて考えます。
上式のプログラムが以下になります。
clear all
x=0:0.1:10;
rho=1.293;
A=1;
c=340; %音速
freq=100; %周波数
w=2*pi*freq;
k=w/c;
t=0;
p=j*w*rho*A./(4*pi*x).*exp(j.*(k*x-w*t));
figure
plot(x,abs(p));
xlabel('X [m]')
ylabel('P [Pa]')
すると下図の結果が得られます。
大学の授業では上記結果を見て
「あー、音は距離に応じて反比例するのね、おもんないね」
で終わってしまうと思います。
そこで下記のようにアニメーションにしてみてください。
clear all
x=0:0.1:10;
rho=1.293;
A=1;
c=340; %音速
freq=100; %周波数
w=2*pi*freq;
k=w/c;
t=(0:0.1:2*pi)/w;
figure(1)
for ii1=1:length(t)
clf
p=j*w*rho*A./(4*pi*x).*exp(j*(k*x-w*t(ii1)));
plot(x,real(p));
xlabel('X [m]')
ylabel('P [Pa]')
ylim([-300 300])
pause(0.1)
% saveas(gcf,[num2str(ii1+100) '.png'],'png')
end
そうすると、波動が伝搬しているようなアニメーションが確認できます。
そこから少しだけ深堀をするために最初の式を見てください。
簡単そうに見える理論式や結果を見ても深堀するとまた違った景色が見えるものです。
この結果はxという座標上の定数xnという位置に、定数の振幅Anが複素数平面上で回転しているという見方もできます。
このような理解ができると、音響工学のプログラミングをするのが飛躍的に簡単になります。
今回はこのへんでGood luck
私の個人ブログもよろしくお願いします。
Discussion