简介这份资料围绕菲涅尔波带片的数值模拟展开面向正在学习光学衍射、信息光学课程或需要完成相关实验作业的高校学生与Matlab初学者帮助解决波带片条纹绘制与半波带数判断这类编程实现问题。包内仅含1个doc文档体积约163KB以图文与代码片段穿插的形式讲解波长、半径、焦距等参数设置并给出奇数波带片与偶数波带片两种绘制的完整思路包括将屏幕划分为1001×1001个点、逐点求半波带数、按奇偶性决定涂黑或透光再借助灰度映射与image函数输出黑白相间的波带片图像。内容预览可见clear、linspace、fix、mod等关键语句的用法说明读者可据此复现模拟结果理解菲涅尔波带片对光波干涉与衍射的作用机制并迁移到更一般的衍射光学仿真中。目前已有859人学习下载适合作为课程实验与自学参考。1. 先把波带片的物理尺度算对再谈 MATLAB 画图一块刻着同心圆环的玻璃片不用任何透镜就能把平行光聚成一个亮点甚至沿着光轴冒出好几个次焦点——这是菲涅尔波带片最反直觉的地方。很多人第一次用 matlab 做这类模拟直接imagesc画一张圈圈图就收工结果传播出来的场根本不聚焦。问题几乎都出在环带半径那一行公式上它同时绑定了波长、焦距和环带序号任何一个量代错单位后面所有 matlab 画图出来的结果都是错的。本文瞄准的是这样一类需求你手里有一个设计波长和焦斑要求需要算环带半径、生成透过率分布再用衍射传播把它送到焦平面看场分布顺便验证效率对不对。适合做光学、太赫兹、X 射线成像的工程师也适合想拿这个例子练手 matlab 图像处理与傅里叶光学的读者。物理先站住代码才不会白跑。2. 菲涅尔波带片的环带半径公式与透过率矩阵生成2.1 从光程差反推每一圈的半径波带片的原理是把波前切成若干半波带。轴上某点到波带片第 n 圈边缘与到中心的距离差为 nλ/2 时相邻波带对该点的贡献正好反相。把几何关系展开成 r² 的二次式就得到这圈的外半径r_n² nλf (nλ/2)²对可见光和毫米波第二项相对第一项常常只有千分之几但 X 射线或长焦场景下不能省。一旦省略最外几圈半径会偏小导致边缘环带宽度算错模拟出的焦斑会明显变宽。所以第一步就是把完整公式写进代码而不是用近似版。顺序上n 从 1 开始计数第 1 圈的“外半径”用的是公式的 n1它的“内半径”为 0。相邻两圈之间就是一条环带奇数环透光、偶数环挡光或者反过来相位恰好差 π效果一样。2.2 用 zone_idx 一次性生成振幅型透过率手工去拼每一圈(Rr_n(k)) (Rr_n(k-1))会写得又长又容易错。更稳的做法是先给每个像素标出它属于第几圈再统一按奇偶赋透过率% ---------- 参数全部 SI 单位---------- lambda 532e-9; % 波长 532 nm f 0.1; % 设计焦距 100 mm n_zones 20; % 环带总数 N 2048; % 采样点数 N x N % ---------- 环带外半径 ---------- n (1:n_zones); r_n sqrt(n*lambda*f (n*lambda/2).^2); % 保留完整公式的第二项 % ---------- 采样网格 ---------- L 2.4 * r_n(end); % 窗口略大于最外环直径 x (-N/2:N/2-1) * (L/N); % 以中心为原点 [X, Y] meshgrid(x, x); R hypot(X, Y); % ---------- 逐环带打标签 ---------- zone_idx zeros(N); r_in [0; r_n(1:end-1)]; % 每圈的内半径 for k 1:n_zones zone_idx(R r_in(k) R r_n(k)) k; end % ---------- 振幅型奇数环透明 ---------- T_amp double(mod(zone_idx, 2) 1); % ---------- 看一眼结构 ---------- figure; imagesc(x*1e3, x*1e3, T_amp); axis image; colormap gray; xlabel(x / mm); ylabel(y / mm); title(Fresnel zone plate, amplitude type);逻辑说明zone_idx把空间结构一次性编码成整数标签之后无论你要做振幅型、相位型还是多焦点的变体只需要换一行对标签取模的表达式即可。参数上L取到 2.4 倍最外环半径是为了留点余量、避免光场在窗口边缘被截断产生额外的衍射条纹如果只看透过率图这个余量可以更小但后面做传播时必须留够。注意所有长度必须统一成米。用毫米代进r_n会让结果差 10³ 量级而图上圈数看起来“还挺像”这才是最坑的地方。2.3 采样率和最外环宽度决定了模拟能不能信真正决定成败的不是 N 取多大而是最外圈宽度与像素尺寸之比。最外一圈的宽度约为 Δr λf / (2 r_n(end))。用上面的参数算一下r_20 ≈ 1.03 mmΔr ≈ 25.8 μm。取 L 2.47 mm、N 2048像素约 1.2 μm一个最外环约 21 个像素——这个采样密度做角谱传播是不会混叠的。如果实际工作里波长或焦距让你算出最外环只有几个像素宽就有两个方向增大 N或者利用波带片自带f/3、f/5次焦点的性质用更少环带、更大最外环去等效模拟。参数取值直观影响波长 λ532 nm决定所有环带绝对半径焦距 f100 mm与 λ 相乘决定半径平方的斜率环带数20越大 NA 越大、焦斑越小、环带越细采样 N2048决定每个最外环能被分成几个像素窗口 L2.4·r_N太小则截断太大则像素变粗3. 角谱法传播把波带片的场送到焦平面3.1 为什么优先选角谱法而不是菲涅尔近似菲涅尔近似的条件是把传播距离 z 与横向尺度的平方作比较满足 z³ ≫ (π/4λ)·[(x−x)²(y−y)²]² 才能用。对一块尺寸只有几个毫米、焦距 100 mm 的波带片这个不等式常常刚好卡在边缘更麻烦的是波带片边缘的环带宽度远小于中心环带宽度对空间高频分量要求很高菲涅尔近似对高频的相位误差会直接体现在焦斑旁瓣上。角谱法把传播拆成两步先用 FFT 把场分解成不同方向的空间频率对每个分量乘上它自己的传播相位再逆变换回来。它不需要任何横向尺度近似唯一前提是采样要能覆盖最高空间频率这正好和上一节的采样条件对上了。3.2 传递函数 H(fx,fy) 在 MATLAB 中的实现传播距离为 z 时频率为 (fx, fy) 的平面波分量对应的传播相位是 exp(i·2πz·√(1/λ² − fx² − fy²))。当 fx² fy² 1/λ² 时是倏逝波指数衰减直接截零即可。% ---------- 频率网格 ---------- fx (-N/2:N/2-1) / L; % 与空间网格对偶 [FX, FY] meshgrid(fx, fx); arg 1/lambda^2 - FX.^2 - FY.^2; arg(arg 0) 0; % 倏逝波截断防止开根号出复数 % ---------- 传播到设计焦平面 z f ---------- z f; H exp(1i*2*pi*z*sqrt(arg)); % 传递函数fftshift 格式 H ifftshift(H); % 转成 fft2 需要的零频在角上的格式 U0 ifftshift(T_amp); % 输入场同样转到零频在角上 Uz fftshift(ifft2(fft2(U0) .* H)); I abs(Uz).^2; % 焦平面强度 figure; imagesc(x*1e3, x*1e3, I); axis image; colormap hot; xlabel(x / mm); ylabel(y / mm); title(焦平面强度分布);逻辑与参数说明ifftshift与fftshift的位置必须对称否则相当于把传递函数整体平移了半个窗口结果会出现奇怪的对称破坏——这是角谱法里最常见的错误。arg 0的分量截零对应忽略倏逝波对波长和特征尺寸相差三个量级以上的波带片完全够用。若传播距离很大还需要检查z·λ·fx_max是否超过 1否则混叠会出现环状伪影直观表现是焦斑外围多出很多不该有的同心条纹。3.3 焦平面到中心点的采样够不够用 FFT 做角谱传播时输出平面的网格和输入完全一致都是 L/N。这意味着焦斑分辨率受限于像素尺寸。理论焦斑半宽大约 0.61·λ/NA用上面参数可估到 NA ≈ r_20/f ≈ 0.0103半宽约 31 μm而像素 1.2 μm焦斑上能摊到二十几个点画图和进一步做定量分析都够用。提示想看轴上强度必须取到 N/21 那个索引MATLAB 从 1 开始计数。其余位置都不是严格中心会让焦点位置测量出现偏移。4. 沿 z 轴扫掠验证焦点位置与效率4.1 用 z 扫掠找真实焦斑峰把 H 里的 z 从 0.6f 一路扫到 1.4f每次取轴上强度就能画出轴上光强随距离的曲线。峰值位置应当落在设计焦距附近偏差一般不超过几个百分点如果偏差很大先怀疑是 r_n 公式里第二项被漏掉或者单位没统一。z_list linspace(0.6*f, 1.4*f, 81); I_axis zeros(size(z_list)); for m 1:numel(z_list) zz z_list(m); H exp(1i*2*pi*zz*sqrt(arg)); Uz fftshift(ifft2(fft2(U0) .* ifftshift(H))); I_axis(m) abs(Uz(N/21, N/21))^2; % 轴上中心点 end [~, idx] max(I_axis); fprintf(设计焦距 %.4f m, 实测峰值 %.4f m\n, f, z_list(idx));逻辑说明sqrt(arg)在循环外算好一次能省不少时间这里为了可读性写在表达式里也无妨若 z 列表更长建议提前把sqrt(arg)存成变量。参数 81 采样点足够看到主峰位置想看次峰可以再加大范围并把区间改到 0.1f~1.5f。4.2 次焦点和振幅型的效率天花板沿轴扫掠还会看到 f/3、f/5 附近出现更高阶的焦斑它们分别对应相位差为 3π、5π 的波带组合。这正是波带片被用于多焦点成像和 X 射线相衬的物理基础也是很多 matlab 图像处理论文里把“多平面重建”当成重点的原因——同一块器件在不同 z 处能聚焦不同相位信息。振幅型波带片的理论聚焦效率只有 1/π² ≈ 10.1%剩下的光要么被挡掉要么被送到次焦点。想验证这个数值把焦斑主瓣径向积分一遍与入射到波带片孔径内的总能量相除就能得到实际效率。如果测出来只有几个百分点通常是被挡光的部分没从总入射能量里扣掉或者焦斑积分半径取得太小主瓣以外的能量漏算了。量理论值模拟关注点主焦斑位置f轴扫掠峰值索引主焦斑效率约 10.1%焦斑能量 / 孔径内总入射能量一级次焦点位置f/3轴扫掠次峰一级次焦点效率约 1/π²·(1/9)与主峰比例是否约 1/9注意把总入射能量取成整个 N×N 矩阵的和是不对的应当只统计落到波带片孔径R r_n(end)内的部分。窗口外的空白是空气不是入射光。4.3 对比菲涅尔近似的差异判断什么时候必须用角谱法做教学演示时可以顺手用菲涅尔近似写一条对照曲线对输入场 fft2 之后乘一个二次相位 exp(iπz(fx²fy²)/λ)再逆变换。两条曲线在中心重合得很好但一旦焦距缩短、环带数变多菲涅尔近似会让第一个旁瓣被压低甚至主峰位置轻微偏移。判定阈值可以简单记作当最外环宽度 Δr 小于 5λ 时菲涅尔近似的相位误差就开始肉眼可见这时候必须回到角谱法。5. 相位型波带片、多焦点结构与参数扫掠技巧把透过率从“0/1”改成“0/π”是波带片模拟里性价比最高的一步。只改一行T_phase exp(1i*pi*double(mod(zone_idx,2) 1)); % 奇数环 π 相移之后所有传播代码保持不变把T_amp换成T_phase即可。理论效率会从约 10.1% 提升到约 40.5%焦斑更亮次焦点被明显压制。做 matlab 图像处理或者成像仿真时相位型的点扩散函数旁瓣更低、对比度更好是常用配置。接下来是参数扫掠。十有八九你要看的不是一组参数而是“不同环带数 / 不同焦距下焦斑半宽怎么变”。最省事的写法是向量化把最内层空间坐标不动只把r_n和T生成放进循环传播还是用同一份频率网格arg和FX、FY省下大量重复计算n_list [10 20 30 40]; for m 1:numel(n_list) rz sqrt((1:n_list(m))*lambda*f ((1:n_list(m))*lambda/2).^2); zone_idx zeros(N); r_in [0; rz(1:end-1)]; for k 1:n_list(m) zone_idx(R r_in(k) R rz(k)) k; end T exp(1i*pi*double(mod(zone_idx,2)1)); Uz fftshift(ifft2(fft2(ifftshift(T)) .* ifftshift(H))); I abs(Uz).^2; % 焦斑半高宽对过中心的一行做归一化后取半高 line I(N/21, :); [~, p] max(line); fwhm sum(line max(line)/2) * (L/N); fprintf(环带数 %2d, 焦斑半高宽 %.2f μm\n, n_list(m), fwhm*1e6); end逻辑说明H和arg在循环外算好一次是扫掠提速的关键fwhm用“超过半高点的像素个数 × 像素尺寸”近似比插值法糙一些但对比较趋势足够了。想更精确可以在半高附近做线性插值代码略长但结果顺滑。最后一条实操经验环带数一多最外圈变细很容易踩到采样不足此时优先增大 N而不是缩小 L。把 L 缩到只剩最外环半径的 1.1 倍看起来像素变细了但窗口截断会把轴上强度压低扫掠出的效率会莫名其妙地低于理论值很容易被误判成“相位型效率不对”。本文还有配套的精品资源点击获取