%菲涅耳全息圖的生成與重現
clear all;clc;
cm=0.01;um=1e-6;
M=512;N=512;                                             %全息圖分辨率
z=0.30;                                                  %菲涅爾傳輸距離
lambda=0.63*um;                                          %波長
k=2*pi/lambda;                                           %波數
hx=0.5*cm;hy=0.5*cm;                                     %全息圖大小
dhx=hx/M;dhy=hy/N;                                       %全息圖各方向增量
x=dhx.*(ones(N,1)*[-M/2:M/2-1]);
y=(dhy.*(ones(N,1)*[-M/2:M/2-1]))';
pic=imread('E256.bmp');
figure;imshow(pic);                                                      %原始圖
f0=zeros(N,1)*[1:M];
f0(129:384,129:384)=pic;                                                  %擴充圖像
dx0=lambda*z/(M*dhx);
dy0=lambda*z/(N*dhy);                                                  %物體采樣間隔
x0=ones(N,1)*[-M/2:M/2-1]*dx0;
y0=[-N/2:N/2-1]'*ones(1,M)*dy0;                                    %Input f0 is in natural order
g=f0.*exp(i*0.5*k*(x0.^2+y0.^2)/z);                                     %First phase factor
G=fftshift(fft2(fftshift(f0)));                                        %Fourier transform
f1=G./(i*lambda*z);
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;                                                     %Baseline for output
yo=2.0*cm;                                                                       %參考光位置
A=1;
ur=A*exp(i*0.5*k*(x1.^2+y1.^2+2*yo*y1)/z);                                    %球面參考光的復波函數
h=A*conj(f1).*exp(j*2*pi*yo*dy1);
H=abs(f1).^2+abs(ur).^2+h+conj(h);              %全息圖
HH=mat2gray(H);                       
figure;imshow(HH);                                                                              %顯示全息圖
%
fh=h;                                         
%a=0.002;b=0.3;
%rp1=pure(rp1,a,b,x1,y1);                                                           %消除零級和共軛像
[f2,dx2,dy2,x2,y2] = fresnelf(fh,N,N,dx1,dy1,z,lambda);
f2=abs(f2).^2;
f2=mat2gray(f2);
figure;imshow(f2);