简介：这份资源面向学习大学物理光学课程、需要动手验证衍射理论的学生与教师，用 MATLAB 编程把抽象的光学衍射现象转化为可视化的光强分布图。内容覆盖两条技术路线：一是衍射积分方法，对单缝衍射、杨氏双缝干涉和黑白光栅衍射逐一给出完整代码与成像结果；二是傅立叶变换方法，借助 fft2、fftshift 等函数模拟圆孔、矩孔、三角孔、单缝、黑白光栅及正弦光栅的夫琅和费衍射，并讨论缺级与干涉因子调制等规律。压缩包共 1 个 doc 文档，约 247KB，以文字讲解配合可运行脚本的形式组织，便于对照复现。目前已有 793 人学习，适合希望用简捷直观方式加深对衍射原理与图像理解、并积累光学仿真代码的读者参考。

1. 光学衍射 MATLAB 模拟：两套方法怎么选、能跑出什么

如果你手头有一份“光学衍射 MATLAB”代码包，打开发现里面是单缝、双缝、光栅、圆孔、矩孔一堆脚本，第一反应大概率是：这些图到底靠不靠谱，参数改一个会不会直接翻车。这份资源解决的正是这个问题——它用两套互相印证的思路把夫琅和费衍射算出来：一套走衍射积分，把孔径拆成离散相干光源做复振幅叠加；另一套走傅立叶变换，直接把屏函数做二维 FFT 取功率谱。两套方法跑同一组孔径，图样能对上，你才敢拿它去讲课、做课设或者验证实验数据。适合光学课程设计的学生、需要快速出衍射图样的工程验证人员，以及想把物理公式落成可调参数脚本的从业者。下面按“先立公式、再跑代码、最后避坑”的顺序拆。

2. 衍射积分法：把孔径拆成离散光源再叠加

衍射积分法的核心不是背公式，而是想清楚“屏幕上某点的复振幅，是孔径上所有面元发出的球面波在该点相干叠加的结果”。资源里单缝被处理成 np 个分立相干光源，每个光源到屏上 P 点有一段光程差，把所有光源的 cos 和 sin 分量分别累加，平方和就是该点光强。这个离散化思路的好处是直观、可调，缝宽、缝数、观察距离都是显式参数；代价是采样数不够时条纹会失真，这一点后面避坑章会细说。

2.1 单缝衍射：离散光源数与观察屏采样

单缝脚本里最容易被忽略的是两个采样量：孔径方向 np 和屏幕方向 ny。np 决定缝被切成多少份，ny 决定屏幕上有多少个采样点。资源给的 np=51、ny=51 是能出图的下限，条纹边缘会有明显锯齿。我一般会把 np 提到 201、ny 提到 401 再出正式图。

clear; lam = 500e-9; % 波长，500 nm a = 1e-3; % 缝宽 1 mm D = 1; % 缝到屏距离 1 m ym = 3*lam*D/a; % 屏幕半宽，取到第3级暗纹附近 ny = 401; % 屏幕采样点数，越大条纹越平滑 ys = linspace(-ym, ym, ny); np = 201; % 缝上离散光源数 yp = linspace(0, a, 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 N = 255; Br = (B/max(B))*N; subplot(1,2,1); image(ym, ys, Br); colormap(gray(N)); axis image; subplot(1,2,2); plot(B, ys); xlabel('相对光强'); ylabel('屏幕坐标 y');

这段代码的逻辑是：外层循环遍历屏幕采样点，内层用向量化方式一次算出所有面元的相位，再分别累加实部和虚部。ym = 3*lam*D/a这个半宽不是随便取的，它对应单缝衍射第 3 级暗纹的位置，取太小看不到完整旁瓣，取太大屏幕大部分是黑的。np和ny要匹配，np 太小相当于孔径欠采样，会出现伪条纹；ny 太小则条纹被“跳采样”，峰谷位置都抓不准。改参数时优先动a和D，ym会跟着自动缩放，这是这份脚本设计得比较顺手的地方。

2.2 杨氏双缝：相位差公式与光强归一化

双缝脚本走的是解析相位差路线，不再离散孔径，而是直接算两缝到屏上点的距离差。资源里r1、r2用平方和开方算，phi = 2*pi*(r2-r1)/lam，光强4*cos(phi/2).^2。这里有个细节：B(i,:)=sum(4*cos(phi/2).^2)对单点求 sum 其实是冗余的，因为 phi 已经是标量，但保留它不影响结果，属于原作者的写法习惯。

clear; lam = 500e-9; a = 2e-3; % 双缝间距 2 mm D = 1; ym = 5*lam*D/a; % 取到第5级干涉条纹附近 n = 401; ys = linspace(-ym, ym, n); B = zeros(n, 1); 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; % 除以理论最大值4做归一化 subplot(1,2,1); image(ym, ys, Br); colormap(gray(N)); axis image; subplot(1,2,2); plot(B, ys); xlabel('相对光强'); ylabel('屏幕坐标 y');

参数上要注意a是缝间距不是缝宽，双缝干涉图样的条纹间距由lam*D/a决定，改a条纹疏密立刻变化。归一化用/4.0是因为两束等强光干涉理论最大光强是单束的 4 倍，不除的话灰度图会整体偏亮、对比度被压掉。ym取 5 倍条纹间距是为了让图上能看到足够多周期，方便数条纹验证公式。这套解析写法比离散积分快得多，但只适用于理想双缝，缝宽带来的衍射包络它没有体现，想要“干涉被衍射调制”的效果得回到积分法或傅立叶法。

2.3 黑白光栅：缝宽因子与多缝因子的乘积

光栅脚本是这份资源里公式最完整的一段，光强写成(sin(alpha)./alpha).^2 .* (sin(N*beta)./sin(beta)).^2，前一项是单缝衍射因子，后一项是多缝干涉因子，两者相乘就得到“干涉条纹被衍射包络调制”的经典图样。资源里N=2其实是双缝特例，改成 10、20 才能看到细锐的主极大和缺级。

clear; lam = 500e-9; N = 10; % 缝数，改成10以上才像光栅 a = 2e-4; % 缝宽 0.2 mm d = 5*a; % 光栅常数，5倍缝宽 D = 5; ym = 2*lam*D/a; n = 1001; ys = linspace(-ym, ym, n); B = zeros(n, 1); 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 B(isnan(B)) = 1; % beta=0 处 0/0，按极限置1 B1 = B/max(B); NC = 255; Br = (B/max(B))*NC; subplot(1,2,1); image(ym, ys, Br); colormap(gray(NC)); axis image; subplot(1,2,2); plot(B1, ys); xlabel('归一化光强'); ylabel('屏幕坐标 y');

这里必须提醒一个数值坑：beta=0时sin(N*beta)/sin(beta)是 0/0，MATLAB 会返回 NaN，中央主极大直接变黑。资源原代码没处理这个点，我加了B(isnan(B))=1兜底，因为该极限等于 N。d=5*a决定了缺级位置，第 5、10 级主极大会被衍射因子压掉，改d/a的比值就能移动缺级，这是验证光栅公式最直接的手段。N越大主极大越锐，但循环里没有向量化，n=1001 时跑起来仍然很快，不用担心性能。

3. 傅立叶变换法：屏函数做 FFT 再取功率谱

傅立叶变换法的物理依据是夫琅和费衍射场复振幅正比于屏函数的傅立叶变换，强度就是功率谱。资源里的流程是：读入一张黑白或灰度孔径图，fft2做二维离散傅立叶变换，abs取模得振幅谱，fftshift把直流分量从四角搬到中心，最后imshow显示。这条路子的优势是孔径形状随便画，圆孔、矩孔、三角孔、字母符号都能直接跑，不用为每种孔径推一遍积分公式；代价是分辨率受输入图像尺寸限制，且要理解 FFT 的坐标映射关系。

3.1 从图像到频谱：fft2、abs、fftshift 三步走

先把最小可运行流程跑通，再谈孔径设计。下面这段把读图、变换、移位、显示串起来，并额外画出中心行剖面，方便和积分法的曲线对照。

clear; a = imread('aperture.bmp'); % 黑白孔径图，背景0孔径1 a = double(a); if ndims(a) == 3 a = rgb2gray(a); % 彩色图先转灰度 end afft = fft2(a); % 二维离散傅立叶变换 aabs = abs(afft); % 取模得振幅谱 aabss = fftshift(aabs); % 直流分量移到中心 aabss = aabss / max(aabss(:)); % 归一化便于显示 figure(1); imshow(a, []); title('孔径屏函数'); figure(2); imshow(aabss, []); colormap(gray); colorbar; title('衍射图样'); figure(3); plot(aabss(round(end/2), :)); title('中心行光强剖面'); figure(4); meshc(aabss); title('频谱三维视图');

逻辑说明：fft2输出的是复数矩阵，直接imshow看不到有意义的东西，必须先abs。fftshift不是可选项，不做的话零频在左上角，图样看起来是“四角亮中间暗”，和实际衍射对不上。归一化那行是我加的，原代码靠imshow(...,[])自动缩放，但画剖面和 mesh 时没有统一量纲，归一化后四张图才对得上。aperture.bmp建议用 256×256 或 512×512 的灰度图，尺寸取 2 的幂，FFT 速度快且频率栅格规整。

3.2 孔径设计：圆孔、矩孔、三角孔与正弦光栅

傅立叶法的真正价值在于孔径可以“画”出来。圆孔用距离判据，矩孔用矩形区域，三角孔用半平面交集，正弦光栅直接按余弦函数生成灰度。下面这段一次生成四种孔径并批量出图，省得反复改文件。

N = 512; [x, y] = meshgrid(linspace(-1, 1, N)); % 圆孔 circle = double(x.^2 + y.^2 <= 0.3^2); % 矩孔 rect = double(abs(x) <= 0.3 & abs(y) <= 0.15); % 三角孔 tri = double(y <= 0.3 & y >= -0.3 & abs(x) <= (0.3 - y)/2); % 正弦光栅 grating= 0.5*(cos(2*pi*8*x) + 1); apertures = {circle, rect, tri, grating}; names = {'圆孔','矩孔','三角孔','正弦光栅'}; for k = 1:4 A = apertures{k}; F = fftshift(abs(fft2(A))); F = F / max(F(:)); figure('Name', names{k}); subplot(1,2,1); imshow(A, []); title(['孔径：' names{k}]); subplot(1,2,2); imshow(F, []); title('衍射图样'); end

参数说明：linspace(-1,1,N)把孔径归一化到 [-1,1] 方便用相对尺寸描述，0.3、0.15这些是归一化半径和半宽，改它们就是改孔径大小。正弦光栅的频率8决定衍射斑间距，周期数越多三个衍射斑分得越开。三角孔用半平面交集写法比逐点判断简洁，但边界会有锯齿，想要平滑边缘可以把判据改成带抗锯齿的软边。四种孔径跑完你会发现圆孔是艾里斑、矩孔是 sinc 十字、三角孔出现六重对称、正弦光栅只有 ±1 级和 0 级三个斑，这些和理论对得上，说明 FFT 路线是可信的。

3.3 多孔阵列与复合光栅：观察调制与缺级

资源后半段提到“两个方孔连线中垂线出现等间距暗纹”“三角形摆放出现蜂巢结构”“两方向光栅叠加”，这些用傅立叶法做特别快，因为孔径就是矩阵加法。下面用两个方孔和双向光栅演示。

N = 512; [x, y] = meshgrid(linspace(-1, 1, N)); % 两个方孔，中心分别在 (-0.4,0) 和 (0.4,0) sq1 = double(abs(x+0.4) <= 0.1 & abs(y) <= 0.1); sq2 = double(abs(x-0.4) <= 0.1 & abs(y) <= 0.1); two_holes = sq1 + sq2; % 双向光栅叠加 g1 = 0.5*(cos(2*pi*10*x) + 1); g2 = 0.5*(cos(2*pi*10*y) + 1); cross_grating = g1 .* g2; for k = 1:2 if k == 1, A = two_holes; t = '双方孔'; else, A = cross_grating; t = '双向光栅'; end F = fftshift(abs(fft2(A))); F = F/max(F(:)); figure('Name', t); subplot(1,2,1); imshow(A, []); title(t); subplot(1,2,2); imshow(F, []); title('衍射图样'); end

双方孔的衍射图样是单孔衍射包络乘以双孔干涉因子，中垂线上出现等间距暗纹，这正是“孔径排布决定干涉调制”的直观体现。双向光栅用两个一维光栅相乘，频谱变成二维点阵，如果两个方向频率不同，点阵会拉成矩形栅格。这里.*不能写成+，乘法对应两个方向透射率叠加，加法会变成两个独立光栅的强度相加，图样完全不同，这是新手最容易写错的一行。

4. 避坑与排查：五个让图样对不上的真实原因

4.1 条纹出现锯齿或伪峰

现象：单缝或光栅剖面曲线毛刺明显，旁瓣位置和理论值差一两个像素。原因：孔径离散采样 np 或屏幕采样 ny 不足，离散光源之间的相位跳变被放大。解决：np 至少取 201，ny 至少取 401，且 ny 要远大于预期条纹数乘以每条纹采样点；改完再看剖面是否光滑。

4.2 中央主极大变黑或出现 NaN

现象：光栅图样中心一个黑点，剖面在 y=0 处断掉。原因：sin(N*beta)/sin(beta)在 beta=0 处是 0/0，MATLAB 返回 NaN。解决：循环后加B(isnan(B)) = N^2或归一化前B(isnan(B)) = 1（归一化后极限为 1），两种写法等价，关键是别让 NaN 进入max和imshow。

4.3 FFT 图样四角亮、中心暗

现象：傅立叶法出的图看起来像“边框发光”，和教材上的中心亮斑相反。原因：忘了fftshift，零频分量还在矩阵左上角。解决：abs之后必须fftshift，顺序不能反；如果先 shift 再 abs，复数移位结果不对。

4.4 读入图像后频谱全是高频噪声

现象：fft2结果几乎均匀，看不到明显衍射结构。原因：输入图是彩色三通道或背景不是 0，imread出来是 uint8 三维矩阵，直接fft2会把通道当空间维度。解决：先rgb2gray再double，并确认孔径区域为 1、背景为 0；用画图工具生成 BMP 时关掉抗锯齿，否则边缘灰度会引入额外高频。

4.5 两种方法结果对不上

现象：同一圆孔，积分法算的艾里斑半径和 FFT 法差一截。原因：积分法里ym和D的单位、FFT 法里像素尺寸和物理尺寸没有建立映射，两者坐标系不是同一套。解决：FFT 法只用于定性看图样结构，定量比对时把像素坐标换算成lam*D/pixel_size再和积分法曲线叠图；或者干脆用积分法做定量、FFT 法做孔径创意验证，别硬凑同一坐标轴。

5. 进阶：把两套方法交叉验证成习惯

真正让这份资源从“能出图”变成“敢用”的，是拿两套方法互相校验。我的做法是固定一组物理参数，先用积分法算单缝中心行剖面，再用 FFT 法对同样缝宽的矩形孔径做变换，把 FFT 结果的像素坐标按lam*D/(N*dx)换算成物理坐标，两条曲线叠在一张图上。第一旁瓣位置误差在 2% 以内，我才认为脚本没写错。下面这段是换算的核心，dx是孔径图一个像素对应的物理尺寸。

% 假设孔径图 N x N，像素尺寸 dx（米），波长 lam，距离 D N = 512; dx = 2e-6; lam = 500e-9; D = 1; F = fftshift(abs(fft2(aperture))); F = F / max(F(:)); fx = (-N/2:N/2-1) / (N*dx); % 空间频率轴 x_phys = fx * lam * D; % 换算到屏幕物理坐标 plot(x_phys, F(N/2+1, :)); % 中心行 xlim([-5e-3, 5e-3]);

参数上dx必须和生成孔径时的归一化坐标对应上，否则换算出来的坐标整体缩放，曲线形状对但横轴是错的。xlim取 ±5 mm 是 1 mm 缝、1 m 距离下前三四个旁瓣的范围，按实际参数调。交叉验证时还要注意积分法用的是小角近似sinphi = y/D，FFT 法是严格变换，大角度下两者本来就会有偏差，比对窗口别开太大。

几个我踩过的边界：孔径图尺寸不是 2 的幂时 FFT 会补零，频谱出现额外调制，看着像“新衍射结构”其实是补零伪影；正弦光栅用cos+1生成时灰度范围是 0 到 2，保存 BMP 前要除以 2 归一化，否则读回来被截断成 0/1 二值光栅，衍射斑从三个变成一堆；光栅脚本里N既是缝数又在别处当灰度级 255 用，变量名冲突过一次，后来我把灰度级统一改成NG才没再翻车。从那以后我每次改完参数都强制走一遍“积分法剖面 + FFT 中心行叠图”的交叉验证，图对不上就先查采样和坐标换算，不急着怀疑物理公式。希望帮到你。

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