简介这份资源面向光学、光电信息与图像处理方向的学习者和研究者提供数字全息图重构的MATLAB实现代码帮助理解从全息图恢复原始光场信息的核心流程。包内共3个m脚本文件压缩包约1KB分别对应基于卷积CONV、Fresnel变换Fresnel_diff与角谱法Angular_spectrum三种重构思路卷积法借助光学传递函数做近场重建计算相对简单Fresnel变换法考虑光源到探测器的传播距离适合中等距离重构角谱法基于傅里叶光学在频域完成位相恢复与逆变换理想条件下可实现无失真重建。三种方法各有适用距离、计算复杂度与精度取舍读者可对照运行脚本比较不同重构结果进而掌握数字全息的基本理论与光学信号处理的实践方法。目前已有1643人学习下载适合作为课程实验、课题入门或算法对比的参考素材。1. 从一张全息图到三种重建结果这个 MATLAB 项目到底在做什么如果你手头有一张数字全息图想把它还原成物体的振幅和相位最直接的问题就是用哪种重构算法卷积法CONV和菲涅耳衍射法Fresnel到底差在哪这个 MATLAB 项目就是冲着这个问题来的——它把数字全息里最常用的三种重构路径打包成可运行的脚本输入一张全息图输出重建后的复振幅分布。适合正在做光学测量、显微成像或者干涉计量方向的研究生和工程师尤其是那些不想从零推导衍射公式、只想先把结果跑出来看看的人。项目本身不依赖任何工具箱之外的第三方库纯 MATLAB 实现拿到手改几个路径参数就能出图。2. 三种重构方法的代码结构CONV、Fresnel 和角谱怎么选2.1 目录组织与入口脚本这个项目的文件组织方式很直白没有搞复杂的包管理。根目录下通常能看到这样几类文件主入口脚本比如main.m或run_reconstruction.m、三个重构函数文件recon_conv.m、recon_fresnel.m、recon_angular.m、以及一个放全息图的data/文件夹。主脚本负责读图、设参数、调函数、画结果重构函数各自独立互不干扰。我一般会先看主脚本里的参数区因为那里决定了整个重建的物理尺度。典型的参数包括波长lambda、像素尺寸pixel_size、重建距离z、以及全息图的尺寸。这些参数不是随便填的它们直接决定了重建平面上的采样间隔和视场大小。% main.m 参数设置区常见写法 lambda 532e-9; % 波长单位米 pixel_size 3.45e-6; % 相机像素尺寸单位米 z 0.15; % 重建距离单位米 hologram imread(data/holo.bmp); hologram double(hologram(:,:,1)); % 取单通道这里lambda和pixel_size决定了衍射计算的尺度z是物平面到全息面的距离。如果z设错了重建出来的像会模糊或者尺寸不对。常见做法是先用一个已知的校准样本跑一遍确认z的符号和大小。2.2 CONV 重构卷积法的实现细节卷积法的核心思想是把菲涅耳衍射积分写成卷积形式然后用 FFT 加速。在离散域里它等价于对全息图乘以一个二次相位因子再做一次卷积。代码里通常这样写function U recon_conv(hologram, lambda, pixel_size, z) [M, N] size(hologram); k 2*pi/lambda; % 空间坐标 x (-N/2:N/2-1) * pixel_size; y (-M/2:M/2-1) * pixel_size; [X, Y] meshgrid(x, y); % 卷积核二次相位因子 h exp(1j*k/(2*z) * (X.^2 Y.^2)); % 频域卷积 H fft2(fftshift(h)); U ifft2(fft2(hologram) .* H); U fftshift(U); end这段代码里h是菲涅耳衍射的脉冲响应fftshift把零频移到中心。卷积法的一个特点是重建平面的采样间隔和全息面一致都是pixel_size。这意味着如果你想要更大的视场卷积法不会自动帮你放大它保持原尺寸。参数z在这里出现在分母上z越小二次相位因子的曲率越大卷积核越“尖锐”。注意卷积法在z很小时会出现混叠因为脉冲响应的相位变化太快离散采样跟不上。这时候要么增大采样率要么换方法。2.3 Fresnel 重构单次 FFT 的轻量路径菲涅耳衍射法也叫单傅里叶变换法是另一种常见实现。它把衍射积分近似成一次傅里叶变换重建平面的采样间隔会随z变化。代码通常长这样function U recon_fresnel(hologram, lambda, pixel_size, z) [M, N] size(hologram); k 2*pi/lambda; x (-N/2:N/2-1) * pixel_size; y (-M/2:M/2-1) * pixel_size; [X, Y] meshgrid(x, y); % 二次相位因子 phase exp(1j*k/(2*z) * (X.^2 Y.^2)); % 单次 FFT U fftshift(fft2(fftshift(hologram .* phase))); end这里只做了一次 FFT计算量比卷积法小。但它的重建采样间隔变成了lambda*z/(N*pixel_size)也就是说z越大重建出来的像素间隔越大视场越大但分辨率越低。这个特性在需要大视场观察时很有用但如果你要精确测量微小结构就得注意分辨率损失。2.4 角谱法第三种重构的补充虽然标题里只提了 CONV 和 Fresnel但很多数字全息项目会带上角谱法作为对比。角谱法在频域里传播不做近轴近似适合短距离或者大角度的情况。代码结构通常是function U recon_angular(hologram, lambda, pixel_size, z) [M, N] size(hologram); k 2*pi/lambda; fx (-N/2:N/2-1) / (N*pixel_size); fy (-M/2:M/2-1) / (M*pixel_size); [FX, FY] meshgrid(fx, fy); % 传递函数 H exp(1j*k*z*sqrt(1 - (lambda*FX).^2 - (lambda*FY).^2)); H(lambda*FX 1 | lambda*FY 1) 0; % 倏逝波截止 U ifft2(fft2(hologram) .* H); end角谱法的传递函数里有个根号项当lambda*FX超过 1 时变成倏逝波代码里直接置零。这个方法在z很小的时候比 Fresnel 更准但计算量稍大。三种方法放在一起你可以根据z的大小和需要的分辨率来选。一般经验是z小、要求精度高用角谱z中等、想快速看结果用 Fresnelz大、想保持采样间隔不变用 CONV。3. 跑通第一个重建从读图到出图的完整操作3.1 环境准备与数据放置这个项目对 MATLAB 版本要求不高R2016b 以上都能跑因为用到的都是基础函数。不需要额外工具箱fft2、ifft2、fftshift都是自带的。唯一要注意的是图像读取如果全息图是彩色相机拍的记得转灰度。把全息图放到data/目录下命名成holo.bmp或者改主脚本里的路径。我一般会先确认图像尺寸因为后面所有坐标计算都依赖[M, N]。% 检查图像尺寸和数据类型 info imfinfo(data/holo.bmp); disp([图像尺寸: , num2str(info.Width), x , num2str(info.Height)]);如果图像太大比如 4000x3000FFT 会吃内存。常见做法是裁剪到 1024x1024 或者 2048x2048 再跑。裁剪的时候注意把感兴趣的区域放在中心因为后面fftshift会假设零频在中心。3.2 参数配置与第一次运行主脚本里通常有一个参数区你需要根据实际光路填三个值波长、像素尺寸、重建距离。波长看激光器标签像素尺寸看相机手册重建距离用卷尺量或者用自动对焦算法找。% 第一次运行建议先用小图测试 lambda 532e-9; pixel_size 3.45e-6; z 0.15; hologram imread(data/holo.bmp); hologram double(hologram(:,:,1)); % 调用三种方法 U_conv recon_conv(hologram, lambda, pixel_size, z); U_fresnel recon_fresnel(hologram, lambda, pixel_size, z); U_angular recon_angular(hologram, lambda, pixel_size, z); % 显示振幅 figure; subplot(1,3,1); imagesc(abs(U_conv)); title(CONV); colormap gray; axis image; subplot(1,3,2); imagesc(abs(U_fresnel)); title(Fresnel); colormap gray; axis image; subplot(1,3,3); imagesc(abs(U_angular)); title(Angular); colormap gray; axis image;跑完这段你会看到三张振幅图。如果三张图都模糊大概率是z不对如果只有某一张模糊那可能是该方法在这个距离下不适用。比如z很小的时候Fresnel 的近似误差会变大角谱和 CONV 更可靠。3.3 结果对比与参数扫描想确认哪个方法更适合你的数据可以做一个z扫描。把z从 0.05 到 0.5 步进 0.01每个距离算一遍看哪个距离下图像最锐利。这个操作在数字全息里叫自动对焦常用评价函数是梯度平方和或者拉普拉斯算子。z_range 0.05:0.01:0.5; sharpness zeros(size(z_range)); for i 1:length(z_range) U recon_conv(hologram, lambda, pixel_size, z_range(i)); sharpness(i) sum(sum(gradient(abs(U)).^2)); end [~, idx] max(sharpness); disp([最佳重建距离: , num2str(z_range(idx)), m]);这段代码用梯度平方和作为清晰度指标峰值对应的z就是最佳距离。注意不同方法得到的最佳z可能略有差异因为它们的近似条件不同。我一般会以角谱法的结果作为基准因为它的近似最少。4. 避坑与排查数字全息重建里最容易翻车的几个点4.1 重建距离符号搞反像面跑到反方向现象重建出来的图像完全模糊怎么调z都不清晰。原因z的符号决定了衍射传播方向如果物平面和全息面的相对位置搞反了相当于往反方向传播。解决把z取负再试一次或者检查光路里物和相机的先后顺序。常见做法是先用一个已知距离的样本标定符号。4.2 像素尺寸填错重建尺度全乱现象重建出来的物体尺寸和实际对不上比如一个 100 微米的球重建出来变成 200 微米。原因pixel_size填的是相机像素尺寸但如果你用了显微物镜或者缩束系统实际物平面的采样间隔是pixel_size / magnification。解决确认光路里有没有放大或缩小把放大倍率算进去。我一般会在代码里单独定义一个mag变量避免直接改pixel_size。4.3 FFT 混叠导致高频伪影现象重建图里出现周期性条纹或者网格状伪影。原因全息图的高频信息超过了采样定理或者fftshift用错了位置。解决先检查全息图本身有没有欠采样如果条纹太密说明载频太高需要调整光路或者用更小的像素尺寸。另外fft2之前记得fftshift之后记得ifftshift顺序错了会引入相位斜坡。4.4 零级项和共轭项没滤掉现象重建图里除了物体像还有一团亮斑或者对称的鬼影。原因全息图里的零级衍射和共轭项没有分离直接做重建会把它们一起传过去。解决在频域里加一个带通滤波器把载频附近的一个边带选出来。常见做法是手动找载频峰值然后用一个圆形窗口滤出来。% 频域滤波示例 H fftshift(fft2(hologram)); [rows, cols] size(H); [X, Y] meshgrid(1:cols, 1:rows); center [cols/2, rows/2]; radius 50; % 根据载频位置调整 mask ((X - center(1)).^2 (Y - center(2)).^2) radius^2; H_filtered H .* mask; hologram_filtered ifft2(ifftshift(H_filtered));这段代码在频域里选了一个圆形区域只保留一个边带。radius需要根据实际载频的展宽来调太小会丢信息太大会把零级带进来。4.5 内存不够导致 FFT 报错现象跑大图的时候 MATLAB 提示内存不足。原因fft2对 4000x4000 的复数矩阵会占用几百 MB加上中间变量很容易爆。解决分块处理或者降采样。我一般会先把图裁到 2048x2048 以内如果非要大视场就用single精度代替double能省一半内存。5. 二次开发与改造把三种重构封装成可复用的工具链5.1 统一接口与参数校验原始项目里三个函数各自独立参数传递靠位置。如果你想把它集成到更大的 pipeline 里建议统一接口加一个method参数来切换。同时加上参数校验避免传错顺序导致玄学错误。function U reconstruct(hologram, lambda, pixel_size, z, method) % 统一入口method 可选 conv、fresnel、angular validateattributes(lambda, {double}, {positive}); validateattributes(pixel_size, {double}, {positive}); validateattributes(z, {double}, {nonzero}); switch lower(method) case conv U recon_conv(hologram, lambda, pixel_size, z); case fresnel U recon_fresnel(hologram, lambda, pixel_size, z); case angular U recon_angular(hologram, lambda, pixel_size, z); otherwise error(不支持的重构方法: %s, method); end end这样改完之后上层代码只需要调reconstruct换方法不用改调用方式。参数校验能挡住大部分低级错误比如把z传成 0 导致除零。5.2 自动对焦与批量处理如果你有一批全息图要处理手动调z不现实。可以把第 3 章里的对焦扫描封装成函数自动返回最佳距离和重建结果。再配合dir遍历文件夹就能批量出结果。function [U_best, z_best] autofocus_reconstruct(hologram, lambda, pixel_size, z_range, method) sharpness zeros(size(z_range)); for i 1:length(z_range) U reconstruct(hologram, lambda, pixel_size, z_range(i), method); sharpness(i) sum(sum(gradient(abs(U)).^2)); end [~, idx] max(sharpness); z_best z_range(idx); U_best reconstruct(hologram, lambda, pixel_size, z_best, method); end这个函数返回最佳距离和对应的复振幅。批量处理的时候把每个全息图的结果存成 MAT 文件或者 PNG方便后续分析。5.3 相位解包裹的衔接数字全息重建出来的相位是包裹在[-pi, pi]之间的直接看会有跳变。如果你要做定量相位测量需要接一个解包裹算法。MATLAB 里可以用unwrap函数做一维解包裹二维的话常见做法是用质量引导或者最小二乘法。我一般会在重建之后加一步相位解包裹再转成高度或者厚度。phase angle(U_best); phase_unwrapped unwrap(unwrap(phase, [], 1), [], 2); % 转成高度反射式 height phase_unwrapped * lambda / (4*pi);这段代码先按行解包裹再按列解包裹适用于连续变化的相位。如果噪声大解包裹会出错这时候需要先做相位滤波比如用均值滤波或者中值滤波。5.4 性能优化GPU 加速与内存复用如果数据量大FFT 是瓶颈。MATLAB 支持gpuArray把全息图传到 GPU 上做 FFT 能快好几倍。但要注意显存限制大图还是会爆。另一个技巧是复用 FFT 结果比如 CONV 和角谱都需要fft2(hologram)可以算一次存起来。if canUseGPU hologram_gpu gpuArray(hologram); U gather(ifft2(fft2(hologram_gpu) .* H_gpu)); else U ifft2(fft2(hologram) .* H); end这段代码先判断有没有 GPU有就传上去算算完再取回来。H_gpu是提前算好的传递函数避免每次重复计算。5.5 从脚本到函数库的改造习惯我刚开始写这类项目的时候喜欢把所有代码堆在一个脚本里改一个参数要翻半天。后来强迫自己把每个重构方法拆成独立函数参数用结构体传主脚本只负责调度和画图。这样改的好处是下次遇到新的重构算法只需要加一个函数文件不用动主逻辑。从那以后我每次拿到新数据都先跑一遍三种方法的对比图确认哪种最适合当前的光路参数再往下做定量分析。希望帮到你。本文还有配套的精品资源点击获取