splxidb.m

微穿孔板是声学领域重要的研究分支。基于此的得到的噪音吸收有重要的应用前景。这是基于matlab实现的共振吸声程序

function y=splxidb(f,xi,yi,zi)
%微穿孔板
%f frequency
%xi,yi,zi innerposition
%fi.m
%intgLx.m
%intgLy.m
%mapot.m
%vk.m
%each mfile above should have the same parameters below:
%Lx
%Ly
%Lz
%st
%wc
%c0
w=2*pi*f;
Lx=0.1;
Ly=0.2;
Lz=1;
st=0;
wc=2*pi*200;%低频范围最高频率为200Hz
c0=344;
PE=1;%外界声压
p0=1.21;
ksi=p0*15.6*10^(-6);
K=2*pi*f/ksi*p0;
a=0.05*10^(-3);%孔半径
dep=0.001;%孔深
Za=8*ksi*dep/(pi*a^4)*(1+(K*a)^2/32)^(1/2)+i*2*pi*f*p0*dep/(pi*a^2)*(1+1/((9+(K*a)^2/2)^(1/2)));%声阻抗率
wc=2*pi*200;%低频范围最高频率为200Hz
kc=wc/c0;
lmax=fix(kc*Lx/pi);
mmax=fix(kc*Ly/pi);
nmax=fix(kc*Lz/pi);
vol=(lmax+1)*(mmax+1)*(nmax+1);%P矩阵的行数
%以下为A矩阵点源
ma=3;%A矩阵x方向的点数
na=5;%A矩阵y方向的点数
numa=ma*na;%A矩阵的总点数
xunit=Lx/ma;%x方向点距
yunit=Ly/na;%y方向点距
A=rand(numa,2);
for j=0:numa-1,
    A(j+1,1)=((fix(j/na))+1)*xunit;
    A(j+1,2)=(j-(fix(j/na))*na+1)*yunit;
end
%以下为B矩阵
B=rand(vol,vol);
for h=1:vol,
    for j=1:vol,
        B(h,j)=0;
        for n=1:numa,
            B(h,j)=B(h,j)+fi(A(n,1),A(n,2),Lz+(A(n,2)-Ly/2)*tan(st),h)*fi(A(n,1),A(n,2),Lz+(A(n,2)-Ly/2)*tan(st),j)/Za;
        end
    end
end
%以下为F矩阵
F=rand(vol,1);
for h=1:vol,
    F(h)=0;
    for n=1:numa,
    F(h)=F(h)+fi(A(n,1),A(n,2),Lz+(A(n,2)-Ly/2)*tan(st),h)*PE*exp(-i*w/c0*distance(A(n,1),A(n,2),Lz+(A(n,2)-Ly/2)*tan(st)))/Za;
    end
end
%以下为V矩阵
V=rand(vol,vol);
for h=1:vol,
    for j=1:vol,
        if h==j
            V(h,j)=vk(w,j);
        else
            V(h,j)=0;
        end
    end
end
%以下为C矩阵
C=rand(vol,vol);
for h=1:vol,
    for j=1:vol,
        randh=mapot(h);
        l1=randh(1);
        m1=randh(2);
        n1=randh(3);
        randj=mapot(j);
        l2=randj(1);
        m2=randj(2);
        n2=randj(3);
        C(h,j)=-tan(st)/16*m2*pi/Ly*(intgLx(l1+l2)+intgLx(l1-l2))*...
               (...
               intgLy(m2+m1,n1+n2)+intgLy(m2+m1,-n1-n2)+intgLy(m2-m1,n1+n2)+intgLy(m2-m1,-n1-n2)...
               +intgLy(m2+m1,n1-n2)+intgLy(m2+m1,-n1+n2)+intgLy(m2-m1,n1-n2)+intgLy(m2-m1,-n1+n2)...
               )+...
               1/16*n2*pi/Lz*(intgLx(l1+l2)+intgLx(l1-l2))*...
               (...
               intgLy(m2+m1,n1+n2)+intgLy(-m1-m2,n2+n1)+intgLy(m1+m2,n2-n1)+intgLy(-m1-m2,n2-n1)...
               +intgLy(m1-m2,n1+n2)+intgLy(m2-m1,n1+n2)+intgLy(m1-m2,n2-n1)+intgLy(m2-m1,n2-n1)...
               );
           C(h,j)=-C(h,j);
       end
   end
   %以下为P矩阵
   LP=(i*p0*w).*B+C+V;
   RP=i*p0*w.*F;
   %x2=x\x1;
   P=LP\RP;%P=rand(vol,1)
   %以下为FI矩阵
   FI=rand(1,vol);
   for h=1:vol,
       FI(h)=fi(xi,yi,zi,h);
   end
   pc=FI*P;
   pc=abs(pc);
   y=20*log10(pc/(2*10^(-5)));