简介：这份资源面向光学、物理及光电专业的学生与教师，以及需要开展衍射实验仿真验证的科研人员，提供用MATLAB编程模拟光学衍射的完整文档。内容围绕衍射积分与傅立叶变换两条主线展开：前者通过计算接收屏光强分布，模拟单缝衍射、杨氏双缝干涉与黑白光栅衍射；后者借助fft2、fftshift等函数对屏函数做功率谱分析，复现圆孔、矩孔、三角孔、单缝、黑白光栅及正弦光栅的夫琅和费衍射图样，并附有可直接运行的代码与光强分布图。资源包为1个doc文档，约247KB，结构紧凑，便于按章节查阅与复现。目前已有791人学习，适合作为课程实验、仿真作业或自学衍射原理的参考材料，帮助读者在编程实践中直观理解干涉条纹受衍射因子调制及缺级等现象。
1. 光学衍射 MATLAB 模拟：从物理公式到灰度图的完整实现路径
如果你正在做光学课程设计，或者需要快速验证一个孔径的夫琅和费衍射图样长什么样，又不想真的去搭光路、调激光器、买光栅，那这份 MATLAB 光学衍射模拟资源就是为你准备的。它用两种完全不同的思路——衍射积分和傅立叶变换——把单缝衍射、杨氏双缝干涉、黑白光栅衍射、圆孔/矩孔/三角孔衍射以及正弦光栅衍射全部跑通，最终输出灰度图和光强分布曲线。代码量不大，但覆盖了从近场到远场、从一维到二维的典型场景。我第一次拿到这份代码时，最直观的感受是：原来课本上那些“光强分布公式”真的可以几行矩阵运算就变成屏幕上的条纹。它适合光学、物理、电子信息类专业的本科生做仿真作业，也适合想用 MATLAB 快速验证衍射规律的工程师。
2. 两种模拟路线：衍射积分与傅立叶变换的架构拆解
2.1 衍射积分方法：把孔径拆成相干点源叠加
衍射积分方法的核心思想很朴素：把单缝或光栅看成由许多个分立相干点光源组成，接收屏上任意一点的光场是所有这些点源在该点产生复振幅的叠加。对于单缝，代码把缝宽a分成np=51个点，每个点对应一个相位alpha=2*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大约能看到中央主极大和两侧第一、第二级暗纹。
杨氏双缝干涉的代码思路类似，但不再对缝内积分，而是直接计算两个缝到屏上某点的光程差，用phi=2*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增大。黑白光栅代码里N=2表示只取两个缝，改成N=10就能看到多缝干涉的细锐主极大。
傅立叶变换法里，衍射图样的尺度由孔径图像的像素尺寸和物理尺寸的对应关系决定。如果你用 500×500 像素表示 1cm×1cm 的孔径，那频谱中心到第一暗环的像素距离就对应衍射角。实际做的时候，我习惯先用小尺寸图像（比如 128×128）快速验证，确认无误后再换成高分辨率图。
提示：
fft2的输出默认零频在左上角，一定要fftshift后再显示，否则你会看到四个对称的“角”，容易误判为高阶衍射。
3.3 结果验证：和理论值对一对
单缝衍射第一暗纹位置理论值y = lam*D/a。用lam=500e-9、a=1e-3、D=1算出来是 0.5mm。跑完代码后，在光强曲线上找第一个零点，应该就在 0.5mm 附近。杨氏双缝干涉条纹间距理论值Δy = lam*D/a，这里a是双缝间距。黑白光栅的缺级条件d/a = k/n，比如d=5a时，第 5、10、15 级主极大缺失。这些验证点能帮你确认代码没有量纲错误或索引错位。
傅立叶变换法的验证更直接：圆孔衍射应该得到艾里斑，中心亮斑加同心暗环；矩孔衍射得到十字形亮纹；三角孔衍射得到六角对称花样。如果跑出来不对称，检查孔径图像是否居中、是否有非零边界。
4. 避坑与排查：衍射模拟中容易翻车的五个细节
4.1 现象：光强曲线出现负值或 NaN
原因：sin(alpha)./alpha在alpha=0时产生0/0。虽然 MATLAB 对sin(0)/0返回NaN，但实际极限是 1。光栅衍射代码里alpha和beta都可能为零。
解决：在计算前加一个极小量eps，或者用逻辑索引单独处理零点。我一般写成alpha = alpha + (alpha==0)*eps;，这样既不影响非零值，又避免 NaN。
4.2 现象：傅立叶变换法得到的衍射图样一片白或一片黑
原因：imshow(aabss, [])里的[]会自动缩放显示范围，但如果图像本身动态范围极大（零频分量比高频高好几个数量级），高频细节会被压扁。另外，如果读入的图像是 RGB 三通道，fft2会对三个通道分别变换，结果维度不对。
解决：读入后先rgb2gray转灰度，再用log(1+aabss)做对数压缩显示，这样高频和低频都能看清。或者手动设置显示范围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的零频位置在(N+1)/2而不是N/2+1，视觉上会偏一格。这是离散傅立叶变换的固有特性，不是 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(1+abs(fft2(a)))再显示，细节立刻清晰。另一个坑是fft2不做归一化，不同尺寸图像的频谱幅度不可比。如果要做定量分析，记得除以numel(a)。
从那以后我每次做傅立叶变换衍射，都强制走一遍“转灰度 → 补零到偶数尺寸 → fft2 → fftshift → log 压缩 → 显示”的流程，再也没出现过一片白或一片黑的情况。希望这些经验能帮到你，少走几个弯路。
本文还有配套的精品资源，点击获取
