这一阵子在做一个光学衍射传播的仿真任务核心需求很直接把一张灰度肖像图当作物平面光场用 MATLAB 算出它经过一段自由空间传播后的衍射强度图并定量解析这些图案的特征。选来选去最终还是落在角谱传输Angular Spectrum Propagation这套方法上。跑通之后我最大的感受是角谱法在工程仿真里实在太顺手了比直接写菲涅尔积分公式省心得多精度和适用范围也更好。所以我把整个实现过程、参数推导、代码细节和调试中踩过的坑都整理出来给做计算光学、数字全息、衍射仿真方向的同学做个参考。这篇内容适合两类人一类是刚接触衍射传播、想快速把菲涅尔衍射跑起来的初学者另一类是已经有仿真经验、但在某个具体环节采样、混叠、边界效应卡住的工程师。全文按照从理论到落地的顺序展开代码可以直接拿走改参数用。1. 为什么用角谱法而不是直接写菲涅尔积分1.1 角谱法的核心思想一句话讲透角谱法的物理图像非常直观。任何标量光场都可以分解成无数个沿不同方向传播的平面波的叠加每个平面波对应一个空间频率组合 (fx, fy)这个空间频率就决定了平面波的传播方向。光场从物平面传播到观察平面本质上就是每个平面波在 z 方向上多走了一段距离多积累了一个相位。换句话说衍射传播在频域里就是一个逐点相乘的操作根本不需要在空间域做卷积。具体实现就三步对输入光场做二维傅里叶变换得到角谱乘上自由空间传递函数再做一次逆傅里叶变换回到空间域。全程是 FFT 操作计算量在 O(N log N) 级别1024×1024 的复数矩阵单次传播在普通笔记本上也就是几十毫秒。这个效率优势在做多距离扫描、参数优化的时候非常明显。角谱法为什么会叫角谱因为平面波的传播方向可以用它与三个坐标轴的夹角来刻画而这些夹角又直接映射到空间频率上。角谱就是按角度展开的频谱这个含义。理解了这一点后面再看传递函数的表达式就不会觉得突兀了。1.2 和直接菲涅尔积分对比适用范围差很多教科书上最常见的菲涅尔衍射公式是那个带二次相位因子的卷积积分U2(x, y) exp(jkz) / (jλz) · ∫∫ U1(x1, y1) · exp{ jk[(x-x1)² (y-y1)²] / 2z } dx1 dy1这个公式当然可以数值实现利用卷积定理也能用 FFT 加速。问题在于它成立的前提是菲涅尔近似条件具体来说要求传播距离 z 满足 z³ 远大于 [(x-x1)² (y-y1)²]² 的某个量级而且当 z 很小的时候二次相位因子 exp(jk(x-x1)²/2z) 在空间域振荡极快离散采样根本跟不上振荡速度结果全是混叠噪声。角谱法没有这个烦恼。它在标量衍射理论框架内是相对严格的从 z0 的紧贴面到远场都能用。实际仿真中我经常要扫描几毫米到两百毫米的距离范围如果用菲涅尔积分公式不同距离段可能还要换不同的近似形式而角谱法从头到尾一套代码通吃只需要改 z 这一个变量。工程维护成本低到可以忽略不计。1.3 角谱法和菲涅尔近似的内在联系角谱法其实可以推导出菲涅尔近似。把传递函数中的平方根做泰勒展开保留到一阶项再回到空间域就能得到菲涅尔衍射的卷积形式。所以在菲涅尔近似成立的区域两种方法的结果高度一致在近场区域角谱法依然有效而菲涅尔近似会失效。这也是我后来坚持用角谱法的根本原因——同样的一个函数从近场到远场全覆盖不需要为不同距离写两套代码。当然角谱法也有自己的适用限制主要是采样条件的约束。这一点很多人会忽略后面用一整节说清楚。2. 仿真参数设计动笔写代码之前先把网格算明白2.1 五个核心参数一个都不能含糊做角谱传播仿真一开始就要定死五个参数波长 λ、像素间距 dxdy、采样点数 M×N、传播距离 z、物面尺寸 Lx×Ly。其中 Lx 和 Ly 不是独立变量它们由采样数和像素间距决定Lx N·dxLy M·dy。所以真正需要手动确定的只有 λ、dx、N 和 z。波长容易理解可见光常用 532nm 绿光或者 633nm 红光。像素间距这个参数在真实系统里往往对应空间光调制器的像素尺寸比如常见的 8μm 或 4μm。采样点数 N 要结合物面大小和计算量来平衡N 太小频谱分辨率不够N 太大内存和计算时间都会涨。我做个人电脑上的仿真习惯用 1024×1024这个规模下复数矩阵占内存约 16MB单次 FFT 很快精度也足够。2.2 角谱法自己的采样约束公式角谱法的传递函数 H(fx, fy) exp(jkz·sqrt(1 - λ²fx² - λ²fy²))在频域是一个快速振荡的复指数函数。既然是数字离散计算就必须保证频域采样间隔足够小能够采到它的振荡细节。频域采样间隔由空间域视场决定Δfx 1/(N·dx)。从奈奎斯特条件的角度推导要求传递函数在相邻两个频域采样点之间的相位变化不超过 π。对传递函数求频率的偏导可以得到它的局部频率近似处理后能得到一个很实用的上限公式zmax ≈ N·dx² / λ这个公式和直接菲涅尔数值积分里的采样条件形式很像但物理来源完全不同——一个是二次相位因子在空间域振荡太快一个是角谱传递函数在频域振荡太快。要注意区分。2.3 实际算例与参数速查表拿我实际用的一组参数举例λ 532nmN 1024dx 8μm。zmax 1024 × (8×10⁻⁶)² / (532×10⁻⁹) 1024 × 64×10⁻¹² / 532×10⁻⁹ ≈ 0.123m也就是说在 1024×1024、8μm 像素的参数下角谱法单次传播最多算到 123mm。想算更远的距离就得把传播过程拆成多段每一段都小于这个上限依次传播。不同参数组合下的单次最大传播距离我整理成了速查表采样数 N像素间距 dx波长 λ单次最大传播距离 zmax5128 μm532 nm61.6 mm10248 μm532 nm123.2 mm20488 μm532 nm246.4 mm10244 μm532 nm30.8 mm注意最后一行像素间距减半之后zmax 缩小到原来的四分之一因为公式里是 dx² 的关系。所以加密像素不是没有代价的——要么牺牲单次传播距离要么增加分段次数。设计仿真参数时这个权衡要心里有数。3. MATLAB代码实现核心传播函数与完整流程3.1 角谱传播函数的标准写法先给出我封装好的核心函数。这个函数我用了很久各种边界情况都处理过直接拿走就能用function U2 angSpectrumProp(U1, lambda, z, dx, dy) % 角谱传播函数 % U1 : 输入光场MxN 复数矩阵 % lambda: 波长m % z : 传播距离m必须为正数 % dx, dy: 采样间距m % U2 : 输出光场MxN 复数矩阵 [M, N] size(U1); % 频域坐标注意与数组维度对应 fx (-N/2 : N/2-1) / (N * dx); fy (-M/2 : M/2-1) / (M * dy); [FX, FY] meshgrid(fx, fy); k 2 * pi / lambda; % 传递函数 kz2 1 - (lambda * FX).^2 - (lambda * FY).^2; valid kz2 0; % 倏逝波区域直接忽略 H zeros(M, N); H(valid) exp(1i * k * z * sqrt(kz2(valid))); % FFT - 乘传递函数 - IFFT U1F fftshift(fft2(ifftshift(U1))); U2F U1F .* H; U2 fftshift(ifft2(ifftshift(U2F))); end有几个细节值得专门说明。第一频域坐标的定义。fx 和 fy 的数组中心是零频但 MATLAB 的 fft2 默认把零频放在左上角所以必须在变换前后用 ifftshift 和 fftshift 做对称的搬移。很多初学者只做一次 fftshift结果相位全乱衍射图完全不对。第二meshgrid 的参数顺序。输入矩阵 U1 是 M×N行对应 y 方向列对应 x 方向。meshgrid(fx, fy) 生成的 FX 会让每一行内部的 x 频率变化FY 会让每一列的 y 频率变化这样乘法才能和矩阵维度对齐。这个顺序写反了结果也是错的而且错得非常隐蔽光看强度图不一定能发现。第三倏逝波的处理。当 (λfx)² (λfy)² 1 时kz 变成纯虚数对应倏逝波。物理上倏逝波在离开物面几个波长之后就衰减到可以忽略所以在自由空间传播的仿真中直接把这一部分传递函数置零是合理的。不置零的话exp(负实数·z) 的数值虽然 MATLAB 也能处理但没有物理意义纯属浪费计算量。3.2 肖像图预处理随机相位这个关键步骤很多人第一次做衍射仿真时会把灰度图 I(x,y) 直接当成振幅也就是 U1 sqrt(I)然后就去传播了。这样跑出来的结果非常干净中心一个亮斑周围是规则的同心环一点也不像实际光学实验里看到的东西。原因在于真实的物体表面是粗糙的光被散射后每个微观面元都带有一个随机相位所以整个物光场的相位是杂乱无章的。纯振幅、零相位的物体只存在于理想化模型里现实中很难找到。正确做法是给物场附加随机相位img imread(portrait.png); img rgb2gray(img); img imresize(img, [1024, 1024]); img double(img) / 255; U0 sqrt(img) .* exp(1i * 2 * pi * rand(1024, 1024));rand 生成 0 到 1 均匀分布的随机数乘上 2π 就是 0 到 2π 的随机相位均匀覆盖整个相位周期。这一步做完之后衍射得到的强度图才会有自然的散斑颗粒感能量分布也更接近真实情况。相位分布用均匀随机就够了。如果想让散斑的尺寸受控可以用空间相关的高斯随机相位但那属于更精细的仿真需求一般场景用不上。3.3 主程序与多距离扫描有了传播函数和预处理后的物场主程序就非常简单了。我一般会扫描一组传播距离观察衍射图案的演化clear; close all; lambda 532e-9; dx 8e-6; dy 8e-6; N 1024; % 读取并预处理肖像图 img imread(portrait.png); img rgb2gray(img); img imresize(img, [N, N]); img double(img) / 255; U0 sqrt(img) .* exp(1i * 2 * pi * rand(N, N)); % 多距离传播 zList [0.005, 0.01, 0.02, 0.05, 0.1]; for iz 1:length(zList) Uz angSpectrumProp(U0, lambda, zList(iz), dx, dy); Iz abs(Uz).^2; Iz Iz / max(Iz(:)); figure; imagesc(Iz); axis image; colormap gray; title(sprintf(z %.0f mm, zList(iz)*1000)); imwrite(uint8(Iz*255), sprintf(prop_%.0fmm.png, zList(iz)*1000)); end这段代码跑起来很快五个距离总共不到两秒。如果想做参数扫描比如扫描波长或者像素尺寸外层再套一层循环即可。3.4 超长距离传播分段处理的标准姿势如果目标距离超过了 zmax比如刚才的例子想到 500mm单次角谱传播就不行了。解决办法是分段传播把长距离切成若干段每段用一次角谱传播逐段迭代function Uz propagateSafely(U0, lambda, zGoal, dx, dy) % 安全分段传播避免超出采样约束 zmax min(size(U0)) * dx^2 / lambda; % 保守估计 nSeg ceil(zGoal / (0.5 * zmax)); % 每段不超过zmax的一半 zStep zGoal / nSeg; U U0; for i 1:nSeg U angSpectrumProp(U, lambda, zStep, dx, dy); end Uz U; end每段步长取 zmax 的一半留出安全余量。分段传播本质上是把多次传播的传递函数相乘而角谱传递函数满足叠加性质H(z1 z2) H(z1) · H(z2)所以分段传播在数学上是精确的不会引入额外误差。这也是角谱法的一个隐性优势——它的传递函数具有半群性质分段毫无心理负担。4. 菲涅尔肖像图的获取与特征解析4.1 不同距离下的图案形态演化把肖像图我用的是一张公开渠道下载的普通灰度人像测试图作为物场按 3.3 的代码扫描五个距离得到的强度图就是菲涅尔肖像图。不同距离下的形态差异非常明显。z 5mm 时衍射效应很弱图案基本还能看出人脸的轮廓只是边缘略微软化像是对焦轻微不准的照片。z 10mm 时边缘开始出现淡淡的振铃纹路图像整体依然可辨认但细节处已经有菲涅尔条纹的雏形。z 20mm 时明显的菲涅尔衍射环带叠加在肖像上明暗交替的条纹开始主导视觉高频细节被明显抹平。z 50mm 时衍射条纹非常丰富肖像的轮廓还在但内部细节已经大面积模糊散斑颗粒感变得很强。z 100mm 时基本只剩下低频包络和粗颗粒散斑人脸轮廓变成一团光晕中的形状原始信息几乎完全被打散。这个演化过程直观地展示了菲涅尔衍射区域从哪个距离开始进入主导地位。仿真结果和光学教材上描述的现象完全吻合。4.2 能量守恒检查验证算对没算对的第一道关卡拿到衍射结果之后第一件事不是看图而是检查能量守恒。自由空间传播不吸收能量理想情况下输出光场总能量应该等于输入光场总能量energy0 sum(abs(U0).^2, all); energyZ sum(abs(Uz).^2, all); ratio energyZ / energy0;正常情况下 ratio 应该非常接近 1。因为倏逝波被置零了会有极微小的能量损失但自由光场占比通常在 99.9% 以上所以 ratio 落在 0.999 附近都算正常。如果发现能量明显下降比如低于 0.95那一定是采样条件出了问题或者频域坐标定义错了导致大部分频谱被置零。这个检查是我调试代码时最先做的自检步骤能拦下绝大多数低级错误。4.3 强度图频谱特征与等效带宽分析解析光场传播的变化不能只看强度图还要看它的空间频谱。我习惯对每个距离的强度图再做一次 FFT分析功率谱分布。这里有一个容易误解的点角谱传播的传递函数只改变频率分量的相位不改变振幅谱所以光场的振幅谱在传播前后逐点不变。但强度图的频谱就不一样了——强度是光场模平方引入了非线性操作会带来新的频率成分所以不能用线性滤波的观点去理解强度图的变化。定量比较不同距离的强度图可以计算功率谱的等效带宽IzF fftshift(fft2(Iz)); P abs(IzF).^2; % 计算功率谱的二阶矩等效带宽 [Ny, Nx] size(P); fxAxis (-Nx/2 : Nx/2-1) / (Nx * dx); fyAxis (-Ny/2 : Ny/2-1) / (Ny * dy); [FXg, FYg] meshgrid(fxAxis, fyAxis); f2 FXg.^2 FYg.^2; bandwidth sqrt(sum(f2 .* P, all) / sum(P, all));随着 z 增大等效带宽逐渐下降说明强度图的高频细节确实在传播中被打散、抹平。这个指标可以把看起来模糊了这种主观感受变成一条定量曲线写报告或者做对比实验的时候特别好用。4.4 逆向解析从菲涅尔肖像图还原原始图像获取菲涅尔肖像图之后最有意思的解析操作是试着还原原图。角谱传播函数天然支持反向传播——只需要把传播距离取负值Urec angSpectrumProp(Uz, lambda, -z, dx, dy); Irec abs(Urec).^2;物理上这相当于把衍射光场沿原路传回去。如果正向传播的代码没有问题还原出来的强度图应该和原始肖像高度接近。我把这个正向传播 反向还原的闭环当作算法正确性的终极验证手段。闭环还原度越高说明正向实现越可靠。不过要注意还原的目标是加过随机相位的物场 U0不是原始灰度图。因为随机相位在正向传播时就已经引入反向传播会把它还原回来此时取振幅平方得到的强度图应该接近原始灰度图的平方。如果直接拿原始灰度图做对比会有系统性的偏差别被这个迷惑了。这套正反向闭环在数字全息、光学加密、相位恢复等领域都有直接应用场景属于角谱法的经典操作。5. 常见问题与调试技巧实录5.1 衍射图中心出现周期性条纹这是最常见的翻车现场。表现是输出强度图中央出现十字亮纹或者周期性的网格条纹非常扎眼。原因在于 FFT 默认把输入当成周期信号处理肖像图左右边缘、上下边缘的亮度如果差异过大就相当于在周期边界上引入了一条接缝衍射时这条接缝会产生十字形的高频成分。解决办法有两个思路。一是给图像周围加一圈黑边让边界处的亮度自然过渡到 0二是做边缘羽化用余弦窗把图像四周渐变到 0。我实践下来黑边方式更简单直接padSize 50; U0(:, 1:padSize) 0; U0(:, end-padSize1:end) 0; U0(1:padSize, :) 0; U0(end-padSize1:end, :) 0;注意这个操作要在加随机相位之后做否则又会在边界引入相位突变。5.2 传递函数出现 NaN 或者无穷值sqrt 的输入为负数时结果就是虚数但如果你用的是 sqrt(k2) 而没有提前过滤 valid 区域当 k2 恰好非常接近 0 时可能产生 NaN。我在 3.1 的代码里用 valid 掩模把倏逝波区域直接置零就是为了规避这个问题。另外要注意 exp 的参数如果出现 NaN结果就是 NaN而且会通过后续的 IFFT 扩散到整个输出场排查起来非常隐蔽。所以传递函数的构建必须一步到位不要分步凑。5.3 传播距离超出单次上限这是角谱法最常见的误用。表现是距离稍长一点输出就变成一团随机噪声完全看不出衍射特征。判断方法很简单算一下 zmax N·dx²/λ如果目标距离超过它就必须分段。我见过不少人在网上问为什么角谱法到 200mm 就崩了十有八九都是没做分段处理。分段的思想在 3.4 已经给出代码。还有一个实践建议分段数宁多勿少每段距离取 zmax 的 1/3 到 1/2不要卡着上限走。虽然理论上等于 zmax 也能算但工程上留余量可以避免各种数值边缘问题。5.4 衍射图中心出现异常亮斑如果输出强度图中心有个特别亮的小圆斑周围是正常的衍射图样通常说明随机相位没有加成功。我遇到过一次是 U0 sqrt(img) 之后下一行代码不小心用了 abs(U0)把相位又抹掉了。还有一次是随机数种子固定了但数组维度写错相位矩阵和振幅矩阵大小不一致广播之后相位全是 0。排查思路很直接显示一下 U0 的相位图如果是一片杂乱的颜色说明相位正常如果全是同一个值说明相位丢了。5.5 换了采样点数结果完全对不上同一个物场1024 和 2048 算出来的衍射图差别很大这是正常现象不是代码 bug。角谱法的结果强烈依赖采样网格和视场大小因为视场 Lx N·dx 决定了频域采样间隔频域采样间隔又决定了传递函数的离散化精度。正确的做法是在物面尺寸 Lx 固定的前提下增大 N相当于减小 dx做收敛性测试观察衍射图样是否稳定。如果图案不再明显变化说明当前分辨率已经收敛。如果始终变化那就要检查是否某个环节的近似条件不满足。这种收敛性验证在正式跑参数扫描之前做一次能避免后面大量无效计算。5.6 关于复振幅显示的一个小技巧最后分享一个显示技巧。角谱传播得到的 Uz 是复数场直接看 abs(Uz) 只能看到振幅信息相位信息全丢了。想看相位分布用 angle(Uz)但相位图通常是密集的随机花纹不好解读。更实用的做法是显示干涉图样——把 Uz 和参考平面波叠加Iinterfere abs(Uz exp(1i * 2 * pi * 500 * X)).^2;其中 X 是空间坐标网格500 是参考光的空间频率这样能看到类似离轴全息图的干涉条纹相位信息以条纹弯曲的形式呈现可视化效果好很多。做全息方向的同学应该会经常用到这个技巧。这套角谱传播仿真方法我后来还扩展到了空间光调制器调制模拟、计算全息图的生成验证等多个场景。核心的 angSpectrumProp 函数始终没变变的只是输入光场和距离参数。回过头来看当初花时间把采样条件、边界处理这些问题想透是非常值得的因为它们直接影响仿真结果的可信度。如果你只是需要快速验证某个想法直接复制 3.1 的传播函数和 3.3 的主程序改改参数就能用。如果要做更严格的科研级仿真建议从参数设计开始认真走一遍第 2 节的推导过程。
