😺
MATLABで学ぶ音響工学 閉空間（直方体）2
前回の記事では閉空間の固有値と固有ベクトル（共振周波数とモードベクトル）の算出方法のプログラムを紹介いたしました。
今回は閉空間内の音圧の算出方法について紹介します。
なお、式の導出過程が難しいので、最終的な式のみを記載します。
（詳しい導出過程は右記参考文献でお願いします→P.A. Nelson, "Active Control of Sound", Academic Press(1993)）
下のプログラムを実行すると、下記結果が得られます。本プログラムの応用については次回の記事で説明したいと思います。
今回はこのへんでGood luck
私の個人ブログもよろしくお願いします。
メインコード（実行コード）
clear all
% addpath('E:\function')
% % % % size of rectangular parallelopiped
L1=1.9;
L2=1.3;
L3=1.2;
V=L1*L2*L3;
% % % % Constant value
co=344;
rho=1.205;
Zeta=0.001;
% % % % modal numbers
n_max=100;
Combi=combinator(n_max,3,'p','r')-1;
% % % % setting max frequency
f_max=2000;
freq=0:1:f_max;
w=2*pi*freq.';
% % % % sampling number
d=0.1;
% % % % xi=(x1,x2,x3)
xi=[0 0 0
0.5 0.5 0.5
1 1 1
1.5 1 1];
% % % % y_secondary_i=(ys1,ys2,ys3)
ysi=[0 0 0
1 0 0];
% % % % function wn
[wn]=angular_frequency(Combi,L1,L2,L3,co);
% % % % function psi psi(modal number , x_number)
[psix_front]=Psi(xi,Combi,L1,L2,L3);
[psiy]=Psi(ysi,Combi,L1,L2,L3);
% % % % sort of wn in order to get n
[wn,I] =sortrows(wn);
% % % % % % % % Select appropriate modes
% for n=1:length(wn)
% if wn(n)<=f_max*2*pi*1.5
% Wn(n,1)=wn(n);
% else
% break
% end
% end
Wn(1:4420,1)=wn(1:4420,1);
Psiy=psiy(I(1:length(Wn)),:);
PsiX=psix_front(I(1:length(Wn)),:);
% % % % % % % % % % % % % % %
G0=cellzeros(size(xi,1),size(ysi,1),length(freq));
% % % % % % % % % % % % % % % %
h = waitbar(0,'Please wait...');
for nx=1:1:size(xi,1)
for ny=1:1:size(ysi,1)
for n=1:length(Wn)
G0{nx,ny}(1:length(freq))=w.*rho*co^2......
./ ( V *( 2*Zeta*Wn(n).*w + i*(w.^2 - Wn(n)^2) ) ).....
*Psiy(n,ny)*PsiX(n,nx).....
+G0{nx,ny}(1:length(freq));
end
end
waitbar(nx/size(xi,1),h)
end
close(h)
figure
plot(freq,abs(G0{2,1}),'b')
hold on
plot(freq,abs(G0{3,1}),'r')
plot(freq,abs(G0{4,1}),'g')
% save('G0.mat','G0')
ちなみに前回の記事に記載したfunctionファイルも使用するので、そのファイルは下記にて参照ください。
Discussion