音圧分布のプロットとその計算
17 views (last 30 days)
Show older comments
kaiho inoue
on 6 Jun 2020
Commented: kaiho inoue
on 8 Jun 2020
現在私は点音源が作る音圧分布をプロットしようとしています。
そこで以下のようなプログラミングを書きました。
% 定数定義
f = 100;
w = 2*pi*f;
lambda = 340/f;
k = 2*pi/lambda;
p = 1.293;
c = 340;
A = 1/k^2*p*c;
r = 1;
theta = (2 * pi)/4;
% グリッド作成
[X,Y,Z] = meshgrid(-3:0.1:3);
soundposition_x = [];
soundposition_y = [];
soundposition_z = [];
%主な原因と思われる場所
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%音源の位置計算
x1 = X - r*cos(theta);
%x1 = X; 本来r*cos(theta)はゼロになるため誤差を含めなけらばこれと同等のはず
y1 = Y - r*sin(theta);
z1 = Z;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%音圧計算
soundsorce1 = (1i*w*p*A/sqrt(x1.^2+y1.^2+z1.^2).*exp(-1i*k*sqrt(x1.^2+y1.^2+z1.^2)));
anser = abs(soundsorce1);
xslice = 0;
yslice = [];
zslice = 0;
% プロット部
slice(X,Y,Z,anser,xslice,yslice,zslice)
title('3次元空間の音圧分布')
xlabel('x[m]')
ylabel('y[m]')
zlabel('z[m]')
c = colorbar;
shading flat
なのですが、音源の位置を
x1 = X;
y1 = Y - r*sin(theta);
z1 = Z;
%プログラミング中ではr = 1,theta = 2*pi/4
とした時は音圧分布がうまくプロットできるのですが、(下にプロットされた図をしめします。)
x1 = X - r*cos(theta);
y1 = Y - r*sin(theta);
z1 = Z;
%プログラミング中ではr = 1,theta = 2*pi/4
とした時に、うまくプロットできません。
本来r*cos(theta)=0なのでどちらのプログラミングも同じ結果が出るはずなのですが、matlab上では誤差の影響なのか6.1232e-17という計算結果になることはわかりました。
しかしこの誤差がどう悪さをしてプロットがうまくできないのかが解明できていません。
多少の誤差が生じたとしても、音源の位置が本来の場所からちょっとずれて音が広がっていく様子が二つ目の図と同じようにプロットされると思ったのですが、うまく点音源の音圧分布がプロットされず困っています。
もし原因がわかる方がいたらご教授よろしくお願いいたします。
※今回音圧の計算に用いたものは以下の式を変形したものです。
また今後比較的楽に音源の数を増やせるようにするためにcos,sinの計算を使用しています。
0 Comments
Accepted Answer
Shunichi Kusano
on 7 Jun 2020
こんにちは。
私の環境だとプロットされる図が逆になりますが、貼り付けの際の間違いでしょうか。
いずれにせよ、結果に違いが生じるのは、
sqrt(x1.^2+y1.^2+z1.^2)
の値を考えると理解できるかと思います。2つの条件で言えば、音源においてゼロかものすごく小さな値になります。ですので、soundsorce1ではゼロ割になるか、ものすごく小さな値で割っているかの違いとなります。ものすごく小さな値で割るとsoundsorce1の値が音源においてものすごく大きくなるので、他の値が全部ゼロに見えてしまう、というからくりです。カラーバーを表示させて色の表示範囲を確認すれば合点がいくのではと思います。
このように本来音源そのものでは数値計算上値が変になってしまうので、別途処理が必要になるかと思います。sqrt(x1.^2+y1.^2+z1.^2)を(sqrt(x1.^2+y1.^2+z1.^2)+eps)にするとか(epsには適当に小さい値を入れる)、音源の位置は別の値に差し替えるとかでしょうか。
More Answers (0)
See Also
Categories
Find more on オーディオとビデオ in Help Center and File Exchange
Products
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!