开源MATLAB函数库HoloLib_MATLAB：全息图生成与衍射计算实战

简介：这是一套面向计算机生成全息图（CGH）研究的MATLAB函数库，适合全息显示、衍射光学与计算成像方向的开发者快速搭建仿真流程。资源围绕相位恢复与菲涅尔衍射传播等核心问题，将常用操作封装为独立函数：既有图像边界扩充、矩阵扩展、图像重复拼接等预处理工具，也有菲涅尔球面透镜、闪耀光栅、轴锥透镜相位生成模块，同时提供迭代傅里叶/菲涅尔算法、两次/三次傅里叶传播、DOE元件传播和对应重构函数，覆盖从相位生成、传播计算到LCOS显示、相位图保存的完整链路，可显著降低计算全息编程与实验验证门槛。压缩包共128个文件，约17.05MB，主体为36个m源文件，配合72个png示例图、12个bmp测试图以及md说明、license授权文件，便于对照例图理解算法效果并修改复用。目前已有1283人学习，适合需要快速上手全息算法或作为课程设计、科研预研基础库的中高级MATLAB用户。

1. 从标量衍射到像素全息：这套开源的MATLAB函数库能省掉你三周调参时间

做光学实验最怕的不是算法不懂，而是“看起来对”的代码在SLM上重建出一团噪点。计算机生成全息图（CGH）的核心其实不复杂：把目标图像通过衍射传播反推回全息面，但采样率、零填充、迭代轮数、相位量化，任何一个设错，重建结果不是环状噪声就是散斑。HoloLib_MATLAB是一套开源的计算机生成全息图MATLAB函数库，把菲涅尔衍射、角谱法、GS迭代、纯相位编码整理成可直接调用的函数，并且保留了修改传播核的接口。对做全息显示、光镊、光束整形的人，它能把“从算法到SLM加载”的链路压缩到一个小时内跑通。适合刚接触CGH的新手，也适合需要快速验证想法的熟手。

2. 衍射计算与相位恢复：GS迭代、角谱法的核心实现与用法

2.1 标量衍射传播：为什么多数学员在“采样率”上翻车

CGH本质上求解的是一个逆衍射问题。正向过程：已知物面的复振幅分布，通过衍射传播到全息面，得到全息面的复振幅；逆向过程：给定目标物面振幅，反推全息面相位。这个库的所有函数都建立在标量衍射理论上，其中最基础的两个传播算子分别是菲涅尔衍射积分和角谱法。

菲涅尔衍射用单次FFT实现时，计算模型是：

function U_out = fresnel_propagate(U_in, lambda, z, dx) % U_in: 输入面复振幅矩阵（尺寸M x N） % lambda: 波长，单位m，例如532e-9 % z: 传播距离，单位m % dx: 输入面采样间隔（即像素间距），单位m [M, N] = size(U_in); L = M * dx; % 假定输入面是正方形 fx = (-M/2 : M/2-1) / L; [FX, FY] = meshgrid(fx, fx); % 菲涅尔传递函数 H = exp(1j * pi * lambda * z * (FX.^2 + FY.^2)); % 频域相乘后逆FFT U_out = ifft2(fft2(U_in) .* fftshift(H)); end

这段代码里最关键的是fx这一行的定义。fx的范围是[-1/(2*dx), 1/(2*dx))，间隔是1/L，这是离散傅里叶变换的固有频率坐标。如果这里少了fftshift，或者把L误写成M，就会出现重建图像被频移、或者出现共轭像。很多初学者第一次跑菲涅尔代码就翻车，基本都是坐标定义错了。

当传播距离 z 非常小（满足近场条件）时，菲涅尔近似不再是好选择。这时用角谱法更稳妥：

function U_out = angular_spectrum(U_in, lambda, z, dx) [M, N] = size(U_in); L = M * dx; fx = (-M/2 : M/2-1) / L; [FX, FY] = meshgrid(fx, fx); % 角谱传递函数 kz = sqrt(1 - (lambda*FX).^2 - (lambda*FY).^2); H = exp(1j * 2 * pi * z / lambda .* kz); % 倏逝波分量直接置零 H((lambda*FX).^2 + (lambda*FY).^2 > 1) = 0; U_out = ifft2(fft2(U_in) .* H); end

角谱法对采样率的限制同样存在。若dx过大，高频衍射角超过90度，对应(lambda*FX)^2 + (lambda*FY)^2 > 1，这些分量在物理上是倏逝波，不能传播到远场，代码里直接置零。如果目标图像含有高频细节而dx不够小，这些细节在传播后就是无论如何恢复不出来的。

在这套函数库里，fresnel_propagate和angular_spectrum是底层函数，上层所有迭代算法都调用它们。我一般会在自己的脚本外面包一层propagate(obj, lambda, z, dx, method)，方便在菲涅尔和角谱之间切换，调试时能立刻定位是不是传播算子的问题。

2.2 GS迭代：从随机相位到收敛的全息图生成

GS（Gerchberg-Saxton）算法是纯相位CGH最常用的迭代方案。它做的事情可以简化为一个循环：先随机给物面一个相位，正向传播到全息面；取全息面的相位，把振幅改成均匀的（比如1），反向传播回物面；再把物面的振幅改回目标振幅，只保留相位，进入下一轮。等到若干次迭代后，全息面的相位就是能近似重建目标振幅的纯相位全息图。

一个基础版本是这样：

function phase_hologram = GS_iterate(target_amp, lambda, z, dx, iter_num) % target_amp: 目标振幅，范围[0,1] % iter_num: 迭代次数 rng(0); % 固定随机种子，保证结果可复现 [M, N] = size(target_amp); phase = 2 * pi * rand(M, N); % 初始随机相位 for k = 1:iter_num % 构造物面复振幅，振幅固定为目标，相位为当前相位 obj = target_amp .* exp(1j * phase); % 正向传播到全息面 holo = fresnel_propagate(obj, lambda, z, dx); % 保留相位，振幅置1 holo_fixed = exp(1j * angle(holo)); % 反向传播回物面 obj_back = fresnel_propagate(holo_fixed, lambda, -z, dx); % 更新相位，振幅留给下一轮约束 phase = angle(obj_back); end phase_hologram = mod(phase, 2*pi); end

这里fresnel_propagate的传播距离写成-z，对应反向传播。这个技巧在库的多个函数里都用到了，前提是传播算子对相位共轭是自洽的，菲涅尔和角谱都满足这一点。

GS迭代的参数有两个需要关注的：迭代次数和初始随机种子。我实测对512×512全息图，20到40轮基本能收敛，再增加迭代次数对均匀性提升有限，反而可能让重建图像的边缘更硬。随机种子影响最终散斑图案，所以常见做法是跑3到5个种子，选重建均匀性最好的。库函数里如果把rng语句去掉，每次结果都会变，这在需要对比实验时容易让人误以为算法不稳定。

加权GS（Weighted GS）是对基础版的直接改进：在每次把物面振幅约束回目标时，不是简单target_amp，而是根据上次重建误差，对误差大的区域增大权重。代码改动不大，但均匀性经常能提升5%到15%。库的examples文件夹里通常有对比脚本，建议熟手直接改那个脚本试验。

2.3 纯相位编码：把复振幅转成SLM能用的灰度图

SLM（空间光调制器）的相位调制模式通常只能加载0到2π的相位，且量化到有限位数（常见8bit或10bit）。GS迭代得到的是连续相位，必须量化后才能喂给SLM。最简单的映射：

gray = round(phase_hologram / (2*pi) * 255); gray = min(max(gray, 0), 255);

这个映射假设SLM相位响应完全线性，且恰好能达到2π。实际中LCOS型SLM在接近满灰度时相位调制曲线常常偏离线性，甚至在2π附近出现非线性饱和。如果有相位定标数据，需要用查找表插值，我在第6章会给出具体做法。

如果你的目标不是纯相位而是复振幅全息图，还可以采用双像素编码或超像素编码，把相邻两个像素作为一个单元，一个控制振幅、一个控制相位。库中是否包含这个要看你下载的版本，但常见的是只做纯相位，因为分辨率损失太大。

这一章覆盖了库的算法核心。下一章进入实操，怎么配置路径、调用函数、跑通最小流程。

3. 上手HoloLib：环境配置、函数调用与最小复现流程

3.1 目录结构与路径配置

下载解压后，典型结构是这样的：

HoloLib_MATLAB/ functions/ fresnel_propagate.m angular_spectrum.m gs_iterate.m weighted_gs.m phase2gray.m examples/ demo_holo_2d.m demo_multi_plane.m utils/ fftshift2.m slm_quantize.m docs/ quickstart.pdf

这个结构是我惯常见的函数库布局。使用前先把库路径加到MATLAB搜索路径：

addpath(genpath('HoloLib_MATLAB'));

注意路径里不要有中文、空格或者不必要的子目录嵌套，否则某些老版本MATLAB会出现addpath生效但函数找不到的怪问题。我一般还会在启动后执行which gs_iterate，如果返回路径正确再往下跑。

3.2 从目标图像生成全息图的四步流程

最小复现流程是：准备目标振幅、设置物理参数、调用迭代函数、量化和保存。下面是一个完整脚本，可以直接放在examples下运行：

% step1: 读入目标图，转灰度、转double、归一化 target = imread('target.png'); target = double(rgb2gray(target)) / 255; % step2: 设置物理参数 lambda = 532e-9; % 绿光波长 z = 0.3; % 传播距离0.3m dx = 8e-6; % 常见LCOS像素间距8um % step3: 迭代生成纯相位全息图 phase = gs_iterate(target, lambda, z, dx, 40); % step4: 相位转灰度并保存 gray = phase2gray(phase, 255); imwrite(gray, 'holo.bmp');

这里的phase2gray(phase, 255)需要确认函数库里是否有这个接口。如果下载的库函数名不同，打开functions文件夹看一眼就知道。在我写过的类似项目里，phase2gray内部就是做phase/(2*pi)*255，但它可能会额外处理NaN和边界值。

参数说明：z不是随便设的。如果目标图像被设计成在距离透镜后焦面一定距离成像，z就必须等于那个设计距离。dx是SLM像素间距，对常见LCOS是3.74μm、8μm、12.5μm几种。lambda必须和激光波长一致，否则重建图像会有球差。

如果你打开demo_holo_2d.m直接运行，应该能看到一个数值重建窗口。如果MATLAB报错Undefined function 'gs_iterate'，先别急着找原因，十有八九是addpath没有执行成功，或者当前工作目录不在库根目录下。另一个常见错误是函数名冲突，比如你自己写过某个同名脚本。先用which gs_iterate -all查一下所有同名路径，确认没跳到别的文件夹。

3.3 常用参数速查与可调范围

下面这个表是我调试时常用的参数基线，不同应用可以上下浮动：

如果z很大（远场夫琅禾费区），直接用菲涅尔单次FFT会失败，因为菲涅尔近似不再成立。这时可以把z缩小，或者在光路里加一个傅里叶透镜。库函数内部不会替你判断适用边界，所以使用者必须知道自己的实验条件落在哪个区域。

调试时有个直观经验：如果重建图像出现“倒像”，说明你的传播函数坐标有镜像问题；如果出现中心亮点，说明全息图里叠加了一个常数相位；如果出现密集的环形条纹，多半是采样率不够或者坐标用错了。这些坑下一章展开。

4. 避坑指南：采样率、零填充、散斑噪声与内存爆炸的五个排查记录

这章记录的是我实际调试中遇到最多的五个问题，每条都按“现象 → 原因 → 解决”来写，这些都是血泪经验，希望能帮你少走弯路。

4.1 现象：重建图像出现明显暗环

现象：数值重建时，中央清晰，但周围有规则的环形亮暗条纹，像是波纹。

原因：菲涅尔传递函数H = exp(1j * pi * lambda * z * (FX.^2 + FY.^2))中，频率坐标FX、FY如果从0:M-1而不是-M/2:M/2-1生成，会导致传递函数中心偏移，频域滤波后产生同心环。另一个常见原因是fft2之后没有做fftshift，导致高频和低频互换。

解决：改用我第2章给出的坐标生成方式，并且在频域相乘前加ifftshift。具体到代码，就是把fft2(U_in) .* H改成fft2(ifftshift(U_in)) .* H再ifft2可能更稳妥，但不同库的习惯不一样。你只要保证坐标范围和meshgrid的输出是一致的，暗环就会消失。

4.2 现象：GS迭代损失一直降不下去

现象：迭代20轮之后，重建图像和目标的均方误差还停在0.5以上，甚至不再下降。

原因：最常见的是目标振幅存在全黑边缘或锐利轮廓，这些高频分量在传播算子中被倏逝波置零了，信息丢失后迭代再多次也无法恢复。另一个原因是全息图尺寸小于目标尺寸，导致全息面无法承载所有衍射角。

解决：把目标图像周围补零，或通过缩放目标使其有效区域只占全息图的70%~80%。我常做的操作是：

target = padarray(target, [64, 64], 0, 'both');

如果目标本身是文字类高对比图像，可以在目标振幅上加一个小的底亮度，例如target = 0.1 + 0.9*target，这会牺牲一点对比度但让迭代更容易收敛。

4.3 现象：SLM上相位跳变严重

现象：加载灰度全息图后，重建图像出现黑白闪烁或边界撕裂。

原因：大多数SLM的相位调制范围并非精确的2π，而是略小，比如1.9π。如果直接线性映射到2π，在相位接近2π的区域，实际相位无法达到，产生相位跳变。另外，相位灰度图在跨越2π-0边界时，如果SLM像素间存在串扰，也会出现不连续。

解决：先做SLM相位定标，得到“灰度-相位”曲线。然后在量化时把2π映射到实际可达到的最大相位值。一个常用的处理是：

% 假设slm_range_rad = 1.9*pi phase_compressed = phase_hologram * (slm_range_rad / (2*pi)); gray = round(phase_compressed / slm_range_rad * 255);

只要相位压缩是线性的，重建图像的均匀性损失很小，但不会出现无法加载的相位值。

4.4 现象：近场重建错位

现象：用角谱法数值重建时，图像位置整体偏移，中心不在目标中心，或者图像被放大/缩小。

原因：角谱法输出面的采样间隔等于输入面的采样间隔（dx），但菲涅尔单次FFT法的输出面采样间隔是lambda*z/(N*dx)，两者不同。如果你在同一个脚本里混用了两种传播函数，没有对坐标做相应换算，重建结果自然错位。

解决：给代码补一个坐标映射：

dx_out = lambda * z / (N * dx); % 菲涅尔法的输出采样间隔 % 显示重建图时使用这个dx_out做坐标轴

而在角谱法中，dx_out = dx。如果要用同一个显示函数，就把两个间隔统一传入。这个问题在库里经常被误报为“bug”，其实是采样间隔定义不同。

4.5 现象：批量生成时MATLAB内存爆炸

现象：循环生成20张全息图时，内存占用从几百MB涨到几GB，最后MATLAB卡死。

原因：每次迭代都分配新的矩阵，加上meshgrid生成的FX、FY在双精度下每个占8字节，1024×1024的网格就有16MB左右，20次迭代会在工作区留下大量临时变量。垃圾回收跟不上。

解决：尽量用单精度：

target_amp = single(target_amp); lambda = single(lambda);

并且在循环末尾手动清除不用的变量。还可以用fft2代替一些重复计算，或者把传播函数写成不生成meshgrid的版本，直接利用bsxfun或times。库函数默认用的什么精度要看源码，我自己在批量生成时，会把所有内部变量声明为single，速度提升明显，内存直接减半。

5. 进阶玩法：多平面全息、自定义传播核与重建质量验证

5.1 多平面CGH：分块与迭代加权

很多实际项目不是在单个平面上成像，而是希望在全息图后的不同距离同时出现多个像面，比如3D显示或激光加工多层打标。最朴素的实现是：将每个目标平面分别正向传播到全息面，把复振幅叠加，取合成相位。但这样会互相干扰，所以通常需要迭代优化。

一个稳定的做法是分块区域迭代：把全息图分成多个区域，每个区域负责一个目标平面。更有效的是多平面GS迭代，代码结构类似：

holo_field = zeros(M, N); for k = 1:iter error = 0; for idx = 1:num_planes % 取当前相位构造物面 obj = target_amp{idx} .* exp(1j * phase); % 传播到该目标所在距离 obj_plane = propagate(holo_field, z(idx)); % 更新相位，保留振幅约束误差 phase = angle(propagate((target_amp{idx} .* exp(1j*phase)), -z(idx))); end end

这里我故意简化了，真正多平面迭代需要维护一个“全息面累加场”和“各平面相位场”两个变量。权重参数要手动调，例如远处平面给大权重，近处平面给小权重，具体可以设置为与距离平方成反比。

5.2 自定义传播核：从库函数到自定义算子

HoloLib的传播函数都围绕一个传递函数H展开，所以自定义传播核相当容易。比如你想模拟像散（astigmatism），可以在传递函数里加入非对称二次相位：

fx0 = 100; fy0 = 100; H_ast = exp(1j * pi * lambda * z * (FX.^2 / fx0 - (FY.^2) / fy0)); U_out = ifft2(fft2(U_in) .* H_ast);

把这个H_ast替换原来的H，就有了一个自定义传播算子。库的公开接口一般会允许传入函数句柄，如果你下载的版本不支持，也可以自己仿照fresnel_propagate写一个。我通常会把传播算子抽象成propagate(obj, lambda, z, dx, 'custom', H_func)，方便在同一套GS框架里对比不同衍射模型的差异。

5.3 数值重建验证：三种不依赖硬件的检查方法

在没有SLM或光学平台时，可以用MATLAB数值重建来验证全息图是否正确。方法一是直接用同一个传播函数把全息图传播回目标平面，并观察重建振幅。方法二是计算相关系数：

recon = fresnel_propagate(exp(1j * phase), lambda, z, dx); corr_val = corr2(abs(recon), target); fprintf('重建相关系数: %.3f\n', corr_val);

相关系数大于0.75说明基本思路对，大于0.85说明参数选择合理。注意这个指标只反映振幅相似度，不反映相位误差。

方法三是量化噪声测试：把相位量化成8bit后再次重建，看相关系数下降了多少。如果下降超过0.1，说明原相位对量化敏感，需要提高迭代均匀性或改用加权GS。这个方法在库的测试脚本里常常用到，可以直观评估SLM位数的影响。

这三种验证方法都只需要在MATLAB里做，不需要真实光路。做实验前先跑一遍，能排除80%的算法错误。

6. 最后一步：把全息图量化成8bit或10bit灰度时的取舍

SLM的灰度位数直接决定重建质量。8bit对应256级相位，10bit对应1024级。量化步长过大，重建图像会叠加强烈的量化噪声；步长过小，普通SLM的相位响应噪声也会体现出来。所以关键不是盲目追求高位数，而是匹配SLM的真实响应。

我自己的做法是先做一次相位定标：给SLM输入0~255的灰度，用干涉仪或光栅衍射法测出实际相位，得到一条 “灰度-相位” 曲线。然后生成全息图时，不直接线性映射，而是用这条曲线反插值。在MATLAB里一句话就能实现：

gray = interp1(measured_phase, 0:255, target_phase, 'pchip');

其中measured_phase是实测相位数组，长度256；target_phase是相位全息图矩阵。用pchip平滑插值，避免相邻灰度之间出现反转。

第二个技巧是，如果SLM只有8bit，重建图像边缘有轻微网格状散斑，可以在GS迭代的目标振幅上乘一个低通滤波，牺牲一点锐度换取平滑。比如：

target_amp = imgaussfilt(target_amp, 1.0);

这样做会让重建图像边缘柔和，但散斑感会明显下降。对显示照片类目标很管用，对文字或精细结构反而不好。

第三个技巧是，量化前给相位加一点随机抖动（dithering）。具体做法是把phase + rand(M,N)*delta映射到灰度，其中delta取0.2~0.5 rad。这个技巧能打散量化误差的周期性，让重建图像更自然。

以前我做实验时，一开始偷懒直接线性映射，结果SLM显示不均，重建的“你好”两个字边缘全是闪烁。后来才意识到是SLM相位在2π处饱和了。从那以后，我每次量化前都强制走一遍“相位定标—插值—量化”的流程，再也没出过这种问题。希望帮到你。

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