简介：这份资源面向学习大学物理光学、需要动手验证衍射与干涉理论的学生和教师，用 MATLAB 编程把抽象的光学现象转化为可视化的光强分布图。内容围绕两条主线展开：一是衍射积分方法，涵盖单缝衍射、杨氏双缝干涉与黑白光栅衍射的完整程序；二是傅立叶变换方法，通过 fft2、fftshift 等函数模拟简单孔径、黑白光栅及正弦光栅的夫琅和费衍射，并讨论缺级现象与孔径尺寸对图样的影响。资源包共 1 个 doc 文档，约 247KB，以文字讲解配合可直接运行的代码片段为主，便于对照复现。目前已有 792 人学习。读者可从中获得两类模拟思路的公式推导与实现细节，理解衍射因子与干涉因子的调制关系，并借助现成脚本快速搭建自己的光学仿真实验，加深对物理原理与图像特征的理解。
1. 光学衍射 MATLAB 模拟：两套方法把物理实验搬进命令行
做光学课程设计那会儿，最头疼的不是公式推不出来，而是实验室的钠灯一开、狭缝一调，眼睛贴着观察屏看了半天，条纹位置和理论值对不上，也不知道是缝宽量错了还是光路没调平。后来我把整个实验搬进 MATLAB，用衍射积分和傅立叶变换两套方法各跑一遍，屏幕上直接出灰度图和光强曲线，参数改一个数字就能看到条纹怎么变——这才算把单缝衍射、杨氏双缝干涉、黑白光栅衍射这些现象真正吃透。这份资源就是围绕这两条技术路线展开的 MATLAB 工程代码集合，适合正在学大学物理光学、要做仿真作业，或者想拿衍射模拟做二次开发练手的人。它不依赖任何工具箱之外的第三方库，纯脚本就能跑，拿到手改几个参数就能出图。
2. 两套模拟方法的技术选型：衍射积分 vs 傅立叶变换
2.1 为什么同一个物理现象要写两遍
单缝衍射、杨氏双缝、光栅衍射这些现象，物理上都是惠更斯-菲涅耳原理的推论，但落到代码里，实现路径可以完全不同。衍射积分方法走的是“暴力合成”路线：把缝面拆成 np 个分立相干光源，对屏上每个采样点，逐个光源算光程差、累加复振幅，最后取模平方得到光强。这条路直观、物理意义清晰，缺点是双重循环，采样点一多就慢。傅立叶变换方法走的是“频域捷径”：夫琅和费衍射场的强度分布等于屏函数的功率谱，直接把孔径图像做二维 FFT、取模、fftshift 移中，就能得到衍射图样。这条路快、能处理任意形状孔径，缺点是对近场（菲涅耳衍射）不适用，而且初学者容易搞不清直流分量为什么跑到四个角上去了。
两套方法放在一起，恰好覆盖了从“理解物理”到“快速出图”的完整需求。我一般建议先跑积分法把单缝和双缝的条纹位置对上理论值，建立信心，再用傅立叶法去玩圆孔、矩孔、三角孔甚至字母形状的衍射花样——后者才是真正能出“好看图”的部分。
2.2 衍射积分法的核心循环结构
以单缝衍射为例，代码的主循环长这样：
clear lam = 500e-9; % 波长，单位米 a = 1e-3; % 缝宽，单位米 D = 1; % 缝到屏的距离，单位米 ym = 3*lam*D/a; % 屏上观察范围半宽，取3倍中央条纹半宽 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); % 右图：光强随屏坐标的曲线这段代码的关键在alpha = 2*pi*yp*sinphi/lam这一行。yp是缝面上每个分立光源的位置，sinphi是屏上该点对应的衍射角正弦，两者相乘就是光程差，再乘2*pi/lam转成相位。sumcos和sumsin分别累加所有光源的复振幅实部和虚部，最后(sumcos^2 + sumsin^2)就是合成复振幅的模平方，即光强。除以np^2是为了归一化，避免光源数一多光强数值爆掉。
参数怎么调：np从 51 加到 201，条纹会变得更平滑，但循环次数翻四倍；ny控制纵向分辨率，太小的话灰度图会出现明显的锯齿；ym取3*lam*D/a是经验值，刚好把中央主极大和两侧第一暗纹包进来，想多看几级次就改成5*或7*。这里有个血泪经验：yp的起点必须是 0 而不是-a/2，因为单缝积分里我们只关心缝内各点到屏点的相对光程差，起点平移不影响最终光强分布，但如果你拿这段代码去改双缝，起点就得对称取了。
2.3 傅立叶变换法的图像处理链路
傅立叶方法的核心链路是“读图 → 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) % 三维网格，看强度起伏fft2出来的是复数矩阵，abs取模得到振幅谱，但此时直流分量（零频）位于矩阵的四个角上，直接imshow会看到中心暗、四角亮的怪图。fftshift的作用是交换第一、三象限和第二、四象限，把零频挪到正中间，这样频谱图才和实际衍射图样的视觉习惯一致。imshow(aabss, [])里的空方括号是自动缩放显示范围，不加的话灰度可能全黑或全白。
这里有个容易翻车的点：imread读进来的 BMP 如果是 RGB 三通道，fft2会直接报错，得先rgb2gray转灰度。另外图像尺寸最好是 2 的幂次（256×256、512×512），FFT 计算效率最高，非 2 的幂虽然也能跑，但速度会明显下降。我一般会在读图后加一句a = imresize(a, [512 512])统一尺寸，省得后面调参时被不同图片的尺寸差异坑到。
3. 目录结构与运行方式：从脚本到可复用函数
3.1 资源里的文件组织
这份资源本质是一组 MATLAB 脚本的集合，按方法分成两大块。衍射积分部分包含单缝、双缝、黑白光栅三个独立脚本，每个脚本都是“清空变量 → 设参数 → 循环算光强 → 双图显示”的固定套路。傅立叶部分包含一个主处理脚本（读图、FFT、四图显示）和若干孔径生成脚本（圆孔、矩孔、三角孔、正弦光栅），孔径脚本负责用矩阵运算或cos函数生成灰度图并保存为 BMP，主处理脚本再读入这些 BMP 做变换。
实际使用时，我建议在本地建一个这样的目录：
diffraction_matlab/ ├── integral/ │ ├── single_slit.m │ ├── double_slit.m │ └── grating.m ├── fourier/ │ ├── aperture_circle.m │ ├── aperture_rect.m │ ├── aperture_triangle.m │ ├── sine_grating.m │ └── fft_process.m └── images/ ├── circle.bmp ├── rect.bmp └── sine.bmpintegral/下的脚本互相独立，直接run就能出图；fourier/下的孔径脚本先生成图片存到images/，fft_process.m里把imread的路径改成对应文件即可。这种拆分的好处是改一个孔径不用动主处理逻辑，主处理逻辑也不用关心孔径是怎么来的。
3.2 运行前的环境检查
MATLAB 版本方面，这些脚本用的都是基础函数（linspace、fft2、fftshift、imshow、meshc），R2014b 之后的版本都能跑，不需要额外工具箱。唯一要注意的是imshow在 Image Processing Toolbox 里，如果装的是纯 MATLAB 没带这个工具箱，可以换成imagesc加colormap(gray)，效果基本一样。
运行顺序上，傅立叶部分必须先跑孔径生成脚本、确认images/下有 BMP 文件，再跑fft_process.m。我见过有同学直接跑主处理脚本，结果imread报“文件不存在”，回头查了半天以为是路径写错了，其实是孔径图还没生成。衍射积分部分没有这个依赖，三个脚本随便挑一个跑都行。
3.3 参数扫描的批量跑法
单次出图只能看一个参数下的条纹，要做对比实验（比如缝宽从 0.5mm 变到 2mm 条纹怎么展宽），手动改参数重跑太慢。我一般会套一层for循环把参数扫一遍：
lam = 500e-9; D = 1; a_list = [0.5e-3, 1e-3, 2e-3]; % 三组缝宽 figure; for k = 1:length(a_list) a = a_list(k); ym = 3*lam*D/a; ys = linspace(-ym, ym, 201); np = 101; yp = linspace(0, a, np); B = zeros(length(ys), 1); for i = 1:length(ys) sinphi = ys(i)/D; alpha = 2*pi*yp*sinphi/lam; B(i) = (sum(cos(alpha))^2 + sum(sin(alpha))^2)/np^2; end subplot(3,1,k) plot(ys*1e3, B) % 横轴转毫米，方便读数 title(sprintf('a = %.1f mm', a*1e3)) xlabel('屏坐标 / mm'); ylabel('相对光强') end注意这里ym是随a变的，因为缝越窄衍射越强、观察范围要越大，固定ym的话窄缝的条纹会跑出画面。横轴乘1e3转成毫米是为了图上读数直观，不然坐标轴上一串10^-4的科学计数法看着累。这种批量跑法在写实验报告做参数对比时特别省事，一次出三张图直接贴进去。
4. 核心模块拆解：从单缝到光栅的公式落地
4.1 杨氏双缝：光程差直接算，不用拆光源
双缝的代码和单缝思路不同，它不拆分立光源，而是直接算两个缝到屏点的距离差：
lam = 500e-9; a = 2e-3; % 双缝间距 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,:) = 4*cos(phi/2)^2; % 双光束干涉光强 end N = 255; Br = (B/4.0)*N; subplot(1,2,1) image(ym, ys, Br); colormap(gray(N)); subplot(1,2,2) plot(B, ys);r1和r2用勾股定理算，a/2是每个缝到中心轴的距离。相位差phi出来后，双光束干涉的光强公式就是4*cos(phi/2)^2，这个 4 是两束光振幅相等时的最大光强倍数。除以 4 归一化再乘 255 映射灰度。这里ym取5*lam*D/a是因为双缝干涉条纹间距是lam*D/a，取 5 倍能看五条亮纹，比单缝的 3 倍宽一些。
有个细节值得说：r1和r2里的ys(i)是标量，但如果你把ys直接传进去做向量化运算（去掉for循环），r1和r2会变成向量，B(i,:)的写法就得改成B = 4*cos(phi/2).^2。向量化版本跑得快，但初学者容易在.^和^上翻车——矩阵幂和逐元素幂搞混了，结果要么报错要么出一堆复数。我建议先用循环版把物理逻辑跑通，确认条纹位置对了，再改向量化提速。
4.2 黑白光栅：两个因子的乘积
光栅衍射的公式比前两个复杂，是单缝衍射因子和多缝干涉因子的乘积：
lam = 500e-9; N = 2; % 缝数，改大能看到更细的干涉条纹 a = 2e-4; % 缝宽 D = 5; d = 5*a; % 光栅常数 ym = 2*lam*D/a; 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; end B1 = B/max(B); NC = 255; Br = (B/max(B))*NC; subplot(1,2,1) image(ym, ys, Br); colormap(gray(NC)); subplot(1,2,2) plot(B1, ys);alpha对应单缝衍射的包络，beta对应多缝干涉的细条纹，两者相乘就得到“干涉条纹被衍射包络调制”的效果。N是缝数，改成 5 或 10 能看到主极大变锐、次极大增多。d = 5*a意味着光栅常数是缝宽的 5 倍，这个比例决定了缺级出现在第 5 级——如果你把d改成3*a，缺级就挪到第 3 级，这是验证缺级现象最直接的办法。
n = 1001比前面几个脚本的 51、101 大得多，因为光栅的干涉条纹很细，采样点少了根本分辨不出来。代价是循环 1001 次，跑起来明显比单缝慢。如果嫌慢，可以把beta和alpha的运算向量化，把for循环去掉，速度能快一个数量级。
4.3 正弦光栅：先造图，再做 FFT
正弦光栅的玩法和其他几个不一样，它先用cos函数生成一张灰度图存成 BMP，再走傅立叶流程：
clear all xm = 10*pi; ys = xm; xs = linspace(-xm, xm, 500); B = cos(xs) + 1; % 正弦透射率，加1保证非负 N = 255; Br = B/2*N; % 归一化到0~255 image(xs, ys, Br); colormap(gray(N)); % 手动保存为 sine.bmp，或用 imwrite 写入cos(xs)+1把透射率从[-1,1]平移到[0,2]，再除以 2 归一化到[0,1]，乘 255 映射灰度。这张图存成 BMP 后丢进fft_process.m，出来的频谱只有三个亮点——零级和正负一级，和理论完全一致。这个例子特别适合用来验证 FFT 流程对不对：如果出来的不是三个点而是一团糊，那八成是fftshift忘了加，或者图像不是严格的灰度。
5. 避坑与排查：跑衍射模拟时最容易翻车的几件事
5.1 灰度图全黑或全白
现象：image或imshow出来的图要么一片黑，要么一片白，看不到条纹。
原因：光强矩阵B没有归一化就直接乘 255，数值远超 255 或者远小于 1，映射到灰度后全被截断。单缝代码里(sumcos^2+sumsin^2)/np^2的归一化如果漏了/np^2，光强会随np线性增长，np=51时数值能到几千。
解决：显示前统一做Br = (B/max(B(:)))*255，用max(B(:))而不是max(B)，后者对矩阵返回的是每列最大值组成的行向量，除出来维度不对。这个坑我踩过不止一次，尤其是从别人代码里复制片段时，归一化那行经常被顺手删掉。
5.2 FFT 频谱中心是暗的、四角是亮的
现象：fft2之后直接imshow(abs(afft))，看到中心暗、四角亮的图，和教科书上的衍射图样对不上。
原因：二维 FFT 的零频分量位于矩阵左上角（MATLAB 的索引习惯），视觉上就是四个角各有一块低频区，中心反而是高频。不fftshift的话，人眼会把四角误认为“中心”。
解决：取模后立刻fftshift，把零频挪到矩阵正中心。顺序不能反——先fftshift再abs也行，但先abs再fftshift更符合“先取振幅谱再移中”的物理直觉。另外fftshift只对abs后的实数矩阵做，别对复数矩阵直接移中再取模，虽然结果一样，但中间变量容易搞混。
5.3 光栅条纹采样不足，细条纹糊成一片
现象：光栅衍射的灰度图上看不到细密的干涉条纹，只有一团模糊的亮区。
原因：n太小。光栅干涉条纹的角间距是lam/d，d=5*a=1e-3时条纹间距约5e-4弧度，屏上对应D*5e-4 = 2.5e-3米。如果n=101、ym又是2*lam*D/a量级，每个条纹只分到一两个采样点，自然糊。
解决：把n加到 1001 以上，或者缩小ym只看中央几级。我一般先n=2001跑一遍确认条纹结构，再降到 501 出图，兼顾速度和清晰度。另外image的插值方式默认是最近邻，采样够了之后可以加imagesc让显示更平滑。
5.4 imread 读 RGB 图导致 fft2 报错
现象：fft2(a)报错“输入必须为二维矩阵”或“不支持 N 维数组”。
原因：BMP 文件是 24 位真彩色，imread读进来是M×N×3的三维数组，fft2只接受二维。
解决：读图后加if ndims(a)==3, a = rgb2gray(a); end，或者干脆在画图工具里存成 8 位灰度 BMP。我习惯在fft_process.m开头统一做一次维度检查，省得换一张图就报一次错。
5.5 循环版和向量化版结果对不上
现象：把for循环改成向量化运算后，条纹位置偏移或者光强分布变形。
原因：向量化时alpha从标量变成向量，sin(alpha)/alpha里的/如果写成/而不是./，MATLAB 会当成矩阵右除，结果完全错误。同理^2要写成.^2。
解决：所有逐元素运算统一用点运算符，.*、./、.^一个都不能漏。改完先拿max(abs(B_vec - B_loop))和循环版对比，差值在1e-10量级才算对。这个对比步骤我强烈建议保留，向量化提速很爽，但错了很难肉眼发现。
6. 二次开发：把脚本改成可配置的仿真函数
原始脚本都是“改参数 → 重跑”的模式，做几次对比实验还行，真要拿来做课程设计或者批量出图，就得改成函数。我的做法是把每个物理模型封装成function，参数用结构体传，返回光强矩阵和坐标轴：
function [I, y] = sim_single_slit(lam, a, D, ny, np) % 单缝衍射光强分布 % lam 波长(m), a 缝宽(m), D 屏距(m), ny 屏采样数, np 缝面光源数 ym = 3*lam*D/a; y = linspace(-ym, ym, ny); yp = linspace(0, a, np); I = zeros(ny, 1); for i = 1:ny alpha = 2*pi*yp*(y(i)/D)/lam; I(i) = (sum(cos(alpha))^2 + sum(sin(alpha))^2)/np^2; end I = I/max(I); end封装之后，参数扫描就变成一行arrayfun或者简单的for循环调函数，不用每次复制粘贴改数字。更进一步，可以把单缝、双缝、光栅统一到一个sim_diffraction(type, params)入口，type用字符串区分，内部switch分发——这样上层做 GUI 或者批量脚本时只认一个接口，扩展新模型也不用动调用方。
验证改造是否正确的办法很朴素：拿封装后的函数跑一组和原始脚本完全相同的参数，把两张光强曲线叠在一起看是否重合。我一般用max(abs(I_new - I_old))卡一个1e-8的阈值，超了就说明封装过程中哪里改错了。傅立叶部分同理，把fft_process.m改成function [S, S_shift] = sim_fft(aperture_path)，返回原始频谱和移中频谱，方便上层做对比显示。
还有一个值得做的改造是给傅立叶流程加一个“孔径生成器”函数，用矩阵运算直接画圆孔、矩孔、三角孔，而不是依赖外部 BMP 文件。比如圆孔就是(X.^2 + Y.^2) <= r^2生成逻辑矩阵再转double，几行代码的事，但省掉了“画图 → 存盘 → 读盘”的来回折腾，批量跑不同孔径时特别顺手。从那以后我每次做仿真类项目，都强制先把“参数 → 模型 → 输出”这条链路函数化再开始调参，不然脚本越复制越多，最后自己都分不清哪个版本是对的。希望帮到你。
本文还有配套的精品资源，点击获取
