请高手帮我看看统计最优声全息算法问题出在哪里，多谢了

苏格

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

mao

粗略地看了下你的程序和结果,感觉没有什么明显错误(由于各人编程习惯不同,别人看得很难受,所以没有细看),从设置的例子来看,你分析的频率为1500Hz,对应的波长约为0.23m,而你的全息面位于0.2m的地方,虽然小于一个波长,理论上没有错误,但是最好的做法是全息面位于1/3波长以内,这样才能得到比较精确的重建结果,所以我觉得你的重建结果在中心位置出现于理论值相差很大的原因上是由于全息面位置稍大了,这样损失了很多的倏逝波成分,重建精度自然不可能很高,建议你仿真时设置全息面的位置小一点,另外可以对一些影响因素如zh,zs,SNR等与重建精度的关系多验证,这样才是仿真的意义所在,不然就做一个例子,发现不了问题,达不到仿真研究的目的,一点小意见,欢迎指正!

苏格

首先十分感谢您的帮助，
介质是水，所以声速为1500，波长为1M，所以全息面距离源面1/5波长，而重建距离是控制在1/10波长的，这样应该得到很好的重建结果。我尝试使用不同重建位置，孔径，阵元间距，信噪比，都无济于事，所以还是确定是程序本身有问题
我的QQ是564242533
能留一下您的QQ或其他联系方式吗？
我想请您多指点我一二，这个程序我快检查一个月了，实在是找不到问题所在，现在马上要使用，十分着急
多谢您了

[ 本帖最后由 苏格 于 2009-7-14 12:47 编辑 ]

苏格

恳请大家帮忙指正，这里现在十万火急，先谢谢大家

苏格

谁来救小弟一命 啊

w89986581

你的A重复定义了，我无法确定你最初定义的A在哪里被调用！

盗草人

首先你的Kz有误，这样和严格按照公式编出来的结果是不一样的
其次ph(Nx,Ny)=exp(-j*Rh(Nx,Ny))./Rh(Nx,Ny);
        ps(Nx,Ny)=exp(-j*Rs(Nx,Ny))./Rs(Nx,Ny);
从这两句话我不知道你定义的是不是声压，因为基本的如频率声速声源波动你的都没有，虽然那些都是常数似乎不影响重建。
再次你的测量点数和测量点之间的间距的确定能不能有效地进行重建，也需要多试探。
另外既然是仿真信噪比的取值也有待商榷；
d=abs(zh-zs);这句话应该不需要abs吧？
仿真研究在于多试。

sc13579

你的kz没有讨论，还有另外一种情况，另外，点声源的辐射公式少了个k