急，请教高人有关统计最优声全息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 编辑 ]