企业官网建设流程全解析

简介：这份资源是面向光学、物理及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，数值上同样成立，但不能和 R=3 mm 混用。原 doc 中p=sqrt(x(m).^2+y(n).^2)得到像素点到中心的距离，k=fix(p^2./(lam.*f))得到该点所在的半波带编号。环半径按 sqrt(k) 增长，所以外圈越来越密，这是波带片和普通同心圆的根本区别。

2.2 参数单位对齐：lam=600e-6、R=3、f=1000 到底意味着什么

原代码里最容易埋雷的是单位。lam=600e-6 看着像“0.0006 米”，但结合 R=3 和 f=1000，它实际是 600e-6 mm，也就是 600 nm。f=1000 是 1000 mm，即 1 m。R=3 是波带片半径 3 mm。三个量都用 mm，p、k、R 的比较才成立。常见错误是把 lam 当米、把 R 当毫米，结果 k 值巨大，整个波带片被压成中心一个小黑点，或者图里只有一圈灰。

用这组参数算一下外边界：k_max = R^2/(lam*f) = 9/0.6 = 15。也就是说半径 3 mm 内包含 15 个半波带。如果程序画出来的黑白环数明显少于 15，或者中心区域尺寸不对，优先回头查单位。

2.3 奇偶波带的透光/遮光策略

菲涅尔波带片分为偶数波带片和奇数波带片。偶数波带片让偶数 k 透光、奇数 k 遮光；奇数波带片反过来。原代码用mod(k,2)==1判断奇数，奇数赋 0 表示涂黑，偶数赋 1 表示透光。背景区域p>R不能和遮光黑混在一起，否则整张图看起来像一个大黑圆，无法判断波带是否算对。常见做法是把背景赋 0.5，用灰色和黑白两档区分开。

k=0 的中心区域按偶数处理，所以偶数波带片中心透光，奇数波带片中心遮光。实际加工时中心点尺寸很小，视觉上不一定明显，但判断逻辑上要一致。

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).^2，y(n)42实际是y(n).^2。

clear; clc; lam = 600e-6; % 波长，单位 mm，600 nm f = 1000; % 焦距，单位 mm R = 3; % 波带片半径，单位 mm N = 1001; % 每边采样点数，奇数保证有中心点 x = linspace(-R, R, N); y = linspace(-R, R, N); I = zeros(N, N); % 初始化 mask：0 黑，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 都用 mm，p 也是 mm；N=1001 时 I 是 1001×1001 的 double 矩阵；colormap(gray(256)) 能显示 0、0.5、1 三档灰度，若只要二值图，可以去掉背景灰。

3.3 向量化版：一行计算 k，两行生成 mask

双重循环在 N=1001 时还能接受，但做参数扫描或提高分辨率时，向量化写法会省很多时间。

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(P>R)=0.5用逻辑索引批量赋值，比循环里逐个判断更简洁。向量化版本没有显式循环，但内存占用更大：N=4001 时单个矩阵约 128 MB，如果机器内存吃紧，可以分块处理。

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 中Ir=I*N; image(lr);大概率是 OCR 误读，按上下文看，实际就是对 mask 做显示，不要把I*N当成物理公式。

4. 复现原代码时最常翻车的点：索引、尺寸不匹配与波带数边界

4.1 x 长度 1010 与 y 长度 1001 的尺寸不匹配

原 doc 中有一行x=linspace(-xm,xm,1010); y=linspace(-ym,ym,1001);，但循环用的是for m=1:1001和for n=1:1001。这意味着 x 只取了前 1001 个点，实际 x 范围不再是 -3 到 3，而是 -3 到约 2.988。波带片在水平方向被轻微压缩，中心点也偏了。这种错误在图上不一定一眼看出，但外环半径会整体偏小。

N = 1001; x = linspace(-R, R, N); y = linspace(-R, R, N);

参数说明：x、y 必须同长度，且 N 取奇数，中心点才正好落在 x=0、y=0。N 取偶数时，中心落在两个像素之间，中心波带会出现不对称。若必须用偶数 N，可以手动把中心点索引附近的像素做对称修正，但更省事的做法是直接用奇数。

4.2 p=R 边界与 k 取整：fix、floor、round 的差别

fix对正数截断小数，等价于 floor；对负数向零取整。p 是非负的，所以fix(p^2/(lam*f))和floor(p^2/(lam*f))结果一致。边界处浮点误差会让 k 差 1，比如商算出来是 3.0000000001 和 2.9999999999，前者 k=3，后者 k=2，奇偶相反，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((N+1)/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))));

说明：中心行取的是(N+1)/2行，对应 y=0。diff(row)~=0找灰度变化的位置，abs(x(edge_idx))得到跳变半径。理论半径由 sqrt(klamf) 算出，前几个波带误差应在半个像素以内。若误差明显偏大，优先检查 x、y 长度是否一致，以及 lam、f、R 的单位是否统一。

5. 从静态 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/mm；lambda*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((N+1)/2, (N+1)/2))^2; end [~, idx_max] = max(I_axis); fprintf('设计焦距 %.1f mm，实测轴上峰值位置 %.1f mm\n', f, z_list(idx_max));

这个扫描能暴露采样是否足够。如果峰值偏离设计焦距超过几个采样间隔，优先检查波长、焦距单位是否统一，以及 mask 的奇偶类型有没有选对。不同焦距对应的外环半径也可以先算表，再决定 N 和 R 是否覆盖完整孔径。

把 z_list 和 I_axis 存成两列 CSV，再画归一化曲线，比直接看 hot 图更容易发现次焦点。

本文还有配套的精品资源，点击获取