请各位检查以下我的程序那里有问题

S0704082

小弟对声全息技术编程不是很了解，请各位给看看程序那里有问题。
    我的点声源频率为1000Hz，点声源大小r0=0.001m，振速2.5m/s，全息面距离点声源0.2m，重建面距离点声源0.01m，接受孔径Lx=2m，Ly=2m，N=32，计算出来点声源辐射到重建面的声压以及全息面上的声压后，然后将其变换到频域，再次有空间域变换到波数域，然后再由波数域变换到空间域，就能够的到重建的点声源以及重建面声压。
    我做的主要工作就是将全息面的声压变换到频域，研究傅里叶变换的误差，以及减小泄露误差，请问如果得到频域准确的声压幅值频率以及相位，那么这个声压如何表示才能进行空间域到波数域的变换？
   哪位大侠有关于声全息方面的程序可以分享以下，再次谢谢！！！
以下是我编写的程序，希望各位指教！！！
clear
f=1000;
c0=344;
p0=1.25;
w=2*pi*f;
k=2*pi*f/340;
lx=2;
ly=2;
N=32;
zh=0.2;
z=0.01;
N=32;
l1=0:N-1;
m1=0:N-1;
r1=0:N-1;
s1=0:N-1;
l=-N/2:N/2-1;
m=-N/2:N/2-1;
r=-N/2:N/2-1;
s=-N/2:N/2-1;
dx=lx/N;
dy=ly/N;
x=l*dx;
y=m*dy;
r11=sqrt(x.^2+y.^2+zh.^2);
r21=sqrt(x.^2+y.^2+z.^2);
r0=0.001;
ua=2.5;
'点声源强度'
Q0=4*pi*r0^2*ua;
t=0:0.001:0.031;
ff1=exp(j*(w*t-k*r11));
ff2=exp(j*(w*t-k*r21));
pa1=k*p0*c0*Q0./(4*pi*r11);
p1=j*pa1.*ff1;
sf=1024;
N0=2048;
td=0:1/sf:(N0-1)/sf;
fd=0:sf/N:(N0-1)/(N0*sf);
P1=fft(p1,N0)/N0;
zj=sqrt(1+(k*r1).^2)./(p0*c0*k*r11);
vra1=pa1.*zj;
pha1=atan(-1./(k*r11));   
vr1=j*vra1.*exp(w*t-k*r11+pha1);  %表面质点速度
I=pa1.^2/(2*p0*c0)  %声强大小
'理论全息面的声压'
pa2=j*k*p0*c0*Q0./(4*pi*r21);
p2=pa2.*ff2;
P2=fft(p2,N0)/N0;
'理论重建面的声压'
%%波数域的参数
kx=r*2*pi/lx;
ky=s*2*pi/ly;
ff=exp(-j*2*pi*(r/N.*l+s/N.*m));
%%全息面的复声压转换到波数
PP1=lx*ly/(2*pi*N^2)*sum(P1'*ff);
%%重建面的声压转换到波数域
PP2=lx*ly/(2*pi*N^2)*sum(P2'*ff);
%%%全息面的格林函数的傅立叶变换和声源面的格林函数
if kx.^2+ky.^2<=k^2
Gd1=exp(j*(z-zh)).*sqrt(k.^2-kx.^2-ky.^2);
  Gd2=exp(j*(zh)).*sqrt(k.^2-kx.^2-ky.^2);
else Gd1=exp(-(z-zh)*sqrt(kx.^2+ky.^2));
    Gd2=exp(-(zh)*sqrt(k^2-kx.^2-ky.^2));
end
k0=kx.^2+ky.^2;
Gd=Gd1;
zp0=PP1.*Gd1;
zp1=PP2./Gd2;
ss0=exp(j*(kx.*x+ky.*y));
%%%重构的全西面声压和重建面的声压
pp1=2*pi/(lx*ly)*sum(zp0'*ss0);
pp2=2*pi/(lx*ly)*sum(zp1'*ss0);
'重建面的声压'
result1=abs(pp1)
figure(1)
plot(result1)
'重建生源面的声压'
result2=abs(pp2)
     figure(2)
     plot(result2)