指向性をもったスピー​カの音圧分布のプロッ​トについて

現在自分はアクティブノイズコントロールについての研究を行おうと思っており、研究をおこなうに当たって、自由音場において適当に配置した無数のスピーカがどのような音場を作るのか（自分もこの分野に対しての知識がまだ浅いためこの言い方が正しいかはわかりませんが）というものをシミュレーション上で見れないかということを模索しています。理想は３次元で音圧分布が可視化できればと思っています。難しければ２次元でも構いません。

現在自分で行ったこととしては、とある書籍で進行波の音圧の複素実効値が以下の式であらわせるというものをみてこの式をもとに計算結果をプロットしてみました。

kが波数、ρが密度（自分が再現しようとしている音場では空気の密度だと思っています)

Aが適当な定数、rが音源からの距離です。

この式を使いプロットを行うことで、点音源が自由音場においてどのような音圧分布を作るかということは可視化することができたのですが、ここから一般的なスピーカのようにある程度の指向性を持たせる方法が分からず困っています。今の使用している式では距離でのみ音圧が決まると思うのですが、そうではなく、スピーカの正面は音圧が高く、後ろ側は音圧が低いというようなものを作成したいと思っています。

MATLAB上でそのような物理的なシミュレーションを行うことは可能でしょうか？

参考までに現在作成しているプログラムを載せておきます

function a = Three_dimention_soundsource_monopole(f,max,interval,min)
% 定数定義
w = 2*pi*f;
lambda = 340/f;
k = 2*pi/lambda;
p = 1.293;
c = 340
A = 1/k^2*p*c

% 計算部
[X,Y,Z] = meshgrid(min:interval:max);
  
x1 = X;
y1 = Y - max/2;
z1 = Z;

x2 = X - max/2;
y2 = Y ;
z2 = Z ;

x3 = X + max/2;
y3 = Y ;
z3 = Z ;

x4 = X ;
y4 = Y + max/2
z4 = Z ;

x5 = X - max/(2*sqrt(2));
y5 = Y - max/(2*sqrt(2));
z5 = Z ;

x6 = X + max/(2*sqrt(2));
y6 = Y + max/(2*sqrt(2));
z6 = Z ;

x7 = X + max/(2*sqrt(2));
y7 = Y - max/(2*sqrt(2));
z7 = Z ;

x8 = X - max/(2*sqrt(2));
y8 = Y + max/(2*sqrt(2));
z8 = Z ;
   

soundsorce1 =  1i*w*p*A/sqrt(x1.^2+y1.^2+z1.^2).*exp(-1i*k*sqrt(x1.^2+y1.^2+z1.^2));
soundsorce2 =  1i*w*p*A/sqrt(x2.^2+y2.^2+z2.^2).*exp(-1i*k*sqrt(x2.^2+y2.^2+z2.^2));
soundsorce3 =  1i*w*p*A/sqrt(x3.^2+y3.^2+z3.^2).*exp(-1i*k*sqrt(x3.^2+y3.^2+z3.^2));
soundsorce4 =  1i*w*p*A/sqrt(x4.^2+y4.^2+z4.^2).*exp(-1i*k*sqrt(x4.^2+y4.^2+z4.^2));
soundsorce5 =  1i*w*p*A/sqrt(x5.^2+y5.^2+z5.^2).*exp(-1i*k*sqrt(x5.^2+y5.^2+z5.^2));
soundsorce6 =  1i*w*p*A/sqrt(x6.^2+y6.^2+z6.^2).*exp(-1i*k*sqrt(x6.^2+y6.^2+z6.^2));
soundsorce7 =  1i*w*p*A/sqrt(x7.^2+y7.^2+z7.^2).*exp(-1i*k*sqrt(x7.^2+y7.^2+z7.^2));
soundsorce8 =  1i*w*p*A/sqrt(x8.^2+y8.^2+z8.^2).*exp(-1i*k*sqrt(x8.^2+y8.^2+z8.^2));

anser = abs(soundsorce1 + soundsorce2 + soundsorce3 + soundsorce4 + soundsorce5 + soundsorce6 + soundsorce7 + soundsorce8);

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;

Answer 1

[Shunichi Kusano]

簡単のため、音源が原点にあるとしまして、音が最大になる向き(スピーカーの正面方向：以下p1)を決めれば、任意の単位ベクトルp2とのなす角は、内積から計算できますので、あとは角度に対して、指向性がどうなるのかという関数を決めてあげれば、指向性を実装できます。

例えばx方向にcos関数で指向性を与える場合、

p1 = [1 0 0];
[X,Y,Z] = meshgrid([-1:0.01:1],[-1:0.01:1],[-1:0.01:1]);
p2 = [X(:),Y(:),Z(:)];
p2 = normalize(p2, 2, 'norm'); % 単位ベクトルに
theta = reshape(acos(p2*p1.'), size(X)); % 各グリッド上で、p1とのなす角
directivity = cos(theta); % 指向性の関数
directivity(directivity < 0) = 0; % 後ろ側はゼロ
% 確認のため表示
slice([-1:0.01:1],[-1:0.01:1],[-1:0.01:1],directivity, [0.3 0.6 1],[0],[0])
shading flat

directivityをsoundsourceに要素ごとに掛ければ指向性を含んだ音圧表現になるかと思います。

[kaiho inoue]

回答ありがとうございます。

お答えいただいたところ大変申し訳ないのですが、2つほど質問をさせていただきます。

1:このreshapeというものは何を目的に行っていますでしょうか？

調べてみたのですがいまいちreshape関数の意味が分からず、このような動作をしている理由が掴めませんでした、、、。

theta = reshape(acos(p2*p1.'), size(X));

2:指向性がどうなるかという関数は自分で適当に決めよいものでしょうか？

もし何か参考にするものなどございましたらご教授いただけると幸いです。

お答えいいただいたところ大変申し訳ございませんがよろしくお願いいたします。

[Shunichi Kusano]

reshape関数は、配列の形を変えるものです。

%例
x = [1:9]
reshape(x,[3,3])

今回は、meshgrid関数から出力したX,Y,Zをp2でX(:),Y(:),Z(:)とすることで、いったん縦ベクトルに変換しています。これは内積の計算を行列計算として効率よく行うための操作です。この結果acos(p2*p1.')は、meshgridで作成したグリッドの要素数分の縦ベクトルになっていますので、これをreshapeで元のグリッドの形にしています。文章だけだとピンときづらいと思います。各コマンドごとに変数の形がどのようになっているのかを確認しながらトレースしてみるとわかりやすいと思います。

次に指向性についてですが、シミュレーションの目的によるのではないかと思います。具体的な製品を想定しているのでしたらその指向性を使うのが一番かと思いますし、ざっくり傾向がわかればいいというものでしたら、代表的な特性にしておけばいいと思います。残念ですが私は音響分野の専門ではないのでどんな関数を使えばいいかのアイデアはありません。ざっと調べたところマイクだとカーディオイドが基本の一つのように見えますが、スピーカーだとどうなんでしょう。

[kaiho inoue]

丁寧なご説明ありがとうございます、、、！！

指向性の関数についてはまた自分でも検討してみたいと思います。

ありがとうございました