马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。
您需要 登录 才可以下载或查看,没有账号?我要加入
x
clc
clear all
close all
%%%%%%%%%% 介质%%%%%%%%
ro=1.29;
c=344;
%%%%%%%%%% 源 %%%%%%%%%
rr=0.0001 %声源半径
f=1500; %频率
w=2*pi*f;
k=w/c;
x0=0;
y0=0;
z0=0;
%%%%%%%%% 测量面 %%%%%%%%
zh=0.1;
d_jianju=0.2;
num=7;
Lx=(num-1)*d_jianju; %测量面长度为Lx
Ly=Lx; %测量面宽度为Ly
n=1:num;
xh=-Lx/2+(n-1)*d_jianju; %一维向量,代表x轴上测量点的坐标
yh=-Ly/2+(n-1)*d_jianju; %一维向量,代表y轴上测量点的坐标
%%%%%%%%% 重建面 %%%%%%%%
zs=0.05;
xs=xh;
ys=yh;
for Nx=1:num
for Ny=1:num
Rh(Nx,Ny)=sqrt((zh-z0)^2+(xh(Nx)-x0)^2+(yh(Ny)-y0)^2); %测量面上各点到原点的距离
Rs(Nx,Ny)=sqrt((zs-z0)^2+(xs(Nx)-x0)^2+(ys(Ny)-y0)^2); %重建面上各点到原点的距离
ph(Nx,Ny)=j*k*ro*c*rr^2*exp(1i*k*Rh(Nx,Ny))./Rh(Nx,Ny);
ps(Nx,Ny)=j*k*ro*c*rr^2*exp(1i*k*Rs(Nx,Ny))./Rs(Nx,Ny);
ph1((Nx-1)*num+Ny)=ph(Nx,Ny); %把二维矩阵排列成一维向量,便于画图展示
ps1((Nx-1)*num+Ny)=ps(Nx,Ny);
end
end
delta_kx=2*pi/(d_jianju*(num-1));
delta_ky=2*pi/(d_jianju*(num-1));
kx=-pi/d_jianju:delta_kx:pi/d_jianju; %对应《近场声全息技术及其应用》p29的m△kx
ky=-pi/d_jianju:delta_ky:pi/d_jianju;
for i1=1:num;
for j1=1:num;
kz(i1,j1)=sqrt(k^2-(kx(i1)).^2-(ky(j1)).^2);
end
end
psp=ifft2(fft2(ph)*exp(-1i*kz*(zh-zs))); %《近场声全息技术及其应用》30页的(2.1.89)
for t=1:num
for y=1:num
psp1((t-1)*num+y)=psp(t,y);
end
end
figure(1)
plot(abs(ps1))
hold on
plot(abs(psp1),'--')
xlabel('y/m')
ylabel('声压幅值/ Pa')
legend('理论值','重建声压')
grid on
参考了论坛里的一个程序,点声源发出声场,用测量面的声压重建重建面的声压。公式参考的是《近场声全息技术及其应用》30页的(2.1.89)
|