模拟点声源求全息面的声压出问题了

S0704082

导师让我用点声源来模拟近场声全息，已知点声源的声压和格林函数编MATLAB程序，请问把点声源放在原点处，在求全息面的声压时积分元面积是源面的面积吗？声压和格林函数在源面的曲面积分相当于二维离散卷积吗？请各位大哥大姐替小弟检查一下此程序的错误。
zh=0.1;
Lx=2;Ly=2;
f=343;c=340;
k=2*pi*f/c;w=2*pi*f;
p=1.293;
dx0=0.03125;dy0=0.03125;
x0=0;y0=0;z0=0;
x=-1:Lx/N:1;y=-1:Lx/N:1;
N=64;fs=10240;N1=1024;
t=0:1/fs:64/fs;
ua=0.1;r0=0.0005;
R=sqrt(x.^2+y.^2+zh.^2);
Q0=4*pi*r0.^2*ua;
L=exp(j*(w*t-k*R))/4*pi*R';
p=j*k*p*c*Q0*L;
K=exp(j*k*R)/2*pi*(R.^3)';
Gd=zh*(1-j*k*R)*K;
PP=sum(conv2(p,Gd)*dx0*dy0);
meshc(PP);

djh713

单从程序上就错误不少吧。另外你编的这全息程序不太对吧。

[ 本帖最后由 djh713 于 2008-4-3 17:24 编辑 ]

S0704082

楼主好，请教一个问题，我如果对NAH进行空间域到波数域的二维傅里叶变换进行校正，由于二位傅里叶变换的肯定存在泄露误差但是对空间声压变换到波数域时，乘以exp(-j*(kx*x+ky*y)),请问一下进行矫正的意义大吗？对全息面时域声压进行频谱校正时，如果用点声源进行模拟，得到全息面的声压就是理论结果啊，只有当用传声器阵列进行测量得到的的声压存在泄漏误差，进行频谱校正，就可以得到比较正确的结果。请问我的理解正确吗？

S0704082

请各位检查一下我的NAH仿真结果有没有错误，感觉幅值不对
点声源频率100HZ，L=2，N=32，全息面距离点声源0.1m，点声源位于（-0.5，0，0）点和（0.5，0，0）点，点声源半径0.001m,质点振速2.5m/s,图一为全息面上的声压和相位，图二为重建源面的声压和相位
clear
f1=300;
f2=300;
c0=340;
p0=1.25;
w1=2*pi*f1;
w2=2*pi*f2;
k1=2*pi*f1/c0;
k2=2*pi*f2/c0;
x1=-0.5;y1=0;
x2=0.5;y2=0;
lx=2;
ly=2;
N=16;
zh=0.1;
m=0:N-1;
n=0:N-1;
dx=lx/N;
dy=ly/N;
x=(m-N/2)*dx;
y=(n-N/2)*dy;
[X,Y]=meshgrid(x,y);
r1=sqrt((X-x1).^2+(Y-y1).^2+zh^2);
r2=sqrt((X-x2).^2+(Y-y2).^2+zh^2);
r0=0.001;
ua=2.5;
Q0=4*pi*r0.^2*ua;
ff1=exp(-j*k1*r1);
ff2=exp(-j*k2*r2);
p1=j*k1*p0*c0*Q0*ff1./(4*pi*r1);
p2=j*k2*p0*c0*Q0*ff2./(4*pi*r2);
p=p1+p2;
figure(1)
subplot(211)
meshc(X,Y,abs(p));
shading interp
colormap
subplot(212)
mesh(X,Y,angle(p));
shading interp
colormap

Y1=fft2(p);
YY=fftshift(Y1);
%mesh(X,Y,YY)
%shading interp
q=0:N-1;
s=0:N-1;
Kx=q*2*pi/lx;
Ky=s*2*pi/ly;
[kx,ky]=meshgrid(Kx,Ky);
rr1=kx.^2+ky.^2-k1.^2;
rr2=kx.^2+ky.^2-k2^2;
R01=sqrt(k1.^2-kx.^2-ky.^2);

Gd1=exp(j*(zh-0)*R01).*(rr1<=0);
     
R01=sqrt(kx.^2+ky.^2-k1.^2);
Gd2=sqrt(-2*pi*(zh-0)*R01).*(rr1>0);
      
Gd=Gd1+Gd2;
Gdd=1./(Gd);
ZZ=YY.*Gdd;
pp=ifft2(ZZ);
p0=pp;
figure(2)
subplot(211)
mesh(X,Y,abs(p0));
shading interp
colormap
subplot(212)
mesh(X,Y,angle(p0)*180/pi)
shading interp
colormap

[ 本帖最后由 S0704082 于 2008-8-27 21:27 编辑 ]

w89986581

呵呵，没有亲自做过就没有发言权。
所以，我提出的建议，仅供参考，或许有助于你的思考。

造成空间-波数域变换或者时间-频率域变换误差的原因有两个：“空间”采样“间隔”不够，“空间”采样“尺度”不够。具体一点来说，就空间-波数域变换而言，空间采样间隔不足，造成高波数域信息丢失，而且能量混迭到低波数分量上了。空间采样位置总是有限大的，其它位置上的信息丢失将造成截断误差。
对于第一种造成误差的因素，是任意一种处理方法都无法弥补的。不过一种情况除外，采样率不足的情况下没有发生混跌，即允许混迭采样的情况。
对于第二种造成误差的因素，可以采取虚拟阵元的方法，即根据测点位置外推测量区域以外的声场。如果外推得准，那么就可以降低误差。对于简单声源或者已知声源某些特征的情况，是可以外推准确的。

以上说的是变换时误差，还有一个就是声场逆变换自身的误差，即信息量误差。实际辐射面表面声波可以分为两类：一类对远场有贡献，一类对远场没有贡献，在辐射面垂直方向上迅速衰减。然而如，测量面与实际辐射面具有一定距离，因此理想情况（测点无限密、空间无限大或者完全包围了辐射面）下测量到的仅为对远场有贡献的波以及对远场无贡献、尚未衰减没（即高于测量通道噪声）的波，衰减那部分波对应的源就无从得知了。随着测量面距离实际辐射面越远，可测量到的第二类波越少。因此，这部分误差是由信息量损失造成的。

[ 本帖最后由 w89986581 于 2008-9-1 13:41 编辑 ]

mao

首先指出楼主的一些错误：

1.你解释为N=32,但程序中N取的是16啊，
2.你在回复djh713 的帖子说：“楼主好”其实是在向自己问好，:)，哈，找个乐了，不要介意。

言归正传，你的程序有不少问题，甚至错误，列举如下：
错误1：
你程序中的格林函数不对，格林函数在辐射圆外的形式与辐射圆内的形式不同，但你在程序中将这两个形式相加（对应于：Gd=Gd1+Gd2;Gdd=1./(Gd);），这是造成你重建结果不对的主要原因，试想公式都错了，结果会对吗？
改正方法：按NAH的步骤重新编写格林的形式
错误2：
你程序中的重建只是机械地按NAH的公式编写，其实NAH中会出现很多误差，正如w89986581 提到的，因此必须对全息面的数据补零后处理，相应地格林函数也有变化才行，否则会出现“Ghost Image”的现象，即卷绕误差。
改正方法：好好参考E.G.Williamas的相关文章，其中有详细的解释。

另外一个建议：
以你程序为准的话，你的分析系统所能达到的最高波数为：pi/(delt)=8*pi=25.13(rad/m)，而你的声源频率对应的波数k=5.54(rad/m),也就是说将有大量的波数成分位于辐射圆之外，
要是你的数据是实际中测量得到的，重建结果将会出现很大的错误，必须在K空间加滤波器才行。

欢迎你继续交流探讨！

花中飞霖

压和格林函数在源面的曲面积分不是二维离散卷积,曲面积分是连续的，程序中的时间T意义不大吧