急，请教高人有关统计最优声全息MATLAB程序的问题
小弟最近在做统计最优方面的问题，用MATLAB进行最简单的单声源仿真，频率1500Hz，介质声速1500M/S，测量面距离全息面0.1M（1/10波长），测量面阵元间隔0.2M（1/5波长），孔径1.2M，7×7个阵元。但重建结果不理想，声源位置重建正确，但是重建幅值和理论值有较大的差距。仿真条件应该没有问题。算法也检查了很多次，实在找不到解决方法，把图和MATLAB程序放在论坛上，请大家帮我找找问题究竟出在哪里了，先谢谢大家了，在线等候。
程序：
clc
clear all
close all
%%%%%%%%%% 介质%%%%%%%%
ro=1000;
c=1500;
%%%%%%%%%% 源 %%%%%%%%%
f=1500;%频率
w=2*pi*f;
k=w/c;
A=1; %距点源一米处幅值
x0=0;
y0=0;
z0=0; %点源位置
%%%%%%%%% 测量面 %%%%%%%%
zh=0.2;%1/5波长
d_jianju=0.2;
num=7;
Lx=(num-1)*d_jianju;
Ly=Lx;
n=1:num;
xh=-Lx/2+(n-1)*d_jianju;
yh=-Ly/2+(n-1)*d_jianju;
%%%%%%%%% 重建面 %%%%%%%%
zs=0.10;%1/10波长
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)=exp(-j*Rh(Nx,Ny))./Rh(Nx,Ny);
ps(Nx,Ny)=exp(-j*Rs(Nx,Ny))./Rs(Nx,Ny);
ph1((Nx-1)*num+Ny)=ph(Nx,Ny);
ps1((Nx-1)*num+Ny)=ps(Nx,Ny);
end
end
%%%%%%%%%%%声压与振速 %%%%%%%%%%%%
for a1=1:num;
for a2=1:num;
Zh(a1,a2)=j*ro*c*k*Rh(a1,a2)*(1-j*k*Rh(a1,a2))/(1+(k*Rh(a1,a2))^2);%??
vh(a1,a2)=ph(a1,a2)./Zh(a1,a2);
vhz(a1,a2)=vh(a1,a2).*(zh-z0)./Rh(a1,a2); %法向阵速
Zs(a1,a2)=j*ro*c*k*Rs(a1,a2)*(1-j*k*Rs(a1,a2))/(1+(k*Rs(a1,a2))^2);
vs(a1,a2)=ps(a1,a2)./Zs(a1,a2);
vsz(a1,a2)=vs(a1,a2).*(zs-z0)./Rs(a1,a2); % 法向阵速
vhz1((a1-1)*num+a2)=vhz(a1,a2);
vsz1((a1-1)*num+a2)=vsz(a1,a2);
end
end
d=abs(zh-zs);
SNR=40;
theta2=(1+1/(2*k*d)^2)*10^(-SNR/10);
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;
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);
for m1=1:num; %测量面空间域
for n1=1:num;
A((j1-1)*num+i1,(m1-1)*num+n1)=exp(-i*(kx(i1).*xh(m1)+ky(j1).*yh(n1)+kz(i1,j1)*zh));
end
end
end
end
for m2=1:num
for n2=1:num
for i2=1:num
for j2=1:num
ALPHA((j2-1)*num+i2)=exp(-i*(kx(i2).*xs(m2)+ky(j2).*ys(n2)+ kz(i2,j2)*zs));
BELTA((j2-1)*num+i2)=exp(-i*(kx(i2).*xs(m2)+ky(j2).*ys(n2)+ kz(i2,j2)*zs)).*(-i*kz(i2,j2)); %p-v
DELTA((j2-1)*num+i2)=i/kz(i2,j2)*exp(-i*(kx(i2).*xs(m2)+ky(j2).*ys(n2)+ kz(i2,j2)*zs)); %v-p
end
end
vsz_v((m2-1)*num+n2)=vhz1*inv(A'*A+theta2*eye(num^2,num^2))*A'*ALPHA.'; %v-v
vsz_p((m2-1)*num+n2)=-1/(i*w*ro)*ph1*inv(A'*A+theta2*eye(num^2,num^2))*A'*BELTA.';%p-v
ps_v((m2-1)*num+n2)=(-i*w*ro)*vhz1*inv(A'*A+theta2*eye(num^2,num^2))*A'*DELTA.';%v-p
ps_p((m2-1)*num+n2)=ph1*inv(A'*A+theta2*eye(num^2,num^2))*A'*ALPHA.'; %p-p
end
end
figure(1)
plot(abs(vsz1))
hold on
plot(abs(vsz_p),'--')
hold on
plot(abs(vsz_v),'+-')
xlabel('y/m')
ylabel('振速幅值/ m/s')
legend('理论值','声压法','振速法')
grid on
figure(2)
plot(abs(ps1))
hold on
plot(abs(ps_p),'--')
hold on
plot(abs(ps_v),'+-')
xlabel('y/m')
ylabel('声压幅值/ Pa')
legend('理论值','声压法','振速法')
grid on
[ 本帖最后由 苏格 于 2009-7-13 14:13 编辑 ]
页:
[1]