菲涅耳计算全息图原代码

这是本人勉强做出来的菲涅耳全息图，请各位大虾帮忙看看，有没出错的地方。
本人还要做一个傅立叶全息图，但matlab用不熟练，实在是搞不出来了，请好心人看看能否帮个忙。

%菲涅耳积分傅立叶变换方法
function[f1,dx1,dy1,x1,y1]=fresnell(f0,M,N,dx0,dy0,z,lambda)
k=2*pi/lambda;
x0=ones(N,1)*[-M/2:M/2-1]*dx0;y0=[-N/2:N/2-1]'*ones(1,M)*dy0;
g=f0.*exp(j*0.5*k*(x0.^2+y0.^2)/z);
G=fftshift(fft2(fftshift(g)));
du=1./(M*dx0);dv=1./(N*dy0);
dx1=lambda*z*du;dy1=lambda*z*dv;
x1=ones(N,1)*[-M/2:M/2-1]*dx1;y1=[-N/2:N/2-1]'*ones(1,M)*dy1;
f1=G.*exp(i*0.5*k*(x1.^2+y1.^2)/z);
f1=f1./(i*lambda*z);


%菲涅耳全息图的生成与再现
clear;clc;
%读取图像
px=imread('1.jpg');
p1=rgb2gray(px);
figure;imshow(p1,[0 255]);
p1=im2double(p1);
%…………………………………………………………
picturesiz=size(p1);
N=picturesiz(1);
M=picturesiz(2);
cm=0.001;um=1e-6;
hx=100*cm;hy=100*cm;
dhx=hx/M;dhy=hy/N;
x=dhx.*(ones(N,1)*[-M/2:M/2-1]);
y=(dhy.*(ones(M,1)*[-N/2:N/2-1]))';
z=0.00002;
w=0.63*um;
k=2*pi/w;
[f,dx1,dy1,x1,y1]=fresnell(p1,M,N,dhx,dhy,z,w);
sit=pi/6;
ur=exp(-j*k*y*sin(sit));
f1=f+ur;
f1=mat2gray(f1);
figure;imshow(f1);
rp1=f1.*conj(ur);
[f2,dx2,dy2,x2,y2]=fresnell(rp1,M,N,dx1,dy1,z,w);
f2=abs(f2).^2;
f2=mat2gray(f2);
figure;imshow(f2);

兄弟，你这是抄来的吧

我仿佛也见过这个程序

见过此程序，有错误