简介这份资源面向光学课程学习者、物理实验教学人员及需要做衍射仿真的理工科学生用 MATLAB 编程复现光学衍射现象解决实验条件受限、抽象公式难以直观理解的问题。内容以衍射积分与傅立叶变换两条技术路线展开前者覆盖单缝衍射、杨氏双缝干涉与黑白光栅衍射后者借助 fft2、fftshift 等函数模拟圆孔、矩孔、三角孔、单缝、黑白光栅及正弦光栅的夫琅和费衍射并给出完整可运行代码与光强分布图。资源包共 1 个 doc 文档约 247KB集中呈现理论推导、程序清单与结果图示便于对照学习。已有 792 人学习下载适合作为课堂演示、课程设计或自学衍射光学的参考材料读者可据此快速搭建仿真脚本、理解衍射图样随孔径与光栅参数变化的规律。1. 光学衍射 MATLAB 模拟两套方法把物理实验搬进命令行做光学课程设计那会儿实验室的导轨、狭缝片和光屏永远不够用一组人调半天也未必能等到稳定的条纹。后来某导师丢过来一份 MATLAB 脚本说“你先在电脑上把条纹跑出来再去动真仪器”。这份资源就是干这个的用衍射积分和傅立叶变换两种方法把单缝衍射、杨氏双缝干涉、黑白光栅衍射、圆孔/矩孔/三角孔衍射、正弦光栅夫琅和费衍射全部用代码复现出来。它适合光学、物理电子、光电信息这类专业的本科生做仿真作业也适合需要快速验证衍射图样规律的工程开发者。整套东西不依赖任何工具箱纯 MATLAB 基础函数就能跑拿到手改几个参数就能看到条纹怎么变。2. 两套模拟路线衍射积分与傅立叶变换的架构分工这份资源的骨架其实很清晰它没有做成一个大而全的 GUI 工程而是按物理建模思路分成两条独立的技术路线每条路线下挂若干具体实验场景。理解这个分工后面改代码才不会乱。2.1 衍射积分路线从光程差累加到屏上光强衍射积分方法的核心思想是把孔径离散成若干个相干点光源对接收屏上每一个采样点累加所有点光源的复振幅贡献最后取模平方得到光强。这条路线的优点是物理图像直观你能清楚看到“缝被切成 np 份”这件事缺点是计算量随采样点数平方增长参数开大了会明显变慢。单缝衍射的脚本是这条路线里最基础的一个先看它的完整结构clear lam 500e-9; % 波长500nm 绿光 a 1e-3; % 缝宽 1mm D 1; % 屏距 1m 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);这段代码里几个参数值得说清楚。ym决定了你观察屏画多宽取3*lam*D/a刚好能包住中央主极大和两侧第一级暗纹改小了条纹会被裁掉。np是缝的离散精度51 个点对单缝够用但如果换成多缝光栅缝数乘以每缝点数会让循环变重这时候要么降np要么把内层循环向量化。B(i,:)那个/np^2是归一化因子保证中央光强接近 1方便后面统一映射灰度。杨氏双缝的脚本换了个思路不再离散缝而是直接算两束光到屏上每点的几何路程差clear lam 500e-9; a 2e-3; % 双缝间距 2mm D 1; ym 5*lam*D/a; % 双缝干涉条纹更密观察范围取 5 倍 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,:) sum(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);这里a是缝间距而不是缝宽ym的系数从 3 变成 5 是因为双缝条纹间距是lam*D/a比单缝中央宽度小观察范围要放大才能看到足够多级。B(i,:) sum(...)里的sum其实对单点没意义原代码保留它可能是为了和向量化写法兼容实际跑的时候可以去掉不影响结果。光栅衍射是这条路线里公式最完整的一个把单缝衍射因子和多缝干涉因子乘在一起clear lam 500e-9; N 2; % 缝数改大能看到更细的干涉峰 a 2e-4; % 缝宽 0.2mm D 5; d 5*a; % 光栅常数5 倍缝宽 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; B1 B/max(B); end 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);N是缝数原代码给 2 其实更像双缝改成 5 或 10 才能看到光栅特有的细密主极大。d 5*a这个比例决定了缺级位置当d/a为整数时第d/a级主极大正好落在单缝衍射暗纹上条纹会“缺”掉这是光栅衍射最值得观察的现象。alpha和beta在sinphi0处会出现 0/0MATLAB 会给出 NaN但image和plot对单个 NaN 不敏感实际图上看不出来严谨点可以加一句alpha(alpha0)eps。2.2 傅立叶变换路线把屏函数当信号做频谱分析第二条路线完全换了个视角夫琅和费衍射的强度分布等于孔径函数的功率谱所以直接把孔径图像读进来做二维 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) % 看中心截面的光强分布 figure(4) meshc(aabss) % 三维网格看频谱形态imread读进来的矩阵如果是 RGB 三通道fft2会报错常见做法是先rgb2gray再二值化。fftshift这一步不能省不做的话直流分量在四个角上图看起来是“四角亮中间暗”和实际衍射图样对不上。imshow(aabss, [])里的空方括号是自动缩放显示范围不加的话高动态范围的频谱会显示成一片黑。meshc那一步是锦上添花用来观察频谱的立体形态正式出图可以不要。正弦光栅的生成是这条路线的配套小工具先造一张余弦灰度图存成 BMP再喂给上面的 FFT 流程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));cos(xs)1保证灰度非负/2*N把范围压到 0~255。存图的时候注意用imwrite存成灰度 BMP不要用截图截图会引入抗锯齿和压缩噪声FFT 出来的频谱会多出不该有的高频分量。正弦光栅的频谱理论上只有三个亮点零级和正负一级如果跑出来一堆杂散斑八成是存图环节出了问题。3. 从零跑通环境准备、脚本组织与出图验证拿到这份资源后最省事的做法不是把代码一段段复制到命令窗口而是按实验场景拆成独立脚本文件每个文件对应一张图跑完一个看一个。下面按实际操作的顺序讲。3.1 环境与目录组织MATLAB 版本没有硬性要求R2014b 之后的都可以因为用到的fft2、fftshift、imshow、meshc都是基础函数不需要 Image Processing Toolbox 之外的额外工具箱。如果imshow报未定义检查一下是不是装的精简版补装 Image Processing Toolbox 即可。目录建议这样组织每个脚本独立可跑互不干扰optical_diffraction/ ├── integral/ │ ├── single_slit.m │ ├── double_slit.m │ └── grating.m ├── fourier/ │ ├── aperture_fft.m │ └── sine_grating_gen.m └── images/ ├── circle.bmp ├── rectangle.bmp └── sine_grating.bmpintegral下放衍射积分的三个脚本fourier下放 FFT 流程和光栅生成脚本images放孔径素材。这样改参数的时候不会误伤别的实验也方便对比两种方法对同一个物理场景的结果。3.2 逐脚本运行与参数调整先跑single_slit.m确认能出图。正常情况下左边是灰度条纹图中间亮两边暗右边是光强曲线中央主极大占大部分能量。如果图是纯黑或纯白检查Br的归一化B/max(B)在B全零时会出 NaN说明循环里alpha算错了回头查lam、a、D的单位是不是统一成了米。再跑double_slit.m应该看到等间距的明暗条纹间距约lam*D/a。把a从2e-3改成1e-3条纹间距应该翻倍这是验证代码正确性最快的方法。如果改了a条纹没变说明你改的是别处的同名变量MATLAB 工作区里残留了旧值养成每个脚本开头clear的习惯就是为了防这个。grating.m跑之前先把N从 2 改成 5否则看到的是双缝而不是光栅。改完应该能看到主极大变细、主极大之间出现次级亮纹。再把d从5*a改成4*a缺级位置会从第 5 级移到第 4 级这是光栅衍射最直观的验证点。FFT 路线先跑sine_grating_gen.m生成sine_grating.bmp存到images目录然后把aperture_fft.m里的imread路径改成这张图。跑完看figure(2)应该只有中心一个亮斑加左右各一个对称亮斑共三个。如果看到一条亮线或者一片噪点回到存图那步查灰度映射。3.3 出图与结果对照两种方法对同一个物理场景的结果可以互相印证。比如单缝衍射积分法给出的是sinc^2形状的光强曲线FFT 法对一张矩形孔径图做变换取中心截面也应该得到sinc^2。把两条曲线画在同一张图上对比如果峰位和零点位置对得上说明两套代码都写对了。% 对比积分法与 FFT 法的单缝光强截面 % B_integral 来自 single_slit.m 的 B % B_fft 来自 aperture_fft.m 中心行取模 figure plot(ys/max(ys), B_integral/max(B_integral), b-, ... linspace(-1,1,length(B_fft)), B_fft/max(B_fft), r--) legend(积分法,FFT法) xlabel(归一化屏坐标) ylabel(归一化光强)横坐标做了归一化处理因为两种方法的采样网格不一样直接画会在不同尺度上。纵坐标也归一化到峰值 1方便比较形状。如果两条曲线在主极大宽度上差很多检查 FFT 法里孔径图的物理尺寸和积分法里的a是否对应FFT 的频谱尺度依赖像素数和像素物理尺寸这一点原代码没有显式处理二次开发时需要补上。提示跑 FFT 路线时孔径图的分辨率不要低于 256×256太小的图频谱分辨率不够条纹会糊成一片。4. 核心代码拆解参数怎么改、现象怎么读这一章把几个关键脚本的参数影响讲透方便你拿到手之后按自己的实验需求调整而不是只会跑默认值。4.1 波长与缝宽对条纹宽度的影响单缝衍射中央主极大的半角宽度约lam/a屏上宽度约2*lam*D/a。这三个参数里lam改大条纹变宽a改大条纹变窄D改大条纹变宽。代码里ym 3*lam*D/a就是跟着这个关系走的你改lam或a之后ym会自动适配不用手动调观察范围。但有个细节ny固定为 51 时如果ym变得很大屏上采样间隔变粗条纹边缘会出现锯齿。常见做法是把ny提到 201 或 501代价是循环变慢。对单缝来说 201 足够对光栅这种条纹细密的场景n给到 1001 是合理的原代码光栅脚本里n1001就是这个原因。4.2 光栅常数与缝数的组合效应光栅衍射图样的两个特征——主极大间距和主极大宽度——分别由d和N控制。主极大位置满足d*sinphi k*lam所以d越大主极大越密主极大宽度反比于NN越大主极大越细。代码里beta pi*d*sinphi/lam控制间距sin(N*beta)/sin(beta)控制锐度。改N的时候注意N增大后sin(N*beta)/sin(beta)在beta接近pi整数倍时会出现数值不稳定分母趋零导致 Inf。实际跑的时候如果看到曲线突然冒尖把beta加一个极小量eps兜底beta pi*d*sinphi/lam eps; % 防止 sin(beta)0 除零这个eps对结果影响可以忽略但能避免图上出现孤立的异常尖峰出图干净很多。4.3 FFT 频谱的物理尺度换算FFT 法最容易翻车的地方是频谱坐标没有物理意义。fft2出来的矩阵索引是像素序号不是角度或屏坐标。要把它和积分法的结果对上需要做换算频率分辨率df 1/(Nx*dx)其中Nx是孔径图横向像素数dx是像素物理尺寸。衍射角正弦sinphi lam*f屏坐标y D*sinphi。原代码没有做这步换算所以figure(3)的横轴是索引而不是物理坐标。二次开发时补上这段两种方法的结果才能真正叠在一张图上比较[Ny, Nx] size(a); dx 1e-5; % 假设像素物理尺寸 10um按实际改 df 1/(Nx*dx); f (-Nx/2:Nx/2-1)*df; % 移频后的频率轴 sinphi lam*f; y_axis D*sinphi; % 对应的屏坐标dx这个值必须和生成孔径图时的实际尺寸一致否则换算出来的屏坐标是错的。如果孔径图是从画图工具截的没有物理尺寸信息那就只能做定性比较不能做定量对照。4.4 灰度映射与显示范围所有脚本最后都有一句Br (B/max(B))*N把光强映射到 0~255。这一步看着简单但max(B)取的是全局最大值如果图里有异常尖峰比如前面说的除零归一化后其他条纹会被压得很暗看起来像只有中心一个点。稳妥的做法是用分位数归一化避开极端值Bmax prctile(B(:), 99.5); % 取 99.5 分位作为归一化上限 Br min(B/Bmax, 1) * N;prctile需要 Statistics Toolbox没有的话用sort手动取分位也行。这样即使有个别异常点整体条纹对比度仍然正常。5. 避坑与排查五个真实会卡住你的问题5.1 图像全黑或全白现象跑完脚本image出来的图要么纯黑要么纯白看不到条纹。原因Br归一化时max(B)为 0 或 Inf导致整个矩阵被压成同一个值。常见触发点是alpha或beta在零点除零产生 NaNmax忽略 NaN 后返回正常值但B里混了 NaNimage显示时整行失效。解决在循环里对alpha、beta加eps兜底或者在归一化前用B(isnan(B)) 0清掉 NaN。跑完先max(B(:))和min(B(:))看一眼确认数值范围正常再画图。5.2 改了参数图没变化现象把a从1e-3改成5e-4条纹宽度纹丝不动。原因MATLAB 工作区残留了上一次运行的变量新脚本里如果没写clear或者你是在命令窗口逐段粘贴而不是整文件运行旧值会覆盖新值。另一种可能是你改的变量名和脚本里实际用的不一致比如脚本里是a你改成了slit_width。解决每个脚本开头强制clear养成整文件F5运行的习惯不要逐段执行。改完参数后在脚本里加一句disp(a)确认生效。5.3 FFT 频谱四角亮中间暗现象aperture_fft.m跑出来figure(2)是四个角亮、中间暗。原因忘了fftshift或者fftshift作用在了错误的矩阵上。fft2的直流分量在矩阵左上角不移频的话显示出来就是四角亮。解决确认aabss fftshift(abs(fft2(a)))这一串顺序正确先fft2再abs再fftshift。如果顺序反了fftshift作用在复数矩阵上虽然也能跑但后续取模的相位关系会乱。5.4 正弦光栅频谱出现多余亮斑现象理论上只有三个衍射斑实际跑出来五六个甚至一片散斑。原因存图环节引入了噪声。用截图工具保存 BMP 会带抗锯齿用 JPEG 中转会有压缩伪影这些都会在频谱里变成高频杂散。另外cos(xs)1如果量化到 8 位灰度时截断太狠也会产生谐波。解决用imwrite(Br, map, sine_grating.bmp)直接写灰度 BMP不要经过任何图像编辑软件。生成时把xs采样点提到 1000 以上减少量化误差。5.5 光栅脚本跑得特别慢现象grating.m里n1001时循环要跑十几秒甚至更久。原因外层 1001 次循环每次内部对alpha、beta做数组运算虽然单次很快但累计起来在低配机器上明显。如果再把n提到 5001时间线性增长。解决把外层循环向量化ys整个作为向量一次性算sinphi ys/D; alpha pi*a*sinphi/lam; beta pi*d*sinphi/lam eps; B (sin(alpha)./alpha).^2 .* (sin(N*beta)./sin(beta)).^2;向量化之后n5001也是毫秒级。原代码用循环可能是教学演示需要实际用的时候没必要逐点算。6. 二次开发把零散脚本收成一个可配置的仿真工具这套代码目前是“一个实验一个脚本”的形态跑通没问题但要做课程设计或者批量出图就有点散。我一般会做一层薄封装把所有参数收进一个结构体主函数根据场景名分发到不同计算分支输出统一的光强矩阵和坐标轴。这样改参数只改一处出图格式也统一。核心思路是把“算什么”和“怎么显示”分开function [I, y] diffraction_sim(cfg) % cfg.scene: single | double | grating | fft % cfg.lam, cfg.a, cfg.d, cfg.D, cfg.N, cfg.n 按场景取用 switch cfg.scene case single [I, y] calc_single_slit(cfg); case double [I, y] calc_double_slit(cfg); case grating [I, y] calc_grating(cfg); case fft [I, y] calc_fft_aperture(cfg); end I I / max(I(:)); % 统一归一化 endcfg里给默认值调用时只覆盖关心的字段。calc_*系列函数就是把前面拆过的脚本主体搬进去去掉clear和绘图部分只返回I和y。绘图统一放到外层想画灰度图就imagesc想画曲线就plot想导出就imwrite不用每个脚本重复写一遍。验证封装是否正确最直接的办法是拿封装后的结果和原脚本逐点对比。选单缝场景cfg用原脚本的默认参数跑出来I和原B做max(abs(I-B))应该在1e-12量级说明只是搬了位置没改逻辑。再换 FFT 场景和aperture_fft.m的中心截面比形状对得上就说明分发逻辑没问题。还有一个值得加的技巧是参数扫描。比如研究缝宽对中央条纹宽度的影响写个循环扫a从0.2e-3到2e-3每次调diffraction_sim把半高全宽记下来画曲线理论上应该反比于a。这种批量验证比单张图有说服力也方便写进报告。a_list linspace(0.2e-3, 2e-3, 20); fwhm zeros(size(a_list)); for k 1:numel(a_list) cfg.a a_list(k); [I, y] diffraction_sim(cfg); idx find(I 0.5, 1, last) - find(I 0.5, 1, first); fwhm(k) y(idx); end plot(a_list, fwhm, o-)这段跑完如果曲线是漂亮的反比关系说明整个封装链路是通的后面加新场景比如椭圆孔、双圆孔只需要补一个calc_*分支不用动主流程。血泪经验是封装的时候千万别顺手“优化”物理公式。我有一次觉得sin(alpha)./alpha在零点不优雅改成了sinc(alpha/pi)结果 MATLAB 的sinc定义是sin(pi*x)/(pi*x)参数差一个pi条纹宽度全错排查了半天才定位到这一行。从那以后我每次封装完都强制拿原脚本跑一遍逐点对比确认数值一致再往下做。希望帮到你。本文还有配套的精品资源点击获取