菲涅尔圆孔衍射MATLAB仿真：从公式到视频的完整实现

简介面向本硕博教研学习提供菲涅尔圆孔衍射的Matlab仿真实现与配套操作录像帮助理解标量衍射理论中圆孔衍射的光场分布计算与可视化方法。资源包共2个文件包含一个可直接运行的Matlab主程序脚本和一个avi格式的操作演示视频压缩包仅126KB轻量易携带。运行需使用Matlab2021a或更高版本执行主脚本即可不建议直接运行子函数文件并注意当前文件夹路径需指向工程所在目录配套录像完整演示了操作流程可有效降低上手门槛。已有1819人学习下载适合正在学习物理光学、信息光学或从事衍射仿真的本科生、研究生及科研人员尤其适用于课程实验、毕业设计与算法复现场景。通过仿真代码与录屏结合读者可掌握菲涅尔圆孔衍射数值模拟的编程思路、参数设置与图像输出方式便于迁移到自建项目中。1. 实验报告需要一张图样我却用菲涅尔圆孔衍射matlab仿真把光学课重新学了一遍第一次交光学实验报告时我对着投影仪拍下的圆孔衍射照片犯了难中心亮斑灰蒙蒙旁边的环勉强能看到三四圈想量个半径比都要拿尺子在屏幕上比划。后来换上菲涅尔圆孔衍射matlab仿真这条路几分钟就能从孔径函数算出观察屏的强度分布参数随便改图样随便存连老师追问的「中心为什么不是永远最亮」都能用代码回答。这个方向做的事其实很具体在菲涅尔近似成立的条件下把圆孔的透过率函数与传输核做卷积或直接算衍射积分得到观察面的复振幅再取模方画成强度图。适合正在做光学课程设计、准备近场衍射实验或者刚接触标量衍射理论、想验证公式推导的人。下文从模型离散化讲到录像脚本全部是能直接跑通的东西。2. 把菲涅尔积分变成矩阵运算模型离散化与算法选型菲涅尔圆孔衍射的物理基础是标量衍射理论中的菲涅尔近似。设孔径平面坐标为 (ξ,η)观察屏坐标为 (x,y)单色平面波垂直入射经过距离 z 后观察屏上的复振幅可以写成U(x,y) e^{jkz} / (jλz) · ∬ U₀(ξ,η) · exp{ jk[(x-ξ)² (y-η)²] / (2z) } dξdη其中 U₀(ξ,η) 是孔径平面的透射函数圆孔就是 ξ²η² ≤ a² 内取 1、外面取 0k2π/λ 是波数。这个积分没有解析闭式解圆孔的情况可以化为洛梅尔函数但那套贝塞尔级数用起来还不如直接数值算所以仿真本质上是在做一件事把连续积分离散化成矩阵运算并且离散化方式直接决定结果对不对、跑得快不快。2.1 菲涅尔近似成立的条件先看菲涅尔数再看距离很多人一上来就套公式忘了菲涅尔近似本身有适用范围。菲涅尔近似要求观察屏上某点与孔径上各点的光程差展开式中高于二次的项可以忽略。工程上更实用的判据是菲涅尔数 Nf a²/(λz)其中 a 是圆孔半径。Nf 是个无量纲数物理直觉也很直接它粗略表示观察屏上可见亮环的数量级Nf 越大图样里包含的细节越多越接近「近场」。菲涅尔数 Nf近似层次图样特征远小于 0.1夫琅禾费区艾里斑包络中心亮斑为主次级环很弱0.1 1过渡区中心开始出现亮暗起伏环数增多1 10典型菲涅尔区中心亮暗交替明显环数约等于 Nf 量级大于 10深近场大量细密环边缘衍射效应强对采样率要求高我一般会在代码里先把 Nf 算出来打屏上再决定后续参数。比如波长 632.8nm、孔半径 1mm、距离 0.5mNf≈3.16属于典型菲涅尔区中心轴上强度应该能明显看到亮暗交替的振荡——这个先验判断能帮你快速发现仿真结果是不是出了问题。2.2 两条数值路线逐点数值积分与 FFT 卷积我为什么选后者菲涅尔积分离散化有两条常见路线。第一条是直接数值积分在孔径平面上划分网格对每个观察点 (x,y) 做二维求和。假设观察屏有 N×N 个点、孔径内有 M×M 个采样点复杂度是 O(N²M²)。N 取 1024 时这个循环基本跑不动等一分钟都算不出一帧强度图。第二条是把积分改写成卷积形式。观察积分式括号里的交叉项展开后exp{ jk(x²y²)/(2z) } · exp{ jk(ξ²η²)/(2z) } · exp{ -jk(xξyη)/z }后两项合起来正是傅里叶变换核的结构。更常见的做法是直接定义菲涅尔传输核 h(x,y)e^{jkz}/(jλz)·exp{jk(x²y²)/(2z)}则 U U₀ ⊗ h这个卷积可以借助 FFT 在频域完成先对 U₀ 和 h 各做一次二维 FFT相乘后再 IFFT 回来复杂度降到 O(N²logN)。数值路线计算量精度特点适用场景逐点数值积分O(N²M²)误差可控直接逼近积分定义小网格校验、单点验证FFT 卷积O(N²logN)快但有抽样误差和周期延拓误差大网格、参数扫描、视频录制FFT 卷积的误差来源同样要心里有数一是采样间隔不够导致的频谱混叠二是 FFT 隐含周期延拓带来的边缘卷绕。这两个坑在第四章专门讲。选型结论很明确——只要观察屏点数超过 256×256就用 FFT 卷积逐点积分只用来做小规模交叉验证。3. 手把手实现圆孔掩模、菲涅尔核与强度分布这一章把完整流程拆成四步初始化参数、构造掩模和核、做 FFT 卷积、画图验证。代码我都按可直接运行的方式给出单位统一用米避免出现「把毫米当成米算错三个数量级」这种低级事故。3.1 波长、孔径与观察距离参数初始化这样做第一步先把物理参数和网格参数分开定义。物理参数决定你算的是哪个物理问题网格参数决定数值精度混在一起改起来容易顾此失彼。% 物理参数单位一律用米 lambda 632.8e-9; % 波长 632.8nm氦氖激光常见谱线 a 1e-3; % 圆孔半径 1mm z 0.5; % 孔到观察屏距离 0.5m k 2*pi/lambda; % 波数 Nf a^2/(lambda*z); % 菲涅尔数约 3.16先算出来备用 % 网格参数 N 2048; % 单边采样点数必须是 2 的幂方便 FFT L 0.012; % 观察屏物理宽度 12mm dx L/N; % 采样间隔约 5.86 微米 x (-N/2:N/2-1)*dx; % 观察屏坐标向量以 0 为中心 [xi, eta] meshgrid(x); % 孔径平面坐标矩阵R2016b 后可用单参数写法这里 L 的取值不是随便拍的。孔径半径 1mm衍射图样的主要能量集中在 λz/a ≈ 0.32mm 的尺度内观察屏取 12mm 是为了覆盖中心振荡区加上若干次级环同时给 FFT 的周期延拓留出缓冲。采样间隔 dx 决定了可分辨的最高空间频率2048 个点分 12mm每个环能分到几十个像素画出来边缘光滑。3.2 构造圆孔掩模与菲涅尔核ifftshift 是个关键细节掩模就是孔径平面上的 0/1 矩阵用距离判断当前网格点是否落在圆孔内。菲涅尔核 h 的定义是传输函数注意它和掩模一样都要做一次 ifftshift——这是初学者最容易忽略的一步。MATLAB 的 fft2 约定频谱原点在矩阵的 (1,1) 位置空间域原点在矩阵中心ifftshift 的作用就是把中心原点到左上角否则卷积结果会整体平移半个周期图样中心对着矩阵角上。% 圆孔掩模半径内为 1外为 0 U0 double(xi.^2 eta.^2 a^2); U0 ifftshift(U0); % 原点从中心挪到 (1,1)配合 fft2 % 构造菲涅尔传输核 [X, Y] meshgrid(x); % 观察屏坐标用于核的表达式 h exp(1i*k*z) / (1i*lambda*z) .* exp(1i*k/(2*z)*(X.^2 Y.^2)); h ifftshift(h); % 同样先把原点移到左上角 % 频域卷积菲涅尔积分等价于掩模与核做卷积 U ifft2(fft2(U0) .* fft2(h)) * dx^2; I abs(U).^2; % 强度 复振幅模方代码最后一行乘 dx² 不是可选项。连续积分里的 dξdη 在离散化后是每个网格小方块的面积而 FFT 卷积默认按单位间隔求和不乘回物理面积中心强度数值就会差上百万倍。这一步我见过太多人漏掉然后怀疑公式写错。验证一下掩模是否正确sum(U0(:))*dx² 应该约等于 πa²也就是圆孔面积。这个检查一秒钟就能做能提前排除坐标偏移和半径单位错误。3.3 画强度图和一维截面先判断整体像不像算完复振幅先用强度图看整体分布再在 y0 的中线上画一维截面方便对照理论上的轴上振荡。% 二维强度图 figure(Color, w); imagesc(x*1e3, x*1e3, I); % 坐标转成毫米显示 axis image; colormap(gray); % 光学衍射图样用灰度显示更直观 xlabel(x / mm); ylabel(y / mm); title(sprintf(z%.3f m, N_f%.2f, z, Nf)); % 一维截面取 y0 那一行 figure(Color, w); plot(x*1e3, I(ceil(N/2), :), LineWidth, 1.2); xlabel(x / mm); ylabel(Intensity); title(y0 截面强度分布);按 3.1 节的参数跑出来后你应该看到中心有一个亮斑周围有几圈暗环亮环交替且中心的强度不是单调衰减——这是菲涅尔区区别于夫琅禾费区的重要特征。如果画出来一幅平滑的高斯亮斑八成是 z 取得太大钻进夫琅禾费区了或者菲涅尔核构造时少了交叉项。3.4 把参数扫描录成操作视频从静态图到动态演示标题既然带了「操作视频」仿真最后一步自然是把参数变化过程录下来。常见做法是用 VideoWriter 把每一帧的 figure 抓成视频帧配合 title 动态显示当前参数。我一般会固定颜色映射范围否则每帧会自动缩放到自己的最大最小值视频里亮度会一闪一闪。% 录制距离 z 从 5cm 扫到 1.5m 的衍射演化视频 v VideoWriter(fresnel_z_scan.mp4, MPEG-4); v.FrameRate 10; % 每秒 10 帧 open(v); zList linspace(0.05, 1.5, 60); % 60 帧覆盖整个距离范围 for idx 1:numel(zList) zk zList(idx); % 每帧重新构造核并做卷积 hk ifftshift(exp(1i*k*zk) / (1i*lambda*zk) .* exp(1i*k/(2*zk)*(X.^2Y.^2))); Uk ifft2(fft2(U0) .* fft2(hk)) * dx^2; Ik abs(Uk).^2; imagesc(x*1e3, x*1e3, Ik); axis image; colormap(gray); clim([0, 0.9*max(Ik(:))]); % 固定颜色映射范围避免逐帧闪烁 title(sprintf(z%.3f m, N_f%.1f, zk, a^2/(lambda*zk))); drawnow; writeVideo(v, getframe(gcf)); % 抓取当前图形窗口 end close(v);注意 clim 是较新版本 MATLAB 的写法如果你的版本较旧会提示未定义改成 caxis 即可。VideoWriter 在 MATLAB 里默认支持 MPEG-4 编码直接生成 mp4绝大多数播放器和报告里嵌入都能用。录制前建议先把 figure 窗口调成固定尺寸避免视频尺寸每帧变化。4. 参数与边界的五个坑菲涅尔圆孔仿真最容易翻车的地方仿真是个「算出来容易算对难」的活。这一章把最常见的五个翻车现场整理出来每条都按「现象 → 原因 → 解决」说清楚每一件都是我踩过或用同样问题帮别人排查过的。4.1 图样出现网格状纹理抽样点数不足导致频谱混叠现象强度图上明暗环之间的过渡带有规律的网格状斜纹像打印机网点而且网格密度随 N 减小而变粗。原因观察屏采样间隔 dx 太大菲涅尔核在频域的高频分量被周期延拓折叠回低频区产生混叠。圆孔边缘的突变和高频衍射分量要求采样间隔满足奈奎斯特条件经验上每个最密条纹至少要有 23 个采样点。解决最简单的是把 N 从 1024 提高到 2048 或 4096观察屏宽度 L 不动。如果内存吃紧可以把 L 缩小到只覆盖主要振荡区比如 4mm 到 6mm但边缘卷绕风险会上升。我一般先跑 N512 快速看趋势确认没问题再开 2048 出正式图。4.2 中心强度数值对不上解析公式忘了乘采样间隔平方现象强度分布形状完全正确但中心极大值比用菲涅尔积分直接算出的理论值大了几十万倍或整体亮度异常刺眼。原因离散卷积的求和少了 dξdη 这一项。FFT 卷积默认网格间隔是 1而连续积分的积分元是面积元物理网格面积是 dx²。不乘回去结果只是形状对、量级全错。解决在 ifft2 之后乘 dx²。另外用 sum(U0(:))*dx² 验证掩模面积如果这个值和 πa² 不一致优先怀疑坐标向量长度或半径单位。4.3 图样边缘出现半圈假条纹FFT 周期延拓引起的卷绕现象观察屏左右边缘各出现一段对称的圆弧状亮带像是把中心图样平移到了角落尤其当 L 取得很小比如只比孔径大 2 倍时更明显。原因FFT 卷积把有限尺寸的观察域当成周期信号处理左侧被截断的信号会从右侧卷回来。本质上就是频域相乘等价于循环卷积而不是线性卷积。解决最有效的办法是扩大观察屏宽度 L让边缘处的强度自然衰减到接近 0卷绕回来的能量就小到看不见。更严格的做法是采用「边缘补零」把 U0 和 h 都放到 2 倍大小的零填充矩阵中再做 FFT 卷积最后裁剪中间区域。代价是 N 翻倍后内存占用变成 4 倍2048 点下已经较大所以实际项目中我优先选扩大 L。4.4 z 很小时输出乱码相位因子振荡速度超出采样率现象距离 z 从 0.5m 改到 0.02m 后图样变成一片随机噪点没有任何环形结构或者环很密但位置忽左忽右。原因菲涅尔核 exp{jk(x²y²)/(2z)} 的空间振荡频率与 1/z 成正比。z 变小时相位随坐标变化极快相邻两个采样点之间的相位差超过 π采样定律被破坏算出来的相位就是混叠后的随机值。解决两招——缩小观察屏宽度 L相位二次项的最大值也跟着变小或者增大 N。一般经验是保证 exp{jkx²/(2z)} 在观察屏边缘处的相位变化不超过 π 的几分之一。更稳妥的方案是改用角谱法数值传播它对大角度分量处理更好但代码复杂度高一些。课程设计场景里我建议直接避开极小 z把最小距离设在 Nf≈20 以内的地方。4.5 录制视频时画面忽明忽暗颜色映射范围被逐帧重置现象录制的 z 扫描视频里前面几帧亮得刺眼中间突然整体变黑亮度像在呼吸一样跳动。原因imagesc 默认按当前帧数据的最小/最大值自动设定颜色映射范围。衍射图样中心亮斑的强度随 z 变化剧烈自动范围每帧都在变像素值相同的物理强度画出来亮度完全不同。解决像 3.4 节代码里那样用 clim([0, c]) 手动固定上限c 可以取整段扫描中最大强度的 90% 或一个固定物理值。这样整段视频亮度标准统一才能直观对比不同距离下图样的相对强弱变化。5. 把菲涅尔数当标尺做参数扫描一个能直接写进报告的分析习惯仿真跑通之后真正拉高报告档次的是「用一条物理主线把参数串起来」。我习惯用菲涅尔数 Nf 而不是原始参数 z 来做扫描轴因为它把波长、孔径、距离三个量压缩成一个无量纲数不同参数组合只要 Nf 相同图样就相似——这就是菲涅尔衍射的相似性。具体做法是固定一个目标菲涅尔数序列比如 Nf 8、4、2、1、0.5对每个 Nf 反推出一组 (z, a)分别仿真后并排画出来。观众一眼就能看出「远离孔时环数变少、中心趋于平滑」的物理过程而不是看到一堆乱七八糟的距离数值。这个思路同样适用于录制视频把 3.4 节的扫描变量从 z 换成 Nftitle 里同时显示 Nf 和实际距离视频的叙事逻辑会清楚很多。配套一个常用的数值自检菲涅尔衍射满足能量守恒不考虑反射损耗时观察屏上所有像素强度之和乘以像素面积应约等于入射到孔径上的总能量 πa²。我通常在算完强度后补一行totalPower sum(I(:)) * dx^2; % 观察屏总功率应与 pi*a^2 同量级如果这个值和 πa² 差超过 5%回去查上面四个坑如果差在 1% 以内这组参数基本可信。我的习惯是先跑 N512 的低精度版本把流程走通再上 2048 或 4096 出正式图最后统一用 Nf 扫描录视频——低配调试永远比高配调试省时间也能让你更早暴露参数设错的问题。这套流程帮我少熬夜改图很多次希望帮到你。本文还有配套的精品资源点击获取