简介这份资源面向工程光学、物理光学课程的学习者与MATLAB仿真入门者聚焦圆孔与矩孔两类孔径的菲涅尔衍射数值模拟帮助读者理解衍射原理并掌握编程实现方法。压缩包内仅含1个doc文档约263KB以图文与代码结合的方式组织内容便于在电脑端阅读与对照练习。文档从惠更斯-菲涅尔原理与基尔霍夫修正出发推导菲涅尔衍射的傅里叶变换表达式再以二维矩阵存储衍射屏与观察屏场分布给出采样点数、屏尺寸与传播距离的选取建议。圆孔部分以He-Ne激光波长632.8nm为例展示半径12mm、20mm、50mm三组衍射图样矩孔部分给出长宽15mm、20mm、30mm的模拟结果并附完整MATLAB程序涵盖衍射屏构造、二次相位因子计算与图像显示。目前已有267人学习适合需要完成课程综合作业或快速复现衍射仿真的读者参考。1. 圆孔矩孔的菲涅尔衍射模拟为什么你算的衍射图样总是一片糊做工程光学的人迟早会撞上菲涅尔衍射这道坎。课程设计要交图样实验要对条纹间距项目里要评估孔径对光路的影响——可当你兴冲冲打开 MATLAB把圆孔或矩孔的参数代进去跑出来的结果要么是一团均匀的亮斑要么是边缘糊成一片的噪声跟教科书上那圈漂亮的同心环、那排清晰的明暗带完全对不上。问题往往不在物理公式而在采样、在坐标、在数值积分的边界处理。这篇笔记就围绕圆孔和矩孔的菲涅尔衍射模拟把从理论到 MATLAB 落地的完整路径拆开讲怎么选衍射积分形式、怎么设采样网格、怎么处理观察屏坐标、怎么避开那些让图样翻车的参数陷阱。适合正在做工程光学课程设计、需要快速拿到可信衍射图样的同学也适合想把衍射计算嵌进自己光路仿真流程的工程师。2. 菲涅尔衍射的数值实现从积分公式到可跑的 MATLAB 代码菲涅尔衍射的核心是基尔霍夫衍射公式在近场近似下的简化形式。对于圆孔和矩孔这类规则孔径直接套用菲涅尔衍射积分就能得到观察屏上的复振幅分布。但要把这个积分变成 MATLAB 里能跑、跑得对、跑得快的代码需要先想清楚三件事用哪种积分形式、怎么离散化、怎么处理坐标原点。2.1 选卷积形式还是直接积分两种菲涅尔衍射计算路径的取舍菲涅尔衍射积分有两种常见写法。一种是直接对孔径平面做二重积分观察屏上每一点都要遍历整个孔径区域另一种是利用菲涅尔衍射的卷积形式把衍射过程看成孔径透射函数与一个二次相位因子的卷积然后用 FFT 加速。直接积分的物理意义最直观代码也最好写但计算量随采样点数平方增长观察屏点一多就慢得让人想砸键盘。卷积形式借助 FFT速度能快好几个数量级但需要小心处理采样间隔和频率坐标的对应关系否则会出现混叠或者相位符号错误。我一般会这样选如果只是做课程设计、观察屏点数在 512×512 以内、孔径尺寸和观察距离都在厘米量级直接用二重积分循环就够了代码短、调试直观不容易在坐标映射上翻车。如果要做参数扫描、要迭代优化孔径形状、或者观察屏点数上千那就必须上 FFT 卷积形式。下面先给一个直接积分的圆孔菲涅尔衍射实现把物理量都摆清楚再讨论 FFT 版本的注意事项。% 圆孔菲涅尔衍射 - 直接二重积分实现 % 波长、孔径半径、观察距离 lambda 632.8e-9; % 氦氖激光波长单位米 a 1e-3; % 圆孔半径1 mm z 0.5; % 观察距离0.5 m k 2*pi/lambda; % 波数 % 孔径平面采样 N_ap 200; % 孔径平面每边采样点数 L_ap 4*a; % 孔径平面采样范围取孔径的4倍 x_ap linspace(-L_ap/2, L_ap/2, N_ap); y_ap x_ap; [X_ap, Y_ap] meshgrid(x_ap, y_ap); % 圆孔透射函数孔内为1孔外为0 aperture double(X_ap.^2 Y_ap.^2 a^2); % 观察屏采样 N_obs 128; % 观察屏每边点数先设小一点看速度 L_obs 5e-3; % 观察屏范围5 mm x_obs linspace(-L_obs/2, L_obs/2, N_obs); y_obs x_obs; [X_obs, Y_obs] meshgrid(x_obs, y_obs); % 预计算孔径平面上的二次相位因子 phase_ap exp(1i*k/(2*z)*(X_ap.^2 Y_ap.^2)); % 二重积分对每个观察点遍历孔径 U_obs zeros(N_obs, N_obs); dx_ap x_ap(2) - x_ap(1); dy_ap y_ap(2) - y_ap(1); for m 1:N_obs for n 1:N_obs % 观察点到孔径上各点的距离引起的相位 phase_obs exp(1i*k/(2*z)*((X_obs(m,n)-X_ap).^2 ... (Y_obs(m,n)-Y_ap).^2)); % 积分核孔径函数 × 孔径二次相位 × 观察点二次相位 integrand aperture .* phase_ap .* phase_obs; U_obs(m,n) sum(integrand(:)) * dx_ap * dy_ap; end end % 归一化并计算强度 I_obs abs(U_obs).^2; I_obs I_obs / max(I_obs(:)); % 显示 figure; imagesc(x_obs*1e3, y_obs*1e3, I_obs); colormap(hot); colorbar; xlabel(x (mm)); ylabel(y (mm)); title(圆孔菲涅尔衍射图样); axis square;这段代码的逻辑很直接先定义波长、孔径半径和观察距离然后在孔径平面和观察屏上分别建立采样网格。圆孔透射函数用X_ap.^2 Y_ap.^2 a^2生成一个二值矩阵孔内为 1孔外为 0。核心循环里对观察屏上每一个点计算它到孔径平面所有采样点的二次相位差乘上孔径函数和孔径自身的二次相位因子求和后再乘采样面积。最后取模平方得到强度归一化后显示。参数上有几个关键点。N_ap决定孔径平面的采样密度太小会导致圆孔边缘变成锯齿状衍射图样出现虚假的高频条纹一般取到孔径直径的 50 到 100 倍采样比较稳妥。L_ap要至少覆盖孔径的 2 到 3 倍让孔径外的零场区域也被采样到否则积分边界会截断相位因子引入误差。L_obs要根据菲涅尔衍射的尺度来估第一暗环半径大约在sqrt(lambda*z)量级对于 632.8 nm、0.5 m 的距离这个值约 0.56 mm所以观察屏取 5 mm 见方是合理的。N_obs在调试阶段可以设小一点比如 64 或 128确认图样形状对了再加大。2.2 矩孔衍射的坐标分离技巧把二重积分拆成两个一维积分矩孔的菲涅尔衍射有一个圆孔没有的好处矩形孔径的透射函数在 x 和 y 方向是可分离的。如果孔径是边长为 2a 和 2b 的矩形透射函数可以写成rect(x/(2a)) * rect(y/(2b))菲涅尔衍射积分里的二次相位因子也是exp(i*k/(2z)*(x^2y^2))的形式可以拆成 x 方向和 y 方向两个因子的乘积。这意味着二重积分可以拆成两个一维积分的乘积计算量从 N^2 降到 2N速度提升非常明显。不过要注意这个分离只在孔径中心与光轴重合、且观察屏坐标原点也在光轴上时才严格成立。如果孔径有横向偏移或者观察屏要算偏心区域分离就不成立了得老老实实做二重积分。下面给一个利用分离性质的矩孔衍射实现同时保留一个可切换的二重积分版本做对照。% 矩孔菲涅尔衍射 - 利用坐标分离的一维积分实现 lambda 632.8e-9; a 0.5e-3; % 矩孔半宽x方向 b 1e-3; % 矩孔半高y方向 z 0.8; % 观察距离 k 2*pi/lambda; % 孔径平面一维采样 N_ap 400; L_ap 6e-3; x_ap linspace(-L_ap/2, L_ap/2, N_ap); dx x_ap(2) - x_ap(1); % 矩孔透射函数的一维形式 rect_x double(abs(x_ap) a); rect_y double(abs(x_ap) b); % 观察屏一维采样 N_obs 256; L_obs 4e-3; x_obs linspace(-L_obs/2, L_obs/2, N_obs); % 预计算孔径平面的二次相位 phase_ap_x exp(1i*k/(2*z)*x_ap.^2); phase_ap_y exp(1i*k/(2*z)*x_ap.^2); % x方向一维积分 U_x zeros(1, N_obs); for m 1:N_obs phase_obs_x exp(1i*k/(2*z)*(x_obs(m) - x_ap).^2); integrand_x rect_x .* phase_ap_x .* phase_obs_x; U_x(m) sum(integrand_x) * dx; end % y方向一维积分 U_y zeros(1, N_obs); for n 1:N_obs phase_obs_y exp(1i*k/(2*z)*(x_obs(n) - x_ap).^2); integrand_y rect_y .* phase_ap_y .* phase_obs_y; U_y(n) sum(integrand_y) * dx; end % 二维复振幅由两个一维结果外积得到 U_obs U_y. * U_x; % 注意维度y方向列向量 × x方向行向量 I_obs abs(U_obs).^2; I_obs I_obs / max(I_obs(:)); figure; imagesc(x_obs*1e3, x_obs*1e3, I_obs); colormap(hot); colorbar; xlabel(x (mm)); ylabel(y (mm)); title(矩孔菲涅尔衍射图样坐标分离法); axis square;这里的关键操作是把二维积分拆成两个独立的一维积分。U_x是 x 方向的一维衍射复振幅U_y是 y 方向的最终的二维复振幅U_obs是它们的外积。注意外积的顺序U_y.是列向量U_x是行向量相乘得到N_obs × N_obs的矩阵对应观察屏上的二维分布。如果写成U_x. * U_y维度就反了图样会转置 90 度这是新手很容易踩的坑。参数方面L_ap要覆盖矩孔的长边这里取 6 mm 是为了让 2 mm 高的矩孔两侧各有足够的零场区域。N_ap取 400 是为了让一维积分的精度足够因为分离法虽然快但每个一维积分的采样密度直接决定最终图样的平滑程度。L_obs取 4 mm对应第一暗纹位置大约在lambda*z/(2b)量级对于 1 mm 半高的矩孔、0.8 m 距离这个值约 0.25 mm所以 4 mm 的观察范围能覆盖十几个暗纹足够看清衍射结构。2.3 采样参数怎么定三个必须算清楚的尺度菲涅尔衍射模拟翻车十有八九是采样参数没算对。有三个尺度必须提前估出来不能靠试。第一个是菲涅尔尺度sqrt(lambda*z)。这个量决定了衍射图样的整体大小。圆孔的第一暗环半径约1.22*lambda*z/(2a)矩孔的第一暗纹位置约lambda*z/(2a)或lambda*z/(2b)。观察屏范围L_obs至少要覆盖 3 到 5 个暗纹间距否则图样会被截断看起来像只有中心亮斑。第二个是孔径平面的采样间隔dx_ap。它必须足够小才能分辨孔径边缘的相位变化。一个经验规则是dx_ap lambda*z/(2*L_ap)这来自菲涅尔衍射的局部空间频率限制。如果dx_ap太大孔径边缘的二次相位因子会欠采样积分结果出现虚假的旁瓣。我一般会让N_ap满足dx_ap小于这个上限的一半留足余量。第三个是观察屏的采样间隔dx_obs。它决定了图样的空间分辨率。如果dx_obs太大细密的衍射条纹会被平均掉看起来一片糊。对于圆孔第一暗环到第二暗环的间距约0.5*lambda*z/(2a)所以dx_obs要小于这个值的三分之一。对于矩孔暗纹间距约lambda*z/(2b)同样要留够采样点。这三个尺度互相制约。L_obs和dx_obs决定了N_obsL_ap和dx_ap决定了N_ap。如果N_obs和N_ap都很大直接积分的循环就会很慢。这时候要么换 FFT 卷积形式要么利用矩孔的坐标分离性质要么就接受较低的分辨率先看趋势。下面这个表把三个尺度的估算公式和典型值列出来方便直接套用。尺度估算公式圆孔示例a1mm, z0.5m矩孔示例a0.5mm, z0.8m菲涅尔尺度sqrt(lambda*z)0.56 mm0.71 mm第一暗纹位置圆孔1.22lambdaz/(2a)矩孔lambda*z/(2a)0.19 mm0.51 mm孔径采样上限lambdaz/(2L_ap)20 μm17 μm观察屏采样上限暗纹间距/30.06 mm0.17 mm注意表中的数值是按氦氖激光 632.8 nm 估算的。如果你用的波长不同所有尺度都要按比例缩放。比如换成 532 nm 的倍频激光菲涅尔尺度会缩小到原来的 0.92 倍左右。3. 从复振幅到强度图相位因子、归一化和显示的那些坑拿到复振幅矩阵只是第一步怎么把它变成一张能看的、物理上正确的强度图还有不少细节要处理。相位因子的符号、归一化的方式、显示时的坐标轴和色标每一个都可能让图样看起来不对劲。3.1 二次相位因子的符号正负号搞反图样会怎样菲涅尔衍射积分里的二次相位因子是exp(±i*k/(2z)*(...))符号取决于你采用的约定是exp(-i*omega*t)还是exp(i*omega*t)。在 MATLAB 里如果符号搞反了强度图样abs(U).^2其实不会变因为正负号只是复共轭的关系模平方是一样的。但如果你要看相位分布或者做干涉叠加符号就至关重要了。我遇到过一种情况把圆孔衍射的复振幅和另一个光路的复振幅叠加结果干涉条纹的弯曲方向跟实验对不上。查了半天才发现是二次相位因子的符号跟实验约定反了。所以建议在代码开头就明确注释你用的约定比如% 采用 exp(-i*omega*t) 约定二次相位因子取正号这样后面调试干涉问题时不会抓瞎。对于纯强度图样符号不影响结果但要注意孔径平面和观察屏的二次相位因子必须用同一个约定。如果孔径平面用了正号观察屏用了负号那相当于在积分核里引入了一个额外的共轭因子强度图样会变成原本的镜像圆孔还是圆孔看不出来但矩孔的长宽方向可能会对调或者偏心孔径的图样会跑到错误的位置。3.2 归一化与显示为什么你的图样中心总是一片白直接算出来的强度矩阵中心点的值可能比边缘大好几个数量级。如果直接用imagesc显示中心会饱和成一片白周围的衍射环完全看不见。常见的做法是除以最大值做归一化但这样边缘的弱条纹还是会被压得很暗。更好的方式是用对数显示或者手动设置色标的范围。% 强度归一化与对数显示 I_obs abs(U_obs).^2; I_norm I_obs / max(I_obs(:)); % 方式一线性归一化适合看主瓣和第一旁瓣 figure; imagesc(x_obs*1e3, x_obs*1e3, I_norm); colormap(hot); colorbar; caxis([0 0.05]); % 只显示最大值的5%以下让弱条纹可见 axis square; % 方式二对数显示适合看高阶衍射环 I_log 10*log10(I_norm 1e-6); % 加一个小量避免log(0) figure; imagesc(x_obs*1e3, x_obs*1e3, I_log); colormap(hot); colorbar; caxis([-40 0]); % 显示40 dB的动态范围 axis square;线性显示时caxis设成[0 0.05]能把弱条纹提亮但主瓣会过曝。对数显示用10*log10把动态范围压缩到 40 dB主瓣和高阶环能同时看清。加1e-6是为了避免强度为零的点取对数变成负无穷。这两种显示方式各有用途线性图适合量暗纹位置对数图适合看整体衍射结构。还有一个容易忽略的点是坐标轴。imagesc默认的坐标是像素索引不是物理坐标。必须把x_obs和y_obs乘上合适的单位比如毫米再传给imagesc否则你量出来的暗纹间距是像素数不是毫米。另外axis square要加上不然圆孔的衍射环会被拉成椭圆。3.3 用 FFT 加速卷积形式的采样条件与混叠避免当观察屏点数超过 512×512或者要做参数扫描时直接积分的循环就太慢了。这时候用 FFT 实现菲涅尔衍射的卷积形式是必然选择。基本思路是把菲涅尔衍射积分写成孔径透射函数与一个二次相位因子的卷积然后利用卷积定理在频域做乘法再反变换。但 FFT 方法有两个硬性条件。第一采样间隔必须满足dx_ap * dx_obs lambda*z/N其中 N 是 FFT 点数。这个条件来自离散傅里叶变换的尺度和菲涅尔衍射的坐标映射关系。如果采样间隔不满足这个等式衍射图样的尺度就会算错暗纹位置对不上。第二孔径平面和观察屏的采样范围要足够大让卷积的周期性延拓不会把相邻周期的图样混叠进来。具体来说L_ap和L_obs都要大于孔径尺寸加上衍射扩展尺度。% 菲涅尔衍射的FFT卷积实现圆孔 lambda 632.8e-9; a 1e-3; z 0.5; k 2*pi/lambda; N 1024; % FFT点数 L_ap 8e-3; % 孔径平面范围 dx L_ap / N; % 采样间隔 L_obs lambda * z / dx; % 由采样条件确定的观察屏范围 x_ap (-N/2:N/2-1) * dx; [X_ap, Y_ap] meshgrid(x_ap, x_ap); % 圆孔透射函数 aperture double(X_ap.^2 Y_ap.^2 a^2); % 菲涅尔卷积核二次相位因子 h exp(1i*k/(2*z)*(X_ap.^2 Y_ap.^2)); % 频域卷积FFT(aperture) .* FFT(h)再IFFT H fft2(fftshift(h)); A fft2(fftshift(aperture)); U_obs fftshift(ifft2(A .* H)); % 强度与显示 I_obs abs(U_obs).^2; I_obs I_obs / max(I_obs(:)); x_obs (-N/2:N/2-1) * (lambda*z/L_ap); figure; imagesc(x_obs*1e3, x_obs*1e3, I_obs); colormap(hot); colorbar; caxis([0 0.05]); axis square; xlabel(x (mm)); ylabel(y (mm)); title(圆孔菲涅尔衍射FFT卷积法);这段代码里dx L_ap/N是孔径平面的采样间隔L_obs lambda*z/dx是由采样条件反推的观察屏范围。fftshift把零频移到中心ifft2之后再用fftshift移回来。卷积核h是二次相位因子注意这里没有取共轭因为菲涅尔衍射的卷积形式里核就是exp(i*k/(2z)*(x^2y^2))符号取决于约定。FFT 方法的速度比直接积分快两个数量级以上但代价是采样条件被锁死。如果你想自由选择观察屏范围就得调整L_ap或N来满足dx_ap * dx_obs lambda*z/N。如果强行用不满足条件的参数图样会出现混叠表现为衍射环上叠加了虚假的周期性条纹或者图样在边缘突然折叠。这是 FFT 方法最常见的翻车点调试时先检查这个等式是否成立。4. 避坑与排查圆孔矩孔衍射模拟的五个血泪教训做衍射模拟这些年踩过的坑比衍射环还多。下面这五条是反复出现、每次都要重新查一遍的典型问题按现象、原因、解决三段式写出来方便对照排查。4.1 图样中心一片白周围全黑现象跑完代码imagesc显示的图样中心是一个巨大的白色亮斑几乎看不到任何衍射环或暗纹。原因强度矩阵的动态范围太大中心主瓣的强度比第一旁瓣高两三个数量级线性显示时色标自动缩放到最大值旁瓣全被压到黑色。这不是计算错误是显示问题。解决用caxis手动限制色标范围比如caxis([0 0.05])只显示最大值的 5% 以下。或者改用对数显示10*log10(I_norm 1e-6)把动态范围压缩到 40 dB。如果对数显示下还是看不到环那才是计算有问题检查孔径采样是否太粗。4.2 圆孔衍射环变成多边形现象本该是同心圆的衍射环边缘出现明显的多边形轮廓采样点越多多边形边数越多但始终不是光滑的圆。原因孔径平面的采样网格是矩形的圆孔透射函数用X^2 Y^2 a^2生成时边缘是阶梯状的。如果N_ap不够大这个阶梯效应就会在衍射图样里显现出来高阶环尤其明显。解决增大N_ap让孔径边缘的阶梯宽度远小于波长尺度的相位变化。经验值是孔径直径上至少有 100 个采样点。如果N_ap已经很大但还有多边形检查L_ap是否太小导致孔径边缘被采样边界截断。另一个技巧是用抗锯齿的孔径函数比如用0.5*(1 - tanh((r-a)/w))代替硬边二值函数w取几个采样间隔能把边缘平滑掉。4.3 矩孔衍射图样长宽方向对调现象矩孔是横向长、纵向短但算出来的衍射图样是横向短、纵向长跟理论预期相反。原因坐标分离法里U_obs U_y. * U_x的外积顺序写反了。U_x对应 x 方向的衍射应该作为行向量U_y对应 y 方向应该作为列向量。如果写成U_x. * U_y得到的矩阵就是转置的图样长宽对调。解决检查外积的顺序。记住U_obs(m,n)对应观察屏上(x_obs(n), y_obs(m))的点所以第一个索引是 y 方向第二个是 x 方向。U_y.是N_obs × 1列向量U_x是1 × N_obs行向量相乘得到N_obs × N_obs矩阵维度正确。如果还是不确定用一个小尺寸的测试用例比如 3×3 的观察屏手动算几个点的值跟代码对比。4.4 FFT 方法算出的暗纹位置偏了现象用 FFT 卷积法算圆孔衍射第一暗环的半径跟理论值1.22*lambda*z/(2a)对不上偏差超过 10%。原因采样条件dx_ap * dx_obs lambda*z/N没有满足。FFT 方法的坐标映射是固定的dx_obs lambda*z/(N*dx_ap)。如果手动指定的L_obs跟这个值不一致暗纹位置就会按比例缩放。解决不要手动指定L_obs而是由dx_ap和N反推L_obs lambda*z/dx_ap。如果需要的L_obs跟这个值差太多就调整L_ap或N。比如想让L_obs大一点就减小dx_ap增大N或减小L_ap。调整后重新检查暗纹位置直到跟理论值吻合。4.5 观察屏边缘出现虚假的亮斑或暗带现象衍射图样的中心区域正常但观察屏的四条边缘出现不该有的亮斑、暗带或者条纹折叠。原因直接积分法里观察屏边缘的点到孔径平面的距离最大二次相位因子的空间频率最高。如果观察屏采样间隔dx_obs不够小这些高频相位会欠采样产生混叠。FFT 方法里这是卷积的周期性延拓导致的相邻周期的图样折叠到了边缘。解决直接积分法要减小dx_obs让每个相位周期至少有 4 到 8 个采样点。估算方法是dx_obs lambda*z/(2*L_ap)因为最大的空间频率出现在孔径边缘。FFT 方法要增大L_ap和L_obs让卷积的周期延拓区域落在观察范围之外。具体来说L_ap要大于孔径尺寸加上lambda*z/L_obs的扩展量。如果边缘问题不影响中心区域的读数也可以直接裁掉边缘 10% 的区域再显示。5. 进阶技巧用相位屏和角谱法把衍射模拟推到极限直接积分和 FFT 卷积能解决大部分圆孔矩孔的菲涅尔衍射问题但如果你要做更复杂的孔径、要算倾斜入射、或者要把衍射嵌入到整个光路的传播中就需要更灵活的工具。这里分享两个我常用的进阶方法相位屏法和角谱法。它们不是替代品而是在不同场景下各有优势的补充。5.1 相位屏法把任意孔径变成可复用的相位调制相位屏法的核心思想是把孔径的透射函数看成一个纯相位调制U_out U_in * exp(i*phi(x,y))。对于圆孔和矩孔这种振幅型孔径phi在孔内为 0孔外为无穷大或者用一个很大的虚部让振幅衰减到零。对于更复杂的孔径比如带相位延迟的透镜、有表面误差的光学元件相位屏法能直接把这些效应乘到光场上不需要重新推导积分公式。在 MATLAB 里相位屏法通常和角谱法配合使用。先把入射光场乘上孔径的相位屏得到孔径后的光场然后用角谱法传播到观察屏。角谱法的传播因子是exp(i*kz*z)其中kz sqrt(k^2 - kx^2 - ky^2)在频域里做乘法再反变换。这个方法对近场和远场都适用而且能处理倾斜传播。% 相位屏 角谱法传播圆孔 lambda 632.8e-9; k 2*pi/lambda; a 1e-3; z 0.5; N 1024; L 8e-3; dx L / N; x (-N/2:N/2-1) * dx; [X, Y] meshgrid(x, x); % 入射平面波 U_in ones(N, N); % 圆孔相位屏孔内相位0孔外用大虚部吸收 r sqrt(X.^2 Y.^2); absorption 1e3; % 孔外吸收系数 phi 1i * absorption * double(r a); U_ap U_in .* exp(phi); % 角谱法传播 fx (-N/2:N/2-1) / L; [FX, FY] meshgrid(fx, fx); kz sqrt(k^2 - (2*pi*FX).^2 - (2*pi*FY).^2); kz real(kz) 1i*abs(imag(kz)); % 倏逝波衰减 H exp(1i * kz * z); U_obs ifft2(ifftshift(fftshift(fft2(U_ap)) .* H)); % 强度显示 I_obs abs(U_obs).^2; I_obs I_obs / max(I_obs(:)); figure; imagesc(x*1e3, x*1e3, I_obs); colormap(hot); colorbar; caxis([0 0.05]); axis square; xlabel(x (mm)); ylabel(y (mm)); title(圆孔菲涅尔衍射相位屏角谱法);这里的关键操作是用exp(1i*absorption*double(ra))生成孔径的相位屏。孔内radouble为 0exp(0)1光场原样通过孔外radouble为 1exp(1i*absorption)是一个模为 1 的相位因子但absorption取 1e3 时exp(1i*1e3)在数值上会剧烈振荡实际效果是让孔外的光场在后续传播中因相位不匹配而迅速衰减。更稳妥的做法是用exp(-absorption*double(ra))让孔外振幅直接衰减到零但那样就不是纯相位屏了。两种方式都可以看你的物理模型需要哪种。角谱法的传播因子H exp(1i*kz*z)里kz在FX^2 FY^2 (1/lambda)^2时变成虚数对应倏逝波。代码里用real(kz) 1i*abs(imag(kz))让倏逝波的虚部为正保证exp(1i*kz*z)在传播后衰减而不是发散。这个细节如果处理错了高频噪声会被放大图样出现网格状伪影。5.2 用角谱法验证直接积分结果交叉验证的实操步骤角谱法和直接积分法在菲涅尔区的结果应该一致但数值实现不同可以作为交叉验证。我一般会这样做先用直接积分算一个 128×128 的小图再用角谱法算同样参数下的图样把两者的中心行强度曲线画在一起对比。如果暗纹位置和相对强度对得上说明两种实现都正确如果对不上分别检查采样条件和相位因子符号。% 交叉验证直接积分 vs 角谱法取中心行强度对比 % 假设已有直接积分结果 I_direct 和角谱法结果 I_angspec % 两者都是 N_obs × N_obs 的归一化强度矩阵 center_row_direct I_direct(round(N_obs/2), :); center_row_angspec I_angspec(round(N_obs/2), :); figure; plot(x_obs*1e3, center_row_direct, b-, LineWidth, 1.5); hold on; plot(x_obs*1e3, center_row_angspec, r--, LineWidth, 1.5); xlabel(x (mm)); ylabel(归一化强度); legend(直接积分, 角谱法); title(圆孔衍射中心行强度对比); grid on;对比时重点看三个地方第一暗环的位置是否一致旁瓣的相对高度是否一致远端的衰减趋势是否一致。如果第一暗环位置差了检查两种方法的坐标映射如果旁瓣高度差了检查孔径采样是否足够细如果远端衰减趋势差了检查 FFT 方法的混叠或者直接积分的截断误差。5.3 从模拟到实验用实测条纹间距反推孔径尺寸模拟的最终目的是跟实验对上。如果你手头有实测的衍射图样可以用暗纹间距反推孔径尺寸验证模拟的可靠性。对于圆孔第一暗环半径r1 1.22*lambda*z/(2a)所以a 1.22*lambda*z/(2*r1)。对于矩孔第一暗纹位置x1 lambda*z/(2a)所以a lambda*z/(2*x1)。实测时要注意观察屏上的图样可能被透镜缩放过或者相机的像素尺寸不是 1:1。我一般会先用一个已知尺寸的孔径做标定算出系统的放大倍率再测未知孔径。另外激光的波长要确认632.8 nm 是氦氖激光的标称值实际可能有零点几纳米的偏差对暗纹位置的影响在千分之一量级一般可以忽略但如果做高精度测量就要修正。提示用实测反推孔径时暗纹位置的读取误差会直接传到孔径尺寸上。第一暗环的半径最好用图像处理的方法自动提取比如沿径向做强度平均找极小值点比手动读图准确得多。最后说一个我自己的习惯每次做新的衍射模拟先不要急着算大图而是用 64×64 的小网格跑一遍把中心行的强度曲线画出来跟理论公式的解析解对比。圆孔和矩孔在菲涅尔区没有简单的解析解但在夫琅禾费区有可以把观察距离设得足够远让菲涅尔数远小于 1这时候直接积分的图样应该趋近于夫琅禾费衍射的解析结果。这个快速检查能在几秒钟内发现符号错误、采样错误或者归一化错误比盯着二维图样猜问题快得多。希望帮到你。本文还有配套的精品资源点击获取