简介这份资源是面向光学、物理及Matlab入门学习者的菲涅尔波带片仿真资料核心解决如何用编程方式复现波带片干涉衍射图像的问题。包内共1个文件为163KB的doc文档正文以文字讲解配合Matlab源码片段展开便于边读边在软件中运行验证。内容从波带片将波面划分为多个半波带、按半径平方与波长焦距关系确定各点半波带数的原理讲起给出完整的参数设置思路波长600nm、半径3mm、焦距1m把屏幕分为1001×1001个点用双重循环逐点计算所在圆半径与半波带数再依据奇偶性判断涂黑或透光并通过灰度映射与image函数绘制结果分别得到黑白相间的偶数波带片与灰底相间的奇数波带片两类图像。文档配有奇偶两种情形的完整代码与运行效果对照可作为光学实验的辅助参考、课程设计或自学练习素材帮助读者理解矩阵运算与数字信号处理在光波仿真中的具体用法。目前已有859人学习。1. 从一张黑白环图说起Matlab 菲涅尔波带片是怎么拼出来的很多人第一次跑菲涅尔波带片仿真都会得到一张看起来“差不多”的同心圆图中心一个黑点往外黑白交替半径越大环越密。但把图发给做光学实验的人对方第一句话往往是外环半径不对焦距也对不上。问题通常不在公式而在坐标单位、采样点数和 k 的取整方式。这份 Matlab 学习资料里的 doc 只给了几十行代码却把三个关键动作都串起来了用半波带数 k 判断奇偶、用 linspace 建立 1001×1001 的像素网格、用灰度映射输出遮光与透光区域。它适合想入门光学仿真的人也适合已经会用 Matlab 画图、但没认真处理过物理量单位的开发者。2. 半波带数 k 与菲涅尔半径公式坐标、单位与奇偶判定2.1 半波带半径 r_k sqrt(k λ f) 从光程差来菲涅尔波带片不是普通同心圆光栅。它把波前分成若干半波带相邻半波带到焦点的光程差为 λ/2奇偶波带在焦点处的相位相反。如果挡住奇数波带、只让偶数波带透光剩余波带在焦点附近同相叠加轴上就会出现亮点。第 k 个半波带外边界半径满足近轴近似% 半波带半径快速估算 lam 600e-6; % 波长单位 mm即 600 nm f 1000; % 焦距单位 mm即 1 m k 1:10; r_k sqrt(k * lam * f); disp(table(k, r_k, VariableNames, {k, r_mm}));这段代码用 lam 和 f 直接算前 10 个半波带的理论半径。lam 和 f 必须统一到 mm算出的 r_k 才是 mm。若把 lam 写成 600e-9、f 写成 1单位变成 m半径也会变成 m数值上同样成立但不能和 R3 mm 混用。原 doc 中psqrt(x(m).^2y(n).^2)得到像素点到中心的距离kfix(p^2./(lam.*f))得到该点所在的半波带编号。环半径按 sqrt(k) 增长所以外圈越来越密这是波带片和普通同心圆的根本区别。2.2 参数单位对齐lam600e-6、R3、f1000 到底意味着什么原代码里最容易埋雷的是单位。lam600e-6 看着像“0.0006 米”但结合 R3 和 f1000它实际是 600e-6 mm也就是 600 nm。f1000 是 1000 mm即 1 m。R3 是波带片半径 3 mm。三个量都用 mmp、k、R 的比较才成立。常见错误是把 lam 当米、把 R 当毫米结果 k 值巨大整个波带片被压成中心一个小黑点或者图里只有一圈灰。变量物理量原代码取值换算到 mm注意点lam波长600e-6600e-6 mm即 600 nm不要当成米f焦距10001000 mm即 1 mR波带片半径33 mm外边界半径N每边采样点数1001无量纲奇数保证中心落在网格点x, y坐标向量linspace(-R,R,N)mm区间对称步长 0.006 mm用这组参数算一下外边界k_max R^2/(lam*f) 9/0.6 15。也就是说半径 3 mm 内包含 15 个半波带。如果程序画出来的黑白环数明显少于 15或者中心区域尺寸不对优先回头查单位。2.3 奇偶波带的透光/遮光策略菲涅尔波带片分为偶数波带片和奇数波带片。偶数波带片让偶数 k 透光、奇数 k 遮光奇数波带片反过来。原代码用mod(k,2)1判断奇数奇数赋 0 表示涂黑偶数赋 1 表示透光。背景区域pR不能和遮光黑混在一起否则整张图看起来像一个大黑圆无法判断波带是否算对。常见做法是把背景赋 0.5用灰色和黑白两档区分开。波带片类型k 为奇数k 为偶数p R偶数波带片0黑遮光1白透光0.5灰奇数波带片1白透光0黑遮光0.5灰k0 的中心区域按偶数处理所以偶数波带片中心透光奇数波带片中心遮光。实际加工时中心点尺寸很小视觉上不一定明显但判断逻辑上要一致。3. 用 1001×1001 网格生成波带片 mask双重循环到向量化3.1 网格构造linspace 与 meshgrid 的差别原代码用二重循环逐点计算这种写法直观也最容易看出每个像素对应的物理坐标。linspace(-R,R,1001)生成 1001 个等间隔点步长是 6/1000 0.006 mm也就是 6 μm。像素索引 m、n 不是物理坐标不能直接拿去算半径。常见做法是先用 linspace 生成 x、y 向量再在循环里用x(m)、y(n)取实际坐标。如果改用meshgrid[X,Y] meshgrid(x,y)会生成两个 1001×1001 的坐标矩阵。X 的每一行相同Y 的每一列相同。矩阵索引和 imagesc 的显示方向需要对齐I(m,n) 的第一个索引是行对应 y第二个索引是列对应 x。原代码中I(m,n)和x(m)、y(n)的搭配是正确的但一旦换成向量化写法就要保证meshgrid(x,y)的顺序不颠倒。3.2 双重循环版脚本与逐行解释先给一个可复现的完整版本修掉了原 doc 里明显的 OCR 错误x(m).A2实际是x(m).^2y(n)42实际是y(n).^2。clear; clc; lam 600e-6; % 波长单位 mm600 nm f 1000; % 焦距单位 mm R 3; % 波带片半径单位 mm N 1001; % 每边采样点数奇数保证有中心点 x linspace(-R, R, N); y linspace(-R, R, N); I zeros(N, N); % 初始化 mask0 黑1 透光0.5 背景灰 for m 1:N for n 1:N p sqrt(x(m)^2 y(n)^2); % 到中心的径向距离 k fix(p^2 / (lam * f)); % 半波带数向下取整 if p R I(m, n) 0.5; % 波带片外的背景 elseif mod(k, 2) 1 I(m, n) 0; % 奇数波带遮光 else I(m, n) 1; % 偶数波带透光 end end end imagesc(x, y, I); colormap(gray(256)); axis image tight; colorbar; title(偶数波带片 mask);逻辑说明外层循环 m 遍历 x 方向内层循环 n 遍历 y 方向每个像素独立算 p 和 k。fix对正数向下取整等价于 floor。p R放在奇偶判断之前避免把波带片外部误判成某个波带。mod(k,2)1对非负整数 k 来说就是奇数判断。参数说明lam、f、R 都用 mmp 也是 mmN1001 时 I 是 1001×1001 的 double 矩阵colormap(gray(256)) 能显示 0、0.5、1 三档灰度若只要二值图可以去掉背景灰。3.3 向量化版一行计算 k两行生成 mask双重循环在 N1001 时还能接受但做参数扫描或提高分辨率时向量化写法会省很多时间。clear; clc; lam 600e-6; f 1000; R 3; N 1001; x linspace(-R, R, N); y linspace(-R, R, N); [X, Y] meshgrid(x, y); P sqrt(X.^2 Y.^2); % 每个像素到中心的距离 K fix(P.^2 / (lam * f)); % 半波带数矩阵 I double(mod(K, 2) 0); % 偶数波带透光奇数遮光 I(P R) 0.5; % 背景灰 imagesc(x, y, I); colormap(gray(256)); axis image tight;参数说明X.^2和Y.^2是逐元素平方不能写成X^2后者是矩阵乘法。double(mod(K,2)0)把逻辑矩阵转成 0/1 数值矩阵方便后续做图像显示或衍射计算。I(PR)0.5用逻辑索引批量赋值比循环里逐个判断更简洁。向量化版本没有显式循环但内存占用更大N4001 时单个矩阵约 128 MB如果机器内存吃紧可以分块处理。写法适用场景优点注意点双重循环教学、调试、小尺寸容易看懂每个像素N 大时速度慢meshgrid 向量化参数扫描、大尺寸代码短、速度快内存占用高分块向量化超大网格平衡速度和内存边界要重叠处理3.4 显示与导出colormap(gray(2))、imagesc 与 imwrite显示时最容易被坑的是灰度映射。二值 mask 用colormap(gray(2))0 映射黑1 映射白。如果矩阵里有 0.5 的背景灰再用 gray(2) 会把 0.5 映射成非预期的中间值所以含背景的图用 gray(256)。导出 PNG 时imwrite 要求整数类型直接写 double 矩阵会出问题。% 二值版只有黑白没有背景灰 I_bw double(mod(K, 2) 0); I_bw(P R) 0; figure; imagesc(x, y, I_bw); colormap(gray(2)); % 两档灰度0 黑1 白 axis image tight; % 导出 PNG先归一化到 0-255 I_export uint8(255 * mat2gray(I_bw)); imwrite(I_export, fresnel_zone_plate.png);说明mat2gray把最小值映射到 0、最大值映射到 1避免 double 矩阵直接写图导致全黑或全白。imagesc(x,y,I_bw)带上 x、y 向量后坐标轴单位就是 mm导出的比例尺才有物理意义。原 doc 中IrI*N; image(lr);大概率是 OCR 误读按上下文看实际就是对 mask 做显示不要把I*N当成物理公式。4. 复现原代码时最常翻车的点索引、尺寸不匹配与波带数边界4.1 x 长度 1010 与 y 长度 1001 的尺寸不匹配原 doc 中有一行xlinspace(-xm,xm,1010); ylinspace(-ym,ym,1001);但循环用的是for m1:1001和for n1:1001。这意味着 x 只取了前 1001 个点实际 x 范围不再是 -3 到 3而是 -3 到约 2.988。波带片在水平方向被轻微压缩中心点也偏了。这种错误在图上不一定一眼看出但外环半径会整体偏小。N 1001; x linspace(-R, R, N); y linspace(-R, R, N);参数说明x、y 必须同长度且 N 取奇数中心点才正好落在 x0、y0。N 取偶数时中心落在两个像素之间中心波带会出现不对称。若必须用偶数 N可以手动把中心点索引附近的像素做对称修正但更省事的做法是直接用奇数。4.2 pR 边界与 k 取整fix、floor、round 的差别fix对正数截断小数等价于 floor对负数向零取整。p 是非负的所以fix(p^2/(lam*f))和floor(p^2/(lam*f))结果一致。边界处浮点误差会让 k 差 1比如商算出来是 3.0000000001 和 2.9999999999前者 k3后者 k2奇偶相反mask 在极细的边界上会翻转。常见做法是加一个极小容差让边界更稳定。K_fix fix(P.^2 / (lam * f)); K_floor floor(P.^2 / (lam * f) 1e-9); diff_idx find(K_fix ~ K_floor); fprintf(fix 与加容差 floor 不同的像素数%d\n, numel(diff_idx));说明容差 1e-9 是经验值针对双精度计算中接近整数的商。不同的像素数通常很少集中在环边界。round会把边界四舍五入到最近整数导致半波带宽度偏移半个环不建议拿来算波带数。物理上波带边界本来就是半波长过渡边界像素对应的遮光或透光对焦点贡献很小所以只要环半径正确少量边界翻转可以接受。4.3 用径向剖面和理论半径快速验证生成 mask 后不要只看图。沿中心行取剖面找灰度跳变位置再和理论半径对比能快速判断参数是否对。row I((N1)/2, :); % 中心行 edge_idx find(diff(row) ~ 0); % 灰度变化位置 r_pixel abs(x(edge_idx)); % 对应半径单位 mm k_theory 1:8; r_theory sqrt(k_theory * lam * f); % 理论半径单位 mm fprintf(理论半径(mm): ); disp(r_theory); fprintf(实测跳变半径(mm): ); disp(r_pixel(1:min(8,numel(r_pixel))));说明中心行取的是(N1)/2行对应 y0。diff(row)~0找灰度变化的位置abs(x(edge_idx))得到跳变半径。理论半径由 sqrt(klamf) 算出前几个波带误差应在半个像素以内。若误差明显偏大优先检查 x、y 长度是否一致以及 lam、f、R 的单位是否统一。现象常见原因修复中心全黑或全白k 奇偶判断写反交换 mod(k,2)1 的赋值外环半径比理论小x 长度 1010只取了前 1001 个点x、y 统一为 1001边界锯齿明显采样点数不足N 提高到 2001 或 4001背景和遮光混在一起没有区分 pR 与 k 奇数背景单独赋 0.5 灰导出图片全黑double 矩阵没有归一化用 mat2gray 后转 uint85. 从静态 mask 到焦点验证菲涅尔波带片的进阶玩法静态图只能看环真正验证波带片是否做对要看焦点。把 mask 当成复振幅用角谱传播到设计焦距再观察中心强度是一个很直接的检查方法。下面这段代码接前面的I_bw透光处为 1遮光处为 0。% 角谱传播验证设计焦距处的聚焦 lambda lam; % mm z f; % 传播距离取设计焦距mm dx x(2) - x(1); % 采样间隔mm N size(I_bw, 1); U0 double(I_bw); % 波带片出射复振幅 k0 2 * pi / lambda; fx (-N/2 : N/2-1) / (N * dx); [FX, FY] meshgrid(fx, fx); H exp(1i * k0 * z * sqrt(1 - (lambda*FX).^2 - (lambda*FY).^2)); H(abs(lambda*FX) 1 | abs(lambda*FY) 1) 0; % 消逝波置零 Uf fftshift(ifft2(fft2(ifftshift(U0)) .* ifftshift(H))); I_focus abs(Uf).^2; imagesc(abs(I_focus)); colormap(hot); axis image tight; title(设计焦距处的焦点强度);说明fx是频率坐标单位 cycles/mmlambda*FX是无量纲方向余弦超过 1 的部分是消逝波直接置零。fftshift和ifftshift用来把零频移到中心。焦点处中心会出现亮斑但菲涅尔波带片有多个焦点轴上强度会在 f 附近出现峰值。若想确认峰值位置扫描 z 从 0.8f 到 1.2f记录中心像素强度。z_list linspace(0.8*f, 1.2*f, 41); I_axis zeros(size(z_list)); for ii 1:numel(z_list) z z_list(ii); H exp(1i * k0 * z * sqrt(1 - (lambda*FX).^2 - (lambda*FY).^2)); H(abs(lambda*FX) 1 | abs(lambda*FY) 1) 0; Uf fftshift(ifft2(fft2(ifftshift(U0)) .* ifftshift(H))); I_axis(ii) abs(Uf((N1)/2, (N1)/2))^2; end [~, idx_max] max(I_axis); fprintf(设计焦距 %.1f mm实测轴上峰值位置 %.1f mm\n, f, z_list(idx_max));这个扫描能暴露采样是否足够。如果峰值偏离设计焦距超过几个采样间隔优先检查波长、焦距单位是否统一以及 mask 的奇偶类型有没有选对。不同焦距对应的外环半径也可以先算表再决定 N 和 R 是否覆盖完整孔径。设计焦距 f / mmlam*f / mm²R3 mm 内半波带数第 1 环半径 / mm5000.3300.54810000.6150.77515000.9100.94920001.271.095把 z_list 和 I_axis 存成两列 CSV再画归一化曲线比直接看 hot 图更容易发现次焦点。本文还有配套的精品资源点击获取 Matlab菲涅尔波带片仿真：半波带mask生成与焦点验证
简介这份资源是面向光学、物理及Matlab入门学习者的菲涅尔波带片仿真资料核心解决如何用编程方式复现波带片干涉衍射图像的问题。包内共1个文件为163KB的doc文档正文以文字讲解配合Matlab源码片段展开便于边读边在软件中运行验证。内容从波带片将波面划分为多个半波带、按半径平方与波长焦距关系确定各点半波带数的原理讲起给出完整的参数设置思路波长600nm、半径3mm、焦距1m把屏幕分为1001×1001个点用双重循环逐点计算所在圆半径与半波带数再依据奇偶性判断涂黑或透光并通过灰度映射与image函数绘制结果分别得到黑白相间的偶数波带片与灰底相间的奇数波带片两类图像。文档配有奇偶两种情形的完整代码与运行效果对照可作为光学实验的辅助参考、课程设计或自学练习素材帮助读者理解矩阵运算与数字信号处理在光波仿真中的具体用法。目前已有859人学习。1. 从一张黑白环图说起Matlab 菲涅尔波带片是怎么拼出来的很多人第一次跑菲涅尔波带片仿真都会得到一张看起来“差不多”的同心圆图中心一个黑点往外黑白交替半径越大环越密。但把图发给做光学实验的人对方第一句话往往是外环半径不对焦距也对不上。问题通常不在公式而在坐标单位、采样点数和 k 的取整方式。这份 Matlab 学习资料里的 doc 只给了几十行代码却把三个关键动作都串起来了用半波带数 k 判断奇偶、用 linspace 建立 1001×1001 的像素网格、用灰度映射输出遮光与透光区域。它适合想入门光学仿真的人也适合已经会用 Matlab 画图、但没认真处理过物理量单位的开发者。2. 半波带数 k 与菲涅尔半径公式坐标、单位与奇偶判定2.1 半波带半径 r_k sqrt(k λ f) 从光程差来菲涅尔波带片不是普通同心圆光栅。它把波前分成若干半波带相邻半波带到焦点的光程差为 λ/2奇偶波带在焦点处的相位相反。如果挡住奇数波带、只让偶数波带透光剩余波带在焦点附近同相叠加轴上就会出现亮点。第 k 个半波带外边界半径满足近轴近似% 半波带半径快速估算 lam 600e-6; % 波长单位 mm即 600 nm f 1000; % 焦距单位 mm即 1 m k 1:10; r_k sqrt(k * lam * f); disp(table(k, r_k, VariableNames, {k, r_mm}));这段代码用 lam 和 f 直接算前 10 个半波带的理论半径。lam 和 f 必须统一到 mm算出的 r_k 才是 mm。若把 lam 写成 600e-9、f 写成 1单位变成 m半径也会变成 m数值上同样成立但不能和 R3 mm 混用。原 doc 中psqrt(x(m).^2y(n).^2)得到像素点到中心的距离kfix(p^2./(lam.*f))得到该点所在的半波带编号。环半径按 sqrt(k) 增长所以外圈越来越密这是波带片和普通同心圆的根本区别。2.2 参数单位对齐lam600e-6、R3、f1000 到底意味着什么原代码里最容易埋雷的是单位。lam600e-6 看着像“0.0006 米”但结合 R3 和 f1000它实际是 600e-6 mm也就是 600 nm。f1000 是 1000 mm即 1 m。R3 是波带片半径 3 mm。三个量都用 mmp、k、R 的比较才成立。常见错误是把 lam 当米、把 R 当毫米结果 k 值巨大整个波带片被压成中心一个小黑点或者图里只有一圈灰。变量物理量原代码取值换算到 mm注意点lam波长600e-6600e-6 mm即 600 nm不要当成米f焦距10001000 mm即 1 mR波带片半径33 mm外边界半径N每边采样点数1001无量纲奇数保证中心落在网格点x, y坐标向量linspace(-R,R,N)mm区间对称步长 0.006 mm用这组参数算一下外边界k_max R^2/(lam*f) 9/0.6 15。也就是说半径 3 mm 内包含 15 个半波带。如果程序画出来的黑白环数明显少于 15或者中心区域尺寸不对优先回头查单位。2.3 奇偶波带的透光/遮光策略菲涅尔波带片分为偶数波带片和奇数波带片。偶数波带片让偶数 k 透光、奇数 k 遮光奇数波带片反过来。原代码用mod(k,2)1判断奇数奇数赋 0 表示涂黑偶数赋 1 表示透光。背景区域pR不能和遮光黑混在一起否则整张图看起来像一个大黑圆无法判断波带是否算对。常见做法是把背景赋 0.5用灰色和黑白两档区分开。波带片类型k 为奇数k 为偶数p R偶数波带片0黑遮光1白透光0.5灰奇数波带片1白透光0黑遮光0.5灰k0 的中心区域按偶数处理所以偶数波带片中心透光奇数波带片中心遮光。实际加工时中心点尺寸很小视觉上不一定明显但判断逻辑上要一致。3. 用 1001×1001 网格生成波带片 mask双重循环到向量化3.1 网格构造linspace 与 meshgrid 的差别原代码用二重循环逐点计算这种写法直观也最容易看出每个像素对应的物理坐标。linspace(-R,R,1001)生成 1001 个等间隔点步长是 6/1000 0.006 mm也就是 6 μm。像素索引 m、n 不是物理坐标不能直接拿去算半径。常见做法是先用 linspace 生成 x、y 向量再在循环里用x(m)、y(n)取实际坐标。如果改用meshgrid[X,Y] meshgrid(x,y)会生成两个 1001×1001 的坐标矩阵。X 的每一行相同Y 的每一列相同。矩阵索引和 imagesc 的显示方向需要对齐I(m,n) 的第一个索引是行对应 y第二个索引是列对应 x。原代码中I(m,n)和x(m)、y(n)的搭配是正确的但一旦换成向量化写法就要保证meshgrid(x,y)的顺序不颠倒。3.2 双重循环版脚本与逐行解释先给一个可复现的完整版本修掉了原 doc 里明显的 OCR 错误x(m).A2实际是x(m).^2y(n)42实际是y(n).^2。clear; clc; lam 600e-6; % 波长单位 mm600 nm f 1000; % 焦距单位 mm R 3; % 波带片半径单位 mm N 1001; % 每边采样点数奇数保证有中心点 x linspace(-R, R, N); y linspace(-R, R, N); I zeros(N, N); % 初始化 mask0 黑1 透光0.5 背景灰 for m 1:N for n 1:N p sqrt(x(m)^2 y(n)^2); % 到中心的径向距离 k fix(p^2 / (lam * f)); % 半波带数向下取整 if p R I(m, n) 0.5; % 波带片外的背景 elseif mod(k, 2) 1 I(m, n) 0; % 奇数波带遮光 else I(m, n) 1; % 偶数波带透光 end end end imagesc(x, y, I); colormap(gray(256)); axis image tight; colorbar; title(偶数波带片 mask);逻辑说明外层循环 m 遍历 x 方向内层循环 n 遍历 y 方向每个像素独立算 p 和 k。fix对正数向下取整等价于 floor。p R放在奇偶判断之前避免把波带片外部误判成某个波带。mod(k,2)1对非负整数 k 来说就是奇数判断。参数说明lam、f、R 都用 mmp 也是 mmN1001 时 I 是 1001×1001 的 double 矩阵colormap(gray(256)) 能显示 0、0.5、1 三档灰度若只要二值图可以去掉背景灰。3.3 向量化版一行计算 k两行生成 mask双重循环在 N1001 时还能接受但做参数扫描或提高分辨率时向量化写法会省很多时间。clear; clc; lam 600e-6; f 1000; R 3; N 1001; x linspace(-R, R, N); y linspace(-R, R, N); [X, Y] meshgrid(x, y); P sqrt(X.^2 Y.^2); % 每个像素到中心的距离 K fix(P.^2 / (lam * f)); % 半波带数矩阵 I double(mod(K, 2) 0); % 偶数波带透光奇数遮光 I(P R) 0.5; % 背景灰 imagesc(x, y, I); colormap(gray(256)); axis image tight;参数说明X.^2和Y.^2是逐元素平方不能写成X^2后者是矩阵乘法。double(mod(K,2)0)把逻辑矩阵转成 0/1 数值矩阵方便后续做图像显示或衍射计算。I(PR)0.5用逻辑索引批量赋值比循环里逐个判断更简洁。向量化版本没有显式循环但内存占用更大N4001 时单个矩阵约 128 MB如果机器内存吃紧可以分块处理。写法适用场景优点注意点双重循环教学、调试、小尺寸容易看懂每个像素N 大时速度慢meshgrid 向量化参数扫描、大尺寸代码短、速度快内存占用高分块向量化超大网格平衡速度和内存边界要重叠处理3.4 显示与导出colormap(gray(2))、imagesc 与 imwrite显示时最容易被坑的是灰度映射。二值 mask 用colormap(gray(2))0 映射黑1 映射白。如果矩阵里有 0.5 的背景灰再用 gray(2) 会把 0.5 映射成非预期的中间值所以含背景的图用 gray(256)。导出 PNG 时imwrite 要求整数类型直接写 double 矩阵会出问题。% 二值版只有黑白没有背景灰 I_bw double(mod(K, 2) 0); I_bw(P R) 0; figure; imagesc(x, y, I_bw); colormap(gray(2)); % 两档灰度0 黑1 白 axis image tight; % 导出 PNG先归一化到 0-255 I_export uint8(255 * mat2gray(I_bw)); imwrite(I_export, fresnel_zone_plate.png);说明mat2gray把最小值映射到 0、最大值映射到 1避免 double 矩阵直接写图导致全黑或全白。imagesc(x,y,I_bw)带上 x、y 向量后坐标轴单位就是 mm导出的比例尺才有物理意义。原 doc 中IrI*N; image(lr);大概率是 OCR 误读按上下文看实际就是对 mask 做显示不要把I*N当成物理公式。4. 复现原代码时最常翻车的点索引、尺寸不匹配与波带数边界4.1 x 长度 1010 与 y 长度 1001 的尺寸不匹配原 doc 中有一行xlinspace(-xm,xm,1010); ylinspace(-ym,ym,1001);但循环用的是for m1:1001和for n1:1001。这意味着 x 只取了前 1001 个点实际 x 范围不再是 -3 到 3而是 -3 到约 2.988。波带片在水平方向被轻微压缩中心点也偏了。这种错误在图上不一定一眼看出但外环半径会整体偏小。N 1001; x linspace(-R, R, N); y linspace(-R, R, N);参数说明x、y 必须同长度且 N 取奇数中心点才正好落在 x0、y0。N 取偶数时中心落在两个像素之间中心波带会出现不对称。若必须用偶数 N可以手动把中心点索引附近的像素做对称修正但更省事的做法是直接用奇数。4.2 pR 边界与 k 取整fix、floor、round 的差别fix对正数截断小数等价于 floor对负数向零取整。p 是非负的所以fix(p^2/(lam*f))和floor(p^2/(lam*f))结果一致。边界处浮点误差会让 k 差 1比如商算出来是 3.0000000001 和 2.9999999999前者 k3后者 k2奇偶相反mask 在极细的边界上会翻转。常见做法是加一个极小容差让边界更稳定。K_fix fix(P.^2 / (lam * f)); K_floor floor(P.^2 / (lam * f) 1e-9); diff_idx find(K_fix ~ K_floor); fprintf(fix 与加容差 floor 不同的像素数%d\n, numel(diff_idx));说明容差 1e-9 是经验值针对双精度计算中接近整数的商。不同的像素数通常很少集中在环边界。round会把边界四舍五入到最近整数导致半波带宽度偏移半个环不建议拿来算波带数。物理上波带边界本来就是半波长过渡边界像素对应的遮光或透光对焦点贡献很小所以只要环半径正确少量边界翻转可以接受。4.3 用径向剖面和理论半径快速验证生成 mask 后不要只看图。沿中心行取剖面找灰度跳变位置再和理论半径对比能快速判断参数是否对。row I((N1)/2, :); % 中心行 edge_idx find(diff(row) ~ 0); % 灰度变化位置 r_pixel abs(x(edge_idx)); % 对应半径单位 mm k_theory 1:8; r_theory sqrt(k_theory * lam * f); % 理论半径单位 mm fprintf(理论半径(mm): ); disp(r_theory); fprintf(实测跳变半径(mm): ); disp(r_pixel(1:min(8,numel(r_pixel))));说明中心行取的是(N1)/2行对应 y0。diff(row)~0找灰度变化的位置abs(x(edge_idx))得到跳变半径。理论半径由 sqrt(klamf) 算出前几个波带误差应在半个像素以内。若误差明显偏大优先检查 x、y 长度是否一致以及 lam、f、R 的单位是否统一。现象常见原因修复中心全黑或全白k 奇偶判断写反交换 mod(k,2)1 的赋值外环半径比理论小x 长度 1010只取了前 1001 个点x、y 统一为 1001边界锯齿明显采样点数不足N 提高到 2001 或 4001背景和遮光混在一起没有区分 pR 与 k 奇数背景单独赋 0.5 灰导出图片全黑double 矩阵没有归一化用 mat2gray 后转 uint85. 从静态 mask 到焦点验证菲涅尔波带片的进阶玩法静态图只能看环真正验证波带片是否做对要看焦点。把 mask 当成复振幅用角谱传播到设计焦距再观察中心强度是一个很直接的检查方法。下面这段代码接前面的I_bw透光处为 1遮光处为 0。% 角谱传播验证设计焦距处的聚焦 lambda lam; % mm z f; % 传播距离取设计焦距mm dx x(2) - x(1); % 采样间隔mm N size(I_bw, 1); U0 double(I_bw); % 波带片出射复振幅 k0 2 * pi / lambda; fx (-N/2 : N/2-1) / (N * dx); [FX, FY] meshgrid(fx, fx); H exp(1i * k0 * z * sqrt(1 - (lambda*FX).^2 - (lambda*FY).^2)); H(abs(lambda*FX) 1 | abs(lambda*FY) 1) 0; % 消逝波置零 Uf fftshift(ifft2(fft2(ifftshift(U0)) .* ifftshift(H))); I_focus abs(Uf).^2; imagesc(abs(I_focus)); colormap(hot); axis image tight; title(设计焦距处的焦点强度);说明fx是频率坐标单位 cycles/mmlambda*FX是无量纲方向余弦超过 1 的部分是消逝波直接置零。fftshift和ifftshift用来把零频移到中心。焦点处中心会出现亮斑但菲涅尔波带片有多个焦点轴上强度会在 f 附近出现峰值。若想确认峰值位置扫描 z 从 0.8f 到 1.2f记录中心像素强度。z_list linspace(0.8*f, 1.2*f, 41); I_axis zeros(size(z_list)); for ii 1:numel(z_list) z z_list(ii); H exp(1i * k0 * z * sqrt(1 - (lambda*FX).^2 - (lambda*FY).^2)); H(abs(lambda*FX) 1 | abs(lambda*FY) 1) 0; Uf fftshift(ifft2(fft2(ifftshift(U0)) .* ifftshift(H))); I_axis(ii) abs(Uf((N1)/2, (N1)/2))^2; end [~, idx_max] max(I_axis); fprintf(设计焦距 %.1f mm实测轴上峰值位置 %.1f mm\n, f, z_list(idx_max));这个扫描能暴露采样是否足够。如果峰值偏离设计焦距超过几个采样间隔优先检查波长、焦距单位是否统一以及 mask 的奇偶类型有没有选对。不同焦距对应的外环半径也可以先算表再决定 N 和 R 是否覆盖完整孔径。设计焦距 f / mmlam*f / mm²R3 mm 内半波带数第 1 环半径 / mm5000.3300.54810000.6150.77515000.9100.94920001.271.095把 z_list 和 I_axis 存成两列 CSV再画归一化曲线比直接看 hot 图更容易发现次焦点。本文还有配套的精品资源点击获取