苏格 发表于 2009-7-13 13:45

请教,有关统计最优声全息的问题,急,在线等

    小弟最近在做统计最优方面的问题,用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 编辑 ]

苏格 发表于 2009-7-13 13:50

利用声压重构声压,利用振速重构声压

利用声压重构声压,利用振速重构声压

parascape 发表于 2009-8-17 17:22

txl 发表于 2017-4-14 11:25

楼主的声压理论公式是不是有错,最优统计近场声全息不用考虑倏逝波的影响么。程序中通过声压求粒子的法向振速看不懂
页: [1]
查看完整版本: 请教,有关统计最优声全息的问题,急,在线等