简介这份资源面向学习大学物理光学课程、需要动手验证衍射理论的学生与教师用 MATLAB 编程把抽象的光学衍射现象转化为可视化的光强分布图。内容覆盖两条技术路线一是衍射积分方法对单缝衍射、杨氏双缝干涉和黑白光栅衍射逐一给出完整代码与成像结果二是傅立叶变换方法借助 fft2、fftshift 等函数模拟圆孔、矩孔、三角孔、单缝、黑白光栅及正弦光栅的夫琅和费衍射并讨论缺级与干涉因子调制等规律。压缩包共 1 个 doc 文档约 247KB以文字讲解配合可运行脚本的形式组织便于对照复现。目前已有 793 人学习适合希望用简捷直观方式加深对衍射原理与图像理解、并积累光学仿真代码的读者参考。1. 光学衍射 MATLAB 模拟两套方法怎么选、能跑出什么如果你手头有一份“光学衍射 MATLAB”代码包打开发现里面是单缝、双缝、光栅、圆孔、矩孔一堆脚本第一反应大概率是这些图到底靠不靠谱参数改一个会不会直接翻车。这份资源解决的正是这个问题——它用两套互相印证的思路把夫琅和费衍射算出来一套走衍射积分把孔径拆成离散相干光源做复振幅叠加另一套走傅立叶变换直接把屏函数做二维 FFT 取功率谱。两套方法跑同一组孔径图样能对上你才敢拿它去讲课、做课设或者验证实验数据。适合光学课程设计的学生、需要快速出衍射图样的工程验证人员以及想把物理公式落成可调参数脚本的从业者。下面按“先立公式、再跑代码、最后避坑”的顺序拆。2. 衍射积分法把孔径拆成离散光源再叠加衍射积分法的核心不是背公式而是想清楚“屏幕上某点的复振幅是孔径上所有面元发出的球面波在该点相干叠加的结果”。资源里单缝被处理成 np 个分立相干光源每个光源到屏上 P 点有一段光程差把所有光源的 cos 和 sin 分量分别累加平方和就是该点光强。这个离散化思路的好处是直观、可调缝宽、缝数、观察距离都是显式参数代价是采样数不够时条纹会失真这一点后面避坑章会细说。2.1 单缝衍射离散光源数与观察屏采样单缝脚本里最容易被忽略的是两个采样量孔径方向 np 和屏幕方向 ny。np 决定缝被切成多少份ny 决定屏幕上有多少个采样点。资源给的 np51、ny51 是能出图的下限条纹边缘会有明显锯齿。我一般会把 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和Dym会跟着自动缩放这是这份脚本设计得比较顺手的地方。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前一项是单缝衍射因子后一项是多缝干涉因子两者相乘就得到“干涉条纹被衍射包络调制”的经典图样。资源里N2其实是双缝特例改成 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; % beta0 处 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);这里必须提醒一个数值坑beta0时sin(N*beta)/sin(beta)是 0/0MATLAB 会返回 NaN中央主极大直接变黑。资源原代码没处理这个点我加了B(isnan(B))1兜底因为该极限等于 N。d5*a决定了缺级位置第 5、10 级主极大会被衍射因子压掉改d/a的比值就能移动缺级这是验证光栅公式最直接的手段。N越大主极大越锐但循环里没有向量化n1001 时跑起来仍然很快不用担心性能。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(x0.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 至少取 201ny 至少取 401且 ny 要远大于预期条纹数乘以每条纹采样点改完再看剖面是否光滑。4.2 中央主极大变黑或出现 NaN现象光栅图样中心一个黑点剖面在 y0 处断掉。原因sin(N*beta)/sin(beta)在 beta0 处是 0/0MATLAB 返回 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结果几乎均匀看不到明显衍射结构。原因输入图是彩色三通道或背景不是 0imread出来是 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/21, :)); % 中心行 xlim([-5e-3, 5e-3]);参数上dx必须和生成孔径时的归一化坐标对应上否则换算出来的坐标整体缩放曲线形状对但横轴是错的。xlim取 ±5 mm 是 1 mm 缝、1 m 距离下前三四个旁瓣的范围按实际参数调。交叉验证时还要注意积分法用的是小角近似sinphi y/DFFT 法是严格变换大角度下两者本来就会有偏差比对窗口别开太大。几个我踩过的边界孔径图尺寸不是 2 的幂时 FFT 会补零频谱出现额外调制看着像“新衍射结构”其实是补零伪影正弦光栅用cos1生成时灰度范围是 0 到 2保存 BMP 前要除以 2 归一化否则读回来被截断成 0/1 二值光栅衍射斑从三个变成一堆光栅脚本里N既是缝数又在别处当灰度级 255 用变量名冲突过一次后来我把灰度级统一改成NG才没再翻车。从那以后我每次改完参数都强制走一遍“积分法剖面 FFT 中心行叠图”的交叉验证图对不上就先查采样和坐标换算不急着怀疑物理公式。希望帮到你。本文还有配套的精品资源点击获取