MATLAB中的角谱法菲涅尔衍射仿真与相位恢复

做光学成像仿真的朋友大概率遇到过类似的命题光路还没搭课题要求先用仿真给出理论结果或者你只是想搞清楚某个衍射现象在距离上到底怎么演化。手里最趁手的工具就是MATLAB。前一阵我就在MATLAB里完整跑通了一个小项目——用角谱传输Angular Spectrum Propagation模拟一束相干光照射一幅灰度肖像让它在自由空间中传播在多个距离上拿到菲涅尔衍射图样再从图样里把肖像信息定量解析出来。整个项目技术点非常集中角谱法原理与实现、菲涅尔近似的适用条件、FFT 调用的正确姿势、衍射结果怎么解读。这篇文章会把完整流程、核心代码、参数设计逻辑以及调试中踩过的坑一并写下来适合正在学衍射仿真、或者要做计算成像初版验证的读者。1. 先说清楚这个项目到底在算什么1.1 角谱传输把光拆成无数“斜着走”的平面波光在自由空间传播最朴素的模型是光线直线走。一旦涉及衍射尺度光线模型就不够用了。更精确的做法是把某个平面上的光场看成一个二维信号对它做傅里叶变换。变换得到的结果里每一个频率分量都对应一束方向确定的平面波——这个频域表示就叫“角谱”。用代码表示就是U0_f fft2(U0);频域里一个点代表一束波矢确定的平面波方向由横纵频率分量共同决定。这束波沿 z 轴传播一段距离后只产生一个与频率相关的相位累积。这就是角谱法比直接做衍射积分快得多的根本原因整个传播过程可以压缩成“两次 FFT、一次频域乘法、一次 IFFT”。512×512 的图像在普通笔记本上毫秒级就能跑完一次传播这让参数扫描、迭代优化这类操作变得很轻松。如果你还在用卷积积分一个个点去算那速度和角谱法完全不是一个量级。1.2 菲涅尔近似什么时候能简化传输函数平面波沿 z 轴传播的完整相位因子是exp(1j * k * z * sqrt(1 - lambda^2 * (fx^2 fy^2)))根号那部分写起来麻烦数值计算时也要小心。当平面波方向偏离光轴不远也就是 fx、fy 都比较小的时候可以对根号做泰勒展开sqrt(1 - x) ≈ 1 - x/2代入后传播函数化简为exp(1j * k * z) .* exp(-1j * pi * lambda * z * (fx^2 fy^2))这就是菲涅尔近似下的角谱传输函数。加了这层近似表达式简洁很多数值计算也稳定绝大多数实验室尺度下的自由空间传播都适用。我后面给的完整代码默认用这个形式但也会在函數里留一个精确形式开关供需要验证近轴条件的时候切换。1.3 “肖像图”在这个项目里的角色项目标题里的“菲涅尔肖像图”本质就是拿一张人物轮廓的灰度图当输入。用肖像图而不是随便一张测试图案图的是它结构丰富、边缘清晰、人眼容易确认。仿真输出的衍射图样里还能不能认出肖像轮廓直接反映仿真链路是否正确、参数是否合理。如果传了几米距离后图案变得完全认不出那往往不是“衍射本来就模糊”而是采样条件出了问题。所以肖像在这里是一个视觉验证载体换成 cameraman、蜘蛛网图、分辨率板原理完全相同。2. 原理拆解传输函数那行复数指数到底从哪来2.1 从亥姆霍兹方程到角谱传播公式标量光场在均匀介质中满足亥姆霍兹方程。对 z 方向做傅里叶变换后可以写出任意传播距离处光场的频域表示U(fx, fy, z) U(fx, fy, 0) · H(fx, fy, z)其中 H(fx, fy, z) 就是所谓的角谱传输函数。它的完整形式是H(fx, fy, z) exp[ j * k * z * sqrt(1 - λ²(fx² fy²)) ]这个式子没有做近轴假设只要频率分量对应的平面波方向不出现倏逝波即根号内不为负它就是精确解。实际编程时如果发现根号内小于零对应的是倏逝波成分数值上直接令该处传输函数为 0 即可因为倏逝波在远超波长的传播距离上早就衰减得无影无踪。MATLAB 里的实现就是一行t 1 - lambda^2 * (FX.^2 FY.^2); H_exact exp(1j * k * z * sqrt(max(0, t)));用 max(0, t) 是为了防止根号内出现负数导致 NaN。2.2 菲涅尔近似下的传输函数推导对根号项做一阶泰勒展开后得到sqrt(1 - α) ≈ 1 - α / 2带回传输函数H exp(j·k·z) · exp(-j·π·λ·z·(fx² fy²))注意这里还有一个整体的 exp(j·k·z)它代表所有平面波共同的平均相位延迟。很多人写程序时会把这个常数因子漏掉。对强度计算来说漏掉不影响但如果后续要做干涉、叠场或者相位定量分析漏掉它就会导致相位误差随距离线性累积。2.3 离散化与采样条件为什么频率坐标要这样构造MATLAB 的 fft2 输出的频率阵列默认是从零频开始的而我们习惯的频域坐标是零频在中心。所以要先构造一套以零频为中心的频率网格dx L / N; fx (-N/2 : N/2-1) / (N * dx); [FX, FY] meshgrid(fx, fx);N·dx 就是仿真区域的边长 Lfx 的单位是 cycles/m对光学来说常用 cycles/mm。频率分辨率 Δf 1/L最大频率 fmax 1/(2dx)这个 fmax 就是离散傅里叶变换的奈奎斯特频率。角谱法在频域里做乘法必须保证传输函数 H 的相位在相邻频点之间的变化不超过 π否则频域采样不足结果会混叠。实际工程里大部分情况下只要你满足z ≤ N·dx² / λ就不会有问题。以 λ633nm、dx17.2μm、N512 为例右边算出来超过两百米实验室尺度里根本不愁。怕的是把像素尺寸设到微米量级、又传几十米那肯定出问题。2.4 角谱法与单次 FFT 菲涅尔积分法的区别很多资料里还有另一种“菲涅尔衍射的 FFT 实现”公式长这样U1 exp(jkz)/(jλz) · exp(jk(x²y²)/(2z)) · FFT{ U0 · exp(jk(x0²y0²)/(2z)) }这种方法只需要一次 FFT看着更省事但它的输出平面采样间隔被输入参数锁死没法自由指定。角谱法则不同输入输出面采用同样的采样网格可以随时改变传播距离而不用重新设计坐标迭代算法里反复正传逆传非常方便。所以做“传播距离扫描”、相位恢复迭代、多次连续传播时角谱法是首选这也是我在项目里坚持用它的原因。3. MATLAB 实操完整代码与关键参数怎么定3.1 参数表从头到尾先定一套合理参数刚开始学仿真的人最容易犯的错是上来直接抄代码参数随便填。角谱法对参数极其敏感尤其是波长、像素尺寸和距离三者的匹配关系。我建议先按下面这套参数跑通再改参数符号数值说明波长λ633 nmHe-Ne 激光常见波长波数k2π/λ代码里直接算采样规模N512偶数保证 FFT 对称性仿真边长L8.8 mm模拟一寸传感器区域像素间距dxL/N约 17.2 μm传播距离z0.001 / 0.01 / 0.05 / 0.2 / 1 m扫描观察衍射演化3.2 核心代码一个可直接复用的角谱传播函数推荐直接把传播过程封装成函数后续无论做距离扫描、逆传播还是 GS 迭代调用都方便function U_out angularSpectrum(U_in, lambda, z, dx) % 角谱法自由空间传播 % U_in : 输入平面复振幅 % lambda : 波长(m) % z : 传播距离(m)可为负值表示反向传播 % dx : 输入面像素间距(m) [Ny, Nx] size(U_in); fx (-Nx/2 : Nx/2-1) / (Nx * dx); fy (-Ny/2 : Ny/2-1) / (Ny * dx); [FX, FY] meshgrid(fx, fy); H exp(1j * 2*pi/lambda * z) .* exp(-1j * pi * lambda * z * (FX.^2 FY.^2)); U_hat fftshift(fft2(U_in)); U_out ifft2(ifftshift(U_hat .* H)); end注意最后一步的顺序fft2 之后 fftshift 把零频移到中心与按中心网格构造的 H 对齐乘法完成后 ifftshift 再把零频移回角落最后 ifft2 回空间域。这三个移位操作必须配对使用漏一个结果就是错的。3.3 三种常见写法的辨析为什么推荐上面这种我见过网上流传的另外两种写法这里对比一下避免踩坑。第一种是全程不 shift直接按 fft2 的原始频率顺序构造 H。技术上没问题但频率网格从 0 出发绕一圈容易把 fx、fy 的坐标顺序写反一旦反了图案会左右上下颠倒还很难排查。第二种是把 fftshift 放在 ifft2 之后U_out ifft2(fftshift(fft2(U_in) .* H));这种写法在绝大多数偶数尺寸情况下也能跑出正确结果因为偶数尺寸下 fftshift 与 ifftshift 等价。但一旦 N 是奇数结果就会整体错位。写代码时养成 H 居中、移位配对的好习惯后面换尺寸、换分辨率都不怕。3.4 输入图像处理灰度图怎么变成初始光场读进来的图像默认是 uint8 灰度要转成 double 并归一化I0 im2double(imread(portrait.png)); if size(I0, 3) 3 I0 rgb2gray(I0); end I0 imresize(I0, [N N]);最关键的物理细节在这里图像灰度通常代表目标的强度透过率而光场的振幅与透过率的平方根成正比。所以初始光场应该取U0 sqrt(I0); % 振幅型掩膜如果你直接把 I0 当振幅用那么重建出来的强度会变成 I0 的平方整张图对比度会被压扁边缘看起来发灰。4. 菲涅尔肖像图的获取三种方式与结果呈现4.1 方式一振幅型肖像直接衍射最简单的获取方式就是把肖像图当作一个振幅物体放进光路然后扫描传播距离Z_list [0.001, 0.01, 0.05, 0.2, 1]; for i 1:numel(Z_list) U1 angularSpectrum(U0, lambda, Z_list(i), dx); subplot(2, 3, i); imshow(abs(U1).^2, []); title(sprintf(z %.3f m, Z_list(i))); end跑出来的结果非常有层次z 很小时输出几乎就是原图的模糊版边缘变粗背景出现淡淡的光晕随着 z 增大图案逐渐失去轮廓变成一套同心环与散斑交织的复杂图样到 z1m 时原本的肖像已经很难辨认但整体强度分布仍旧保留着原图的低频信息。这时你会直观感受到什么叫“衍射随距离演化”。保存这几张图其实已经是整个项目里最有说服力的结果。4.2 方式二GS 相位恢复生成“隐形肖像”比直接衍射更有意思的是用角度谱法配合 GS 算法设计一块纯相位掩膜。这快掩膜看起来像一片随机噪声但在目标距离上经它调制的光场会自动汇聚成一张清晰的肖像图。这是计算成像里很经典的一个演示。核心迭代思路非常简单在光源面放一个纯相位物体正向传播到目标面把振幅替换成目标肖像的振幅、保留相位再反向传播回光源面把光源面振幅归一化为 1、保留相位。重复这个循环target_amp sqrt(I0); % 目标平面振幅 src exp(1j * 2*pi * rand(N, N)); % 随机初始相位 for iter 1:200 img_plane angularSpectrum(src, lambda, z, dx); img_plane target_amp .* exp(1j * angle(img_plane)); src angularSpectrum(img_plane, lambda, -z, dx); src exp(1j * angle(src)); end phase_mask angle(src); % 这就是纯相位全息图迭代到 200 次左右目标面的重建强度就已经相当接近原图了。GS 算法在 MATLAB 里跑 512×512 大概几十秒配合角谱法速度完全够用。验证时把 phase_mask 当作初始相位重新正向传播U_rec angularSpectrum(exp(1j * phase_mask), lambda, z, dx); imshow(abs(U_rec).^2, []);这时你会发现光源面看来毫无规律目标面却“凭空”出现了肖像。做课程展示或者研究汇报这个效果相当有冲击力。4.3 方式三逆传播解析把衍射图样变回肖像角谱传输函数满足一个很好的性质H(-z) conj(H(z))。所以想从衍射图样反推原始光场只需要用负距离再传播一次U_forward angularSpectrum(U0, lambda, z, dx); % 正向衍射 U_back angularSpectrum(U_forward, lambda, -z, dx); % 逆传播恢复恢复出来的强度基本就是原图但由于数值舍入、采样误差和边界截断直接恢复的图像在边缘处会有一圈振铃。这时可以对频谱做一些平滑或者在输入面边缘加一个软边窗函数。这三种方式覆盖了“生成衍射图样—设计光场—恢复原图”三种诉求也正好对应很多实际项目的三种工作模式正向仿真、逆向设计、数据恢复。三个模式共用同一个 angularSpectrum 函数代码维护成本也低。5. 结果解析从衍射图样中提取信息5.1 强度图、相位图、频谱图三种视角分别看什么拿到传播后的复振幅第一反应是画强度图 abs(U1).^2这没错。但强度图丢掉了相位信息常见的手法是把三样东西放在一起看figure; subplot(1, 3, 1); imshow(abs(U1).^2, []); title(强度图); subplot(1, 3, 2); imshow(angle(U1), []); title(相位图); subplot(1, 3, 3); imagesc(log(1 abs(fftshift(fft2(U1))))); axis image; title(频域幅度谱);相位图在随机散射场景下通常是杂乱无章的但如果传播距离很近相位分布还能看出原图的边缘轮廓说明光场还没有充分“搅拌均匀”。频域幅度谱则体现图案的空间频率成分原图高频占比低衍射后高频逐渐衰减频谱中心越来越亮。这三张图合在一起能快速判断仿真是否合理。5.2 定量指标用 PSNR 和 SSIM 量化恢复质量光靠眼睛判断不够需要量化。逆传播恢复原图后对比恢复强度与原图像素值I_rec abs(U_back).^2; I_rec I_rec / max(I_rec(:)); % 归一化 I0_norm I0 / max(I0(:)); % 需要 Image Processing Toolbox psnr_val psnr(I_rec, I0_norm); ssim_val ssim(I_rec, I0_norm); fprintf(PSNR %.2f dB, SSIM %.4f\n, psnr_val, ssim_val);通常 PSNR 在 25 dB 以上、SSIM 在 0.9 以上就说明恢复质量不错。如果 PSNR 很低先检查是不是忘了归一化或者初始图像取振幅时忘了开根号。5.3 频谱分析与衍射分辨率估算角谱法的另一个价值是可以用来验证系统的理论分辨率。对于一条宽度为 D 的入射孔径衍射极限给出的最小可分辨尺寸大约是Δ λ·z / D比如 λ633nm、z1m、D8.8mm算出来 Δ≈72μm约为 4 个像素。这个数字指导意义很大如果仿真里像素尺寸设得太小低于这个值那么你看到的一切细节都是数值噪声不是真实的物理信息。在频域里也可以做同样的分析频谱最高有效频率受限于 λ·z/D 对应的空间频率超过这个频率的频谱分量即使初始存在传播后也会迅速衰减。手动滤波时可以用这个频率做截止点避免引入无意义的数值振荡。6. 调试记录与常见问题速查6.1 输出图案出现环形条纹或边缘振铃这是我在第一次跑距离扫描时最常遇到的现象图案周围出现一圈圈规则环形条纹。原因通常是输入面孔径边缘太锐利也就是图像四边从 1 突变到 0 的硬边在频域引入大量高频分量。解决方法是给输入图添加软边窗比如用余弦边缘平滑一下边缘过渡win tukeywin(N, 0.1) * tukeywin(N, 0.1); U0 U0 .* win;加了软边之后振铃明显减弱。代价是有效视场变小一点但对大多数定性演示影响不大。6.2 结果整体发糊或完全看不出原图如果逆传播恢复出来的图模糊多半是像素尺寸和距离设置不匹配导致频域出现欠采样。另一个常见错误是 fftshift 和 ifftshift 放错位置这会直接让图案做中心翻转人眼看来就是“乱七八糟”。排查方法很简单把传播距离设为零也就是让角谱函数乘一个全 1 的传输函数如果输出不等于输入说明移位或者网格构造有问题。6.3 大距离、高像素仿真内存爆炸512×512 的双精度矩阵运算几个 GB 内存绰绰有余但跑到 4096×4096 就不是回事了。频域乘法中间变量占用多MATLAB 会临时复制好几份矩阵。我的处理方法是尽量用单精度U0 single(U0); lambda single(lambda);单精度带来的误差对强度图显示几乎不可感知内存却能省一半。如果矩阵实在太大还可以把二维 FFT 拆成两次一维 FFT分块处理不过这类优化一般用不上。6.4 常见问题速查表现象可能原因排查与修正图案环形条纹边缘硬边衍射加 tukeywin 软边窗恢复图模糊采样不足或距离过大检查 z ≤ N·dx²/λ图案上下左右颠倒fftshift 与 ifftshift 不配对统一使用推荐写法出现 NaN精确 H 根号内为负用 max(0, t) 保护强度对比较暗振幅忘了开根号U0 sqrt(I0)运行很慢double 矩阵过大改为 single 精度7. 写在最后的一点体会角谱传输这个工具我第一次用的时候只是把它当成一个“FFT 三连”的数学黑箱代码跑通、图出来就以为完事了。直到后来做逆传播恢复才发现那个共轭传输函数背后有大量细节值得琢磨采样条件是不是满足、边缘效应怎么抑制、归一化怎么处理每一步都直接影响结果物理上是否可信。这个项目里最小的感悟反而是最想告诉大家的光学仿真不是去背公式而是要把每一步操作和物理含义对号入座。拿着这套角谱函数你后续可以把它扩展成多平面叠加、涡旋光束传播、大气湍流相位屏模拟甚至光学神经网络的正向模型。从一个简单的菲涅尔肖像图出发能延伸出去的方向远比想象中多。