MATLAB计算全息函数库：衍射传播与相位恢复实现

简介：面向计算全息（CGH）研究的 MATLAB 函数库，适合光学实验人员与算法开发者快速搭建全息图生成、传播与重建流程，借助封装函数降低编程门槛。资源共 128 个文件：36 个 .m 源文件涵盖图像边界扩充、矩阵扩展、轴锥透镜/球面透镜/闪耀光栅相位生成，迭代傅里叶与迭代菲涅尔算法，以及基于 DFFT/TFFT/DOE 的多种传播计算，还提供相位图保存、LCOS 显示等工具；72 张 .png 示例图和 12 张 .bmp 测试图可用于对比输出结果，另附 license 和说明文档，压缩包约 17.05MB。已有 1283 人学习。通过该库可获得一套完整可运行的 CGH 计算工具，既可直接调用顶层函数生成指定相位图，也能依据源码注释跟踪每一步数值计算，理解傅里叶变换、菲涅尔衍射等核心原理；内置初始化脚本可自动配置 MATLAB 搜索路径，便于在科研或教学中直接二次开发与验证。

1. 计算机生成全息图，为什么你需要一个整理好的MATLAB函数库

手头散落一堆命名混乱的.m脚本，文件名从GS_v2_final.m到Fresnel_3.m，每个脚本里都复制了一遍衍射传播代码——这是不少人接触计算机生成全息图（CGH）时的真实状态。这套HoloLib_MATLAB函数库的价值，就是把 CGH 生产线上最常用的模块整理成可复用函数：标量衍射传播、GS 迭代相位恢复、相位编码与量化、模拟重建评估，全部用 MATLAB 实现，不用碰 C++ 或 CUDA 也能跑出可用的纯相位全息图。适合刚开始做全息显示、全息光镊或投影系统的从业者，也适合手里有实验脚本但不敢改参数的熟手——把核心算法抽出来做成函数，参数边界一次调明白。

2. 从衍射传播到相位恢复：函数库的模块划分与核心函数拆解

做 CGH 的完整链路不是“输入图片、输出相位”这么簡單。你需要把振幅图像编码成相位分布，再让这个相位分布经过真实光学传播后在目标平面重建出预期光场。这个链路里至少涉及三件事：光场怎么传播、相位怎么迭代恢复、相位图怎么适配显示器件。函数库把这三件事拆成了独立模块，下面逐个拆。

2.1 模块一：衍射传播算子，CGH 的地基

光场从全息面传播到重建面，MATLAB 实现里最常见是三种：菲涅尔衍射、角谱法（ASM）、傅里叶变换法。函数库对这三种都有封装，命名通常类似holo_fresnel_prop、holo_asm_prop、holo_ft_prop。三者的适用场景差异很大，我一般按这张表选：

一个常见误用是把菲涅尔单次FFT当成万能函数，距离设得很小时直接失真。原因在于单次FFT形式的菲涅尔传播隐含了一个约束：输出平面采样间隔与输入采样间隔不同，且与波长、距离、像素数绑定。距离过小时，输出采样间隔被压得极小，重建区域随之变小，边缘出现混叠。函数库里holo_fresnel_prop的入参设计了z、lambda、dx、N四个量，目的就是让你显式意识到采样间隔不是随便给的。

2.2 模块二：GS 迭代与相位恢复的核心实现

拿到目标振幅图后，相位不可能直接逆衍射解出来——这是个欠定问题。GS 算法的思路是反复在物面与频谱面之间来回投影：正向传播后保留振幅、替换相位；反向传播后同样保留振幅、替换相位，直到两个平面的振幅都逼近目标。函数库里对应函数通常是holo_gs_phase，典型调用结构如下：

% 目标振幅: 256x256 灰度图，范围 [0,1] target_amp = double(imread('target.png')) / 255; target_amp = imresize(target_amp, [256 256]); % 初始化随机相位 rng(42); phase_init = exp(1j * 2 * pi * rand(256)); % 调用 GS 迭代，20 次迭代，菲涅尔模式传播 [phase_out, err_history] = holo_gs_phase(target_amp, ... 'method', 'fresnel', 'iterations', 20, ... 'lambda', 532e-9, 'z', 0.2, 'dx', 8e-6, ... 'init_phase', phase_init);

这段代码的逻辑并不复杂：目标振幅归一化后，随机相位作为起始猜测，holo_gs_phase内部做正向传播、振幅约束替换、反向传播、振幅约束替换，迭代过程中把每个周期末重建面与目标之间的均方误差记录到err_history。参数上最值得留意的是init_phase——相同目标、相同迭代次数，换一个随机种子结果差别很大，这就是 GS 的局部最优问题。函数库保留这个入参，就是为了让你能对比不同初始化的重建质量。

2.3 模块三：相位编码与量化，从浮点相位到 SLM 可用的灰度图

GS 迭代输出的相位是连续浮点数，范围通常是[-pi, pi]或[0, 2*pi]，但空间光调制器（SLM）或数字全息打印机接受的是 8bit 灰度图。这个转换不是简单的im2uint8——要考虑相位包裹边界、灰度反转、γ 校正。函数库里的holo_slm_encode通常处理三件事：把浮点相位归一化到[0, 1]，按 SLM 的灰度-相位响应曲线做映射，再叠加可选的闪耀光栅偏移。

phase_wrapped = angle(phase_out); % 包裹到 [-pi, pi] phase_norm = (phase_wrapped + pi) / (2*pi); % 归一化到 [0, 1] slm_gray = uint8(round(phase_norm * 255)); % 默认线性映射到 8bit % 若 SLM 实测灰度-相位曲线非线性，传入查找表 load('slm_lut.mat', 'lut'); % lut 是 256x1 的相位值向量 phase_mapped = interp1(lut, 0:255, phase_norm(:), 'linear'); slm_gray = reshape(uint8(round(phase_mapped * 255)), size(phase_norm));

这里的关键坑在于：把包裹后的angle直接归一化到灰度，假设了 SLM 工作在[0, 2*pi]线性相位模式。实际很多 SLM 只覆盖[0, 2*pi]的一部分，或者灰度与相位之间有明显的非线性，这时候不做 LUT 校正，重建图像会出现背景噪声被放大、对比度下降的问题。函数库把这个环节独立成模块，而不是塞进 GS 迭代里，就是为了让你在换用不同 SLM 时只改编码函数，不动相位恢复算法。

3. 跑通第一个 CGH：环境配置、GS 迭代与最小复现路径

拿到函数库后第一件事不是读文档，是把它跑起来。这个库纯 MATLAB 实现，不依赖第三方工具箱，但我建议在动手前先确认 MATLAB 版本和并行计算选项，然后按下面路径走通一次完整流程。

3.1 环境准备：路径设置与依赖检查

% 把整个 HoloLib_MATLAB 仓库加入 MATLAB 路径 addpath(genpath('./HoloLib_MATLAB')); % 确认核心函数可见 assert(exist('holo_gs_phase', 'file') == 2, '核心函数缺失，检查仓库完整性'); assert(exist('holo_fresnel_prop', 'file') == 2, '菲涅尔传播函数缺失'); % 检查是否有 GPU（可选，不强制） try gpuDeviceCount; use_gpu = gpuDeviceCount > 0; catch use_gpu = false; end fprintf('GPU 可用: %d\n', use_gpu);

addpath(genpath(...))会把仓库下所有子目录一次性加入搜索路径，省去逐个目录添加的麻烦。exist(..., 'file') == 2用来验证核心函数是否真的在路径里——不少人是 clone 下来后发现子模块没拉全，函数找不到。GPU 检查那段是可选优化，CGH 迭代里 FFT 占大头，如果你的数据是 512×512 以上，GPU 能提速 3～5 倍，但函数库默认走 CPU，use_gpu这个logical变量后面可以传给传播函数。

3.2 最小复现：从目标图到相位图只需五步

这里给一个最简流程，不做任何花哨优化，保证你照着跑能出结果。目标图像用函数库自带的示例图，或者自己生成一张简单图案。

% 第一步：生成目标振幅（100x100 的圆环图案） N = 256; [x, y] = meshgrid(linspace(-1, 1, N)); target = double((x.^2 + y.^2) > 0.3 & (x.^2 + y.^2) < 0.4); % 第二步：调用 GS 迭代 phase = holo_gs_phase(target, ... 'method', 'fresnel', 'iterations', 30, ... 'lambda', 532e-9, 'z', 0.3, 'dx', 10e-6, ... 'verbose', true); % 第三步：模拟重建，验证相位是否正确 recon = holo_fresnel_prop(exp(1j * phase), ... 'lambda', 532e-9, 'z', 0.3, 'dx', 10e-6); % 第四步：取重建振幅并归一化 recon_amp = abs(recon); recon_amp = recon_amp / max(recon_amp(:)); % 第五步：对比目标与重建 corr_val = corr2(target, recon_amp); fprintf('重建相关系数: %.4f\n', corr_val);

这段代码是把整个管线的最小闭环串起来了。第五步的corr2是衡量重建质量最直观的指标，取值 0 到 1，通常在 0.8 以上就能看出图像轮廓。注意第一步生成的target是逻辑矩阵，直接转成double后值域[0,1]，GS 内部对振幅的假设就是归一化的，传入非归一化数据会导致迭代不稳定。

3.3 重建结果怎么看：相关系数与视觉检查配合

只盯corr2数值容易误判。我见过某开发者的重建相关系数 0.85，图像却明显有中心光斑污染，原因在于corr2对全局灰度趋势敏感，对局部细节损失不敏感。所以模拟重建后一定要同时做三件事：看强度分布、看剖面线、看误差图。

figure('Position', [100 100 1200 400]); subplot(1,3,1); imshow(target); title('目标振幅'); subplot(1,3,2); imshow(recon_amp); title('模拟重建振幅'); subplot(1,3,3); imagesc(abs(recon_amp - target)); colorbar; title('绝对误差');

误差图的用处比相关系数大得多。如果误差集中在图像边缘，说明采样间隔偏大或传播距离与像素数不匹配；如果误差呈现周期性条纹，大概率是 FFT 混叠或零填充不足。函数库在这里的价值是让你快速对比不同参数组合下的误差分布，而不是只给一个标量得分。

4. 把参数调明白：采样间隔、迭代次数与 SLM 适配的边界

GS 迭代跑通容易，调好很难。这一章把三个最常出问题的参数讲透：采样间隔与传播距离的匹配、迭代次数的收敛判据、SLM 硬件参数怎么换算进算法。

4.1 采样间隔与传播距离：一个公式决定重建视场

菲涅尔单次 FFT 传播里，输出平面的采样间隔dx_out由lambda * z / (N * dx_in)决定。这意味着你在算法里设的dx不是随便填的——它和重建图像的实际物理尺寸直接挂钩。

% 参数换算示例 lambda = 532e-9; % 波长，单位米 z = 0.3; % 传播距离，单位米 N = 256; % 像素数 dx_in = 10e-6; % 输入平面采样间隔（对应SLM像素间距） dx_out = lambda * z / (N * dx_in); fprintf('重建平面采样间隔: %.2f um\n', dx_out * 1e6); fprintf('重建平面物理尺寸: %.2f mm\n', dx_out * N * 1e3);

这段换算公式是 CGH 里最容易被忽略的。假设 SLM 像素间距是 10μm，波长 532nm，传播 0.3m，256 像素的重建面物理尺寸约 4.1mm——你的目标图像必须在这个尺寸内才能完整重建。如果你把目标图放大到 8mm 宽，直接用单次 FFT 形式算，边缘区域必然混叠。这时候要么增大N（零填充），要么换用角谱法。函数库中holo_fresnel_prop的入参把lambda、z、dx、N全部暴露出来，目的就是逼你在跑之前把这个公式算清楚。

4.2 迭代次数：不是越多越好，要盯收敛曲线

GS 迭代常见默认值是 20～50 次。但迭代次数和图像内容强相关：目标图案细节多（比如人像照片），50 次往往不够；图案简单（比如圆环、字符），15 次就能收敛。函数库的holo_gs_phase返回了err_history，你应该每次跑完都画一下收敛曲线。

[phase_out, err_history] = holo_gs_phase(target, ... 'method', 'fresnel', 'iterations', 50, ... 'lambda', 532e-9, 'z', 0.3, 'dx', 10e-6, ... 'verbose', false); figure; plot(1:length(err_history), err_history, 'o-'); xlabel('迭代次数'); ylabel('归一化均方误差'); grid on;

观察收敛曲线的经验：前 10 次误差下降最快，20 次以后曲线趋于平缓，如果 30 次后还有明显抖动，说明迭代进入振荡区。振荡的常见原因是目标振幅动态范围太大，GS 无法同时满足亮区和暗区的约束。我一般会在 GS 前对目标振幅做个压缩：target_compressed = target.^0.5，把高亮区域权重降低，改善重建均匀性。这个操作原理和加权 GS 类似，但实现成本为零。

4.3 SLM 适配：像素间距、灰度量化与闪耀光栅

算法侧调好了，硬件侧还有三个参数必须对得上。SLM 像素间距直接决定物理重建尺寸，8bit 灰度意味着相位只有 256 级量化，闪耀光栅是让重建图像偏离零级衍射斑的手法。

灰度量化这件事经常被新手忽略。连续相位被四舍五入到 256 级时，等价于给相位叠加了一个随机量化噪声，这个噪声在重建面表现为散射背景。函数库里的holo_slm_encode有一个可选的'quantize_bit'参数，设成 8 就是模拟 8bit SLM 的量化效果，设成Inf则跳过量化——方便你对比量化前后重建质量的差距。我习惯的做法是先用Inf跑出一版理想结果，再切到 8bit 看损失，如果损失超过 10%，就要考虑换 10bit SLM 或改用混合灰度编码。

5. 避坑指南：五个高频翻车场景与排查路径

CGH 的坑不少，大多是参数边界没摸清导致的。这五个翻车场景是我和几个同事在项目里反复踩过的，按“现象→原因→解决”写清楚。

5.1 重建图像中心偏移或上下颠倒

现象：模拟重建时目标图案出现在视场角落，或者上下左右翻转。

原因：FFT 的频谱中心与数组索引零点错位。MATLAB 的fft2输出低频在数组四角，如果不做ifftshift，传播结果的空间位置就整体偏移。翻转问题则是坐标定义不一致——菲涅尔传播函数里输入输出坐标轴方向反了。

解决：在每个传播函数内部统一做fftshift/ifftshift配对。如果你自己写传播代码，记住这条铁律：正向fft2前ifftshift，fft2后fftshift；逆向反过来。函数库里已经处理好的版本不存在这个问题，但如果你把phase单独导出来用别家代码做重建，一定要检查对方是否做了同样的 shift。

5.2 重建图像对比度差、背景出现光斑

现象：暗区不够暗，亮区不够亮，整个重建像蒙了一层雾。

原因：GS 迭代对振幅动态范围大的目标收敛困难，暗区能量泄漏到背景。另一个常见原因是目标振幅里有低频渐变（比如从中心到边缘亮度递减），GS 对这种平滑分布容易陷入局部最优。

解决：对目标振幅做非线性压缩，常见做法是target^0.5或target^0.3，压缩越狠暗区表现越好但亮区细节损失越大。这个参数要按你的图像内容试，没有通用值。函数库的holo_gs_phase接受'target_compression'参数，传0.5就是做平方根压缩，不用你手动改目标图。

5.3 换了一台 SLM，重建质量急剧下降

现象：同一份相位图，在甲 SLM 上效果不错，换到乙 SLM 上背景噪声暴涨。

原因：两台 SLM 的灰度-相位响应曲线不同，甚至灰度反转方向都可能相反。把[0,1]线性映射到灰度，只适用于理想线性器件，实际器件都有程度不等的非线性。

解决：用函数库的 LUT 接口。先测出当前 SLM 的灰度-相位曲线，生成 256×1 的查找表，让holo_slm_encode在编码阶段做非线性补偿。没有实测条件的话，至少确认灰度 0 对应0相位还是2*pi相位，方向反了重建图像会出现背景增强、目标消失的现象。

5.4 传播距离增大后，重建图像尺寸对不上

现象：算法里设z=0.5，重建出来的图像物理尺寸和手算对不上，差了两三倍。

原因：菲涅尔单次 FFT 的输出采样间隔dx_out = lambda*z/(N*dx_in)随距离线增大，图像尺寸也随之增大，但如果z超过某个临界值，输出像素数不变但采样间隔过大，图像直接超出计算窗口，边缘被截断。

解决：把N翻倍做零填充，同时保持dx_in不变。这样dx_out减半，视场尺寸不变但分辨率提高。代价是 FFT 计算量涨了四倍。函数库的传播函数里设了'padding'参数，0是默认不填充，1表示填充到两倍尺寸。我一般遇到距离超过 0.5m 的情况就直接开填充。

5.5 GS 迭代 30 次后误差不再下降

现象：收敛曲线前 10 次快速下降，然后平台期，误差恒定在某个值不再动。

原因：GS 是确定性迭代，遇到局部最优就卡住。局部最优的产生通常与目标振幅的对称性有关——对称图案（圆、方块）特别容易让 GS 陷入对称分布的相位解。

解决：换随机种子重新初始化，或者改用加权 GS。函数库里有'method', 'weighted_gs'选项，核心思想是对重建误差大的区域加大权重，强制迭代过程优先优化这些区域。经验值是把权重设成误差的 1.5 次方，效果最明显。如果你不想改方法，最简单的手段是换rng种子跑三次，取重建误差最小的一次。

6. 用模拟重建验证全息图：收敛曲线与评价指标的实用技巧

相位图生成后，不要直接上光路，先在 MATLAB 里把重建质量验证一遍。这个习惯能帮你省掉大量调光路的时间。我每次拿到新的目标图像，都会强制跑一遍下面这套验证流程。

第一件事是迭代次数扫描。固定其他参数不变，跑[5, 10, 20, 50, 100]五组迭代，把每组最终的重建相关系数画成曲线。这条曲线的价值在于告诉你“再多迭代还有没有意义”——如果 20 次到 100 次相关系数只涨了 0.01，说明 20 次就是性价比拐点，以后批量跑实验就固定用 20 次。如果曲线还在明显上升，说明你的图像复杂度超过了当前迭代策略的承载能力，这时候不是继续加迭代次数，而是应该换加权 GS 或改初始相位。

第二件事是量化影响评估。用holo_slm_encode的'quantize_bit'分别设Inf和8跑两次重建，算出两次结果的相关系数差。这个差值如果超过 0.05，你的实验结果会受到 SLM 灰度量化的明显干扰，需要在相位编码环节加抖动处理。

第三件事是 LUT 验证。把 SLM 实测灰度-相位曲线输入编码函数，对比线性映射和 LUT 校正后的重建误差图。我经历过一次：线性映射下重建还能看出图案轮廓，但边缘有规律性条纹；LUT 校正后条纹消失。这个差异用相关系数只能看出一丝差别，看误差图才非常明显。从那以后我每次换 SLM 都强制走一遍 LUT 验证流程，再没有出现过“换设备后效果全变”的翻车。希望这套验证路径帮你在上光路前就把问题拦下来。

本文还有配套的精品资源，点击获取