★ FEATURED ARTICLE 简介这份资源面向光学、光电信息与图像处理方向的学习者和研究者提供数字全息图重构的MATLAB实现代码帮助理解从全息图中恢复原始光场信息的核心算法。包内共3个m脚本文件压缩包约1KB分别对应基于卷积的重构、基于Fresnel变换的重构以及角谱法三种典型方案覆盖近场、中等距离到远场不同传播条件的重建需求。代码结构简洁便于逐行阅读与对比实验读者可借此掌握卷积核与光学传递函数的关系、Fresnel积分公式的离散实现以及频域位相恢复后逆傅里叶变换的完整流程并可根据重构距离、计算复杂度和精度要求灵活选择方法。目前已有1643人学习下载适合作为数字全息入门实践与光学信号处理课程实验的参考素材。1. 数字全息重构的三种打开方式一份 MATLAB 代码能跑出什么如果你手头有一张全息图想从干涉条纹里把物体复现出来大概率会撞上同一个问题用哪种重构算法卷积法、菲涅耳衍射法、还是角谱法这三种方法在教科书里各占一节但真到写代码的时候边界条件怎么设、补零补多少、重构距离怎么选全是玄学。这份 MATLAB 代码把 CONV卷积法和 Fresnel菲涅耳衍射法两种重构路径都写成了可直接运行的脚本配合数字全息的典型实验参数能让你在半小时内看到重构结果而不是花三天调一个采样率。它适合刚接触数字全息的研究生、需要快速验证光路设计的工程师以及想把这套流程移植到 Python 或 C 的开发者。代码不依赖任何工具箱纯 MATLAB 基础函数实现拿到就能跑。2. 先搞清楚三种重构方法在算什么CONV、Fresnel 与角谱的物理账2.1 从全息图到复振幅重构到底在解什么方程数字全息的核心就一句话用计算机模拟光波从全息面传播到像面的过程。全息图记录的是物光波和参考光波的干涉条纹强度分布可以写成I(x,y) |O|^2 |R|^2 O*R O R*其中 O 是物光复振幅R 是参考光复振幅。前两项是直流分量后两项才是真正携带物体信息的干涉项。重构的第一步永远是滤波——把 1 级或 -1 级频谱单独拎出来去掉直流和共轭项。这一步在代码里通常用傅里叶变换加带通滤波实现滤波窗口的位置和大小直接决定重构质量。很多新手在这里翻车窗口开太大直流分量漏进来重构面上会出现一条亮带窗口开太小高频信息被切掉边缘变模糊。我一般会先把频谱图画出来用imagesc看一眼 1 级和直流分量的间距再手动定窗口参数。滤波之后得到的是全息面上的物光复振幅接下来才是真正的传播计算。传播距离 d 的符号决定是正向重构还是反向重构这个符号搞反了出来的像就是散焦的。代码里通常用d的正负来区分但不同文献的约定不一样我的习惯是物在记录面之前用负距离物在记录面之后用正距离跑之前先用一个已知距离的模拟全息图验证符号。2.2 CONV 卷积法把衍射积分当成卷积来算卷积法的出发点是菲涅耳-基尔霍夫衍射积分。当传播距离满足近场条件时衍射积分可以写成物光复振幅与一个二次相位因子的卷积% CONV 重构核心卷积法 % U0: 全息面复振幅, lambda: 波长, d: 传播距离 % dx: 采样间隔, N: 矩阵尺寸 k 2*pi/lambda; % 构造卷积核二次相位因子 [x, y] meshgrid((-N/2:N/2-1)*dx, (-N/2:N/2-1)*dx); h exp(1i*k/(2*d) * (x.^2 y.^2)); % 频域卷积FFT - 乘传递函数 - IFFT H fft2(fftshift(h)); U0_fft fft2(U0); U_recon ifftshift(ifft2(U0_fft .* H));这段代码的关键在h的构造。exp(1i*k/(2*d)*(x.^2y.^2))是菲涅耳近似下的脉冲响应本质是一个二次相位因子。注意这里用了fftshift把零频移到中心再和U0的频谱相乘。如果不做fftshift卷积核的相位中心就对不上重构结果会整体偏移。参数说明lambda一般取 532nm 或 632.8nm看你的激光器d是重构距离单位米dx是相机像素尺寸常见 3.45μm 或 5.86μmN是矩阵边长通常取 1024 或 2048。d不能太小否则菲涅耳近似失效这时候得换角谱法。卷积法的优点是计算速度快一次 FFT 加一次 IFFT 就完事。缺点是采样间隔在重构面上保持不变如果原始全息图的采样率不够重构面的分辨率就上不去。另外卷积核的尺寸和全息图一样大边缘会出现循环卷积带来的伪影常见做法是补零到 2N×2N 再算算完裁回 N×N。2.3 Fresnel 菲涅耳衍射法单次 FFT 的轻量路径菲涅耳衍射法和卷积法在数学上等价但实现路径不同。它把衍射积分里的二次相位因子拆成两部分一部分留在空域一部分转到频域最终只需要一次 FFT% Fresnel 重构核心单次 FFT 法 % U0: 全息面复振幅, lambda: 波长, d: 传播距离 % dx: 采样间隔, N: 矩阵尺寸 k 2*pi/lambda; % 空域二次相位因子 [x, y] meshgrid((-N/2:N/2-1)*dx, (-N/2:N/2-1)*dx); phase_screen exp(1i*k/(2*d) * (x.^2 y.^2)); % 频域二次相位因子 [fx, fy] meshgrid((-N/2:N/2-1)/(N*dx), (-N/2:N/2-1)/(N*dx)); phase_freq exp(-1i*pi*lambda*d * (fx.^2 fy.^2)); % 单次 FFT 传播 U_recon phase_screen .* ifft2(fft2(U0 .* phase_screen) .* phase_freq);这段代码比卷积法多了一个空域相位因子phase_screen它在 FFT 前后各乘一次。phase_freq是频域的二次相位因子对应菲涅耳衍射的传递函数。最终U_recon就是重构面上的复振幅。参数说明和卷积法一致但要注意Fresnel 法的重构面采样间隔会随距离变化计算公式是dx_recon lambda*d/(N*dx)。这意味着如果d很小重构面的像素间隔会非常小可能出现欠采样如果d很大重构面会变大但分辨率下降。我一般会先算一下dx_recon和原始dx对比如果差太多就换卷积法或角谱法。Fresnel 法的优点是实现简单一次 FFT 搞定适合实时重构。缺点是适用范围比卷积法窄近场和远场都不太准。代码里通常会加一个判断如果d N*dx^2/lambda就用卷积法否则用 Fresnel 法。这个阈值来自菲涅耳近似的成立条件。2.4 角谱法为什么没写进这份代码适用边界与替代方案角谱法Angular Spectrum Method是三种方法里最精确的它不做菲涅耳近似直接对平面波谱做相位延迟。传递函数是H(fx, fy) exp(1i*2*pi*d*sqrt(1/lambda^2 - fx^2 - fy^2))当fx^2 fy^2 1/lambda^2时根号里是负数对应倏逝波通常直接置零。角谱法的采样间隔在重构面上保持不变和卷积法一样但精度更高适合大角度衍射和近场重构。这份代码没写角谱法我猜是因为标题里只提了 CONV 和 Fresnel。但如果你要重构的物体特征尺寸接近波长或者传播距离小于几个波长角谱法是唯一靠谱的选择。补上它也不难把phase_freq换成上面的H就行其他结构不变。我一般会在代码里留一个method参数用switch-case切换三种方法方便对比。3. 把代码跑起来从全息图输入到重构结果输出的完整链路3.1 环境准备与文件结构MATLAB 版本和依赖检查这份代码对 MATLAB 版本要求不高R2016a 以上都能跑因为只用了fft2、ifft2、fftshift、meshgrid这些基础函数。不需要 Image Processing Toolbox也不需要 Signal Processing Toolbox。如果你用的是 Octave大部分函数也兼容但fftshift的行为在旧版本 Octave 里可能不一致建议用 MATLAB 或新版 Octave。文件结构通常是这样的digital_holography_recon/ ├── main.m % 主脚本读图、滤波、重构、显示 ├── recon_conv.m % 卷积法重构函数 ├── recon_fresnel.m % 菲涅耳法重构函数 ├── filter_spectrum.m % 频谱滤波函数 ├── data/ │ └── hologram.mat % 示例全息图数据 └── utils/ ├── load_hologram.m % 数据加载 └── show_result.m % 结果显示如果你的代码包里没有data文件夹可以用模拟全息图代替。生成模拟全息图的代码一般长这样% 生成模拟全息图一个矩形孔径 N 1024; lambda 532e-9; dx 3.45e-6; d 0.1; [x, y] meshgrid((-N/2:N/2-1)*dx, (-N/2:N/2-1)*dx); aperture double(abs(x) 1e-3 abs(y) 1e-3); U_obj aperture .* exp(1i*pi/(lambda*d)*(x.^2 y.^2)); U_ref exp(1i*2*pi/lambda * x * sin(pi/6)); I abs(U_obj U_ref).^2;这段代码生成一个矩形孔径的全息图U_ref是倾斜平面参考光I是干涉强度。跑完main.m后重构面上应该能看到一个清晰的矩形。如果看不到先检查d的符号和lambda的单位。3.2 主脚本逐段拆解读图、滤波、重构、显示主脚本main.m一般分四段。第一段读图和预处理% 第一段读入全息图并归一化 load(data/hologram.mat); % 变量名通常是 I I double(I); I I / max(I(:)); % 归一化到 [0,1]归一化是为了后续滤波时阈值好设。如果全息图是 uint8 格式double转换后范围是 0-255不归一化的话频谱幅度差异太大滤波窗口不好定。第二段是频谱滤波% 第二段傅里叶变换与带通滤波 H fftshift(fft2(I)); figure; imagesc(log(abs(H)1)); colormap gray; title(频谱); % 手动选择 1 级窗口 cx 512; cy 512; % 窗口中心根据频谱图调整 r 80; % 窗口半径 [xx, yy] meshgrid(1:N, 1:N); mask ((xx-cx).^2 (yy-cy).^2) r^2; H_filtered H .* mask; U0 ifft2(ifftshift(H_filtered));cx、cy、r这三个参数需要根据频谱图手动调。我一般会先跑一次imagesc(log(abs(H)1))看 1 级亮斑的位置然后估一个中心坐标和半径。如果亮斑是椭圆可以把mask改成椭圆方程。第三段是重构% 第三段调用重构函数 method conv; % 或 fresnel if strcmp(method, conv) U_recon recon_conv(U0, lambda, d, dx, N); else U_recon recon_fresnel(U0, lambda, d, dx, N); end第四段是显示% 第四段显示强度与相位 figure; subplot(1,2,1); imagesc(abs(U_recon)); colormap gray; title(强度); subplot(1,2,2); imagesc(angle(U_recon)); colormap jet; title(相位);强度图看物体轮廓相位图看厚度变化。如果相位图全是噪声说明滤波没滤干净或者d不对。3.3 参数怎么调波长、距离、像素尺寸的联动关系这三个参数不是独立的它们通过菲涅耳数N_F dx^2/(lambda*d)耦合在一起。N_F决定了重构面上能分辨的最小特征。如果N_F太大重构面采样不足细节丢失如果N_F太小重构面过大计算量浪费。我一般会先固定lambda和dx这两个由实验设备决定。然后根据物体尺寸估一个d算一下N_F看是否在 0.1 到 10 之间。如果不在调整d或换方法。比如lambda532nm、dx3.45μm、d0.1m时N_F (3.45e-6)^2/(532e-9*0.1) ≈ 0.224在合理范围内。另一个容易忽略的参数是补零尺寸。卷积法里如果直接对 N×N 做 FFT循环卷积会污染边缘。常见做法是补零到 2N×2NN_pad 2*N; U0_pad zeros(N_pad, N_pad); U0_pad(1:N, 1:N) U0; % 后续 FFT 用 N_pad算完再裁回中心 N×N。这样边缘伪影会少很多代价是计算量增加四倍。如果N1024N_pad2048现代电脑跑一次也就几秒可以接受。3.4 结果验证怎么判断重构是对的重构结果对不对不能只看图好不好看。我一般用三个判据第一重构面上物体的边缘是否锐利。如果边缘模糊可能是滤波窗口太小或者d不对。可以试着手动微调d看边缘什么时候最锐。第二相位图是否连续。如果相位出现 2π 跳变说明有包裹需要解包裹算法。这份代码没写解包裹但 MATLAB 有unwrap函数可以用。第三用模拟数据验证。生成一个已知物体的全息图重构后和原始物体对比。如果吻合说明参数对了。这一步很多人跳过结果用实验数据调半天调不出来。4. 避坑与排查数字全息重构里最容易翻车的五个地方4.1 频谱滤波窗口选错重构面出现亮带或模糊现象重构强度图上有一条横向或纵向亮带或者物体边缘像蒙了一层雾。原因滤波窗口包含了直流分量或共轭项。直流分量在频谱中心1 级和 -1 级对称分布在两侧。如果窗口中心没对准 1 级或者半径太大就会漏进直流。解决先画频谱图用imagesc(log(abs(H)1))看清楚 1 级的位置。窗口中心取 1 级亮斑的质心半径取亮斑半径的 1.2 倍左右。如果亮斑和直流有重叠说明参考光角度太小需要重新记录全息图。4.2 传播距离符号搞反重构像散焦现象重构面上物体轮廓模糊调d的绝对值也没用。原因d的符号约定和代码里的公式不匹配。有些代码用exp(-1i*k/(2*d))有些用exp(1i*k/(2*d))符号反了就是反向传播。解决用模拟数据验证符号。生成一个已知距离的全息图分别用正负d重构看哪个清晰。记住这个约定以后都按这个来。我一般会在代码注释里写清楚d 0表示物在记录面之后。4.3 补零不足导致循环卷积伪影现象重构面边缘出现物体重影或者物体被截断。原因FFT 默认是循环卷积不补零的话卷积核的尾部会绕回来污染边缘。解决补零到至少 2N×2N。如果物体靠近边缘补到 3N×3N 更保险。补零后记得裁回中心 N×N否则显示范围不对。4.4 采样间隔不匹配重构面出现混叠现象重构面上出现周期性条纹或者物体细节变成摩尔纹。原因Fresnel 法的重构面采样间隔dx_recon lambda*d/(N*dx)和原始dx不匹配导致欠采样。解决算一下dx_recon如果比dx小很多说明重构面被过度采样可以减小N或增大d。如果比dx大很多说明欠采样需要换卷积法或角谱法。我一般会把这个判断写进代码自动选方法。4.5 相位包裹没解厚度信息读不出来现象相位图上有明显的 2π 跳变像梯田一样。原因相位值超过[-π, π]范围被angle函数截断了。解决用unwrap函数解包裹。二维解包裹可以用unwrap(unwrap(angle(U_recon), [], 1), [], 2)但效果一般。更好的做法是用quality-guided或least-squares解包裹算法MATLAB 社区有开源实现。如果只是看趋势不追求定量可以忽略包裹直接看相位梯度。5. 二次开发与移植从 MATLAB 到 Python 的改造思路5.1 把重构函数抽成独立模块方便替换算法这份代码的recon_conv.m和recon_fresnel.m是独立函数这是好的设计。二次开发时我建议再加一个recon_angular.m把角谱法也做成独立函数然后在main.m里用method参数切换。这样对比三种方法时只需要改一行。函数签名统一成function U_recon recon_xxx(U0, lambda, d, dx, N)输入输出一致方便用cell数组批量跑methods {conv, fresnel, angular}; for i 1:length(methods) U_recon{i} feval([recon_ methods{i}], U0, lambda, d, dx, N); end5.2 移植到 PythonNumPy 和 SciPy 的对应关系MATLAB 的fft2、ifft2、fftshift、ifftshift在 NumPy 里都有对应MATLABNumPyfft2(A)np.fft.fft2(A)ifft2(A)np.fft.ifft2(A)fftshift(A)np.fft.fftshift(A)ifftshift(A)np.fft.ifftshift(A)meshgrid(x,y)np.meshgrid(x,y)注意 NumPy 的fft2默认对最后两个轴做变换如果输入是二维数组行为一致。fftshift和ifftshift对偶数尺寸的数组是逆操作对奇数尺寸有细微差别移植时保持尺寸一致就行。Python 版的卷积法重构import numpy as np def recon_conv(U0, wavelength, d, dx, N): k 2 * np.pi / wavelength x np.arange(-N//2, N//2) * dx y np.arange(-N//2, N//2) * dx xx, yy np.meshgrid(x, y) h np.exp(1j * k / (2 * d) * (xx**2 yy**2)) H np.fft.fft2(np.fft.fftshift(h)) U0_fft np.fft.fft2(U0) U_recon np.fft.ifftshift(np.fft.ifft2(U0_fft * H)) return U_recon这段代码和 MATLAB 版逐行对应跑出来的结果应该一致。如果不一样检查np.arange的边界和 MATLAB 的(-N/2:N/2-1)是否对齐。5.3 性能优化用 GPU 加速和预计算传递函数如果N2048或更大CPU 跑一次 FFT 可能要几秒。二次开发时可以考虑用 GPU 加速。MATLAB 有gpuArray把U0和H转到 GPU 上fft2会自动调用 cuFFTU0_gpu gpuArray(U0); H_gpu gpuArray(H); U_recon gather(ifftshift(ifft2(fft2(U0_gpu) .* H_gpu)));Python 里可以用cupy替代numpy接口几乎一样。另一个优化是预计算传递函数。如果lambda、d、dx、N不变H只算一次存下来复用。批量重构时能省不少时间。5.4 从重构到定量相位解包裹和相位解调的衔接重构出来的相位是包裹的要得到定量厚度需要解包裹和相位解调。解包裹可以用skimage.restoration.unwrap_phasePython或unwrapMATLAB。解调是把相位转成高度h lambda * phase_unwrapped / (2 * pi * (n - 1))其中n是物体折射率。这一步需要知道物体材料否则只能得到相对高度。我一般会在重构函数后面加一个postprocess函数把解包裹和解调串起来。这样主脚本只需要调一次输出就是高度图。从那以后我每次跑重构都会先用模拟数据验证一遍参数再上实验数据。这个习惯帮我省了很多来回调参的时间。希望帮到你。本文还有配套的精品资源点击获取
阅读完成 · 觉得有帮助？
