基于S-FFT的矩形孔衍射计算
光波是一种电磁波,它满足麦克斯韦方程组,将光波表示成球面波的形式,根据推导
[
1
−
4
]
[1-4]
[1−4],可以获得标量场下光波的菲涅尔衍射积分的表达式如下
U
(
x
,
y
)
=
exp
(
j
k
d
)
j
λ
d
∬
−
∞
∞
U
0
(
x
0
,
y
0
)
exp
(
j
k
2
d
[
(
x
−
x
0
)
2
+
(
y
−
y
0
)
2
]
)
d
x
d
y
U(x,y)=\frac{\exp(jkd)}{j \lambda d} \iint_{-\infty}^{\infty}U_0(x_0,y_0)\exp(\frac{jk}{2d}[(x-x_0)^2+(y-y_0)^2])dxdy
U(x,y)=jλdexp(jkd)∬−∞∞U0(x0,y0)exp(2djk[(x−x0)2+(y−y0)2])dxdy
上述式子中
x
0
y
0
,
x
y
x_0y_0,xy
x0y0,xy分别为衍射面和观察面的平面坐标系,
U
0
(
x
,
y
)
U_0(x,y)
U0(x,y)为物平面函数,
k
k
k为波数,
λ
\lambda
λ为光波波长。将上式改写为离散的傅里叶变换的形式为
U
(
p
Δ
x
,
q
Δ
y
)
=
exp
(
j
k
d
)
j
λ
d
exp
[
j
k
2
d
(
p
Δ
x
2
)
+
q
Δ
y
2
]
×
ϝ
[
U
0
(
m
Δ
x
0
,
n
Δ
y
0
)
×
exp
(
j
k
2
d
(
m
Δ
x
0
)
2
,
(
n
Δ
y
0
)
2
)
]
p
Δ
x
λ
d
,
p
Δ
y
λ
d
U(p\Delta x,q\Delta y) = \frac{\exp(jkd)}{j\lambda d}\exp[\frac{jk}{2d}(p\Delta x^2)+q\Delta y^2]\\\times \digamma[U_0(m\Delta x_0,n\Delta y_0)\times\exp(\frac{jk}{2d}(m\Delta x_0)^2,(n\Delta y_0)^2)]_{\frac{p\Delta x}{\lambda d},\frac{p\Delta y}{\lambda d}}
U(pΔx,qΔy)=jλdexp(jkd)exp[2djk(pΔx2)+qΔy2]×ϝ[U0(mΔx0,nΔy0)×exp(2djk(mΔx0)2,(nΔy0)2)]λdpΔx,λdpΔy
其中
ϝ
\digamma
ϝ表示傅里叶变换,设衍射平面尺寸为
L
0
L_0
L0,观察平面为
L
L
L,采样数为
N
N
N,则两者的关系可以由下式表示
L
0
=
N
λ
d
L
L_0=\frac{N\lambda d}{L}
L0=LNλd
根据以上的公式及说明,编写矩形孔的衍射计算小程序(MATLAB)具体如下:
r=512,c=r; %抽样数
a=zeros(r,c); %预设衍射面
a(r/2-r/4:r/2+r/4,c/2-c/4:c/2+c/4)=1; %生成衍射孔
lamda=6328*10^(-10);k=2*pi/lamda; %波长、波数
L0=5*10^(-3); %衍射面尺寸,单位:米
d=0.1; %衍射距离,单位:米
x0=linspace(-L0/2,L0/2,c); %衍射面x轴坐标
y0=linspace(-L0/2,L0/2,r); %衍射面y轴坐标
[x0,y0]=meshgrid(x0,y0); %衍射面二维坐标网格
L=r*lamda*d/L0, %观察面的尺寸,单位:米
x=linspace(-L/2,L/2,c); %观察面x轴坐标
y=linspace(-L/2,L/2,r); %观察面y轴坐标
[x,y]=meshgrid(x,y); %观察面二维坐标网格
%下面开始用式(3-4)计算衍射积分
F0=exp(j*k*d)/(j*lamda*d)*exp(j*k/2/d*(x.^2+y.^2)); % 赋值exp(ikd)/(iλd)exp[ik(x2+y2) /2d]
F=exp(j*k/2/d*(x0.^2+y0.^2)); %赋值exp[ik (x02+y02) /2d]
a=a.*F; %赋值U0(x0,y0)exp[ik (x02+y02) /2d]
Ff=fftshift(fft2(a)); %完成U0(x0,y0)exp[ik (x02+y02) /2d]的傅里叶变换
Fuf=F0.*Ff; %得到观察面上的光场分布U (x,y)
I=Fuf.*conj(Fuf); %计算观察面上的光强分布
figure,
subplot(1,2,1),imshow(I,[0,max(max(I))]),colormap(gray)
subplot(1,2,2),imshow(max(max(I))-I,[]),colormap(gray)
参考书目
[1] 古德曼 ,傅里叶光学导论[M],北京:电子工业出版社,2006,第3版
[2]李俊昌,衍射计算及数字全息[M],北京:科学出版社,2014
[3] 钱晓凡,信息光学数字实验室[M],北京:科学出版社,2014
[4]陈家壁,苏显渝.光学信息技术原理及应用[M].北京:高等教育出版社,2001