☰
MATLAB光学衍射模拟:衍射积分与傅里叶变换的双路径实现
2026/10/15 10:07:00 网站建设 项目流程

简介:一份基于MATLAB的光学衍射模拟学习文档,面向物理专业学生、光学实验爱好者以及需要借助MATLAB进行科学计算与仿真的人员。文档以衍射积分与傅立叶变换两条技术路线为核心:衍射积分部分从单缝出发,逐步扩展到杨氏双缝干涉和黑白光栅衍射,给出了逐点计算光强并绘图显示的完整程序;傅立叶变换部分则通过对屏函数进行二维FFT与频谱位移,模拟圆孔、矩孔、三角孔、黑白光栅及正弦光栅等夫琅和费衍射图样,并讨论了干涉条纹受衍射因子调制及缺级现象。所有代码均可直接运行或修改参数,用于观察不同缝宽、缝数、光栅常数下的光强分布变化,能够帮助读者将抽象的物理公式转化为直观图像,加深对衍射、干涉等概念的掌握。资源包共1个doc文件,大小247KB,内含摘要、关键词、引言、程序清单与图例说明,结构完整便于查阅。该资源已有791人学习浏览,适合作为光学实验课程辅助材料或MATLAB模拟入门参考。

1. 光学衍射 MATLAB 模拟:这门课设为什么值得认真做一遍

做光学实验最怕的就是光学平台没调平、狭缝没对准,辛辛苦苦测得一组数据,回来一画图还和理论差得老远。后来我换了一条路:直接用 MATLAB 把单缝衍射、杨氏双缝、黑白光栅这些经典实验全部在电脑里模拟一遍,效果反而比实验室拍到的照片更干净,也更容易看出物理规律。这套光学衍射模拟代码的核心是两条技术路线:一条是衍射积分法,把波前切开成无数点源再叠加;另一条是傅里叶变换法,直接用 fft2 求屏函数的功率谱。本文就把这两条路从原理讲到参数调节,再补上几个我实际踩过的坑,适合正在做物理仿真课设、或者想把衍射图样快速可视化的人。

2. 光学衍射模拟的两条技术路线:衍射积分还是傅里叶变换

2.1 衍射积分法:把波前拆成点源的叠加

单缝衍射的典型处理方式,是把缝看成 np 个小光源的集合。每个小光源在接收屏上某一点贡献一个复振幅,最后把这些复振幅矢量叠加起来,得到总光强。这个思路对应惠更斯-菲涅耳原理的离散化:每个缝元都是次级子波源,接收屏上任意一点的振动是所有这些子波到达该点的叠加。

模拟里最核心的一段是光程差的计算:

lam = 500e-9; % 波长,单位米,这里是绿光 a = 1e-3; % 缝宽,单位米,即 1mm D = 1; % 缝到接收屏的距离 ym = 3*lam*D/a; % 屏幕显示范围,取中央明纹附近三个次级大 ny = 51; % 屏上采样点数 ys = linspace(-ym,ym,ny); % 屏上 y 坐标序列 np = 51; % 将缝分成 51 个点源 yp = linspace(0,a,np); % 缝内点源位置 for i = 1:ny sinphi = ys(i)/D; % 小角度近似下 sinφ ≈ y/D alpha = 2*pi*yp*sinphi/lam; % 每个点源相对缝中心的相位差 sumcos = sum(cos(alpha)); sumsin = sum(sin(alpha)); B(i,:) = (sumcos^2 + sumsin^2)/np^2; % 归一化光强 end

这段代码里 sinphi = ys(i)/D 是一个小角度近似,实际的 sinφ 是 y 除以斜距离,但当 D 远大于 y 时直接除以 D 误差很小。alpha 是每个点源在屏上某点产生的相位差,cos 和 sin 分别求实部和虚部叠加,最后平方相加得到强度。除以 np^2 是为了让最大光强不超过 1,方便后续映射到灰度值。

2.2 傅里叶变换法:夫琅和费衍射就是屏函数的功率谱

衍射积分在物理上直观,但算起来麻烦。对于夫琅和费衍射,有更快的办法:衍射屏的远场光强分布,恰好等于屏函数傅里叶变换的模平方,也就是功率谱。MATLAB 里用 fft2 做二维离散傅里叶变换,一条命令就代替了几层循环。

a = imread('diffraction_aperture.bmp'); % 读入黑白掩膜图 afft = fft2(a); % 二维 FFT,结果是一个复数矩阵 aabs = abs(afft); % 取模,得到振幅谱 aabss = fftshift(aabs); % 把零频移到画面中心 imshow(aabss,[]); % 显示衍射图样

这里有几个关键点。fft2 得到的频谱里,直流分量分布在四角,fftshift 的作用是把矩阵的第 1、3 象限和第 2、4 象限交换,让零频居中,这样看到的图案才和实际光屏上的分布一致。取模函数 abs 是因为频谱是复数,人眼只能看强度(振幅的平方)。实际模拟中如果需要更接近真实光强,可以用 aabs.^2 而不是 aabs。

2.3 两条路线怎么选:看你要物理过程还是要普适性

衍射积分法适合那些能写出解析表达式的目标:单缝、双缝、光栅,因为这些结构天然可以切成点源或线源,而且能直观看到每个参数(缝宽、缝数、光栅常数)对条纹形态的影响。缺点是计算量随采样点数增加很快,遇到任意形状孔径就要重写几何关系。

傅里叶变换法则几乎不关心孔径形状,只要是平面波照射、夫琅和费条件成立,把孔径图塞进去就能出结果。圆孔、矩孔、三角孔、甚至一行字母,都能一次性搞定。缺点是中间过程像一个黑匣子:输入图片、输出花样,物理图像不直接,而且必须注意采样率和补零的问题,否则会出现很假的混叠条纹。

实际做课设推荐两条路都写:衍射积分法用来验证明文公式,傅里叶法用来应付形状复杂、写不出公式的孔径。两个方法在夫琅和费区给出的是同一物理结果,正好可以互相验证。

3. 衍射积分法实战:单缝、双缝、光栅的代码与调参

3.1 单缝衍射:51 个点源怎么合成一条光强曲线

单缝是最基础的验证。缝宽 a=1mm,波长 500nm,屏距 1m 时,中央明纹半宽度大约是 lamD/a = 0.5mm。上面代码把显示范围设为 ym = 3lam*D/a,也就是能看到中央极大两侧各三个次级大,这个范围选得很聪明:太小看不清条纹结构,太大次极大被压缩到几个像素。

一个常被忽略的细节是 ny 和 np 的取值。np 是缝内点源数,如果 np 太小,点源之间相位差离散化严重,模拟出的旁瓣位置会偏离理论值;ny 是屏上采样数,如果 ny 太小,画出来的光强曲线是折线,次级大可能被漏掉。经验值:ny 取 201 以上,np 取 101 以上,曲线就平滑了。原版 51 是为了演示速度快,自己复现时建议加大。

3.2 杨氏双缝:用 r1、r2 直接算相位差

杨氏双缝和单缝最大的区别是:不把缝拆成连续点源,而是直接算两个缝中心到屏上某点的光程差。代码里 r1、r2 是两点到接收屏上第 i 个采样点的距离:

lam = 500e-9; a = 2e-3; % 双缝间距,注意这里 a 的物理含义是缝间距 D = 1; ym = 5*lam*D/a; % 条纹更密,需要更大的屏幕范围才能看到多根条纹 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

这里 a 是两缝间距,不是缝宽。很多初学者把单缝代码里的 a 直接搬过来,结果发现条纹间距完全不对。双缝干涉的条纹间距是 Δy = lamD/a,a 越大条纹越密,所以 ym 要取大一些,代码里取的是 5lamD/a,可以覆盖 5 个条纹周期。光强公式 4cos(phi/2).^2 来自双光束干涉的叠加,其中每一束本身是等振幅的,合成强度按余弦平方变化。

3.3 黑白光栅:缝间干涉因子与衍射因子相乘

光栅衍射的公式是 I = (sinα/α)^2 * (sin(Nβ)/sinβ)^2,前一项是单缝衍射因子,后一项是 N 个缝之间的干涉因子。程序里对应:

alpha = pi*a*sinphi/lam; % 单缝的相位参数 beta = pi*d*sinphi/lam; % 缝间相位参数 B(i,:) = (sin(alpha)./alpha).^2 .* (sin(N*beta)./sin(beta)).^2;

注意 N 在这里是缝数,不是灰度级。原代码里 N=2 只画了两条缝,如果你想观察缺级现象,把 N 改成 8 或 16,再配合调整 a 和 d 的比例就能看到。缺级条件是 d/a 等于整数倍时,某些主极大恰好落在单缝衍射的暗纹位置。比如 d=5a 时,±5、±10 级主极大消失。这是光学课设最喜欢考的一个观察点,建议改参数重新跑一遍。

beta 接近 π 的整数倍时 sin(beta) 会趋于零,公式出现 0/0,此时需要做极限处理。MATLAB 里一个简单的办法是在分母上加一个极小量 eps,或者直接用 sinc 函数的近似表达式替代,否则图上会出现突兀的 NaN 空白。

4. 傅里叶变换法实战:从一张黑白图到衍射花样

4.1 标准流程:imread、fft2、abs、fftshift、imshow

傅里叶方法的核心是把衍射屏图像当作一个二维函数,直接吃掉它的透射率分布。输入可以是程序生成的矩阵,也可以是 Windows 画图工具存的黑白 BMP。标准流程如下:

clear; close all; a = imread('aperture.bmp'); % 读入掩膜图 a = double(a); % 转 double,否则 fft2 按整数算精度不够 afft = fft2(a); aabs = abs(afft); % 振幅谱 aabss = fftshift(aabs); % 零频居中 imshow(aabss,[]); % 自动拉伸灰度显示

imread 读进来通常是 uint8 类型,范围 0-255。fft2 对 uint8 也能算,但绝对值动态范围很大,中央极大可能比旁边亮几个数量级,直接显示会只剩一个亮点。所以转 double 是第一步。显示时 imshow 第二个参数传空数组 [],意思是把矩阵最小值映射到黑、最大值映射到白,这样旁瓣才能勉强看出来。如果仍然太暗,可以取 log(1+aabss) 做对数变换,把低亮度的细节拉出来。

4.2 圆孔、矩孔和三角孔:几何形状决定衍射斑的分布特征

圆孔的夫琅和费衍射是艾里斑,中心亮圆加同心圆环;矩孔衍射是十字交叉条纹,两个方向的条纹间距分别由矩孔的长宽决定;三角孔的衍射花样更复杂,会出现沿孔边界法线方向的条纹扩展。这些用手写衍射积分验算非常费劲,但傅里叶法只需要换一张输入图。

实际操作中,孔径图建议用程序生成而不是画图工具手绘,因为手绘的边缘锯齿会引入额外的高频噪声,衍射图上出现不该有的星芒。我一般用矩阵逻辑索引生成几何形状,比如圆孔就是:

[x,y] = meshgrid(-1:0.01:1); r = sqrt(x.^2 + y.^2); aperture = double(r <= 0.3); % 半径 0.3 的圆孔

这样生成的孔边缘是理想阶跃,不会混入鼠标画的锯齿。孔径尺寸和矩阵分辨率要一起考虑:如果圆孔半径只占几个像素,FFT 出来的艾里斑会严重变形;经验是让孔径直径至少占矩阵边长的 1/10。

4.3 正弦光栅:为什么输出只有三个衍射斑

正文里给出的正弦光栅生成办法是先画余弦条纹再保存成 BMP,再走 FFT 流程。其实完全不用绕这一圈,直接在 MATLAB 里生成灰度图矩阵:

xm = 10*pi; xs = linspace(-xm, xm, 500); B = cos(xs) + 1; % 余弦条纹,值域 0~2 N = 255; Br = B/2*N; % 归一化到 0~255

画出来是一排黑白渐变的条纹。对这个图像做 fft2 并取模平方,会出现三个斑:中心斑对应直流分量,左右两个斑对应余弦的 ±1 级。原因很直接:余弦函数在频域就是两个冲激,加了一个直流偏置才让灰度非负。所以正弦光栅的衍射级次只有 -1、0、+1 三级,不会出现黑白光栅那种多级谱。这个例子特别适合说明为什么光栅的“正弦”形态决定了谱的稀疏性,也解释了为什么实际光谱仪更喜欢正弦光栅来减少高级次重叠。

4.4 多孔阵列与多方向光栅的调制效应

正文还提到方孔阵列、三角形排布和斜向光栅。这类实验用傅里叶法做几乎零成本:先生成掩膜矩阵,再用同样的 fft2 流程。方孔阵列的衍射花样里,每个孔自身的衍射因子和阵列的干涉因子相乘,结果是在方孔衍射的十字条纹基础上,又叠加上阵列周期对应的细密斑點。两个方向光栅叠加时,如果两族条纹不平行,频谱上会出现沿各自法线方向展开的两组衍射斑,交叉处还会有复合的强度调制。

5. 模拟过程中的高频踩坑:五个最容易翻车的点

5.1 图像读入后没转 double,衍射图全是一个白点

现象:输入图片跑完 fft2,imshow 出来除了中心一个亮点,周围全是黑的。原因:uint8 图像做 FFT 后,直流分量巨大,其他频率分量相对值太小,被 imshow 的自动拉伸压没了。解决:读图后先用 double() 转换,必要时对频谱取对数 log(1+abs(fftshift(fft2(a)))) 再显示。我后来凡是处理图像输入,第一步固定转 double,绝不给 uint8 进 FFT 的机会。

5.2 光栅公式出现 NaN,整条强度曲线断掉

现象:黑白光栅的光强曲线在某个位置突然归零或消失,plot 里出现缺口。原因:beta 等于 π 的整数倍时,sin(beta) 为 0,程序算出 0/0 得到 NaN。解决:在分母加一个极小量,例如 sin(beta) + 1e-10;或者判断 beta 接近 k*pi 时直接用极限值替换。血泪经验:凡是公式里出现分式的仿真,先检查分母零点,别等图画出来再找原因。

5.3 杨氏双缝把双缝间距当成缝宽,条纹数少了一半

现象:双缝程序跑出来条纹间距明显偏大,和理论不符。原因:把单缝代码中的 a 直接粘过来,a 在双缝里应该是缝间距而不是缝宽,二者差一个数量级。解决:先确认物理量定义。双缝的条纹间距由间距 a 决定,缝宽只影响单缝衍射包络。建议在代码注释里写明 a 的物理含义,这两个量混用的概率非常高。

5.4 频谱中心没移到画面中央,图案怎么都不对

现象:衍射图样的中央极大出现在图像角落,整体图案像被切了一刀。原因:fft2 的结果默认把零频放在矩阵左上角,必须做 fftshift 才能移到中心。解决:每次 FFT 后用 fftshift 然后再显示。检查是否处理过的最快方法:看图案对称中心是否在画面中心,如果对称但偏角,就是少了 fftshift 这一步。

5.5 图样有奇怪的周期条纹,不是预期的衍射斑

现象:输出图像出现规则的斜纹或栅格状伪影,和理论衍射差得很远。原因:输入图像的边界处理不当,或者图像尺寸不是 2 的幂时 FFT 出现泄漏;更常见的是手绘图像边缘有抗锯齿渐变,引入了高频分量。解决:用程序生成二值掩膜,保证边缘是硬边;输入矩阵尺寸可以补零到 512×512 或 1024×1024,减少环效应。另外,显示前先做一次 median 滤波滤掉孤立噪点,衍射图样会干净得多。

6. 验证技巧:用峰值位置与条纹间距反推模拟是否正确

写完全部代码,第一件事不是欣赏漂亮的衍射图,而是量化验证。衍射模拟最怕的是“图看着像但数值全错”。我的习惯做法是提取光强曲线的峰值坐标,和理论公式比对。

对单缝衍射,暗纹位置满足 asinφ = mλ,在小角度近似下就是 y = mλD/a。跑完程序后,用 findpeaks 或直接找局部极小值,把前几个暗纹的位置打印出来:

[~, idx] = findpeaks(-B); % 找光强的局部极小 y_dark = ys(idx(1:min(5,length(idx)))); % 前几个暗纹位置 y_theory = (1:length(y_dark))' * lam * D / a; disp([y_dark, y_theory]);

如果误差超过 2%,优先检查 ym 范围和采样点数。ym 取得太大,条纹被压缩到曲线末端,峰值定位误差变大;ym 取得太小,暗纹超出显示范围,根本找不到。

双缝干涉的验证相对简单,条纹间距理论值是 Δy = lam*D/a,直接测量相邻亮纹间距和理论值对比。光栅则验证缺级:如果 d/a = 5,观察第 ±5 级主极大是否消失,这个现象比数值更直观,一旦看到缺级,说明衍射因子和干涉因子的耦合写对了。

还有一招我经常用:把衍射积分法跑出来的单缝结果和傅里叶方法跑出来的结果叠在同一张图里。两条曲线在夫琅和费区应该几乎重合,如果出现系统性偏移,往往是小角度近似被突破了;如果形状不同,多半是某个参数单位不一致,例如波长写成了纳米(500)而不是米(500e-9)。

这个检查习惯帮我避免了很多“看起来对、实际错”的尴尬。从那以后,我每次写完模拟代码都会强制走一遍峰值对比流程,先验证物理量再画图,宁可多花五分钟,也不愿意交作业时被一眼看出数值对不上。希望这篇文章里拆的代码和踩坑记录,能帮你的衍射模拟一次跑通。

本文还有配套的精品资源,点击获取

需要专业的网站建设服务?

联系我们获取免费的网站建设咨询和方案报价,让我们帮助您实现业务目标

立即咨询