简介这份资源面向光学、物理及光电专业的学生与教师以及需要开展衍射实验仿真验证的科研人员提供用MATLAB编程模拟光学衍射的完整文档。内容围绕衍射积分与傅立叶变换两条主线展开前者通过计算接收屏光强分布模拟单缝衍射、杨氏双缝干涉与黑白光栅衍射后者借助fft2、fftshift等函数对屏函数做功率谱分析复现圆孔、矩孔、三角孔、单缝、黑白光栅及正弦光栅的夫琅和费衍射图样并附有可直接运行的代码与光强分布图。资源包为1个doc文档约247KB结构紧凑便于按章节查阅与复现。目前已有791人学习适合作为课程实验、仿真作业或自学衍射原理的参考材料帮助读者在编程实践中直观理解干涉条纹受衍射因子调制及缺级等现象。1. 光学衍射 MATLAB 模拟从物理公式到灰度图的完整实现路径如果你正在做光学课程设计或者需要快速验证一个孔径的夫琅和费衍射图样长什么样又不想真的去搭光路、调激光器、买光栅那这份 MATLAB 光学衍射模拟资源就是为你准备的。它用两种完全不同的思路——衍射积分和傅立叶变换——把单缝衍射、杨氏双缝干涉、黑白光栅衍射、圆孔/矩孔/三角孔衍射以及正弦光栅衍射全部跑通最终输出灰度图和光强分布曲线。代码量不大但覆盖了从近场到远场、从一维到二维的典型场景。我第一次拿到这份代码时最直观的感受是原来课本上那些“光强分布公式”真的可以几行矩阵运算就变成屏幕上的条纹。它适合光学、物理、电子信息类专业的本科生做仿真作业也适合想用 MATLAB 快速验证衍射规律的工程师。2. 两种模拟路线衍射积分与傅立叶变换的架构拆解2.1 衍射积分方法把孔径拆成相干点源叠加衍射积分方法的核心思想很朴素把单缝或光栅看成由许多个分立相干点光源组成接收屏上任意一点的光场是所有这些点源在该点产生复振幅的叠加。对于单缝代码把缝宽a分成np51个点每个点对应一个相位alpha2*pi*yp*sinphi/lam然后分别累加余弦和正弦分量最后取模平方得到光强。这种“离散化求和”的做法本质上就是基尔霍夫衍射公式的数值近似。% 单缝衍射衍射积分法 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; % 近似小角度 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; % 映射到 0~255 灰度 subplot(1,2,1) image(ym, ys, Br); % 左图灰度条纹 colormap(gray(N)); subplot(1,2,2) plot(B, ys); % 右图光强随位置变化曲线这段代码里lam、a、D是三个最关键的物理参数。改a会直接改变衍射条纹的宽度——缝越窄衍射越明显改D则影响观察屏上的尺度。np控制离散精度一般取 51 到 101 之间就足够平滑再大对结果改善有限但计算量线性增长。ym的取值决定了你看到几级条纹取3*lam*D/a大约能看到中央主极大和两侧第一、第二级暗纹。杨氏双缝干涉的代码思路类似但不再对缝内积分而是直接计算两个缝到屏上某点的光程差用phi2*pi*(r2-r1)/lam得到相位差光强按4*cos(phi/2)^2分布。这里有个容易忽略的点双缝代码里B(i,:)sum(4*cos(phi/2).^2)其实是对单点求值sum在这里没有实际聚合作用但结果正确。黑白光栅衍射则引入了多缝干涉因子(sin(N*beta)/sin(beta))^2和单缝衍射因子(sin(alpha)/alpha)^2的乘积能清晰看到缺级现象——当衍射因子为零时干涉主极大消失。2.2 傅立叶变换方法屏函数的功率谱就是衍射场傅立叶变换方法走的是另一条路夫琅和费衍射场的复振幅分布正比于衍射屏透射率函数的傅立叶变换。MATLAB 的fft2可以直接对图像矩阵做二维离散傅立叶变换取模后得到振幅谱再用fftshift把零频分量移到中心最后imshow显示。这条路线的好处是你不需要推导任何衍射积分公式只要把孔径形状画成黑白图像读进来做 FFT 就行。% 傅立叶变换法读取孔径图像计算衍射图样 clear a imread(aperture.bmp); % 黑白孔径图白色透光、黑色挡光 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) % 三维网格观察强度分布这里的关键操作是fftshift。不做这一步零频分量在矩阵四角显示的频谱是“分裂”的看起来不直观。abs取模得到的是振幅谱如果想看光强分布应该取模的平方。imread读入的图像可以是自己用画图工具生成的黑白 BMP也可以用代码动态生成——比如正弦光栅那段就是用cos(xs)1生成灰度渐变条纹再保存成 BMP 的。两种方法各有适用场景衍射积分法物理意义明确适合理解公式和参数影响傅立叶变换法通用性强任意形状孔径只要画出黑白图就能算适合快速探索复杂孔径的衍射花样。3. 从零跑通环境配置、参数调优与结果验证3.1 环境准备与文件组织这份代码对 MATLAB 版本要求不高R2015a 之后的版本都能跑。不需要额外工具箱fft2、fftshift、imread、imshow都是基础函数。我一般会把不同实验的脚本分开存成独立.m文件比如single_slit.m、double_slit.m、grating.m、fft_diffraction.m避免变量互相污染。傅立叶变换法需要一个孔径图像文件建议放在脚本同目录下用相对路径读取这样换电脑不用改路径。% 生成正弦光栅图像并保存供傅立叶变换法使用 clear all xm 10*pi; ys xm; xs linspace(-xm, xm, 500); B cos(xs) 1; % 正弦透过率范围 0~2 N 255; Br B/2*N; % 归一化到 0~255 image(xs, ys, Br); colormap(gray(N)); imwrite(uint8(Br), sine_grating.bmp); % 保存为灰度 BMP这段代码生成的正弦光栅其傅立叶变换理论上只有三个衍射斑零级和正负一级跑出来如果看到更多级次说明图像量化或边界处理引入了高频分量。保存时用uint8转换避免imwrite对 double 类型做额外缩放。3.2 参数调优改什么、看什么衍射积分法的三个核心参数是波长lam、特征尺寸a缝宽或光栅常数、接收距离D。它们共同决定衍射角sinphi ≈ lam/a。如果你想看到更宽的衍射条纹把a减小一个数量级想看到更密集的干涉条纹把光栅常数d增大。黑白光栅代码里N2表示只取两个缝改成N10就能看到多缝干涉的细锐主极大。傅立叶变换法里衍射图样的尺度由孔径图像的像素尺寸和物理尺寸的对应关系决定。如果你用 500×500 像素表示 1cm×1cm 的孔径那频谱中心到第一暗环的像素距离就对应衍射角。实际做的时候我习惯先用小尺寸图像比如 128×128快速验证确认无误后再换成高分辨率图。提示fft2的输出默认零频在左上角一定要fftshift后再显示否则你会看到四个对称的“角”容易误判为高阶衍射。3.3 结果验证和理论值对一对单缝衍射第一暗纹位置理论值y lam*D/a。用lam500e-9、a1e-3、D1算出来是 0.5mm。跑完代码后在光强曲线上找第一个零点应该就在 0.5mm 附近。杨氏双缝干涉条纹间距理论值Δy lam*D/a这里a是双缝间距。黑白光栅的缺级条件d/a k/n比如d5a时第 5、10、15 级主极大缺失。这些验证点能帮你确认代码没有量纲错误或索引错位。傅立叶变换法的验证更直接圆孔衍射应该得到艾里斑中心亮斑加同心暗环矩孔衍射得到十字形亮纹三角孔衍射得到六角对称花样。如果跑出来不对称检查孔径图像是否居中、是否有非零边界。4. 避坑与排查衍射模拟中容易翻车的五个细节4.1 现象光强曲线出现负值或 NaN原因sin(alpha)./alpha在alpha0时产生0/0。虽然 MATLAB 对sin(0)/0返回NaN但实际极限是 1。光栅衍射代码里alpha和beta都可能为零。解决在计算前加一个极小量eps或者用逻辑索引单独处理零点。我一般写成alpha alpha (alpha0)*eps;这样既不影响非零值又避免 NaN。4.2 现象傅立叶变换法得到的衍射图样一片白或一片黑原因imshow(aabss, [])里的[]会自动缩放显示范围但如果图像本身动态范围极大零频分量比高频高好几个数量级高频细节会被压扁。另外如果读入的图像是 RGB 三通道fft2会对三个通道分别变换结果维度不对。解决读入后先rgb2gray转灰度再用log(1aabss)做对数压缩显示这样高频和低频都能看清。或者手动设置显示范围imshow(aabss, [0, max(aabss(:))/100])。4.3 现象衍射积分法跑出来的条纹不对称原因linspace(0, a, np)从 0 开始而不是从-a/2开始导致点源分布整体偏移了半个缝宽。对于单缝这只会引入一个整体相位因子不影响光强但对于双缝或多缝如果两个缝的采样起点不一致干涉条纹会偏移。解决把缝内采样改成linspace(-a/2, a/2, np)保证对称。双缝代码里直接用解析光程差公式不存在这个问题。4.4 现象fftshift之后图像中心不在正中间原因当图像尺寸为奇数时fftshift的零频位置在(N1)/2而不是N/21视觉上会偏一格。这是离散傅立叶变换的固有特性不是 bug。解决如果对中心位置敏感把图像尺寸补零到偶数。或者接受这一格偏移在后续分析时用ceil(size/2)定位中心。4.5 现象正弦光栅衍射出现多余级次原因生成正弦光栅时用了cos(xs)1但xs的范围是-10*pi到10*pi两端不连续相当于给正弦波加了矩形窗。矩形窗的傅立叶变换是 sinc 函数会引入额外高频。解决把范围取成2*pi的整数倍并且首尾值相等。或者用hann窗做渐变过渡。我一般直接取xs linspace(-10*pi, 10*pi, 500)然后检查首尾是否接近不接近就微调点数。5. 二次开发从固定孔径到参数化批量仿真5.1 把脚本改造成函数支持批量参数扫描原始代码都是脚本变量直接写在开头。如果你想研究“缝宽从 0.1mm 到 2mm 变化时衍射条纹如何演化”手动改参数再跑几十次太慢。我一般会把核心计算封装成函数输入物理参数输出光强分布矩阵。function [B, ys] single_slit_intensity(lam, a, D, ny, np) % 单缝衍射光强分布衍射积分法 % lam: 波长(m), a: 缝宽(m), D: 接收距离(m) % ny: 屏上采样点数, np: 缝内点源数 % 返回 B: 光强分布(ny×1), ys: 屏上坐标 ym 3*lam*D/a; ys linspace(-ym, ym, ny); yp linspace(-a/2, a/2, np); % 对称采样 B zeros(ny, 1); for i 1:ny sinphi ys(i)/D; alpha 2*pi*yp*sinphi/lam; sumcos sum(cos(alpha)); sumsin sum(sin(alpha)); B(i) (sumcos^2 sumsin^2)/np^2; end end封装之后批量扫描就变成两层循环外层遍历a内层调用函数把结果存到元胞数组或三维矩阵里。用imagesc把不同a对应的条纹拼成一张图能直观看到衍射角随缝宽减小而增大的趋势。5.2 用 App Designer 做一个交互式衍射模拟器如果你需要给不熟悉 MATLAB 的人演示可以拖几个滑块控件波长、缝宽、屏距离实时刷新衍射图样。核心回调函数就是上面那个single_slit_intensity把滑块的值传进去算完用imagesc或plot更新坐标轴。傅立叶变换法那边可以做一个“画孔径”的面板让用户用鼠标画任意形状然后实时 FFT 显示衍射花样。这个改造的工作量主要在界面布局计算部分几乎不用动。5.3 扩展到二维任意孔径从图像到衍射的通用流程傅立叶变换法最大的优势是通用性。你不需要为每种孔径推导公式只要生成对应的黑白图像。我常用的流程是先用poly2mask或insertShape生成孔径掩膜再做fft2。比如要模拟三角形孔用poly2mask指定三个顶点坐标要模拟字母“A”的衍射直接读入文字图像二值化。唯一需要注意的是图像尺寸要足够大避免边界截断引入伪影。一般取 1024×1024 或 2048×2048计算量在现代笔记本上完全可接受。5.4 一个我踩过的坑归一化方式影响对比度早期我做傅立叶变换衍射时直接imshow(abs(fft2(a)), [])结果图样对比度很低暗纹几乎看不见。后来发现是因为零频分量太大自动缩放把高频压没了。改成log(1abs(fft2(a)))再显示细节立刻清晰。另一个坑是fft2不做归一化不同尺寸图像的频谱幅度不可比。如果要做定量分析记得除以numel(a)。从那以后我每次做傅立叶变换衍射都强制走一遍“转灰度 → 补零到偶数尺寸 → fft2 → fftshift → log 压缩 → 显示”的流程再也没出现过一片白或一片黑的情况。希望这些经验能帮到你少走几个弯路。本文还有配套的精品资源点击获取