简介这份资源面向光学、物理及光电专业的学生与教师提供用MATLAB模拟光学衍射实验的完整文档帮助读者通过编程直观理解衍射与干涉的物理原理。压缩包内共1个doc文件约247KB内容以文字说明配合可运行代码为主涵盖衍射积分与傅立叶变换两种模拟思路。资源具体展开单缝衍射、杨氏双缝干涉、黑白光栅衍射的衍射积分实现并给出圆孔、矩孔、三角孔、单缝、黑白光栅及正弦光栅夫琅和费衍射的傅立叶变换模拟方法每部分均附有可直接运行的MATLAB程序与光强分布图说明。读者可据此复现衍射图样、观察干涉条纹受衍射因子调制及缺级现象并借助fft2、fftshift等函数掌握频谱分析方法。目前已有791人学习适合作为课程实验辅助或自学参考。1. 光学衍射 MATLAB 模拟两种方法怎么选、能跑出什么图做光学课程设计或者物理实验仿真时很多人第一反应是找现成的光学仿真软件但那些工具要么收费、要么黑匣子参数改不动。这份资源走的是另一条路用 MATLAB 直接写代码从衍射积分和傅立叶变换两个角度把单缝衍射、杨氏双缝干涉、黑白光栅衍射、圆孔/矩孔/三角孔衍射以及正弦光栅夫琅和费衍射全部模拟出来。它的价值不在于代码有多长而在于把「公式怎么变成矩阵运算」这件事拆得很透——衍射积分法让你看清每个点光源的相干叠加过程傅立叶变换法则让你理解夫琅和费衍射为什么等价于屏函数的功率谱。适合正在学光学、需要交仿真作业、或者想用 MATLAB 验证衍射理论的人。下面按「先跑通再优化」的顺序把两种方法的实现细节和踩坑点讲清楚。2. 衍射积分法从单缝到光栅的逐点叠加实现2.1 为什么先学衍射积分而不是直接上 FFT衍射积分法的核心思路是把孔径分割成大量分立相干光源接收屏上任意一点的光场是所有这些子光源传播到该点的复振幅叠加。这个思路和物理教材里的惠更斯-菲涅耳原理完全对应代码写出来就是两层循环加三角函数求和没有黑盒操作。对于刚接触光学仿真的人来说这种「笨办法」反而最直观——你能看到 np 从 11 改到 101 时衍射条纹怎么从毛刺变得平滑也能理解为什么采样点数不够时高频条纹会失真。傅立叶变换法虽然代码更短、跑得更快但它跳过了中间物理过程直接给出结果。如果你连单缝衍射的光强分布公式都没推过直接看 FFT 的频谱图很容易一头雾水。所以我的建议是先用衍射积分法把单缝和双缝跑通建立「光强分布长什么样」的直觉再切到傅立叶方法去处理复杂孔径。2.2 单缝衍射采样点数 np 和 ny 怎么定单缝衍射的代码逻辑是把缝宽 a 分成 np 个点光源接收屏范围从 -ym 到 ym 取 ny 个采样点对每个屏上点计算所有子光源的相位差然后求 cos 和 sin 的累加和再平方。关键参数有三个波长 lam、缝宽 a、屏距 D。ym 的取值由 3lamD/a 决定这个系数 3 保证能看到中央主极大和两侧第一级暗纹。% 单缝衍射——衍射积分法 clear lam 500e-9; % 波长 500nm a 1e-3; % 缝宽 1mm D 1; % 屏距 1m 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; % 小角度近似 sin(phi) ≈ 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 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);这段代码里 np 和 ny 都取 51 是原文的默认值实际跑的时候你会发现 ny51 画出来的光强曲线有轻微锯齿。把 ny 提到 201 或 501曲线会平滑很多代价是循环变慢——但单缝这个规模完全感觉不到。np 的影响更隐蔽np 太小比如 11时子光源之间的相位差采样不足中央主极大的峰值会偏低旁瓣位置也会偏移。我一般把 np 设到 101 以上确保相位覆盖足够密。还有一个容易翻车的地方image(ym, ys, Br)这行里 ym 是标量ys 是向量image 函数会把 x 轴范围设成 [1, ym] 而不是 [-ym, ym]。如果你想让横轴对称显示得改成image([-ym, ym], ys, Br)或者直接用imagesc。原文代码能跑出图但坐标轴标注和实际物理量对不上这是血泪经验。2.3 杨氏双缝干涉相位差公式的符号陷阱杨氏双缝的代码比单缝短但有一个符号问题很容易搞反。两个缝到屏上 P 点的距离 r1 和 r2 分别用勾股定理算相位差 phi 2pi(r2-r1)/lam。注意这里 r2-r1 的顺序决定了条纹的左右对称性——如果写成 r1-r2干涉条纹会整体翻转虽然光强分布形状不变但如果你要叠加其他光场或者做相位分析符号错了后面全错。% 杨氏双缝干涉——衍射积分法 clear lam 500e-9; a 2e-3; % 双缝间距 2mm 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 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);这里B(i,:) sum(4*cos(phi/2).^2)其实有点冗余——phi 是标量cos(phi/2)^2 也是标量sum 加不加结果一样。原文这么写可能是为了和单缝的向量化风格保持一致。真正要注意的是 ym 的系数单缝用的是 3双缝用的是 5这是因为双缝干涉条纹间距由 lam*D/a 决定a 越小条纹越密需要更大的接收屏范围才能看到足够多的条纹。如果你把 a 改成 0.5e-3ym 不变的话屏上会挤满条纹看起来就是一片灰。2.4 黑白光栅衍射缺级现象的参数条件光栅衍射的公式比前两个复杂光强由单缝衍射因子和多缝干涉因子相乘得到I I0 * (sin(alpha)/alpha)^2 * (sin(Nbeta)/sin(beta))^2。其中 alpha piasinphi/lam 对应单缝衍射beta pid*sinphi/lam 对应多缝干涉N 是缝数d 是光栅常数。% 黑白光栅衍射——衍射积分法 clear lam 500e-9; N 2; % 缝数改成 5 或 10 看主极大变锐 a 2e-4; % 缝宽 0.2mm D 5; % 屏距 5m d 5*a; % 光栅常数 5 倍缝宽 ym 2*lam*D/a; xs ym; n 1001; % 光栅条纹密采样点要够 ys linspace(-ym, ym, n); 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);这段代码里最值得玩味的是 d 5a 这个设定。当 d/a 5 时干涉主极大的第 5 级、第 10 级……刚好落在单缝衍射的暗纹位置出现缺级。你把 d 改成 3a缺级就出现在第 3 级、第 6 级。这是验证光栅理论最直接的方式——改一个参数图上的亮线就少几条。另外 N2 时其实还是双缝干涉把 N 改成 10 或 20主极大明显变锐次极大被压下去这才是光栅该有的样子。采样点数 n 这里取 1001 是有道理的光栅条纹比单缝密得多n 太小的话高频干涉条纹会混叠图上会出现莫尔条纹一样的假信号。我试过 n201主极大位置还对但次极大已经糊成一片了。3. 傅立叶变换法把衍射屏当矩阵处理3.1 为什么夫琅和费衍射等于屏函数的功率谱傅立叶变换法的物理基础是夫琅和费衍射场的复振幅分布正比于衍射屏透射率函数的傅立叶变换光强就是它的模平方。这意味着你不需要逐点算相位差直接把屏函数扔给 fft2 就行。MATLAB 的 fft2 输出是复数矩阵直流分量在四个角上用 fftshift 把它挪到中心再用 abs 取模、imshow 显示就能看到衍射图样。这个方法的优势在于处理复杂孔径——圆孔、矩孔、三角孔、字母、甚至手绘图案只要你能生成对应的二值或灰度矩阵fft2 都能算。缺点是分辨率受矩阵尺寸限制而且边界效应需要额外处理。3.2 从图像文件到频谱图imread 和 fft2 的衔接原文的傅立叶变换代码读的是一张外部 BMP 图像路径写死在E:\1.bmp。实际用的时候你得先准备一张黑白孔径图或者用 MATLAB 直接生成矩阵。下面这段代码把两种方式都覆盖了% 傅立叶变换法——从图像文件读入 clear a imread(aperture.bmp); % 替换成你的孔径图路径 if size(a,3) 3 a rgb2gray(a); % 如果是彩色图转灰度 end a double(a); % 转 double 避免整数溢出 figure(1) imshow(a, []) title(衍射屏) afft fft2(a); % 二维离散傅立叶变换 aabs abs(afft); % 取模得到振幅谱 aabss fftshift(aabs); % 直流分量移到中心 figure(2) imshow(aabss, []) colormap(gray); colorbar title(衍射图样频谱) figure(3) plot(aabss) colormap(gray); figure(4) meshc(aabss) % 三维网格看光强分布 maxx1 max(max(aabss));这里有几个细节原文没展开但实际会卡住人。第一imread读进来的图像如果是 uint8 类型直接做 fft2 会因为整数运算丢失精度必须先double(a)。第二如果图像是 RGB 三通道fft2 会报错或者对每个通道分别变换得先转灰度。第三imshow(aabss, [])里的空方括号是让 MATLAB 自动缩放显示范围不加的话高频弱信号全黑。第四fftshift只对振幅谱做位移相位谱如果也要看得单独处理。3.3 圆孔、矩孔、三角孔用矩阵生成孔径不依赖外部图片的话可以直接用矩阵运算生成孔径。圆孔就是判断每个像素到中心的距离是否小于半径矩孔就是行列范围判断三角孔稍微麻烦一点但也就是几个线性不等式。% 生成圆孔孔径并做傅立叶变换 N 512; [x, y] meshgrid(linspace(-1, 1, N)); r sqrt(x.^2 y.^2); aperture double(r 0.3); % 半径 0.3 的圆孔 figure(1) imshow(aperture, []) title(圆孔) afft fft2(aperture); aabss fftshift(abs(afft)); figure(2) imshow(aabss, []) colormap(gray); title(圆孔衍射图样)把r 0.3改成abs(x) 0.3 abs(y) 0.2就是矩孔改成y -0.2 y 0.2 x -0.3 x 0.3 (y x 0.2) (y -x - 0.2)就是三角孔。矩阵尺寸 N 建议至少 512256 的话圆孔边缘会有明显锯齿衍射图样上会出现十字形伪影——这是矩形采样网格带来的不是物理效应。3.4 正弦光栅先造屏再衍射正弦光栅的透射率是余弦函数原文的做法是先画一张灰度图保存成 BMP再读进来做 FFT。其实可以一步到位% 正弦光栅——直接生成屏函数并做傅立叶变换 clear xm 10*pi; N 512; xs linspace(-xm, xm, N); [X, Y] meshgrid(xs, xs); B cos(X) 1; % 正弦光栅透射率 Br B/2 * 255; figure(1) imshow(Br, []) colormap(gray); title(正弦光栅) afft fft2(Br); aabss fftshift(abs(afft)); figure(2) imshow(aabss, []) colormap(gray); title(正弦光栅衍射图样)理论上正弦光栅的夫琅和费衍射只有三个衍射斑0 级和 ±1 级。跑出来如果看到更多亮斑通常是两个原因一是灰度量化到 255 级引入了高次谐波二是矩阵边界不连续导致频谱泄漏。把 N 加大到 1024 或者给屏函数加一个平滑窗比如高斯窗可以缓解。4. 避坑与排查跑不出图、图不对、跑得慢怎么办4.1 现象程序报错「索引超出矩阵维度」原因衍射积分法里B(i,:)在循环前没有预分配MATLAB 每次循环都在动态扩维如果 ny 和 np 设得不对或者循环变量 i 和数组维度对不上就会报这个错。更隐蔽的情况是ys和yp长度不一致时alpha的维度会出问题。解决在循环前加B zeros(ny, 1);预分配。另外检查linspace的第三个参数是否和循环上限一致。单缝代码里B(i,:)的第二个维度是 1但写B(i,:)也没错MATLAB 会自动处理。4.2 现象衍射图样全黑或者全白原因imshow或image的显示范围没设对。傅立叶变换后的频谱直流分量比高频分量大好几个数量级如果直接imshow(aabss)不加空方括号MATLAB 默认按 0-255 显示高频部分全黑。衍射积分法里如果Br没有归一化到 0-255image也会显示异常。解决imshow(aabss, [])让 MATLAB 自动缩放。或者手动做对数变换imshow(log(1 aabss), [])这样弱信号也能看见。衍射积分法里确保Br (B/max(B))*N这行没漏。4.3 现象光栅衍射图上看不到缺级原因d/a 的比值没设对。缺级条件是 d/a m/n其中 m 是干涉级次n 是衍射暗纹级次。原文用 d 5*a缺级出现在 5、10、15……如果你把 a 改了但忘了同步改 d缺级位置就跑了。另外 N 太小比如 N2时干涉主极大很宽缺级看起来不明显。解决固定 d/a 整数比如 3 或 5。把 N 提到 10 以上主极大变锐后缺级一目了然。如果还看不到检查ym范围是否覆盖到了缺级对应的角度。4.4 现象傅立叶变换法跑大矩阵时内存爆了原因fft2 对 N×N 矩阵的复杂度是 O(N²logN)N2048 时复数矩阵占 2048×2048×16 字节 ≈ 64MB加上中间变量和 fftshift 的副本内存占用翻倍。如果同时开多个 figure 且没关MATLAB 的图形缓存也会吃内存。解决用clear清理不用的变量close all关掉旧图窗。N 不需要太大512 或 1024 对大多数衍射模拟够用了。如果确实需要高分辨率可以用single类型代替double内存减半精度对显示来说足够。4.5 现象正弦光栅衍射出现多余亮斑原因灰度量化误差和边界不连续。Br B/2 * 255把连续余弦量化到 256 级引入了高次谐波。矩阵边界处光栅相位不连续FFT 的周期性假设导致频谱泄漏。解决把 N 从 512 提到 1024量化误差的相对影响变小。或者在屏函数两端加一个平滑窗比如B (cos(X) 1) .* exp(-X.^2/50)让边界处透射率趋近于零。这样衍射斑会稍微变宽但杂散亮斑明显减少。5. 进阶玩法用 FFT 方法批量扫描参数并验证物理规律跑通基本代码之后真正有意思的是用傅立叶变换法做参数扫描。衍射积分法每改一个参数都要重跑循环而 FFT 方法可以把孔径生成、变换、取模封装成函数然后在一个循环里批量处理不同孔径尺寸、不同光栅常数、不同缝数把结果拼成一张对比图。我一般会写一个diffract_fft函数输入是孔径矩阵输出是归一化频谱然后主脚本里改参数调用它。function spec diffract_fft(aperture) % 输入孔径矩阵double 类型 % 输出归一化振幅谱已 fftshift afft fft2(aperture); aabs abs(afft); aabss fftshift(aabs); spec aabss / max(aabss(:)); end用这个函数可以快速验证几个规律。第一圆孔衍射的艾里斑半径与孔径半径成反比——把r 0.3改成r 0.15中心亮斑明显变大。第二矩孔衍射的中央亮斑在 x 和 y 方向的宽度分别与孔径在对应方向的尺寸成反比。第三两个方孔干涉会在连线中垂线方向出现等间距暗纹三个方孔按三角形摆放会出现类似蜂巢的结构。这些在原文里都提到了但用函数化写法可以几行代码就扫完。还有一个我常用的技巧把meshc的三维网格图旋转一下视角从侧面看光强分布的峰谷比二维灰度图更容易判断主极大和次极大的比值。对于光栅衍射三维图能清楚看到主极大像一排尖峰次极大像矮丘缺级位置直接就是平的。从那以后我每次做光学仿真都强制先跑一遍单缝衍射验证代码路径没问题再切到复杂孔径。这个习惯帮我省了很多调试时间——单缝是最简单的测试用例如果单缝都跑不对后面的光栅和傅立叶变换肯定有问题。希望帮到你。本文还有配套的精品资源点击获取