简介一套基于MATLAB的数字全息仿真代码包面向光学成像、计算机图像处理方向的研究者和学生可用于计算全息图的制作与再现实验。包内共18个文件以6个.m脚本为主体覆盖菲涅尔衍射全息图生成、在线全息重建、并行处理等典型算法同时附有2个jpg、2个bmp、1个tif、1个png共6张示例图像包括经典cameraman测试图以及docx、wps说明文档便于对照代码理解每一环节的输出效果。算法代码从全息记录面的干涉场计算到利用傅里叶变换进行波前重建形成了完整的数字全息处理流程资源中还包含多个变体和备份脚本可观察不同参数设置对再现图像质量的影响适合在修改衍射距离、波长等条件后进一步分析。整个压缩包约515KB体量虽小但信息完整已有690人学习下载通过这份资料不仅能跑通全息图生成到再现的完整示例还能借鉴代码结构设计自己的仿真实验为深入全息三维显示、防伪识别等应用打下基础是入门数字全息编程和验证再现算法的实用参考。1. 数字全息再现zip 文件包里藏的科学比压缩算法更值钱数字全息.zip_cryni1_全息_全息再现_全息图再现这个文件名乍看像随手存的压缩包实际是一个含噪的采集结果cryni1像是采集批次或样品编号全息出现了三次说明这批数据是同一个样品在不同再现条件下的反复处理。数字全息要解决的核心问题不是存储格式而是从一张灰度干涉图全息图里恢复出物光波的振幅和相位进而算出微米甚至纳米级的三维形貌。这个过程的专业说法叫全息再现最常用的算法是菲涅尔衍射积分和角谱法。读这篇文章的人要么手头正握着类似命名的 .zip要么在用显微镜做相衬观察却不知道如何从干涉条纹里挖出定量数据要么在排查为什么明明拍了图却什么都没有重建出来。下面从头到尾讲一遍数字全息图的读取、预处理、数值衍射复现和相位提取代码直接用 Python 写参数全部给到能跑通为止。2. 数字全息图的读入与预处理把干涉条纹变成可计算的复数场2.1 先认清手上有什么参考光、物光和全息图的角色任何一个数字全息记录系统离不开三样东西物光 O、参考光 R 和记录面上的强度分布 H。CCD 或 CMOS 传感器只能记录光强所以干涉条纹被写成H(x, y) |O|² |R|² O·R* O*·R前两项是直流项零级像的源头第三项是物光与共轭参考光的干涉项包含原始物光信息第四项是孪生像。全息再现的目标就是把第三项从 H 里分离出来并让它在数值上重新传播回物平面。解压数字全息.zip之后通常会看到 R、O、H 三张图或直接看到一张带斜条纹的灰度图。cryni1这种后缀一般对应采集时的激光波长或曝光参数标识读数据前先确认像素位深常见 8bit、16bit再确认图像尺寸。位深错了后面的强度归一化和频谱滤波全都白做。常见做法是先用 PIL 或 imageio 读入看最小最大值判断是否需要除以 409512bit或 6553516bit。2.2 减背景与归一化直流项必须尽可能压平直接对原始全息图做再现DC 项会在再现像中心形成一个巨大亮斑把细节全淹没。标准的预处理分三步暗场减除、平场校正、整体归一化。暗场不开激光时采集减掉传感器的固定噪声平场无样品时采集参考光用来校正照明不均匀。import numpy as np from imageio.v2 import imread from scipy.fft import fft2, ifft2, fftshift hologram imread(cryni1_H.tif).astype(np.float64) dark imread(cryni1_dark.tif).astype(np.float64) flat imread(cryni1_flat.tif).astype(np.float64) # 结构先减暗场再除以平场最后做零均值归一化 H (hologram - dark) / (flat - dark 1e-6) H - H.mean() H / H.std()这段代码里的flat - dark对应有效照明分布1e-6是防止除零。归一化到零均值、单位标准差之后DC 项被整体压低后续频谱里三个峰零级、正一级、负一级对比度会明显提升。需要说明的是平场校正不是必需的——如果激光照明本身非常均匀跳过这一步也能出结果但加了之后频谱边缘的环状伪影会少很多。2.3 频谱滤波从干涉图里把物光波单独剥出来全息再现的本质是滤波在空间频谱域把 1 级或 -1 级的频谱切出来移到中心再做一次逆傅里叶变换得到一个复数场。这个复数场就是传感器平面上重建出的物光波分布幅值对应光强相角对应光程差。F fftshift(fft2(H)) rows, cols H.shape mask np.zeros_like(F, dtypenp.float64) cy, cx rows // 2, cols // 2 # 参数说明peak_y、peak_x 是 1 级频谱峰的中心位置需从频谱强度图里目视读出 peak_y, peak_x 178, 223 r 36 # 滤波半径按条纹载频的 2 倍宽度给 yy, xx np.ogrid[:rows, :cols] mask ((yy - peak_y) ** 2 (xx - peak_x) ** 2) r ** 2 F_filtered F * mask complex_field ifft2(fftshift(F_filtered))peak 位置和半径 r 是这段代码里最依赖人工判断的两个参数。半径太小会切掉高频信息形貌边缘变钝半径太大则把零级残留或孪生像的一部分圈进来再现相位图上出现规律性波纹。经验做法是先打印np.abs(fftshift(F))的对数强度图用 matplotlib 直接看峰的位置再用10 * np.log10(np.abs(F).max()) - 20 dB作为阈值自动找峰的范围。需要提醒的是不要试图用一个固定半径跑完所有样品——光学平台的震动状态、样品散射强弱都会让峰的形状变化滤波器半径跟着样品走才稳定。提示频谱峰不在整像素点上时掩模边缘会产生截断伪影。可靠做法是加一个高斯渐变边沿而非硬圆形掩模半径外 5 个像素内让权重从 1 平滑降到 0。3. 数字全息再现的数值算法菲涅尔卷积、角谱法怎么选、怎么算3.1 为什么实际处理很少用直接菲涅尔变换拿到传感器面上的复数场 u(x, y) 之后全息再现要解决的是让光场继续往前传的问题。严格解是瑞利-索末菲衍射积分计算量太大近轴近似下退化为菲涅尔衍射积分U(ξ, η) exp(jkz) / (jλz) · exp[jk(ξ²η²)/2z] · ∫∫ u(x, y) · exp[jk(x²y²)/2z] · exp[-j2π(xξyη) / (λz)] dxdy直接按这个式子离散化输出的像素尺寸会随 z 改变而且为了满足采样率z 被限制在很小的范围内。这导致一个实际问题当你对同一张全息图试 z50mm 和 z55mm 时两幅再现像的横向分辨率都变了没法直接对比。更难受的是当 z 太小比如显微配置下只有几毫米菲涅尔近似的采样条件根本不成立。所以现在主流的全息再现实现里算法库普遍用卷积法或角谱法代替直接菲涅尔变换。3.2 卷积法实现保持像素尺寸不变的核心技巧卷积法守住了这么一条性质把衍射积分写成冲激响应与物场的卷积再利用卷积定理用两次 FFT 完成计算输出平面像素尺寸和输入完全一致。全息再现的卷积形式对应传递函数G(fx, fy) exp[jkz·sqrt(1 - (λfx)² - (λfy)²)]这个式子看着比菲涅尔积分更简洁却保留了像面尺寸不变、任意 z 都能算的两个优势。代码里用scipy.signal.fftconvolve直接算或者手写 FFT 相乘两种方式等价。from scipy.fft import fft2, ifft2, fftshift def angular_spectrum(u0, wavelength, pixel_size, z): 角谱法全息再现像素尺寸不随 z 改变的衍射传播器。 参数 ---- u0 : 2D ndarray传感器平面复振幅频谱滤波后 wavelength : float激光波长单位米 pixel_size : float传感器像素间距单位米 z : float光源/样品到传感器的再现距离单位米 rows, cols u0.shape fx np.fft.fftfreq(cols, dpixel_size) fy np.fft.fftfreq(rows, dpixel_size) FX, FY np.meshgrid(fx, fy) # 频域传递函数大于 1/λ 的部分 evanescent 波直接置 0 term 1 - (wavelength * FX) ** 2 - (wavelength * FY) ** 2 term[term 0] 0 H np.exp(1j * 2 * np.pi * z / wavelength * np.sqrt(term)) U fft2(u0) * H result ifft2(U) return result # 典型参数氮氖激光 632.8nm传感器像素 3.45μm再现距离 85mm reconstructed angular_spectrum(complex_field, wavelength632.8e-9, pixel_size3.45e-6, z0.085)参数选型上有三个坑值得展开。第一wavelength * max(fx)若接近 1则sqrt(term)的斜率很大数值上对 z 极其敏感微小的 z 误差会被放大成明显的离焦模糊所以 z 的给定精度至少到毫米级最好用自动对焦算法去搜。第二网格点fx的排布是[-fs/2, fs/2)零频在正中间而ifft2默认零频在左上角所以注意是否需要ifftshift不匹配会导致输出像发生半周期平移。第三当 z 很小而图像尺寸很大时角谱法的频域采样间隔 δf1/(N·Δ) 可能大于 1/λ 的两倍此时存在混叠请把原始全息图先裁到目标感兴趣区域再传播而不是全幅传播。3.3 再现像面参数的计算放大率、像距和最小分辨尺寸常见实验配置是物光经过显微物镜放大后再与参考光干涉这时再现距离 z 与物镜焦距、筒长、管镜焦距有关。全息再现里最常用的标定公式是M z2 / z1其中 z1 是物镜前焦面到样品的距离近似等于物镜焦距z2 是等效像距。数字全息系统的横向分辨率不取决于物镜数值孔径吗既对也不对NA 决定能收集多少衍射级次但像素尺寸决定了频带能不能完整记录下来。Δx_sample λ·z / (N·Δx_sensor)给出的是无显微放大下的样品面极限分辨尺寸。写代码时建议把上面这个公式保留在注释里因为一旦遇到为什么我重建的微球直径偏大 5%这类问题第一反应应该是去算系统放大率而不是去调滤波器半径。4. 相位提取与解包裹2π 跳变才是全息再现的最后一关4.1 为什么距离信息藏在相位里却读不出来再现得到的复数场result每一项是 a bj相位直接np.angle(result)就行。但 arctan 的输出被截断在 (-π, π] 之间而真实光程差往往大于一个波长导致相位图上一圈一圈的条纹相邻条纹之间差 2π。数字全息再现之后拿到的是包裹相位wrapped phase要变成连续的高度分布必须做相位解包裹phase unwrapping。一维信号可以沿路径逐点累加 2π 偏移二维图像却不行——噪声点会给路径积分引入 2π 整数倍的全局误差这就是二维解包裹算法存在的理由。4.2 质量图引导的路径解包裹抗噪最稳的通用实现最简单的逐行解包裹在条纹密集区比如微球边缘会整行漂移产生胶皮状的伪影。我一般优先用质量图引导的洪水填充算法原理是先计算每个像素的相位导数方差作为可靠性指标从最可靠的像素出发按队列依次展开邻域。def phase_unwrap_quality(phase, reliability_maskNone): 基于质量图(二阶差分)的路径跟踪解包裹。 phase : 2D ndarray包裹相位范围(-π, π] reliability_mask : 可选0~1 的权重图用于屏蔽坏区域 # 相位导数的二阶差分 质量的负指标 dx np.diff(phase, axis1) dy np.diff(phase, axis0) # 处理 2π 跳变后求梯度 dx np.angle(np.exp(1j * dx)) dy np.angle(np.exp(1j * dy)) quality -(np.abs(np.diff(dx, axis0))[1:, :] np.abs(np.diff(dy, axis1))[:, 1:]) # 堆的起点是质量最高的像素逐点弹出并扩展邻域 import heapq rows, cols phase.shape unwrapped np.zeros_like(phase) visited np.zeros_like(phase, dtypebool) heap [] cy, cx np.unravel_index(np.argmax(quality), quality.shape) # 第一行、最后一行、第一列、最后一列的像素质量不可靠跳过做质量图的逻辑里处理 visited[cy, cx] True unwrapped[cy, cx] phase[cy, cx] heapq.heappush(heap, (0, cy, cx)) while heap: _, y, x heapq.heappop(heap) for dy_, dx_ in ((1, 0), (-1, 0), (0, 1), (0, -1)): ny, nx y dy_, x dx_ if 0 ny rows and 0 nx cols and not visited[ny, nx]: delta phase[ny, nx] - phase[y, x] delta - 2 * np.pi * np.round(delta / (2 * np.pi)) unwrapped[ny, nx] unwrapped[y, x] delta visited[ny, nx] True # 质量图的索引比相位图小1要小心对齐 q quality[ny - 1, nx - 1] if ny 0 and nx 0 else -999 heapq.heappush(heap, (-q, ny, nx)) return unwrapped unwrap_phase phase_unwrap_quality(np.angle(reconstructed))这段代码里的delta - 2π * round(delta / 2π)是核心它保证两个像素之间的相位差被调整到最短路径上无论相位图上这一格是 3.1 还是 -3.1相对关系都不会偏。质量图用二阶差分反映局部条纹密度密度越大质量越低扩张就绕过这些区域。参数quality与phase的尺寸差一像素代码里用if ny 0 and nx 0兜底但如果输入图像较大更可靠的做法是把质量图先np.pad到与相位图同尺寸。4.3 从解包裹相位到高度载波相位的扣除方法和标定公式全息再现得到相位差之后还不能直接乘系数。干涉图里的条纹还包含参考光和物光夹角带来的线性载波相位tilt phase这一项不扣掉解出来的高度会叠加一个斜坡面。常见的扣除办法是拿无样品的纯平场全息图做同样的再现与解包裹得到载波相位分布再用样品的解包裹相位减去它。# 假设 flat_unwrapped 是平场全息图用同样流程解包裹后的相位 tilt flat_unwrapped # 线性斜坡但含低频偏移 corrected_phase unwrap_phase - tilt # 高度换算氮氖激光反射式测量光程差为 2h wavelength 632.8e-9 height corrected_phase * wavelength / (4 * np.pi)这里有坑透射式还是反射式换算系数差一倍。透射式测量样品折射率未知时分不清是厚度还是折射率变化反射式最干净公式里除以 4π 因为光走了来回。corrected_phase要做趋势去除再去算高度但趋势去除时不要用多项式拟合整个面——样品占全场 30% 以下时应该只选四角背景区域来拟合平面否则把样品的真实弯曲一起拟合掉了。提示解包裹相位里出现孤立的±2π斑块通常不是物理信号而是传感器坏点或掩模边缘的频谱截断伪影。在质量图里把这些像素权重设 0比解包裹后再做中值滤波更有效。5. 数字全息再现的验证手段与现场检查技巧5.1 三步验证平面镜标定、已知高度台阶、重复性拿到一个全息再现系统第一步不是测量未知样品而是用光学平晶或镀膜硅片当样品记录。理想平面在解包裹后的相位图上应该是平坦的残余起伏就是系统像差。记录平面镜全息图做完整链路滤波、角谱传播、解包裹输出残余相位峰谷值数字全息系统的可用横向分辨率就是峰谷值对应高度的两到三倍。第二步用刻蚀台阶或标准高度块高度已知到纳米级验证高度比例系数这时就会暴露波长校准问题激光器标称 632.8nm 但实际可能是 632.99nm校准系数直接乘到最终高度上。第三步是重复性同一位置连续拍十张解包裹后计算每张相位的标准差。这个值告诉你系统的测量噪声是 0.5nm 还是 5nm也直接决定你论文里误差棒该画多宽。这里有一项值得做把十张全息图先平均再处理和平场平均一样系统随机噪声会按 1/√N 下降。5.2 离焦判据用锐度曲线而不是肉眼自动对焦搜索 z 时很多新手用强度图最锐作为判据但相位图对焦更敏感。我常用的做法是定义相位梯度能量或直接用高频傅里叶能量作为对焦评价函数。具体命令如下def focus_metric(recon): grad np.gradient(np.angle(recon)) return np.sum(grad[0] ** 2 grad[1] ** 2) z_best None score_best -1 for z in np.linspace(0.080, 0.090, 50): r angular_spectrum(complex_field, 632.8e-9, 3.45e-6, z) s focus_metric(r) if s score_best: score_best, z_best s, znp.gradient在 10 像素边缘会放大噪声所以曲线在近焦点附近会变尖扫描步长取 0.2mm 足够。值得补充一点用相位梯度作测度时样品边缘的相位跳变贡献很多评分如果样品本身是阶梯状焦点位置会自动偏向使阶梯最陡的地方——这时要和幅度图的对比度做交叉验证。5.3 快速故障排查表当你什么都重建不出来最后一条经验把最常见的失败现象与对应处置集中在一起方便现场翻。现象原因处置再现像只有亮斑没有细节滤波器中心放在零级没找到 1 级峰检查频谱强度图确认峰位重新滤波相位图出现规则同心环滤波半径包含零级残影缩小掩模半径或中心加圆形阻塞线性斜坡怎么都扣不掉平场全息图和样品全息图采集时刻载波不一致重新采集平场数据确保与样品同批次深沟槽附近出现沿扫描方向拉丝解包裹路径穿过坏像素质量图里坏区权重设 0重新处理z 扫描曲线多峰样品表面有强反射双像或相干噪声先对原始全息图做 3×3 中值滤波再处理磁盘上写完这么多记得把每一步处理后的复数场和数据参数同时保存成.npz而不是只存最终 PNG。因为调整滤波半径或 z 参数的代价是重新跑一遍后端有了中间结果调试周期从两小时缩短到两分钟。这个习惯比任何算法技巧都更能决定全息再现项目能不能按期交付。本文还有配套的精品资源点击获取