MATLAB--傅立叶变换--和衍射积分-模拟光学.doc

MATLAB编程用两种方法模拟光学实验

摘要：
利用MATLAB软件编程实现了用衍射积分的方法对单缝衍射、杨氏双缝干涉、黑白光栅衍射的计算机模拟；以及用傅立叶变换方法对简单孔径衍射、黑白光栅及正弦光栅夫琅和费衍射的模拟。
关键词：
MATLAB；衍射积分；傅立叶变换；计算机模拟
引言：
美国Mathworks公司推出的MATLAB，是一种集数值计算、符号预算、可视化建模、仿真和图形处理等多种功能于一体的优秀图形化软件。本文介绍了通过MATLAB软件编程实现用衍射积分和傅立叶变换实现夫琅和费衍射计算机模拟的方法。
计算机模拟为衍射实验的验证提供一条简捷、直观的途径。从而加深了对物理原理、概念和图像的理解。
正文：
大学教学课程中引入计算机模拟技术正日益受到重视，与Basic、C和Fortran相比，用MATLAB软件做光学试验的模拟，只需要用数学方式表达和描述，省去了大量繁琐的编程过程。下面来介绍利用MATLAB进行光学模拟的两种方法。
（一）衍射积分方法：
该方法首先是由衍射积分算出接收屏上的光强分布，然后根据该分布调制色彩作图，从而得到衍射图案。
1． 单缝衍射。
把单缝看作是np个分立的相干光源，屏幕上任意一点复振幅为np个光源照射结果的合成，对每个光源，光程差Δ=ypsinΦ,sinΦ=ys/D,光强I=I0(Σcosα)2+(Σsinα)2，其中α=2Δ/λ=πypys/λD
编写程序如下，得到图1
clear
lam=500e-9;
a=1e-3;D=1;
ym=3*lam*D/a;
ny=51;
ys=linspace(-ym,ym,ny);
np=51;
yp=linspace(0,a,np);
for i=1:ny
sinphi=ys(i)/D;
图1 单缝衍射的光强分布
alpha=2*pi*yp*sinphi/lam;
sumcos=sum(cos(alpha));
sumsin=sum(sin(alpha));
B(i,:)=(sumcos^2+sumsin^2)/np^2;
end
N=255;
Br=(B/max(B))*N;
subplot(1,2,1)
image(ym,ys,Br);
colormap(gray(N));
subplot(1,2,2)
plot(B,ys);
2． 杨氏双缝干涉
两相干光源到接收屏上P点距离r1=(D2+(y-a/2)2)1/2, r2=(D2+(y+a/2)2)1/2,相位差
Φ=2π（r2-r1）/λ,光强I=4I0cos2（Φ/2）
编写程序如下，得到图2
clear
lam=500e-9
a=2e-3;D=1;
ym=5*lam*D/a;xs=ym;
n=101;ys=linspace(-ym,ym,n);
for i=1:n
r1=sqrt((ys(i)-a/2).^2+D^2);
r2=sqrt((ys(i)+a/2).^2+D^2);
phi=2*pi*(r2-r1)./lam;
B(i,:)=sum(4*cos(phi/2).^2);
end
图2 杨氏双缝干涉的光强分布
N=255;
Br=(B/4.0)*N
subplot(1,2,1)
image(xs,ys,Br);
colormap(gray(N));
subplot(1,2,2)
plot(B,ys)
3． 光栅衍射
公式：I=I0（sinα/α）2（sin(λβ)/sinβ）2
α=(πa/λ)sinΦ
β=(πd/λ)sinΦ
编写程序如下：得到图3
clear
lam=500e-9;N=2;
a=2e-4;D=5;d=5*a;
ym=2*lam*D/a;xs=ym;
n=1001;
ys=linspace(-ym,ym,n);
图3 黑白光栅衍射光强分布
for i=1:n
sinphi=ys(i)/D;
alpha=pi*a*sinphi/lam;
beta=pi*d*sinphi/lam;
B(i,:)=(sin(alpha)./alpha).^2.*(sin(N*beta)./sin(beta)).^2;
B1=B/max(B);
end
NC=255;
Br=(B/max(B))*NC;
subplot(1,2,1)
image(xs,ys,Br);
colormap(gray(NC));
subplot(1,2,2)
plot(B1,ys);
（二）傅立叶变换方法：
在傅立叶变换光学中我们知道夫琅和费衍射场的强度分布就等于屏函数的功率谱。因此我们可以直接将光屏进行傅立叶变换，再处理得到衍射图样。
物体图像的生成可以直接由矩阵运算生成，也可利用Windows下的画图工具，生成一幅黑白图像，并调用命令函数imread（）输入图像，输入的图像是一个巨大的二维矩阵，利用MATLAB函数库中的fft2（）命令对该矩阵进行二位离散傅立叶变换，得到图像的频谱，该频谱是一个复数矩阵，然后用取模函数abs（）对该复数矩阵取模，得到振幅谱矩阵，利用函数fftshift（）对取模后的矩阵进行频谱位移，这是因为变换后的二维矩阵的直流分量位于图像的周边角，该函数交换矩阵的1、3象限和2、4象限，使直流分量移到频谱中心，从而使FFT频谱可视效果与实际图像相吻合。最后利用imshow（）函数将图像显示出来。
编写程序如下：

clear
a=imread('E:\1.bmp');
grid on
figure(1)
imshow(a,[])
afft=fft2(a);
aabs=abs(afft);
aabss=fftshift(aabs);
figure(2)
imshow(aabss,[])
colormap(gray);
colorbar
figure(3)
plot(aabss)
colormap(gray);
figure(4)
meshc(aabss)
maxx1=max(max(aabss));

输入黑白或灰度的衍射屏图像，得到输出的衍射图样和光强分布。
1． 简单孔径的衍射
(a)圆孔衍射：可以看到光强分布随孔径大小的变化

图4 圆孔的光强分布

(b)矩孔衍射：

图5 矩孔衍射的光强分布

（c）三角孔衍射：

图6 三角孔衍射的光强分布

2．单缝衍射：

图7 单缝衍射的光强分布

3．黑白光栅衍射：

缝宽不变，光栅常数变大：

缝宽变小，光栅常数不变：

图8 黑白光栅衍射的光强分布

可以明显地观察到干涉条纹受衍射因子的调制，以及缺级现象
4．正弦光栅衍射：
我们先设计程序制得一张正弦光栅，用正（余）弦函数作图如下：
clear all
xm=10*pi;ys=xm;
xs=linspace(-xm,xm,500);
B=cos(xs)+1;
N=255;
Br=B/2*N
image(xs,ys,Br);
colormap(gray(N));
将该图片保存成灰度BMP文件。调用程序得到衍射图样，只有三个衍射斑，与理论一致。

图9 正弦光栅衍射的光强分布

（三）利用傅立叶方法还可以得到其他有趣的衍射花样，并研究其规律性

每个方孔几何尺寸相同
可以看出，基本衍射图样还是
方孔衍射但是在任意两个方
孔连线的中垂线方向会出现
等间距暗纹，三角形摆放则会
出现类似蜂巢的结构。

两个方向上的光栅衍射叠加，
如果线条不平行则会在每个
中垂线方向都向外发展。

可以看到，
圆孔衍射的图样受到
椭圆孔衍射因子的调制

一些字母和符号也可以
产生很有趣的衍射图样：

-----------------------------------------------------------------------------
参考文献：
[1]：胡守信、李柏年，《基于MATLAB的数学实验》，科学出版社
[2]：张铮、杨文平、石博强、李海鹏，《MATLAB程序设计与实例应用》，中国铁道出版社
[3] 姚启钧. 光学教程[M] . 北京：高等教育出版社，2002.
[4] 陈锺贤. 计算物理学. 哈尔滨：哈尔滨工业大学出版社，2001.
[5] 游璞，于国萍. 光学. 北京：高等教育出版社，2003
[6] 程守洙等. 普通物理学. 北京：高等教育出版社，1999.
[7] 潭浩强. C程序设计. 北京：清化大学出版社，2001.

附录2 具体程序设计
#include<graphics.h>
#include<stdlib.h>
#include<math.h>
#include<stdio.h>

N=512;
disp('衍射孔径类型 1.圆孔 2.单缝 3.方孔')
kind=input('please input 衍射孔径类型：');% 输入衍射孔径类型
while kind~=1&kind~=2&kind~=3
disp('超出选择范围，请重新输入衍射孔径类型');
kind=input('please input 衍射孔径类型：');% 输入衍射孔径类型
end
switch(kind)
case 1
r=input('please input 衍射圆孔半径(mm)：');% 输入衍射圆孔的半径
I=zeros(N,N);
[m,n]=meshgrid(linspace(-N/16,N/16-1,N));
D=(m.^2+n.^2).^(1/2);
I(find(D<=r))=1;
subplot(1,2,1),imshow(I);
title('生成的衍射圆孔');

case 2
a=input('please input 衍射缝宽：');% 输入衍射单缝的宽度
b=1000;% 单缝的长度
I=zeros(N,N);
[m,n]=meshgrid(linspace(-N/4,N/4,N));
I(-a<m&m<a&-b<n&n<b)=1;
subplot(1,2,1);imshow(I);
title('生成的衍射单缝');

case 3
a=input('please input 方孔边长：');% 输入方孔边长
I=zeros(N,N);
[m,n]=meshgrid(linspace(-N/4,N/4,N));
I(-a/2<m&m<a/2&-a/2<n&n<a/2)=1;
subplot(1,2,1),imshow(I);
title('生成的方孔');
otherwise kind=input('please input 衍射孔径类型：');% 输入衍射孔径类型
end

% 夫琅禾费衍射的实现过程
L=500;
[x,y]=meshgrid(linspace(-L/2,L/2,N));
lamda_1=input('please input 衍射波长(nm)：');% 输入衍射波长;
lamda=lamda_1/1e6
k=2*pi/lamda;
z=input('please input 衍射屏距离衍射孔的距离(mm)：');% 衍射屏距离衍射孔的距离
h=exp(1j*k*z)*exp((1j*k*(x.^2+y.^2))/(2*z))/(1j*lamda*z);%脉冲相应
H =fftshift(fft2(h));%传递函数
B=fftshift(fft2(I));%孔频谱
G=fftshift(ifft2(H.*B));
subplot(1,2,2),imshow(log(1+abs(G)),[]);
title('衍射后的图样');
figure
meshz(x,y,abs(G));
title('夫琅禾费衍射强度分布')

实验输入：衍射孔径类型 1.圆孔 2.单缝 3.方孔
please input 衍射孔径类型：1
please input 衍射圆孔半径(mm)：3
please input 衍射波长(nm)：632
lamda =6.3200e-04
please input 衍射屏距离衍射孔的距离(mm)：1000000