简介一个基于MATLAB的全息成像仿真代码包面向光学工程、信息光学及相关方向的学生与科研人员旨在通过数值模拟揭示菲涅尔衍射、干涉记录与数值再现的完整物理过程。压缩包内包含3个.m源文件整体仅2KB代码结构紧凑覆盖了从光源参数设定、物光波前构建、干涉图样计算到衍射重建的关键步骤便于使用者直接运行和修改参数。该资源目前已吸引385人学习下载适合作为全息成像课程实验或科研预研的参考实现。借助fft2、ifft2、fftshift等频域处理函数以及卷积运算代码能够清晰展示全息图频谱分布如何影响重建像并帮助读者掌握MATLAB在波动光学仿真中的典型用法例如利用meshgrid生成二维光场、通过卷积模拟传播过程等。通过调整波长、距离、孔径等参数可直观观察重建像质量的变化进一步理解菲涅尔近似及全息记录条件的适用边界。1. 全息成像仿真先把菲涅尔衍射变成可运行的数值实验全息成像仿真这个标题十有八九是冲着菲涅尔衍射公式来的要在一个二维复振幅数组上生成全息图再把它数值重建回原始物面整个过程不碰激光器、不碰光学平台全靠 FFT 把公式翻译成 NumPy。这类实验的目的很直接就是把「记录振幅与相位 数值重建」这一闭环跑通。适合谁光学、电子、信号处理背景的学生或者刚转计算成像、需要从原理过渡到可复现代码的工程师。核心参数就四个波长 λ、采样间隔 Δx、衍射距离 z以及参考光角度 θ。把它们定对剩下的事情其实就是几段可复现的矩阵运算。这篇文章按我实际做过的流程来写代码可以直接抄参数和坑都标在明处。2. 菲涅尔衍射数值化先选对传播模型再动手写衍射函数菲涅尔衍射的近似积分长这样U(ξ, η) e^{jkz} / (jλz) ∬ U₀(x, y) · exp[jπ((ξ−x)² (η−y)²) / (λz)] dxdy这个式子的意思是物面每个点发出的球面波在传播距离 z 后相位被二次项 (ξ−x)² (η−y)² 调制。看上去要双重积分但它的数学结构是卷积所以计算上只有两条路可走在空域构造脉冲响应做两次 FFT或者在频域构造传递函数做一次 FFT。工程上绝大多数菲涅尔仿真代码都选第二条路因为它快、稳、好写。2.1 从菲涅尔积分到单次 FFT传递函数法为什么是默认选择菲涅尔衍射的频域传递函数是H(fx, fy) e^{jkz} · exp[−jπλz(fx² fy²)]物场 U₀ 做一次 FFT乘上 H再做一次 IFFT就得到传播 z 之后的复振幅 U。常见做法是把这个函数封装成独立模块后续生成全息图和重建都复用同一份代码。下面这个函数是带限角谱法的菲涅尔近似版我一般把它当默认工具用import numpy as np def fresnel_TF(u0, z, lam632.8e-9, dx5e-6, band_limitTrue): 菲涅尔衍射传递函数法单次 FFT。 u0 : 输入复振幅二维 numpy 数组 z : 传播距离单位 m lam : 波长单位 m dx : 采样间隔单位 m band_limit : 为 True 时对传递函数做带限避免近场混叠 k 2 * np.pi / lam ny, nx u0.shape fx np.fft.fftfreq(nx, ddx) # x 方向频率 fy np.fft.fftfreq(ny, ddx) # y 方向频率采样间隔默认相同 FX, FY np.meshgrid(fx, fy) # 菲涅尔近似的传递函数第一项是整体相位不影响振幅图 H np.exp(1j * k * z) * np.exp(-1j * np.pi * lam * z * (FX**2 FY**2)) if band_limit: # 只保留角谱法里可传播的分量防止近场时高频噪声被放大 fmax 1.0 / (np.sqrt(4 * z**2 / (ny * dx)**2 1) * lam) H * (np.sqrt(FX**2 FY**2) fmax) U np.fft.ifft2(np.fft.fft2(u0) * H) return U频域网格用np.fft.fftfreq生成返回的是从负 Nyquist 到正 Nyquist 的频率值单位是 1/m。meshgrid之后 FX 和 FY 形成二维频域坐标这是 FFT 类仿真里最容易写错的一行很多人忘了频率轴单位导致结果左右翻转。带限条件里那个ny*dx是全息面宽度通频带下限由 z 和面宽共同决定这是带限角谱法BL-ASM的经典做法比我之前用固定半径截断要稳得多。这个函数对近场和远场都适用只有一个代价它在频域做了截断近场传播时会损失部分高频信息。但对于全息图生成和重建这个场景信号带宽远小于截止频率实际影响可以忽略。以我常用的参数为例λ632.8nmΔx5μmN2048z0.3m截止频率在 26300 1/m 左右而物面的有效信息带宽只有几千 1/m留了充足余量。2.2 空域脉冲响应法两次 FFT 的对照帮我抓出过不少 bug不是所有场景都适合传递函数法。当 z 特别小小到接近 Δx 量级时带限角谱的截止频率会降得很低导致结果锐度明显变差。这个时候空域脉冲响应法更直观因为它的核函数长在空域物理意义清楚菲涅尔衍射的脉冲响应是 h(x, y) e^{jkz} / (jλz) · exp[jπ(x² y²) / (λz)]把它和物场做循环卷积即可。我通常用下面这个函数作为第一个函数的交叉验证def fresnel_IR(u0, z, lam632.8e-9, dx5e-6): 菲涅尔衍射空域脉冲响应法两次 FFT。 适用于采样条件 z N*dx^2/lambda 的场景 近场小于该距离时结果会出现环形伪影。 k 2 * np.pi / lam ny, nx u0.shape y (np.arange(ny) - ny // 2) * dx # 空域坐标要对称不能从 0 开始 x (np.arange(nx) - nx // 2) * dx Y, X np.meshgrid(y, x) h np.exp(1j * k * z) / (1j * lam * z) * np.exp(1j * np.pi * (X**2 Y**2) / (lam * z)) Hf np.fft.fft2(np.fft.ifftshift(h)) # 原点移到数组中心后再 FFT U np.fft.ifft2(np.fft.fft2(u0) * Hf) return U * dx * dx # 卷积积分里的微元 dx*dy注意两个细节。第一空域坐标要用np.arange(ny) - ny // 2也就是让原点落在数组正中心而不是从 0 开始否则核函数的相位中心会错位结果是整体平移一个像素且相位全乱。第二循环卷积做完之后要乘一个 dx²这是把离散求和还原成连续积分的缩放因子。这两个错误都不会让图像变黑只是重建像位置偏移或者振幅整体缩放排查起来特别费时间。我用这两个函数做过一致性测试同一个输入物场、同一组参数两种方法输出的振幅差别在 1e-6 量级相位差别在 1e-4 量级。这个一致性给了我很大信心。但要注意空域脉冲响应法有个前置条件z 必须大于等于 N·Δx²/λ否则核函数的高频相位在相邻像素间变化超过 πFFT 循环卷积会把高频折叠成低频形成一圈一圈的环状伪影。第 4 章会专门讲这个坑这里先记下。3. 从物面到全息面离轴干涉、记录与数值重建衍射仿真只是第一步全息成像仿真真正的主菜是把物光与参考光干涉记录强度全息图再数值重建出原始物场。这一节从参考光设计开始一路做到三幅像分离代码跑通后你会看到零级项、实像和孪生像在频域里的分布。3.1 参考光设计离轴角先定上限再算分离度离轴全息的参考光在仿真里是一个平面波R(x, y) exp(j2π(fx·x fy·y))。fx 和 fy 是参考光在频域的位置它和入射角的关系是 fx sinθ/λ。记录下来的全息图强度包含三项物光自干涉项 |U|²、参考光强度 |R|²、交叉干涉项 U·R* 和 U*·R。前两项都落在零频附近第三项和第四项分别出现在 fx 和 -fx 附近。要让重建像和零级光分离参考光频率必须离原点足够远。参考光的频率上限由采样率决定。奈奎斯特条件是 fx 1/(2Δx)工程上我会再留一倍余量取 fx ≤ 1/(4Δx)对应条纹周期至少 4 个像素。下限则由物体带宽决定物体空间尺寸越小频谱越宽参考光频率要大于物体最高频率加上零级项半径一般取几十个频域网格以上。下面是一个参数计算表方便直接套用参数数值说明波长 λ632.8 nm氦氖激光采样间隔 Δx5 μm模拟 CMOS 像元采样点数 N2048单边点数面宽 L10.24 mmN·Δx衍射距离 z0.3 m物面到全息面参考光频率 fx12600 1/m对应偏轴角约 0.008 rad频率网格间隔 Δf97.7 1/m1/L参考光频移129 格fx/Δf在这个配置下零级项占中心约几十格实像频移 129 格孪生像在 258 格三者完全分得开。如果参考光频率太小比如只有 20 格零级项和实像在频域就会重叠重建出来一团亮斑糊住像如果频率太大比如超过 100000 1/m全息图上的条纹周期小于 2 个像素直接混叠成摩尔纹。这两头我在第 4 章都会给出具体翻车现象。3.2 记录与重建三步闭环从物面复振幅到频域滤波取实像物面我用一个组合图形一个圆孔加一个矩形孔放在视场中心偏下一点比纯字母更可控不需要依赖字体渲染。代码如下def make_object(N2048, dx5e-6): 生成测试用物面振幅半径 80 像素的圆孔 40x40 像素方孔。 y (np.arange(N) - N // 2) * dx x (np.arange(N) - N // 2) * dx Y, X np.meshgrid(y, x) obj np.zeros((N, N), dtypenp.complex128) circle (X**2 Y**2) (80 * dx)**2 rect_x (np.abs(X) 40 * dx) (np.abs(Y - 200 * dx) 20 * dx) obj[circle | rect_x] 1.0 return obj物面大小是 10.24mm圆孔直径只有 0.8mm占比很小。这个比例是刻意的物体尺寸小频谱就宽参考光分离窗口的时候不容易和零级项打架。如果你把物体铺满整个视场重建时边缘会和零级光晕混在一起肉眼看着像蒙了一层雾。记录和重建的完整流程如下这是整个仿真最核心的一段lam 632.8e-9 dx 5e-6 N 2048 z 0.3 fx 12600.0 # 参考光 x 方向空间频率1/m fy 0.0 # 1. 物面复振幅 - 菲涅尔传播 - 全息面复振幅 U0 make_object(N, dx) UH fresnel_TF(U0, z, lam, dx, band_limitTrue) # 2. 参考光与干涉记录强度全息图 Yc (np.arange(N) - N // 2).reshape(-1, 1) Xc (np.arange(N) - N // 2).reshape(1, -1) R np.exp(1j * 2 * np.pi * (fx * Xc * dx fy * Yc * dx)) I_H np.abs(UH R)**2 # 记录全息图强度 # 3. 重建全息图乘参考光 - 反向传播 z - 频域窗口提取实像 E I_H * R U_back fresnel_TF(E, -z, lam, dx, band_limitTrue) # 反向传播 F np.fft.fftshift(np.fft.fft2(U_back)) yy, xx np.mgrid[0:N, 0:N] mask ((xx - N//2)**2 (yy - N//2)**2) (40)**2 # 0 频附近半径 40 格 U_clean np.fft.ifft2(np.fft.ifftshift(F * mask))第一步的正向传播用了第 2 章的fresnel_TF把物面复振幅传到全息面。第二步构造参考光时Xc和Yc的单位是像素乘以 dx 后才变成米这个单位换算漏掉的话参考光频率会被放大 5μm 倍全息图条纹密到完全无法识别。第三步最关键E I_H * R之后实像分量 U_H·|R|² 落在零频零级项被搬到 fx孪生像被搬到 2fx。所以频域滤波的窗口应该放在原点而不是 fx 处。这一步我第一次做反了取 fx 附近的窗口结果只滤出零级光的模糊亮斑重建像怎么调都出不来。后来把频谱图打出来才明白实像就在原点是我自己把频移方向搞反了。重建结果可以从三个数判断重建像的振幅最大值对应圆孔中心圆孔边缘锐利方孔四角清晰背景噪声在 1e-4 量级实像质心与物面原点的偏移不超过 1 个像素。这三个条件都满足这个全息成像仿真闭环就算真正跑通了。如果只有肉眼看着像说明你可能只是碰巧调出了一个能看的参数换成别的物体马上翻车。4. 菲涅尔全息仿真避坑采样、距离与角度三座山这一类仿真翻车九成翻在三个地方传播距离 z 不满足采样条件、参考光角度超了奈奎斯特、零级项没滤干净。剩下的翻车大多来自数据类型和内存管理。下面五条是我自己踩过、也在帮别人看代码时反复见到的坑每条按现象、原因、解决的顺序写清楚。4.1 距离 z 太小重建像四周出现一圈圈环状波纹现象把 z 从 0.3m 改成 0.02m重建像主体还在但周围多出一圈一圈同心圆环像水波纹一样而且怎么调整滤波半径都消不掉。原因菲涅尔空域脉冲响应在 z 很小时相位变化极快。采样条件要求 z ≥ N·Δx²/λ代入 N2048、Δx5μm、λ632.8nmz_min 0.081m。z0.02m 远小于这个值脉冲响应相邻像素的相位差超过 π循环卷积把本该折叠的高频分量映射成了低频于是出现环状条纹。这是典型的频谱混叠不是算法写错。解决先算 z_min 再定距离。如果因为实验场景限制必须用小 z就把采样点数 N 降下来或者改用带限角谱法。带限角谱法的band_limitTrue能缓解一部分但缓解不了太多它只是把超高频通道关掉信息本身已经丢了。最可靠的做法是保持 z ≥ N·Δx²/λ宁可让物体在视场里小一点。4.2 参考光角度调大后全息图出现斜向摩尔纹现象为了把三个像拉得更开把参考光频率从 12600 1/m 调到 90000 1/m。全息图看起来像有条纹但重建像变成一格一格的斜纹完全看不出物体形状。原因90000 1/m 对应的条纹周期约 11μm而采样间隔是 5μm每个周期不到 2.3 个像素已经逼近奈奎斯特极限。频域里实像和孪生像的频谱发生重叠混叠重建出来自然是一团乱纹。90000 1/m 对应的 θ arcsin(fx·λ) arcsin(0.057) ≈ 0.057 rad而极限角度是 λ/(2Δx) 0.063 rad已经贴着边了。解决参考光频率不要超过奈奎斯特的一半即 fx ≤ 1/(4Δx)。对应本组参数是 50000 1/m我一般取 1000020000 1/m 这个区间既保证分离度又留足采样余量。按这个区间算出来的条纹周期都大于 30μm约 6 个像素以上重建质量稳定得多。4.3 重建像中心一团亮斑物体像被糊在光晕里现象物体重建出来了边缘也清晰但中心有一个很亮的圆斑物体像正好压在上面对比度差到几乎看不清。原因这是零级项没滤干净。零级项包含 |U|² 和 |R|² 两部分它们在全息图里能量最大重建后集中在原点附近尺寸比物体像大得多。如果物体本身放在视场中心零级光斑就会直接盖住实像。解决两步走。第一滤波器半径要大于物体带宽但小于参考光频移量。以本组参数为例物体半宽约 80 像素对应频谱半径约二三十格滤波半径取 40 格就够同时零级频移到 129 格完全不会被窗口圈进来。第二如果物体必须放中心、零级光又实在太强可以在记录时做三步相移把 |U|² 和 |R|² 消掉只剩干涉项。常见做法是采集相位差 120° 的三张全息图加减组合后得到纯干涉场仿真里只需要把参考光相位写成参数循环三次代价是时间翻三倍但效果立竿见影。4.4 全息图存成 PNG 再读回来重建全是噪点现象仿真跑完把强度全息图I_H用plt.imsave()存成 PNG下次直接读图重建结果背景噪声大像被砂纸磨过一样。原因全息图的干涉条纹动态范围很大暗条纹和亮条纹差几个数量级而 PNG 只有 8bit 量化256 个灰度级根本装不下完整的条纹细节。量化误差在重建时会被反向传播的高通特性放大变成高频噪点。更隐蔽的问题是很多同学保存前做了归一化把最大值压到 255暗条纹直接变成 0相当于把干涉信息截断了。解决仿真链路里不要存图全程用 numpy 的复数数组传递数据必须持久化时用np.savez_compressed()保存浮点复数场或者保存I_H.astype(np.float32)读取后用双线性插值恢复尺寸再重建。如果只是要一张示意图存 PNG 无所谓但别拿它当重建输入。这是我交过学费的一条后来所有全息数据都改成 npz 存档再没出过这类问题。4.5 N 从 1024 改到 4096 后内存溢出结果还不可复现现象为了看得更清楚把 N 改成 4096程序报 MemoryError。调回 1024 后重建结果每次都有一点不同特别是背景噪声。原因N4096 时单个复数数组占 4096²×16 字节 268MBfresnel_TF里中间变量至少有三个这样的数组加上全息图和重建场峰值内存轻松超过 1.5GB。至于每次结果不同是因为全息图强度、重建场里有随机初始化的部分或者用了某些依赖随机数的操作没有固定随机种子。FFT 本身是确定的但程序里只要有一个np.random调用背景噪声就会变。解决N 一般取 2048 就足够教学演示再大就分块计算或者用np.empty预先分配内存减少临时数组。同时固定np.random.seed(0)并把所有中间结果用np.allclose()做回归测试。我现在的习惯是每次跑完把重建像的峰值位置和背景均值打出来如果两次运行这两个数不一致一定是有随机源混进来了优先排查而不是继续调参。5. 像质验证与调参顺序让重建结果可量化、可复现重建像看着清晰并不等于正确。我会用三个量化手段验证结果任何一个不达标都会优先怀疑参数而不是算法。第一个是闭合性测试把物面场正向传播 z再反向传播 z理论上应该回到物面本身。用归一化均方误差衡量U0 make_object(N, dx) U_forward fresnel_TF(U0, z, lam, dx, band_limitTrue) U_roundtrip fresnel_TF(U_forward, -z, lam, dx, band_limitTrue) err np.linalg.norm(np.abs(U_roundtrip) - np.abs(U0)) / np.linalg.norm(np.abs(U0)) print(fround-trip NRMSE {err:.2e})闭合性误差在 1e-5 量级说明传播函数本身没问题误差大于 1e-2 就要检查距离和带限条件。第二个是质心偏移校验离轴全息中实像中心相对零级中心的偏移量理论值是 z·tanθ用重建像振幅做质心加权实测偏移和理论值的差应该小于一个像素。这个校验能一次性抓出参考光频率单位写错、频移方向写反、滤波窗口选错三类 bug。第三个是参数扫描最具参考价值。把 z 从 0.02m 扫到 0.6m对每个 z 做完整记录和重建计算重建像与原物的相关系数你会看到结果分三段z 小于 0.08m 时相关系数骤降这是采样条件失效0.1m 到 0.5m 之间是平缓平台超过 0.6m 后物体在频域衰减过大相关系数缓慢下降。这个曲线让我彻底明白菲涅尔仿真的 z 不是随便给一个就行它有一个由 Δx、N、λ 共同决定的「甜蜜区间」。我现在的习惯是把 λ、Δx、N、z、fx 五个参数写成常量块放在脚本最顶部每次只动一个变量跑完先看闭合性误差再看质心偏移最后才看图像。遇到重建糊了先回归测试确认代码没改坏再依次排查采样条件、参考光角度和零级滤除。这个顺序帮我少走了很多弯路也把「调参数看运气」变成了可解释的流程。这套方法不只是应付课设放到真实的计算全息、数字全息显微镜项目里排查思路完全一样。希望帮到你。本文还有配套的精品资源点击获取
菲涅尔全息成像仿真：从FFT原理到数值重建完整指南
简介一个基于MATLAB的全息成像仿真代码包面向光学工程、信息光学及相关方向的学生与科研人员旨在通过数值模拟揭示菲涅尔衍射、干涉记录与数值再现的完整物理过程。压缩包内包含3个.m源文件整体仅2KB代码结构紧凑覆盖了从光源参数设定、物光波前构建、干涉图样计算到衍射重建的关键步骤便于使用者直接运行和修改参数。该资源目前已吸引385人学习下载适合作为全息成像课程实验或科研预研的参考实现。借助fft2、ifft2、fftshift等频域处理函数以及卷积运算代码能够清晰展示全息图频谱分布如何影响重建像并帮助读者掌握MATLAB在波动光学仿真中的典型用法例如利用meshgrid生成二维光场、通过卷积模拟传播过程等。通过调整波长、距离、孔径等参数可直观观察重建像质量的变化进一步理解菲涅尔近似及全息记录条件的适用边界。1. 全息成像仿真先把菲涅尔衍射变成可运行的数值实验全息成像仿真这个标题十有八九是冲着菲涅尔衍射公式来的要在一个二维复振幅数组上生成全息图再把它数值重建回原始物面整个过程不碰激光器、不碰光学平台全靠 FFT 把公式翻译成 NumPy。这类实验的目的很直接就是把「记录振幅与相位 数值重建」这一闭环跑通。适合谁光学、电子、信号处理背景的学生或者刚转计算成像、需要从原理过渡到可复现代码的工程师。核心参数就四个波长 λ、采样间隔 Δx、衍射距离 z以及参考光角度 θ。把它们定对剩下的事情其实就是几段可复现的矩阵运算。这篇文章按我实际做过的流程来写代码可以直接抄参数和坑都标在明处。2. 菲涅尔衍射数值化先选对传播模型再动手写衍射函数菲涅尔衍射的近似积分长这样U(ξ, η) e^{jkz} / (jλz) ∬ U₀(x, y) · exp[jπ((ξ−x)² (η−y)²) / (λz)] dxdy这个式子的意思是物面每个点发出的球面波在传播距离 z 后相位被二次项 (ξ−x)² (η−y)² 调制。看上去要双重积分但它的数学结构是卷积所以计算上只有两条路可走在空域构造脉冲响应做两次 FFT或者在频域构造传递函数做一次 FFT。工程上绝大多数菲涅尔仿真代码都选第二条路因为它快、稳、好写。2.1 从菲涅尔积分到单次 FFT传递函数法为什么是默认选择菲涅尔衍射的频域传递函数是H(fx, fy) e^{jkz} · exp[−jπλz(fx² fy²)]物场 U₀ 做一次 FFT乘上 H再做一次 IFFT就得到传播 z 之后的复振幅 U。常见做法是把这个函数封装成独立模块后续生成全息图和重建都复用同一份代码。下面这个函数是带限角谱法的菲涅尔近似版我一般把它当默认工具用import numpy as np def fresnel_TF(u0, z, lam632.8e-9, dx5e-6, band_limitTrue): 菲涅尔衍射传递函数法单次 FFT。 u0 : 输入复振幅二维 numpy 数组 z : 传播距离单位 m lam : 波长单位 m dx : 采样间隔单位 m band_limit : 为 True 时对传递函数做带限避免近场混叠 k 2 * np.pi / lam ny, nx u0.shape fx np.fft.fftfreq(nx, ddx) # x 方向频率 fy np.fft.fftfreq(ny, ddx) # y 方向频率采样间隔默认相同 FX, FY np.meshgrid(fx, fy) # 菲涅尔近似的传递函数第一项是整体相位不影响振幅图 H np.exp(1j * k * z) * np.exp(-1j * np.pi * lam * z * (FX**2 FY**2)) if band_limit: # 只保留角谱法里可传播的分量防止近场时高频噪声被放大 fmax 1.0 / (np.sqrt(4 * z**2 / (ny * dx)**2 1) * lam) H * (np.sqrt(FX**2 FY**2) fmax) U np.fft.ifft2(np.fft.fft2(u0) * H) return U频域网格用np.fft.fftfreq生成返回的是从负 Nyquist 到正 Nyquist 的频率值单位是 1/m。meshgrid之后 FX 和 FY 形成二维频域坐标这是 FFT 类仿真里最容易写错的一行很多人忘了频率轴单位导致结果左右翻转。带限条件里那个ny*dx是全息面宽度通频带下限由 z 和面宽共同决定这是带限角谱法BL-ASM的经典做法比我之前用固定半径截断要稳得多。这个函数对近场和远场都适用只有一个代价它在频域做了截断近场传播时会损失部分高频信息。但对于全息图生成和重建这个场景信号带宽远小于截止频率实际影响可以忽略。以我常用的参数为例λ632.8nmΔx5μmN2048z0.3m截止频率在 26300 1/m 左右而物面的有效信息带宽只有几千 1/m留了充足余量。2.2 空域脉冲响应法两次 FFT 的对照帮我抓出过不少 bug不是所有场景都适合传递函数法。当 z 特别小小到接近 Δx 量级时带限角谱的截止频率会降得很低导致结果锐度明显变差。这个时候空域脉冲响应法更直观因为它的核函数长在空域物理意义清楚菲涅尔衍射的脉冲响应是 h(x, y) e^{jkz} / (jλz) · exp[jπ(x² y²) / (λz)]把它和物场做循环卷积即可。我通常用下面这个函数作为第一个函数的交叉验证def fresnel_IR(u0, z, lam632.8e-9, dx5e-6): 菲涅尔衍射空域脉冲响应法两次 FFT。 适用于采样条件 z N*dx^2/lambda 的场景 近场小于该距离时结果会出现环形伪影。 k 2 * np.pi / lam ny, nx u0.shape y (np.arange(ny) - ny // 2) * dx # 空域坐标要对称不能从 0 开始 x (np.arange(nx) - nx // 2) * dx Y, X np.meshgrid(y, x) h np.exp(1j * k * z) / (1j * lam * z) * np.exp(1j * np.pi * (X**2 Y**2) / (lam * z)) Hf np.fft.fft2(np.fft.ifftshift(h)) # 原点移到数组中心后再 FFT U np.fft.ifft2(np.fft.fft2(u0) * Hf) return U * dx * dx # 卷积积分里的微元 dx*dy注意两个细节。第一空域坐标要用np.arange(ny) - ny // 2也就是让原点落在数组正中心而不是从 0 开始否则核函数的相位中心会错位结果是整体平移一个像素且相位全乱。第二循环卷积做完之后要乘一个 dx²这是把离散求和还原成连续积分的缩放因子。这两个错误都不会让图像变黑只是重建像位置偏移或者振幅整体缩放排查起来特别费时间。我用这两个函数做过一致性测试同一个输入物场、同一组参数两种方法输出的振幅差别在 1e-6 量级相位差别在 1e-4 量级。这个一致性给了我很大信心。但要注意空域脉冲响应法有个前置条件z 必须大于等于 N·Δx²/λ否则核函数的高频相位在相邻像素间变化超过 πFFT 循环卷积会把高频折叠成低频形成一圈一圈的环状伪影。第 4 章会专门讲这个坑这里先记下。3. 从物面到全息面离轴干涉、记录与数值重建衍射仿真只是第一步全息成像仿真真正的主菜是把物光与参考光干涉记录强度全息图再数值重建出原始物场。这一节从参考光设计开始一路做到三幅像分离代码跑通后你会看到零级项、实像和孪生像在频域里的分布。3.1 参考光设计离轴角先定上限再算分离度离轴全息的参考光在仿真里是一个平面波R(x, y) exp(j2π(fx·x fy·y))。fx 和 fy 是参考光在频域的位置它和入射角的关系是 fx sinθ/λ。记录下来的全息图强度包含三项物光自干涉项 |U|²、参考光强度 |R|²、交叉干涉项 U·R* 和 U*·R。前两项都落在零频附近第三项和第四项分别出现在 fx 和 -fx 附近。要让重建像和零级光分离参考光频率必须离原点足够远。参考光的频率上限由采样率决定。奈奎斯特条件是 fx 1/(2Δx)工程上我会再留一倍余量取 fx ≤ 1/(4Δx)对应条纹周期至少 4 个像素。下限则由物体带宽决定物体空间尺寸越小频谱越宽参考光频率要大于物体最高频率加上零级项半径一般取几十个频域网格以上。下面是一个参数计算表方便直接套用参数数值说明波长 λ632.8 nm氦氖激光采样间隔 Δx5 μm模拟 CMOS 像元采样点数 N2048单边点数面宽 L10.24 mmN·Δx衍射距离 z0.3 m物面到全息面参考光频率 fx12600 1/m对应偏轴角约 0.008 rad频率网格间隔 Δf97.7 1/m1/L参考光频移129 格fx/Δf在这个配置下零级项占中心约几十格实像频移 129 格孪生像在 258 格三者完全分得开。如果参考光频率太小比如只有 20 格零级项和实像在频域就会重叠重建出来一团亮斑糊住像如果频率太大比如超过 100000 1/m全息图上的条纹周期小于 2 个像素直接混叠成摩尔纹。这两头我在第 4 章都会给出具体翻车现象。3.2 记录与重建三步闭环从物面复振幅到频域滤波取实像物面我用一个组合图形一个圆孔加一个矩形孔放在视场中心偏下一点比纯字母更可控不需要依赖字体渲染。代码如下def make_object(N2048, dx5e-6): 生成测试用物面振幅半径 80 像素的圆孔 40x40 像素方孔。 y (np.arange(N) - N // 2) * dx x (np.arange(N) - N // 2) * dx Y, X np.meshgrid(y, x) obj np.zeros((N, N), dtypenp.complex128) circle (X**2 Y**2) (80 * dx)**2 rect_x (np.abs(X) 40 * dx) (np.abs(Y - 200 * dx) 20 * dx) obj[circle | rect_x] 1.0 return obj物面大小是 10.24mm圆孔直径只有 0.8mm占比很小。这个比例是刻意的物体尺寸小频谱就宽参考光分离窗口的时候不容易和零级项打架。如果你把物体铺满整个视场重建时边缘会和零级光晕混在一起肉眼看着像蒙了一层雾。记录和重建的完整流程如下这是整个仿真最核心的一段lam 632.8e-9 dx 5e-6 N 2048 z 0.3 fx 12600.0 # 参考光 x 方向空间频率1/m fy 0.0 # 1. 物面复振幅 - 菲涅尔传播 - 全息面复振幅 U0 make_object(N, dx) UH fresnel_TF(U0, z, lam, dx, band_limitTrue) # 2. 参考光与干涉记录强度全息图 Yc (np.arange(N) - N // 2).reshape(-1, 1) Xc (np.arange(N) - N // 2).reshape(1, -1) R np.exp(1j * 2 * np.pi * (fx * Xc * dx fy * Yc * dx)) I_H np.abs(UH R)**2 # 记录全息图强度 # 3. 重建全息图乘参考光 - 反向传播 z - 频域窗口提取实像 E I_H * R U_back fresnel_TF(E, -z, lam, dx, band_limitTrue) # 反向传播 F np.fft.fftshift(np.fft.fft2(U_back)) yy, xx np.mgrid[0:N, 0:N] mask ((xx - N//2)**2 (yy - N//2)**2) (40)**2 # 0 频附近半径 40 格 U_clean np.fft.ifft2(np.fft.ifftshift(F * mask))第一步的正向传播用了第 2 章的fresnel_TF把物面复振幅传到全息面。第二步构造参考光时Xc和Yc的单位是像素乘以 dx 后才变成米这个单位换算漏掉的话参考光频率会被放大 5μm 倍全息图条纹密到完全无法识别。第三步最关键E I_H * R之后实像分量 U_H·|R|² 落在零频零级项被搬到 fx孪生像被搬到 2fx。所以频域滤波的窗口应该放在原点而不是 fx 处。这一步我第一次做反了取 fx 附近的窗口结果只滤出零级光的模糊亮斑重建像怎么调都出不来。后来把频谱图打出来才明白实像就在原点是我自己把频移方向搞反了。重建结果可以从三个数判断重建像的振幅最大值对应圆孔中心圆孔边缘锐利方孔四角清晰背景噪声在 1e-4 量级实像质心与物面原点的偏移不超过 1 个像素。这三个条件都满足这个全息成像仿真闭环就算真正跑通了。如果只有肉眼看着像说明你可能只是碰巧调出了一个能看的参数换成别的物体马上翻车。4. 菲涅尔全息仿真避坑采样、距离与角度三座山这一类仿真翻车九成翻在三个地方传播距离 z 不满足采样条件、参考光角度超了奈奎斯特、零级项没滤干净。剩下的翻车大多来自数据类型和内存管理。下面五条是我自己踩过、也在帮别人看代码时反复见到的坑每条按现象、原因、解决的顺序写清楚。4.1 距离 z 太小重建像四周出现一圈圈环状波纹现象把 z 从 0.3m 改成 0.02m重建像主体还在但周围多出一圈一圈同心圆环像水波纹一样而且怎么调整滤波半径都消不掉。原因菲涅尔空域脉冲响应在 z 很小时相位变化极快。采样条件要求 z ≥ N·Δx²/λ代入 N2048、Δx5μm、λ632.8nmz_min 0.081m。z0.02m 远小于这个值脉冲响应相邻像素的相位差超过 π循环卷积把本该折叠的高频分量映射成了低频于是出现环状条纹。这是典型的频谱混叠不是算法写错。解决先算 z_min 再定距离。如果因为实验场景限制必须用小 z就把采样点数 N 降下来或者改用带限角谱法。带限角谱法的band_limitTrue能缓解一部分但缓解不了太多它只是把超高频通道关掉信息本身已经丢了。最可靠的做法是保持 z ≥ N·Δx²/λ宁可让物体在视场里小一点。4.2 参考光角度调大后全息图出现斜向摩尔纹现象为了把三个像拉得更开把参考光频率从 12600 1/m 调到 90000 1/m。全息图看起来像有条纹但重建像变成一格一格的斜纹完全看不出物体形状。原因90000 1/m 对应的条纹周期约 11μm而采样间隔是 5μm每个周期不到 2.3 个像素已经逼近奈奎斯特极限。频域里实像和孪生像的频谱发生重叠混叠重建出来自然是一团乱纹。90000 1/m 对应的 θ arcsin(fx·λ) arcsin(0.057) ≈ 0.057 rad而极限角度是 λ/(2Δx) 0.063 rad已经贴着边了。解决参考光频率不要超过奈奎斯特的一半即 fx ≤ 1/(4Δx)。对应本组参数是 50000 1/m我一般取 1000020000 1/m 这个区间既保证分离度又留足采样余量。按这个区间算出来的条纹周期都大于 30μm约 6 个像素以上重建质量稳定得多。4.3 重建像中心一团亮斑物体像被糊在光晕里现象物体重建出来了边缘也清晰但中心有一个很亮的圆斑物体像正好压在上面对比度差到几乎看不清。原因这是零级项没滤干净。零级项包含 |U|² 和 |R|² 两部分它们在全息图里能量最大重建后集中在原点附近尺寸比物体像大得多。如果物体本身放在视场中心零级光斑就会直接盖住实像。解决两步走。第一滤波器半径要大于物体带宽但小于参考光频移量。以本组参数为例物体半宽约 80 像素对应频谱半径约二三十格滤波半径取 40 格就够同时零级频移到 129 格完全不会被窗口圈进来。第二如果物体必须放中心、零级光又实在太强可以在记录时做三步相移把 |U|² 和 |R|² 消掉只剩干涉项。常见做法是采集相位差 120° 的三张全息图加减组合后得到纯干涉场仿真里只需要把参考光相位写成参数循环三次代价是时间翻三倍但效果立竿见影。4.4 全息图存成 PNG 再读回来重建全是噪点现象仿真跑完把强度全息图I_H用plt.imsave()存成 PNG下次直接读图重建结果背景噪声大像被砂纸磨过一样。原因全息图的干涉条纹动态范围很大暗条纹和亮条纹差几个数量级而 PNG 只有 8bit 量化256 个灰度级根本装不下完整的条纹细节。量化误差在重建时会被反向传播的高通特性放大变成高频噪点。更隐蔽的问题是很多同学保存前做了归一化把最大值压到 255暗条纹直接变成 0相当于把干涉信息截断了。解决仿真链路里不要存图全程用 numpy 的复数数组传递数据必须持久化时用np.savez_compressed()保存浮点复数场或者保存I_H.astype(np.float32)读取后用双线性插值恢复尺寸再重建。如果只是要一张示意图存 PNG 无所谓但别拿它当重建输入。这是我交过学费的一条后来所有全息数据都改成 npz 存档再没出过这类问题。4.5 N 从 1024 改到 4096 后内存溢出结果还不可复现现象为了看得更清楚把 N 改成 4096程序报 MemoryError。调回 1024 后重建结果每次都有一点不同特别是背景噪声。原因N4096 时单个复数数组占 4096²×16 字节 268MBfresnel_TF里中间变量至少有三个这样的数组加上全息图和重建场峰值内存轻松超过 1.5GB。至于每次结果不同是因为全息图强度、重建场里有随机初始化的部分或者用了某些依赖随机数的操作没有固定随机种子。FFT 本身是确定的但程序里只要有一个np.random调用背景噪声就会变。解决N 一般取 2048 就足够教学演示再大就分块计算或者用np.empty预先分配内存减少临时数组。同时固定np.random.seed(0)并把所有中间结果用np.allclose()做回归测试。我现在的习惯是每次跑完把重建像的峰值位置和背景均值打出来如果两次运行这两个数不一致一定是有随机源混进来了优先排查而不是继续调参。5. 像质验证与调参顺序让重建结果可量化、可复现重建像看着清晰并不等于正确。我会用三个量化手段验证结果任何一个不达标都会优先怀疑参数而不是算法。第一个是闭合性测试把物面场正向传播 z再反向传播 z理论上应该回到物面本身。用归一化均方误差衡量U0 make_object(N, dx) U_forward fresnel_TF(U0, z, lam, dx, band_limitTrue) U_roundtrip fresnel_TF(U_forward, -z, lam, dx, band_limitTrue) err np.linalg.norm(np.abs(U_roundtrip) - np.abs(U0)) / np.linalg.norm(np.abs(U0)) print(fround-trip NRMSE {err:.2e})闭合性误差在 1e-5 量级说明传播函数本身没问题误差大于 1e-2 就要检查距离和带限条件。第二个是质心偏移校验离轴全息中实像中心相对零级中心的偏移量理论值是 z·tanθ用重建像振幅做质心加权实测偏移和理论值的差应该小于一个像素。这个校验能一次性抓出参考光频率单位写错、频移方向写反、滤波窗口选错三类 bug。第三个是参数扫描最具参考价值。把 z 从 0.02m 扫到 0.6m对每个 z 做完整记录和重建计算重建像与原物的相关系数你会看到结果分三段z 小于 0.08m 时相关系数骤降这是采样条件失效0.1m 到 0.5m 之间是平缓平台超过 0.6m 后物体在频域衰减过大相关系数缓慢下降。这个曲线让我彻底明白菲涅尔仿真的 z 不是随便给一个就行它有一个由 Δx、N、λ 共同决定的「甜蜜区间」。我现在的习惯是把 λ、Δx、N、z、fx 五个参数写成常量块放在脚本最顶部每次只动一个变量跑完先看闭合性误差再看质心偏移最后才看图像。遇到重建糊了先回归测试确认代码没改坏再依次排查采样条件、参考光角度和零级滤除。这个顺序帮我少走了很多弯路也把「调参数看运气」变成了可解释的流程。这套方法不只是应付课设放到真实的计算全息、数字全息显微镜项目里排查思路完全一样。希望帮到你。本文还有配套的精品资源点击获取
